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

SUBROUTINE RK4(R,Y,A,K1,N1,C,RM,RHOM,MUM,LAMBDAM,GZEROM,N,
1 PI,G,H,IR,RMIN,ANGS,COR,WES)

C
C This subroutine completes the last three steps of a fourth-order
C Runge-Kutta integration, given the derivatives K1(6,6) at
C the starting radius R.
C

IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION RM(100),RHOM(100),GZEROM(100),C(100,100),
1 Y(6,6),YY(6,6),A(6,6)
DOUBLE PRECISION MUM(100),LAMBDAM(100),MU,LAMBDA,
1 K1(6,6),K2(6,6),K3(6,6),K4(6,6)
C Set half stepsize, one sixth stepsize.
HH=0.5D0*H
H6=H/6.D0
C Set number of solutions for shell.
NS=6
C If region is inner core, reset number of solutions to three.
IF(IR.EQ.1)NS=3
C Complete remaining steps of Runge--Kutta integration.
DO 10 I=1,6
DO 11 J=1,NS
YY(I,J)=Y(I,J)+HH*K1(I,J)
11 CONTINUE
10 CONTINUE
C Increment radius and solve for new derivatives.
R=R+HH
CALL YPRIME(R,YY,A,K2,N1,C,RM,RHOM,MUM,LAMBDAM,GZEROM,N,PI,G,
1 IR,RMIN,ANGS,COR,WES)
DO 12 I=1,6
DO 13 J=1,NS
YY(I,J)=Y(I,J)+HH*K2(I,J)
13 CONTINUE
12 CONTINUE
C Find new derivatives.
CALL YPRIME(R,YY,A,K3,N1,C,RM,RHOM,MUM,LAMBDAM,GZEROM,N,PI,G,
1 IR,RMIN,ANGS,COR,WES)
DO 14 I=1,6
DO 15 J=1,NS
YY(I,J)=Y(I,J)+H*K3(I,J)
15 CONTINUE
14 CONTINUE
C Increment radius and solve for new derivatives.
R=R+HH
CALL YPRIME(R,YY,A,K4,N1,C,RM,RHOM,MUM,LAMBDAM,GZEROM,N,PI,G,
1 IR,RMIN,ANGS,COR,WES)
DO 16 I=1,6
DO 17 J=1,NS
Y(I,J)=Y(I,J)+H6*(K1(I,J)+2.D0*(K2(I,J)+K3(I,J))+K4(I,J))

17 CONTINUE

16 CONTINUE

RETURN
END

Calculation of the fundamental solutions requires specification of the Newto-
nian constant of gravitation, G , as well as the mean rate of rotation of the Earth, .

A recommended value for the former has been given by the Committee on Data for
Science and Technology (CODATA) for 2006 as G

Ω

10 − 11 m 3 kg − 1 s − 2
(Mohr et al ., 2008, pp. 688-691). The mean rate of rotation of the Earth is specified

= 6.67428 ×