46 lines
1.0 KiB
Fortran
46 lines
1.0 KiB
Fortran
SUBROUTINE PSOLVE
|
|
C =================
|
|
C
|
|
C "formal" solution of the second-order equation for the total
|
|
C pressure - d**2 P/d m**2 = Q/DENS;
|
|
C with known density
|
|
C the resulting tridiagonal system is solved by the standard
|
|
C elimination
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
DIMENSION D(MDEPTH),ANU(MDEPTH)
|
|
C
|
|
C forward elimination
|
|
C
|
|
ID=1
|
|
B=1.D0
|
|
VL=PTOTAL(ID)
|
|
D(ID)=0.
|
|
ANU(ID)=VL/B
|
|
DO ID=2,ND-1
|
|
DMD=HALF*(DM(ID+1)-DM(ID-1))
|
|
A=UN/DMD/(DM(ID)-DM(ID-1))
|
|
C=UN/DMD/(DM(ID+1)-DM(ID))
|
|
B=A+C-A*D(ID-1)
|
|
VL=QGRAV/DENS(ID)
|
|
D(ID)=C/B
|
|
ANU(ID)=(VL+A*ANU(ID-1))/B
|
|
END DO
|
|
ID=ND
|
|
A=TWO/(DM(ID)-DM(ID-1))**2
|
|
B=A-A*D(ID-1)
|
|
VL=QGRAV/DENS(ID)
|
|
ANU(ID)=(VL+A*ANU(ID-1))/B
|
|
PTOTAL(ND)=ANU(ND)
|
|
C
|
|
C backsubstitution
|
|
C
|
|
DO IID=1,ND-1
|
|
ID=ND-IID
|
|
PTOTAL(ID)=ANU(ID)+D(ID)*PTOTAL(ID+1)
|
|
END DO
|
|
RETURN
|
|
END
|