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

DO 38 I=1,6
DO 39 J=1,3
RM1(I,J)=XVS*YSCAL(I,J)
RM2(I,J)=XV4*YSCAL2(I,J)
39 CONTINUE
38 CONTINUE
C Write out terms in power series expansions.
DO 40 J=1,3
WRITE(2,41)J
41 FORMAT(/,5X,'Fundamental solution number ',I3)
WRITE(2,42)
42 FORMAT(/,6X,'first term',10X,'second term',9X,'third term')
WRITE(2,43)(Y(I,J),RM1(I,J),RM2(I,J),I=1,6)
43 FORMAT(5X,D15.8,5X,D15.8,5X,D15.8)
40 CONTINUE
C Set last three columns of arrays Y(6,6) and YP(6,6) to zero
C for inner core.
DO 44 I=1,6
DO 45 J=4,6
Y(I,J)=0.D0
YP(I,J)=0.D0
45 CONTINUE
44 CONTINUE
C Put inner core solutions in array YIC(6,3)
C and free solutions array FS(6,3).
DO 46 I=1,6
DO 47 J=1,3
YIC(I,J)=Y(I,J)
FS(I,J)=Y(I,J)
47 CONTINUE
46 CONTINUE
C Write out starting values of free solutions.
RAD=XV*1.D-3
WRITE(3,48)RAD,(FS(I,1),I=1,6)
WRITE(4,48)RAD,(FS(I,2),I=1,6)
WRITE(7,48)RAD,(FS(I,3),I=1,6)
48 FORMAT(F7.2,6D14.6)
C Begin Runge-Kutta integration for the inner core.
C Find maximum ratio of second coefficient to first term
C in power series expansions.
CALL REL(ERRMAX,YSCAL,YIC,N)
C Set initial stepsize so that fourth-order method would have
C relative error bound EPS.
H=(EPS/ERRMAX)**(0.2D0)
C Find solutions and derivatives at minimum radius RMIN
C by power series expansions.
DO 49 I=1,6
DO 50 J=1,3
C Find solutions.
Y(I,J)=Y(I,J)+XVS*YSCAL(I,J)+XV4*YSCAL2(I,J)
C Find derivatives.
YP(I,J)=2.D0*XV*YSCAL(I,J)+4.D0*XV3*YSCAL2(I,J)
50 CONTINUE
49 CONTINUE
C Set inner core index.
IR=1
51 CONTINUE
C Calculate derivatives at beginning of current step.
CALL YPRIME(XV,Y,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,
1 IR,RMIN,ANGS,COR,WES)