SUBROUTINE PZEVLD
C     =================
C
C     Auxiliary procedure called from RESOLV
C     determination of the total and gas pressures, and logarithmic
C     gradient of pressure
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ARRAY1.FOR'
      COMMON/PRSAUX/VSND2(MDEPTH),HG1,HR1,RR1
      COMMON/DEPTDR/DDM(MDEPTH),DDP(MDEPTH),DD0(MDEPTH),
     *              DDMIN(MDEPTH),DDPLU(MDEPTH),DDA(MDEPTH),
     *              DDC(MDEPTH),DDB(MDEPTH)
      common/grdpra/GRD(MDEPTH),pra(mdepth),pgs0(mdepth),ANTP(MDEPTH)
      common/ifpzpa/ifpzev
      dimension dpp(mdepth),zd1(mdepth),zd2(mdepth),zd3(mdepth),
     *          ZOLD(MDEPTH)
C
      if(ifryb.gt.0.and.ifpzev.eq.0) return
      iheitr=0
    5 continue
      iheitr=iheitr+1
C
C     geometrical distance from the central plane - z
C
      if(ihecor.ge.0) then
         ZD(ND)=ZND
         DO IID=1,ND-1
            ID=ND-IID
            if(iheitr.eq.1) then
            ZD(ID)=ZD(ID+1)+HALF*(DM(ID+1)-DM(ID))*(UN/DENS(ID+1)+
     *             UN/DENS(ID))
            ZOLD(ID)=ZD(ID)
            end if
         END DO
       else
         zd1(1)=-ddc(1)/ddb(1)
         zd2(1)=-dens1(1)/ddb(1)
         do id=2,nd-1
            x=un/(ddb(id)-dda(id)*zd1(id-1))
            zd1(id)=-x*ddc(id)
            zd2(id)=-x*(dens1(id)-dda(id)*zd2(id-1))
         end do
         zd(nd)=znd
         do id=nd-1,1,-1
            zd(id)=zd1(id)*zd(id+1)+zd2(id)
         end do
      end if
C
C     total pressure, gas pressure, and sound speed
C
      DO ID=1,ND
         PTURB=HALF*DENS(ID)*VTURB(ID)*VTURB(ID)
         PGSC=(DENS(ID)/WMM(ID)+ELEC(ID))*BOLK*TEMP(ID)
         PGS(ID)=PGSC
         PTOTL0=PGS(ID)+PRADT(ID)+PTURB
         PTOTAL(ID)=PTOTL0
         VSND2(ID)=PTOTAL(ID)/DENS(ID)
      END DO
      ID=1
      HG1=SQRT(2.*PGS(1)/DENS(1)/QGRAV)
      HR1=PRD0/QGRAV
      RR1=HR1/HG1
C         
C     recalculate the z-distance 
C
      IJ1=1           
      if(icompt.gt.0.and.icombc.gt.0.and.ijex(1).gt.0) IJ1=2
      DO ID=1,ND
         GRP=0.
         FLEX=0.
         IF(NFREQE.GT.0.or.ifryb.eq.0) THEN
            DO IJ=IJ1,NFREQE
               RAD0(IJ)=RADEX(IJ,ID)
               FK0(IJ)=FAKEX(IJ,ID)
               ABSO0(IJ)=ABSOEX(IJ,ID)
               IJT=IJFR(IJ)
               WD0=W(IJT)
               FLUXW=FH(IJT)*RAD0(IJ)-HEXTRD(IJT)
               IF(.NOT.LSKIP(ID,IJT)) THEN
                  IF(ID.EQ.1) THEN
                     GRP=GRP+WD0*FLUXW*ABSO0(IJ)
                   ELSE
                     RADM(IJ)=RADEX(IJ,ID-1)
                     FKM(IJ)=FAKEX(IJ,ID-1)
                     ABSOM(IJ)=ABSOEX(IJ,ID-1)
                     FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ)
                     GRP=GRP+WD0*FRD
                  END IF
               END IF
            END DO
            GRD(ID)=GRP+FPRD(ID)
         END IF
         IF(ID.EQ.1) THEN
            GRV=QGRAV*ZD(ID)
          ELSE
            GRV=QGRAV*(ZD(ID)+ZD(ID-1))*HALF
            dpt=(ptotal(id)-ptotal(id-1))/ddm(id)
            dpr=(pradt(id)-pradt(id-1))/ddm(id)
            dpg=(pgs(id)-pgs(id-1))/ddm(id)
            dpr1=grd(id)/(dm(id)-dm(id-1))*4.19168946e-10
            err=10.
            if(grv.ne.0.) err=(dpg+dpr-grv)/grv
            dpp(id)=(dpg+dpr1)/qgrav
         end if
      end do
c
      if(iter.le.iabs(ifz0)) then 
         zd1(1)=dpp(2)
         zd1(nd)=znd
         do id=2,nd-1
            zd1(id)=(dpp(id)+dpp(id+1))*half
         end do
c
         zd2(nd)=znd
         zd3(nd)=znd
         izdiv=0
         nzdiv=0
         do id=nd-1,1,-1
            zd2(id)=2.*dpp(id+1)-zd2(id+1)
            zd3(id)=dpp(id)*ddmin(id)+dpp(id+1)*ddplu(id)
            if(zd2(id).le.zd2(id+1)) nzdiv=nzdiv+1
            if(nzdiv.eq.1) izdiv=id
         end do
c
         if(ihecor.ge.0) then
            do id=1,nd
               zd(id)=zd2(id)
               if(id.le.izdiv) zd(id)=zd1(id)
            end do
         else
            do id=1,nd
               zd(id)=zd3(id)
            end do
         end if
c
c        recalculate densities (if required)
c
         if(ihecor.gt.0) then
            do id=nd-1,1,-1
               x=HALF*(DM(ID+1)-DM(ID))
               xne=elec(id)/dens(id)
               dens(id)=x*dens(id+1)/((zd(id)-zd(id+1))*dens(id+1)-x)
               dens1(id)=un/dens(id)
               elec(id)=xne*dens(id)
            end do
          else if(ihecor.lt.-1) then
            do id=nd-1,1,-1
               xne=elec(id)/dens(id)
               if(id.gt.1) then
                  dens1(id)=zd(id-1)*dda(id)-zd(id)*ddb(id)-
     *                      zd(id+1)*ddc(id)
                else
                  dens1(id)=-zd(id)*ddb(id)-zd(id+1)*ddc(id)
               end if
               dens(id)=un/dens1(id)
               elec(id)=xne*dens(id)
            end do
         end if
c
         dzmx=0.
         do id=1,nd-1
            dzmx=max(dzmx,abs((zd(id)-zold(id))/zd(id)))
            zold(id)=zd(id)
         end do
         if(iheitr.ge.5.or.dzmx.lt.1.d-3) go to 25
         go to 5
   25    continue
c
      end if
      RETURN
      END