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

The quadratures required for the calculation of the mean density, and the solu-
tion of the Clairaut equation, are performed by the subroutine QUAD. The Earth
model is interpolated by cubic splines with calls to the subroutines SPMAT and
INTPL, described in Section 1.6. The solution of the second-order equation for κ
(5.192) requires both the complementary function and the particular integral. These
are found by the fourth-order Runge-Kutta scheme, described in Section 3.6, by
the subroutine RK. The subroutine DYDR provides the derivatives required by the
Runge-Kutta integration scheme. The starting values of the particular integral and
its derivative at the geocentre require the derivative of the density there. This is
found by cubic spline interpolation by a call to the subroutine DERIV, described in
Section 1.6. The number of steps between Earth model points, for both the quad-
ratures and the Runge-Kutta integrations, is specified by the integer variable NIS,
input at the request of the programme. As well as output to the screen, the file
figure.dat is generated by the programme.

C PROGRAMME FIGURE.FOR

C
C FIGURE.FOR calculates the mean density, the 'flattening', f, and
C departure from ellipticity, kappa, of the internal level surfaces for
C an input Earth model, under the hydrostatic assumption. Related
C second- and fourth-degree coefficients of the gravitational potential,
C J2 and J4, are also calculated.
C

IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Q(100),R(100),RHO(100),GZERO(100),RI(300),RHOI(300),
1 GZEROI(300),ENAME(10),NM(4),NI(4),NK(4),RHOB(100),
2 B(98,198), C(100,100),D(100,100,4),RHOBI(300),ETA(100),ETAI(300),
3 E(100), EI(300),W(6),X(6),YYC(2,100),YYCI(2,300),YYP(2,100),
4 YYPI(2,300),FI(300),YYI(2,300),RHS(100)
DOUBLE PRECISION MU(100),LAMBDA(100),MUI(300),LAMBDAI(300),
1 INT(100,4),INTI(5,5),INITYC(2,5),INITYP(2,5),KAPPA,J2,J4
CHARACTER*20 EMODEL
C Enter Earth model file name.
WRITE(6,10)
10 FORMAT(1X,'Type in Earth model file name.')
READ(5,11)EMODEL
11 FORMAT(A20)
C Open Earth model file.
OPEN(UNIT=1,FILE=EMODEL,STATUS='OLD')
C Open figure file.
OPEN(UNIT=2,FILE='figure.dat',STATUS='UNKNOWN')
C Open inner core data file.
OPEN(UNIT=3,FILE='innerc.dat',STATUS='UNKNOWN')
C Set value of pi.
PI=3.141592653589793D0
C Set value of universal constant of gravitation (CODATA 2006).
G=6.67428D-11
C Set angular frequency of Earth's rotation (WGS84).
WE=7.292115D-5
C Set maximum dimensions for interpolation.
M1=100
M2=198
M3=98
C Read in and write out Earth model.