Geology Reference
In-Depth Information
C Integrate solution forward one full step.
CALL RK4(XV,Y1,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,H,IR,
1 RMIN,ANGS,COR,WES)
C Integrate solution forward two half-steps.
C Halve step.
HH=0.5D0*H
C Reset radius.
XV=XV-H
CALL RK4(XV,Y2,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,HH,IR,
1 RMIN,ANGS,COR,WES)
CALL YPRIME(XV,Y2,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,IR,
1 RMIN,ANGS,COR,WES)
CALL RK4(XV,Y2,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,HH,IR,
1 RMIN,ANGS,COR,WES)
C Find error matrix and maximum relative error.
DO 49 I=1,6
DO 50 J=1,3
DELTA(I,J)=Y1(I,J)-Y2(I,J)
50 CONTINUE
49 CONTINUE
C Copy Y2(6,6) into inner core solution array YIC(6,3).
DO 51 I=1,6
DO 52 J=1,3
YIC(I,J)=Y2(I,J)
52 CONTINUE
51 CONTINUE
C Find maximum relative error.
CALL REL(ERRMAX,DELTA,YIC,N)
C Scale maximum relative error by specified relative error tolerance.
ERRMAX=ERRMAX/EPS
C Test if step was successful.
IF(ERRMAX.LT.1.D0) GO TO 53
C Step too large, reduce stepsize and repeat.
XV=XV-H
H=0.9D0*H*(ERRMAX**(-0.25D0))
GO TO 46
53 CONTINUE
C Step successful, accept solution and improve truncation error.
DO 54 J=1,6
DO 55 K=1,3
Y(J,K)=Y2(J,K)+DELTA(J,K)/15.D0
55 CONTINUE
54 CONTINUE
C Increase stepsize.
H=0.9D0*H*(ERRMAX**(-0.20D0))
C Test if solution is finished.
XVT=XV+H
IF(XVT.LT.R(N1)) GO TO 46
C Find remaining distance to the ICB.
H=R(N1)-XV
C Compute fraction of radius.
RATIO=H/R(N1)
C Continue integration only if remaining distance is more than
C 1.D-9 of radius.
IF(RATIO.GT.1.D-9) GO TO 46
C Convert solutions back to y-variables.
C Set value of alpha for first two solutions.
ANM2=AN-2.D0
C Convert three solutions.
DO 56 J=1,3
C Set value of alpha for third solution.
Search WWH ::
Custom Search