SUBROUTINE NEWDMT
C     =================
C
C     New m-scale, calculated as that corresponding to the new
C     grid better representing temperature variations
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      DIMENSION DM0(MDEPTH),DM11(MDEPTH),DENS0(MDEPTH),ZD0(MDEPTH),
     *          T0(MDEPTH),T1(MDEPTH),ELEC0(MDEPTH),PT0(MDEPTH),
     *          ABRS0(MDEPTH),ABPL0(MDEPTH)
      COMMON/PRSAUX/VSND2(MDEPTH),HG1,HR1,RR1
      COMMON/FACTRS/GAMJ(MDEPTH),GAMH,FAK0
      COMMON/FLXAUX/T4,PGAS,PRAD,PGM,PRADM,ITGMAX,ITGMX0
C
      DO ID=1,ND
         DM0(ID)=LOG10(DM(ID))
         T0(ID)=LOG10(TAUROS(ID))
         ELEC0(ID)=ELEC(ID)
         DENS0(ID)=DENS(ID)
         PT0(ID)=PTOTAL(ID)
         ZD0(ID)=ZD(ID)
         ABRS0(ID)=ABROSD(ID)
         ABPL0(ID)=ABPLAD(ID)
      END DO
      ND1=ND-1
      CALL GRIDP(DM0,T0,DM11,T1,ND1)
      DM11(ND)=DM0(ND)
      T1(ND)=T0(ND)
      DO ID=1,ND
         DM(ID)=EXP(2.3025851*DM11(ID))
         TAUROS(ID)=EXP(2.3025851*T1(ID))
      END DO
      CALL INTERP(DM0,ELEC0,DM11,ELEC,ND,ND,2,0,1)
      CALL INTERP(DM0,DENS0,DM11,DENS,ND,ND,2,0,1)
      CALL INTERP(DM0,PT0,DM11,PTOTAL,ND,ND,2,0,1)
      CALL INTERP(DM0,ZD0,DM11,ZD,ND,ND,2,0,0)
      CALL INTERP(DM0,ABRS0,DM11,ABROSD,ND,ND,2,0,1)
      CALL INTERP(DM0,ABPL0,DM11,ABPLAD,ND,ND,2,0,1)
      DO ID=1,ND
         VSND2(ID)=PTOTAL(ID)/DENS(ID)
      END DO
C
C     New Rosseland opacity and functions theta and tauthe
C
      AMUV0=DMVISC**(ZETA0+UN)
      AMUV1=UN-AMUV0
      DO ID=1,ND
         IF(DM(ID).LE.DMVISC*DM(ND)) THEN
            VISCD(ID)=(UN-FRACTV)*(ZETA1+UN)/
     *                DMVISC**(ZETA1+UN)*(DM(ID)/DM(ND))**ZETA1
            THETA(ID)=(UN-FRACTV)*(DM(ID)/DMVISC/DM(ND))**(ZETA1+UN)
          ELSE
            VISCD(ID)=FRACTV*(ZETA0+UN)/AMUV1*
     *                (DM(ID)/DM(ND))**ZETA0
            THETA(ID)=(UN-FRACTV)+FRACTV*((DM(ID)/DM(ND))**(ZETA0+UN)-
     *                 AMUV0)/AMUV1
         END IF
         GAMJ(ID)=UN
         IF(ID.EQ.1) THEN
            TAUROS(ID)=DM(ID)*ABROSD(ID)
            TAUTHE(ID)=TAUROS(ID)*THETA(ID)/(ZETA1+TWO)
          ANEREL=ELEC(ID)/(DENS(ID)/WMM(ID)+ELEC(ID))
         ELSE
            DDM=DM(ID)-DM(ID-1)
            TAUROS(ID)=TAUROS(ID-1)+DDM*HALF*(ABROSD(ID-1)+ABROSD(ID))
            ZETAD=ZETA0
            IF(DM(ID).LE.DMVISC*DM(ND)) ZETAD=ZETA1
          A0=(ABROSD(ID-1)*DM(ID)-ABROSD(ID)*DM(ID-1))/DDM/
     *         (ZETAD+TWO)
          A1=(ABROSD(ID)-ABROSD(ID-1))/DDM/(ZETAD+3.D0)
            TAUTHE(ID)=TAUTHE(ID-1)+
     *              A0*(THETA(ID)*DM(ID)-THETA(ID-1)*DM(ID-1))+
     *              A1*(THETA(ID)*DM(ID)**2-THETA(ID-1)*DM(ID-1)**2)
         END IF
         TAUR=TAUROS(ID)
         CALL TEMPER(ID,TAUR,1)
      END DO
C
C     Next step - simultaneous solution of the hydrostatic
C     equilibrium and the z-m relation
C
      if(nconit.ge.0) CALL HESOLV
C
C     New temperature and mean opacities for the current density
C     and pressure
C
      DO ID=1,ND
         TAUR=TAUROS(ID)
         CALL TEMPER(ID,TAUR,1)
      END DO
C
C     Once again - simultaneous solution of the hydrostatic
C     equilibrium and the z-m relation
C
      if(nconit.ge.0) CALL HESOLV
C
      IF(IPRING.GE.1) THEN
         WRITE(6,601)
         DO ID=1,ND
            WRITE(6,602) ID,DM(ID),TAUROS(ID),
     *                   TEMP(ID),ELEC(ID),PTOTAL(ID),
     *                   ZD(ID),ABROSD(ID),ABPLAD(ID)
         END DO
      END IF
C
  601 FORMAT(1H1,' NEW DEPTH GRID ESTABLISHED, NEW MODEL:'/
     *           ' --------------------------------------'/
     * ' ID    DM     TAUROSS   TEMP     NE       P',
     * 8X,'ZD     ROSS.MEAN  PLANCK'/)
  602 FORMAT(1H ,I3,1P2D9.2,0PF8.0,1P3D9.2,2X,2D9.2)
C
      RETURN
      END