Geology Reference
In-Depth Information
WRITE(6,69)
69 FORMAT(11X,'(km)',9X,'(gm/cc)',8X,'(gm/cc)',4X/)
C Scale and write out values.
DO 70 I=1,NT
RAD=RI(I)/1000.D0
DEN=RHOI(I)/1000.D0
DENB=RHOBI(I)/1000.D0
RECIPF=1/FI(I)
KAPPA=YYI(1,I)
WRITE(2,71)RAD,DEN,DENB,RECIPF,KAPPA
WRITE(6,71)RAD,DEN,DENB,RECIPF,KAPPA
71 FORMAT(1X,F15.1,3F15.3,0PD15.3)
70 CONTINUE
C Calculate values of second- and fourth-degree gravitational coefficients,
C J2, J4, and write out values.
J2=-FRD/3.D0+2.D0*FI(NT)/3.D0+2.D0*FRD*FI(NT)/21.D0
1 -FI(NT)*FI(NT)/3.D0+8.D0*YYI(1,NT)/21.D0
J4=4.D0*FRD*FI(NT)/7.D0-4.D0*FI(NT)*FI(NT)/5.D0
1 -32.D0*YYI(1,NT)/35.D0
WRITE(6,72)J2
WRITE(2,72)J2
72
FORMAT(//5X,'J2=',4PD14.7/)
WRITE(6,73)J4
WRITE(2,73)J4
73
FORMAT(//5X,'J4=',1PD11.4/)
END
C
SUBROUTINE QUAD(INT,EX,Q,NIS,W,X,N1,C,R,M,M1)
C
C This subroutine performs the quadrature integral of a function, Q, of
C mean radius, from the centre of the Earth through the inner core for
C M=1, or the outer core for M=2, or the mantle for M=3, or the crust
C for M=4. The result is returned in INT. The integration uses six-point
C Gaussian integration which is eleventh-order accurate (exact for
C polynomials up to and including degree eleven). Interpolation of Q to
C the Gaussian integration points is performed by the spline
C interpolator subroutine INTPL. The number of integration steps between
C Earth model points is specified by NIS.
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Q(M1),R(M1),C(M1,M1),W(6),X(6),RI(6),QI(6)
DOUBLE PRECISION INT(M1,4)
N1M1=N1-1
DO 10 I=1,N1M1
IP1=I+1
C Find half-increment in radius.
H=(R(IP1)-R(I))/(2.D0*DFLOAT(NIS))
C Set sum to zero.
SUM=0.D0
C Integrate over subinterval.
DO 11 J=1,NIS
JM1=J-1
C Find mean radius.
RM=R(I)+2.D0*DFLOAT(JM1)*H+H
C Find radii for Gaussian integration.
DO 12 K=1,6
RI(K)=H*X(K)+RM
C Interpolate integrand at Gaussian abscissae.
CALL INTPL(RI(K),QI(K),N1,C,R,Q,M1)
C Accumulate sum.
SUM=SUM+W(K)*QI(K)*((RI(K))**EX)
Search WWH ::
Custom Search