SUBROUTINE DMEVAL
C     =================
C
C     Auxiliary procedure for RESOLV  - for disks
C     recomputation of the m-scale in the case where z-scale is the
C     basic scale
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ITERAT.FOR'
      INCLUDE 'ARRAY1.FOR'
      dimension dma(mdepth),dmb(mdepth)
C
C     total pressure and gas pressure
C
      DO ID=1,ND
         PTURB=HALF*DENS(ID)*VTURB(ID)*VTURB(ID)
         PGS0=(DENS(ID)/WMM(ID)+ELEC(ID))*BOLK*TEMP(ID)
         PGS(ID)=PGS0
         PTOTL0=PGS(ID)+PRADT(ID)+PTURB
         PTOTAL(ID)=PTOTL0
      END DO
c
C     mass at the first depth point
C
      ID=1
      GRD=0.
      DO IJ=1,NFREQE
         IJT=IJFR(IJ)
         FLUXW=W(IJT)*FH(IJT)*RADEX(IJ,ID)
         GRD=GRD+FLUXW*ABSOEX(IJ,ID)
      END DO
      HG1=SQRT(TWO*PGS(1)/DENS(1)/QGRAV)
      HR1=PCK/QGRAV*(GRD+FPRD(1))/DENS(1)
      if(iter.eq.1) pgas0=pgs(1)
      X=(ZD(1)-HR1)/HG1
      IF(X.LT.3.) THEN
         IF(X.LT.0.) X=0.
         F1=8.86226925D-1*EXP(X*X)*ERFCX(X)
       ELSE
         F1=HALF*(UN-HALF/X/X)/X
      END IF
C
      DMa(1)=HG1*DENS(1)*F1
      DMb(1)=DMa(1)
c
c     recompute the DM scale
C
      write(6,600)
  600 format(/'   ID          ZD         DM(old)  DMA     DMB'/)
      DO ID=1,ND
         if(id.gt.1) then
         dmb(id)=dm(id)
         DMA(ID)=DMA(ID-1)-(ZD(ID)-ZD(ID-1))*TWO/
     *          (UN/DENS(ID)+UN/DENS(ID-1))
         DMb(ID)=DMb(ID-1)-(ZD(ID)-ZD(ID-1))*
     *          (DENS(ID)+DENS(ID-1))*half
         end if
         write(6,601) id,zd(id),dm(id),dma(id),dmb(id)
c        DM(ID)=DMB(ID)
         DM(ID)=DMa(ID)
  601 format(i5,1p4e12.4)
      END DO
      DMTOT=DM(ND)
      EDISC=SIG4P*TEFF**4/DMTOT
      RETURN
      END