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