Earth deformations - Earth Dynamics: Deformations and Oscillations of the Rotating Earth - page 259

For n ≥ 1, completion of the power series expansions requires the solution of the

6 × 6 linear systems (3.309), (3.311), (3.312), (3.314) and (3.315), with coe cient

matrix A , or its modification A . The subroutine MATRIX calculates the coe cient

matrix A , allowing the solution of these linear systems with calls elsewhere to the
subroutine LINSOL, described in Section 1.5.

SUBROUTINE MATRIX(A,RHO,LAMBDA,MU,PI,G,AN,WE)

C
C This subroutine calculates the 6X6 matrix A required to find
C higher terms in the power series expansions of the three
C fundamental solutions at the geocentre.
C

IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION A(6,6),RHO(100)
DOUBLE PRECISION LAMBDA(100),MU(100)
BETA=1.D0/(LAMBDA(1)+2.D0*MU(1))
DELTA=2.D0*MU(1)*(3.D0*LAMBDA(1)+2.D0*MU(1))*BETA
GAMMA=(4.D0*PI*G*RHO(1)-2.D0*WE*WE)/3.D0
AN1=AN*(AN+1.D0)
EPSLN=4.D0*AN1*MU(1)*(LAMBDA(1)+MU(1))*BETA-2.D0*MU(1)
A(1,1)=AN+3.D0+2.D0*LAMBDA(1)*BETA
A(1,2)=-BETA
A(1,3)=-AN1*LAMBDA(1)*BETA
A(1,4)=0.D0
A(1,5)=0.D0
A(1,6)=0.D0
A(2,1)=-2.D0*DELTA
A(2,2)=AN+2.D0+4.D0*MU(1)*BETA
A(2,3)=AN1*DELTA
A(2,4)=-AN1
A(2,5)=0.D0
A(2,6)=0.D0
A(3,1)=1.D0
A(3,2)=0.D0
A(3,3)=AN+2.D0
A(3,4)=-1.D0/MU(1)
A(3,5)=0.D0
A(3,6)=0.D0
A(4,1)=DELTA
A(4,2)=LAMBDA(1)*BETA
A(4,3)=-EPSLN
A(4,4)=AN+5.D0
A(4,5)=0.D0
A(4,6)=0.D0
A(5,1)=-3.D0*GAMMA-2.D0*WE*WE
A(5,2)=0.D0
A(5,3)=0.D0
A(5,4)=0.D0
A(5,5)=AN+4.D0
A(5,6)=-1.D0
A(6,1)=0.D0
A(6,2)=0.D0
A(6,3)=AN1*(3.D0*GAMMA+2.D0*WE*WE)
A(6,4)=0.D0
A(6,5)=-AN1
A(6,6)=AN+5.D0
RETURN
END