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

RAD=R(I)
C Set cumulative sums to zero.
SUM(1)=0.D0
SUM(2)=0.D0
DO 11 J=1,NIS
C Determine current Y vector.
Y(1)=YY(1,I)+SUM(1)
Y(2)=YY(2,I)+SUM(2)
CALL DYDR(YP,RAD,Y,IC,RHS,RHO,RHOB,E,N1,C,R,M1)
C Find vector K1.
K1(1)=H*YP(1)
K1(2)=H*YP(2)
C Increment RAD and Y.
RAD=RAD+H/2.D0
Y(1)=Y(1)+K1(1)/2.D0
Y(2)=Y(2)+K1(2)/2.D0
CALL DYDR(YP,RAD,Y,IC,RHS,RHO,RHOB,E,N1,C,R,M1)
C Find vector K2.
K2(1)=H*YP(1)
K2(2)=H*YP(2)
C Increment Y.
Y(1)=Y(1)+K2(1)/2.D0
Y(2)=Y(2)+K2(2)/2.D0
CALL DYDR(YP,RAD,Y,IC,RHS,RHO,RHOB,E,N1,C,R,M1)
C Find vector K3.
K3(1)=H*YP(1)
K3(2)=H*YP(2)
C Increment RAD and Y.
RAD=RAD+H/2.D0
Y(1)=Y(1)+K3(1)
Y(2)=Y(2)+K3(2)
CALL DYDR(YP,RAD,Y,IC,RHS,RHO,RHOB,E,N1,C,R,M1)
C Find vector K4.
K4(1)=H*YP(1)
K4(2)=H*YP(2)
C Accumulate sums.
SUM(1)=SUM(1)+(K1(1)+2.D0*K2(1)+2.D0*K3(1)+K4(1))/6.D0
SUM(2)=SUM(2)+(K1(2)+2.D0*K2(2)+2.D0*K3(2)+K4(2))/6.D0
11 CONTINUE
C Increment YY.
YY(1,IP1)=YY(1,I)+SUM(1)
YY(2,IP1)=YY(2,I)+SUM(2)

10 CONTINUE

RETURN
END

To illustrate the output from the programme FIGURE.FOR, we show the results
for Earth model Cal8. The radius, density, mean density, reciprocal of the flatten-
ing, 1/ f , and the departure of the internal equipotentials from ellipsoidal shape, κ,
are listed for each point of the Earth model. To ensure convergence, the number
of integration steps per Earth model subinterval is chosen as 100. The results are
showninTable5.3.
The second- and fourth-degree coe cients of the external gravitational potential,

J 2 and J 4 , are also computed. For Earth model Cal8 they are,
J 2 = 10 − 6 1071.3219 × , (5.211)

10 − 6 J 4 =−

×

2.9605 . (5.212)