Geology Reference
In-Depth Information
C Initialize function values at beginning of integration step.
DO 52 I=1,6
DO 53 J=1,6
Y1(I,J)=Y(I,J)
Y2(I,J)=Y(I,J)
53 CONTINUE
52 CONTINUE
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 54 I=1,6
DO 55 J=1,3
DELTA(I,J)=Y1(I,J)-Y2(I,J)
55 CONTINUE
54 CONTINUE
C Copy Y2(6,6) into inner core solution array YIC(6,3).
DO 56 I=1,6
DO 57 J=1,3
YIC(I,J)=Y2(I,J)
57 CONTINUE
56 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 58
C Step too large, reduce stepsize and repeat.
XV=XV-H
H=0.9D0*H*(ERRMAX**(-0.25D0))
GO TO 51
58 CONTINUE
C Step successful, accept solution and improve truncation error.
DO 59 I=1,6
DO 60 J=1,3
Y(I,J)=Y2(I,J)+DELTA(I,J)/15.D0
FS(I,J)=Y(I,J)
60 CONTINUE
59 CONTINUE
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)
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 51
Search WWH ::
Custom Search