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

IF(I.GT.3)J=J+1
CM1(1,I)=Y(2,J)
CM1(2,I)=Y(4,J)
CM1(3,I)=Y(5,J)+D*Y(6,J)/(AN+1.D0)
CM1(4,I)=0.D0
CM1(5,I)=0.D0
79 CONTINUE
CM1(4,1)=GZEROM
CM1(4,2)=-1.D0/RHOBM
CM1(4,4)=-1.D0
CM1(5,1)=4.D0*PI*G*RHOBM
CM1(5,2)=0.D0
CM1(5,5)=1.D0
C Set coefficient matrices equal.
DO 80 I=1,5
DO 81 J=1,5
CM2(I,J)=CM1(I,J)
CM3(I,J)=CM2(I,J)
81 CONTINUE
80 CONTINUE
C Set values of constant vectors.
BV1(1)=0.D0
BV1(2)=0.D0
BV1(3)=0.D0
BV1(4)=1.D0
BV1(5)=0.D0
DO 82 I=1,4
BV2(I)=0.D0
82 CONTINUE
BV2(5)=-1.D0
C Store first constant vector.
DO 83 I=1,5
BV3(I)=BV1(I)
83 CONTINUE
C Solve for linear combination coefficients.
CALL LINSOL(CM1,BV1,5,CAUGS,DET,5,11)
CALL LINSOL(CM2,BV2,5,CAUGS,DET,5,11)
C Alter coefficient matrix to apply subseismic condition.
CM3(5,1)=GZEROM
CM3(5,2)=-1.D0/RHOBM
CM3(5,5)=0.D0
C Solve for subseismic coefficients.
CALL LINSOL(CM3,BV3,5,CAUGS,DET,5,11)
C Calculate potential condition parameter BN.
BN=(AN+1.D0)/(4.D0*PI*G*RHOBM*BV3(1)+BV3(5))
GO TO 84
C Construct coefficient matrix for N=0.
C Reduce Y-matrix to 4x4.
78

CONTINUE
DO 85 I=1,4
K=I
IF(K.GT.2) K=K+2
DO 86 J=1,4
L=J
IF(L.GT.2) L=L+2
Y(I,J)=Y(K,L)

86 CONTINUE

85 CONTINUE

DO 87 I=1,4
CM01(1,I)=Y(2,I)
CM01(2,I)=Y(3,I)+D*Y(4,I)/(AN+1.D0)