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

X(5)=-X(2)
X(6)=-X(1)
C Initialize starting values of the five quadratures at the geocentre.
INTI(1,1)=0.D0
INTI(2,1)=0.D0
INTI(3,1)=0.D0
INTI(4,1)=0.D0
INTI(5,1)=0.D0
DO 22 M=1,4
MP1=M+1
C Set up interpolation for region M.
N1=NM(M)
N2=2*N1-2
N3=N1-2
C Put radius and density for region M in active locations.
DO 23 I=1,N1
J=NK(M)+I
R(I)=RI(J)
RHO(I)=RHOI(J)
23 CONTINUE
C Construct interpolation matrix.
CALL SPMAT(N1,N2,N3,C,R,B,M1,M2,M3)
C Store interpolation matrices.
DO 24 I=1,N1
DO 25 J=1,N1
D(I,J,M)=C(I,J)
25 CONTINUE
24 CONTINUE
C Begin quadrature for mean density.
C Set exponent in integrand.
EX=2.D0
C Initialize starting value of integral INT.
INT(1,M)=INTI(1,M)
C Form integrand function Q.
DO 26 I=1,N1
Q(I)=RHO(I)
26 CONTINUE
C Perform quadrature.
CALL QUAD(INT,EX,Q,NIS,W,X,N1,C,R,M,M1)
C Store starting value for next region.
INTI(1,MP1)=INT(N1,M)
C Form mean density and store in RHOBI.
DO 27 I=1,N1
IF(I.EQ.1.AND.M.EQ.1)GO TO 28
RHOB(I)=3.D0*INT(I,M)/(R(I)**3)
GO TO 29
28 RHOB(I)=RHO(I)
C Store mean density.
29 J=NK(M)+I
RHOBI(J)=RHOB(I)
27 CONTINUE
C Perform first estimate of Radau variable ETA.
C Set exponent in integrand.
EX=4.D0
C Initialize starting value of integral INT.
INT(1,M)=INTI(2,M)
C Form integrand function Q.
DO 30 I=1,N1
Q(I)=RHOB(I)
30 CONTINUE
C Perform quadrature.