Earth’s figure and gravitation - Earth Dynamics: Deformations and Oscillations of the Rotating Earth - page 383

C Set geocentre values.
RI(1)=R(1)
RHOI(1)=RHO(1)
RHOBI(1)=RHOB(1)
FI(1)=F(1)
C Construct new equally spaced model.
DO 16 I=1,NSM1
IP1=I+1
RI(IP1)=H*DFLOAT(I)
CALL INTPL(RI(IP1),RHOI(IP1),NM,C,R,RHO,M1)
CALL INTPL(RI(IP1),RHOBI(IP1),NM,C,R,RHOB,M1)
CALL INTPL(RI(IP1),FI(IP1),NM,C,R,F,M1)
CALL DERIV(RI(IP1),FPI(IP1),NM,C,R,F,M1)
16 CONTINUE
C Begin integration for the calculation of the difference between
C the axial and equatorial moments of inertia of the inner core.
C Set weights and abscissae for six-point Gaussian integration.
W(1)=0.171324492379170D0
W(2)=0.360761573048139D0
W(3)=0.467913934572691D0
W(4)=W(3)
W(5)=W(2)
W(6)=W(1)
X(1)=-0.9324695142031521D0
X(2)=-0.6612093864662649D0
X(3)=-0.2386191860831970D0
X(4)=-X(3)
X(5)=-X(2)
X(6)=-X(1)
C Set exponent in integrand.
EX=4.D0
C Initialize starting value of integral INT.
INT(1,1)=0.D0
C Form integrand function Q.
DO 17 I=1,NS
Q(I)=RHOI(I)*(5.D0*FI(I)+RI(I)*FPI(I))
17 CONTINUE
C Perform quadrature.
CALL QUAD(INT,EX,Q,1,W,X,NS,C,RI,1,M1)
C Calculate difference between axial and equatorial moments of inertia.
DO 18 I=1,NS
CMA(I)=8.D0*PI*INT(I,1)/15.D0
18 CONTINUE
C Calculate enclosed mass as a function of radius.
DO 19 I=1,NS
MASS(I)=4.D0*PI*(RI(I)**3)*RHOBI(I)/3.D0
19 CONTINUE
C Calculate mass of the inner core as though its density was that
C at the bottom of the outer core.
MASSIP=4.D0*PI*(RI(NS)**3)*RHOBOC/3.D0
C Find coefficient of gravity restoring torque.
GAMMA=2.D0*CMA(NS)*(G*MASS(NS)*(2.D0*FI(NS)+RI(NS)*
1 FPI(NS))/(RI(NS)**3)+3.D0*G*MASSIP*FI(NS)/(RI(NS)**3)
2 +5.D0*(OMEGA**2)/2.D0)/5.D0
C Write out radius, density, 1/f, f', C-A, and mass for inner core.
C Write out headings.
WRITE(2,20)
WRITE(6,20)

20 FORMAT(//10X,'Radius',8X,'Density',5X,'1/f',5X,'df/dr ',5X,

1 'C-A',7X,'Mass')
WRITE(2,21)