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

C Read in and write out headers.
C Earth model name.
READ(1,12)(ENAME(I),I=1,10)
12 FORMAT(10A7)
WRITE(2,13)(ENAME(I),I=1,10)
WRITE(3,12)(ENAME(I),I=1,10)
WRITE(6,13)(ENAME(I),I=1,10)
13 FORMAT(11X,10A7)
READ(1,14)NN
14 FORMAT(I10)
C NN is not used.
C Number of model points (NI is not used).
READ(1,15)(NM(I),NI(I),I=1,4)
15 FORMAT(8I10)
WRITE(2,16)
WRITE(6,16)
16 FORMAT(//5X,'Number of model points for inner core, outer core, ',
1 'mantle and crust:'/)
WRITE(2,15)(NM(I),I=1,4)
WRITE(6,15)(NM(I),I=1,4)
C Enter number of steps in each region for integrations.
WRITE(6,17)
17 FORMAT(//1X,'Select number of steps per subinterval for',
1 ' integrations.')
READ(5,*)NIS
WRITE(2,18)NIS
18 FORMAT(/5X,'Number of steps per subinterval for integrations =',
1 I3)
C Initialize model point count.
K=0
C Read in Earth model.
DO 19 M=1,4
N1=NM(M)
READ(1,20)(R(I),RHO(I),LAMBDA(I),MU(I),GZERO(I),I=1,N1)
20 FORMAT(1X,F10.1,F10.2,F10.1,F10.1,F10.1)
DO 21 I=1,N1
C Scale Earth model to SI values and store.
J=K+I
RI(J)=R(I)*1.D3
RHOI(J)=RHO(I)*1.D3
LAMBDAI(J)=LAMBDA(I)*1.D8
MUI(J)=MU(I)*1.D8
GZEROI(J)=GZERO(I)*1.D-2
21 CONTINUE
C Record number of model points up to end of previous region.
NK(M)=K
K=K+N1
19 CONTINUE
C Begin quadratures for solution of the Clairaut equation by
C the Radau transformation.
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)