Geology Reference
In-Depth Information
Y(5,3)=2.D0*PI*G*RHO(1)
Y(6,3)=0.D0
YSCAL(1,3)=-RHO(1)*(4.D0*GAMMA+ANGS+2.D0*WES)/
1 (10.D0*(LAMBDA(1)+2.D0*MU(1)))
YSCAL(2,3)=(5.D0*LAMBDA(1)+6.D0*MU(1))*YSCAL(1,3)
YSCAL(3,3)=0.D0
YSCAL(4,3)=0.D0
YSCAL(5,3)=PI*G*RHO(1)*YSCAL(1,3)
YSCAL(6,3)=0.D0
YSCAL2(1,3)=-RHO(1)*(4.D0*GAMMA+ANGS+2.D0*WES)*YSCAL(1,3)/
1 (4.D0*(7.D0*LAMBDA(1)+14.D0*MU(1)))
YSCAL2(2,3)=(7.D0*LAMBDA(1)+10.D0*MU(1))*YSCAL2(1,3)
YSCAL2(3,3)=0.D0
YSCAL2(4,3)=0.D0
YSCAL2(5,3)=2.D0*PI*G*RHO(1)*YSCAL2(1,3)/3.D0
YSCAL2(6,3)=0.D0
37 CONTINUE
C Set last three columns of array Y(6,6) to zero for inner core.
DO 40 I=1,6
DO 41 J=4,6
Y(I,J)=0.D0
41 CONTINUE
40 CONTINUE
C Put inner core solutions in array YIC(6,3).
DO 42 I=1,6
DO 43 J=1,3
YIC(I,J)=Y(I,J)
43 CONTINUE
42 CONTINUE
C Begin Runge-Kutta integration for the inner core.
C Find maximum ratio of second to first coefficient in power series
C expansion at the geocentre.
CALL REL(ERRMAX,YSCAL,YIC,N)
C Set initial stepsize so that fifth-order method would have relative
C error bound EPS.
H=(EPS/ERRMAX)**(0.2D0)
C Find solutions and derivatives at minimum radius RMIN by power series
C expansions.
XV=RMIN
XVS=XV*XV
XV3=XV*XVS
XV4=XVS*XVS
C Find solutions.
DO 44 I=1,6
DO 45 J=1,3
Y(I,J)=Y(I,J)+XVS*YSCAL(I,J)+XV4*YSCAL2(I,J)
45 CONTINUE
44 CONTINUE
C Set inner core index.
IR=1
46 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)
C Initialize function values at beginning of integration step.
DO 47 I=1,6
DO 48 J=1,6
Y1(I,J)=Y(I,J)
Y2(I,J)=Y(I,J)
48
CONTINUE
47
CONTINUE
Search WWH ::
Custom Search