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

72 CONTINUE
C Prepare for spline interpolation.
CALL SPMAT(N1,N2,N3,C,R,B,M1,M2,M3)
C Set initial values for fundamental solutions.
DO 73 I=1,6
Y(I,I)=1.D0
DO 74 J=1,6
IF(I.NE.J) Y(I,J)=0.D0
74 CONTINUE
73 CONTINUE
C Begin Runge-Kutta integration for the mantle.
M=NI(3)
AM=DFLOAT(M)
C Find stepsize.
H=(R(N1)-R(1))/AM
C Set initial radius.
XV=R(1)
C Set integration region counter.
IR=3
DO 75 I=1,M
CALL YPRIME(XV,Y,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,
1 IR,RMIN,ANGS,COR,WES)
CALL RK4(XV,Y,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,H,IR,
1 RMIN,ANGS,COR,WES)
75 CONTINUE
C Set up interpolation for the crust.
N1=NM(4)
N2=2*N1-2
N3=N1-2
C Put crustal values in active locations.
DO 76 I=1,N1
J=NK(4)+I
R(I)=RI(J)
RHO(I)=RHOI(J)
MU(I)=MUI(J)
LAMBDA(I)=LAMBDAI(J)
GZERO(I)=GZEROI(J)
76 CONTINUE
C Prepare for spline interpolation.
CALL SPMAT(N1,N2,N3,C,R,B,M1,M2,M3)
C Begin Runge-Kutta integration for the crust.
M=NI(4)
AM=DFLOAT(M)
C Find stepsize.
H=(R(N1)-R(1))/AM
C Set initial radius.
XV=R(1)
C Set integration region counter.
IR=4
DO 77 I=1,M
CALL YPRIME(XV,Y,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,
1 IR,RMIN,ANGS,COR,WES)
CALL RK4(XV,Y,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,H,IR,
1 RMIN,ANGS,COR,WES)
77 CONTINUE
C Begin construction of coefficient matrix.
AN=DFLOAT(N)
C Test if N=0.
IF(N.EQ.0) GO TO 78
DO 79 I=1,5
J=I