The subseismic equation and boundary conditions - Earth Dynamics: Deformations and Oscillations of the Rotating Earth - page 432

15 FORMAT(I10)
WRITE(2,15)NN
C Number of model points and number of integration steps.
READ(1,16)(NM(I),NI(I),I=1,4)
16 FORMAT(8I10)
WRITE(6,17)
17 FORMAT(//16X,'Number of model points and number of integration',
1 ' steps.'/)
WRITE(6,16)(NM(I),NI(I),I=1,4)
WRITE(6,18)
18 FORMAT(//5X,'Radius',6X,'Rho',4X,'Lambda',7X,'Mu',6X,'Gzero')
WRITE(6,19)
19 FORMAT(6X,'(km)',5X,'(gm/cc)',2X,'(kbars)',4X,'(kbars)',2X,
1 '(cm/sec/sec)'/)
C Read in Earth model, inner core, outer core, mantle and crust.
K=0
DO 20 M=1,4
N1=NM(M)
READ(1,21)(R(I),RHO(I),LAMBDA(I),MU(I),GZERO(I),I=1,N1)
21 FORMAT(1X,F10.1,F10.2,F10.1,F10.1,F10.1)
DO 22 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
C Write Earth model out to screen.
WRITE(6,21)R(I),RHO(I),LAMBDA(I),MU(I),GZERO(I)
22 CONTINUE
NK(M)=K
K=K+N1
20 CONTINUE
C Find density above, lambda above, gravity at the inner core boundary.
RHOAC=RHOI(NK(2)+1)
GZEROC=GZEROI(NK(2)+1)
C Find density below, lambda below, gravity at the core-mantle boundary.
RHOBM=RHOI(NK(3))
GZEROM=GZEROI(NK(3))
C Store radii of inner core and core-mantle boundaries
C and surface radius.
RIOB=RI(NK(2))
RCMB=RI(NK(3))
D=RI(NK(4)+NM(4))
C Write out radii of ICB and CMB.
WRITE(2,23) RIOB,RCMB
C Write out gravity at the ICB and the CMB.
WRITE(2,23) GZEROC,GZEROM
23 FORMAT(2D23.15)
C Begin calculation of Love numbers for the inner core.
C Set up interpolation for the inner core.
N1=NM(1)
N2=2*N1-2
N3=N1-2
C Put inner core values in active locations.
DO 24 I=1,N1
J=NK(1)+I
R(I)=RI(J)
RHO(I)=RHOI(J)
MU(I)=MUI(J)