Geology Reference
In-Depth Information
LAMBDA(I)=LAMBDAI(J)
GZERO(I)=GZEROI(J)
22 CONTINUE
C Construct interpolation matrix.
CALL SPMAT(N1,N2,N3,C,R,B,M1,M2,M3)
C Begin integration by Simpson's rule.
WRITE(6,23)
23 FORMAT(//9X,'radius',9X,'exponent',7X,'f factor')
EINT=0.D0
F(N1)=1.D0
RKM=R(N1)*1.D-3
WRITE(2,14)N1
WRITE(2,24)R(N1),F(N1),RHO(N1)
24 FORMAT(3D23.15)
WRITE(6,25)RKM,EINT,F(N1)
25 FORMAT(/1X,F15.1,2F15.3)
N1M1=N1-1
DO 26 I=1,N1M1
C Interpolate rho, g, lambda to midpoint of interval.
J=N1-I
JP1=J+1
C Find stepsize for integration.
H=0.5D0*(R(JP1)-R(J))
C Find midpoint of interval.
RM=R(J)+H
CALL INTPL(RM,RHOM,N1,C,R,RHO,M1)
CALL INTPL(RM,GMEAN,N1,C,R,GZERO,M1)
CALL INTPL(RM,LAMBDAM,N1,C,R,LAMBDA,M1)
C Apply Simpson's rule.
F1=RHO(J)*GZERO(J)/LAMBDA(J)
F2=RHOM*GMEAN/LAMBDAM
F3=RHO(JP1)*GZERO(JP1)/LAMBDA(JP1)
EINT=EINT+H*(F1+4.D0*F2+F3)/3.D0
F(J)=DEXP(EINT)
RKM=R(J)*1.D-3
C Write out results.
WRITE(2,24)R(J),F(J),RHO(J)
WRITE(6,25)RKM,EINT,F(J)
26 CONTINUE
C Interpolate decompression factor onto integration points.
C Find number of integration steps and interval size.
M=NI(2)
AM=DFLOAT(M)
H=(R(N1)-R(1))/AM
RKM=R(1)*1.D-3
EINT=DLOG(F(1))
FI=F(1)
WRITE(2,25)RKM,EINT,FI
DO 27 I=1,M
AI=DFLOAT(I)
XV=R(1)+AI*H
CALL INTPL(XV,FI,N1,C,R,F,M1)
EINT=DLOG(FI)
RKM=XV*1.D-3
WRITE(2,25)RKM,EINT,FI
27
CONTINUE
END
The programme DECOMP can be used to calculate the decompression factor for
the Earth models Cal8, 1066A, PREM and Core11, listed in Appendix C, simply
Search WWH ::
Custom Search