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

YA(2)=RHO(1)*(-(4.D0*GAMMA+ANGS+2.D0*WES)*YSCAL(1,2)
1 +(AN*(AN+1.D0)*GAMMA+COR)*YSCAL(3,2)-YSCAL(6,2))
YA(3)=0.D0
YA(4)=RHO(1)*((GAMMA+COR/(AN*(AN+1.D0)))*YSCAL(1,2)
1 -(ANGS-COR/(AN*(AN+1.D0)))*YSCAL(3,2)-YSCAL(5,2))
YA(5)=0.D0
YA(6)=0.D0
C Copy stored matrix AS into matrix A.
DO 33 I=1,6
DO 34 J=1,6
A(I,J)=AS(I,J)
34 CONTINUE
33 CONTINUE
C Solve for third terms in second solution.
CALL LINSOL(A,YA,6,CAUG,DET,6,13)
YSCAL2(1,2)=YA(1)
YSCAL2(2,2)=YA(2)
YSCAL2(3,2)=YA(3)
YSCAL2(4,2)=YA(4)
YSCAL2(5,2)=YA(5)
YSCAL2(6,2)=YA(6)
C Set right-hand side for calculation of third terms in third solution.
YA(1)=0.D0
YA(2)=RHO(1)*(-(4.D0*GAMMA+ANGS+2.D0*WES)*YSCAL(1,3)
1 +(AN*(AN+1.D0)*GAMMA+COR)*YSCAL(3,3)-YSCAL(6,3))
YA(3)=0.D0
YA(4)=RHO(1)*((GAMMA+COR/(AN*(AN+1.D0)))*YSCAL(1,3)
1 -(ANGS-COR/(AN*(AN+1.D0)))*YSCAL(3,3)-YSCAL(5,3))
YA(5)=0.D0
YA(6)=0.D0
C Copy stored matrix AS into matrix A and add 2.0 to diagonal of A.
DO 35 I=1,6
DO 36 J=1,6
A(I,J)=AS(I,J)
IF(I.EQ.J)A(I,J)=A(I,J)+2.D0
36 CONTINUE
35 CONTINUE
C Solve for third terms in third solution.
CALL LINSOL(A,YA,6,CAUG,DET,6,13)
YSCAL2(1,3)=YA(1)
YSCAL2(2,3)=YA(2)
YSCAL2(3,3)=YA(3)
YSCAL2(4,3)=YA(4)
YSCAL2(5,3)=YA(5)
YSCAL2(6,3)=YA(6)
GO TO 37
27 CONTINUE
C Begin power series expansion of fundamental solution for degree N=0.
C Set first two solutions to zero.
DO 38 I=1,6
DO 39 J=1,2
Y(I,J)=0.D0
YSCAL(I,J)=0.D0
YSCAL2(I,J)=0.D0
39 CONTINUE
38 CONTINUE
C Set values for first three terms of fundamental solution.
Y(1,3)=1.D0
Y(2,3)=3.D0*LAMBDA(1)+2.D0*MU(1)
Y(3,3)=0.D0
Y(4,3)=0.D0