Geology Reference
In-Depth Information
The matrix
A
is derived from
A
by the element modifications
1
A
12
=
r
A
12
,
A
21
=
A
23
=
rA
23
,
A
26
=
rA
21
,
rA
26
,
1
A
34
=
r
A
34
,
A
41
=
A
43
=
rA
43
,
A
45
=
r
2
A
45
,
rA
41
,
(3.308)
1
1
A
51
=
r
A
51
,
A
56
=
r
A
56
,
A
65
=
rA
65
.
The numerical propagation of the fundamental solutions by the Runge-Kutta
scheme (3.281) requires the repeated calculation of their derivatives as the product
of the matrix
A
(
r
) with the propagator matrix
Y
(
r
,ρ), as in (3.299). The sub-
routine YPRIME performs this operation, after a call elsewhere to the subroutine
SPMAT to enable cubic spline interpolation of a specified Earth model through
calls to the subroutine INTPL. The subroutines SPMAT and INTPL are described
in Section 1.6.
SUBROUTINE YPRIME(R,Y,A,YP,N1,C,RM,RHOM,MUM,LAMBDAM,GZEROM,N,
1 PI,G,IR,RMIN,ANGS,COR,WES)
C
C This subroutine finds the derivatives, YP(6,6), of the six
C fundamental solutions represented by the propagator matrix, Y(6,6),
C in the mantle and crust (shell), at an arbitrary radius R
C for degree N>0. For N=0, the third and fourth solution vectors
C can be ignored. In the inner core, it finds the derivatives of
C the fundamental solutions in terms of the z-variables, which are
C transformations of the y-variables.
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION RM(100),RHOM(100),GZEROM(100),C(100,100),
1 Y(6,6),YP(6,6),A(6,6)
DOUBLE PRECISION MUM(100),LAMBDAM(100),MU,LAMBDA
C Zero fill derivative matrix YP(I,J) and coefficient matrix A(I,J).
DO 10 I=1,6
DO 11 J=1,6
YP(I,J)=0.D0
A(I,J)=0.D0
11 CONTINUE
10 CONTINUE
C Set maximum number of points for interpolation.
M1=100
C If R is less than or equal to RMIN, take properties to be those
C at the geocentre.
IF(R.LE.RMIN) GO TO 12
C Interpolate Earth properties at radius R.
CALL INTPL(R,RHO,N1,C,RM,RHOM,M1)
CALL INTPL(R,MU,N1,C,RM,MUM,M1)
CALL INTPL(R,LAMBDA,N1,C,RM,LAMBDAM,M1)
CALL INTPL(R,GZERO,N1,C,RM,GZEROM,M1)
GO TO 13
C Assign geocentric Earth properties.
12
RHO=RHOM(1)
MU=MUM(1)
LAMBDA=LAMBDAM(1)
GZERO=(4.D0*PI*G*RHO-2.D0*WES)*R/3.D0
Search WWH ::
Custom Search