SUBROUTINE CONTMD
C     =================
C
C     Auxiliary procedure for LTEGRD
C     Determination of temperature in convectively unstable layers
C     for disks.
C     This is done by solving the energy balance equation
C     F(rad)+F(conv)=F(mech), which yields a cubic equation for
C     the logarithmic temperature gradient
C
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      COMMON ESEMAT(MLEVEL,MLEVEL),BESE(MLEVEL),
     *       DEPTH(MDEPTH),DEPTH0(MDEPTH),TAU(MDEPTH),TAU0(MDEPTH),
     *       TEMP0(MDEPTH),ELEC0(MDEPTH),DENS0(MDEPTH),DM0(MDEPTH)
      DIMENSION DELTR(MDEPTH),TEMPR(MDEPTH),ICON0(MDEPTH)
      COMMON/CUBCON/A,B,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
      COMMON/PRSAUX/VSND2(MDEPTH),HG1,HR1,RR1
      PARAMETER  (ERRT=1.D-3)
C
C     First, store the temperature(rad) and gradient Delta(rad) -
C     quantities for the  purely raditive equilibrium model
C
      T4=TEFF**4
      FLXTO0=SIG4P*T4
      DPRAD=1.891204931D-15*T4
      if(ifprad.eq.0) dprad=0.
      PRAD0=DPRAD/1.732D0
      DO ID=1,ND
         TEMPR(ID)=TEMP(ID)
         IF(ID.EQ.1) THEN
            DELTR(ID)=0.
          ELSE
            DELTR(ID)=
     *      (TEMP(ID)-TEMP(ID-1))/(PTOTAL(ID)-PTOTAL(ID-1))*
     *      (PTOTAL(ID)+PTOTAL(ID-1))/(TEMP(ID)+TEMP(ID-1))
         END IF
      END DO
      ICONIT=0
C
C     ------------------------------------------------------
C     Global iteration loop for calculating convective model
C     ------------------------------------------------------
C
   20 ICONIT=ICONIT+1
      ICONBE=0
      HR1=FLXTO0*PCK*ABROSD(1)/QGRAV
      CHANTM=0.
      DOID=1,ND
         T=TEMP(ID)
         PTOT=PTOTAL(ID)
         PGAS=PGS(ID)
         PTURB=HALF*DENS(ID)*VTURB(ID)*VTURB(ID)
         PRAD=PRADT(ID)
         FLXTOT=FLXTO0*(UN-THETA(ID))
         GRAVD=ZD(ID)*QGRAV
         ICON0(ID)=0
C
         IF(ID.EQ.1) GO TO 40
         J=0
         IF(ICONIT.EQ.1) T=T-TEMPR(ID-1)+TEMP(ID-1)
         TM=TEMP(ID-1)
         IF(T.LT.0.) T=TM
         PGM=PGS(ID-1)
         PTOTM=PTOTAL(ID-1)
         PT0=HALF*(PTOT+PTOTM)
         DELR=DELTR(ID)
C
C        Inner iteration loop for determining temperature in the
C        conectively unstable layers
C
   30    J=J+1
         TOLD=T
         T0=HALF*(T+TM)
         PG0=HALF*(PGAS+PGM)
         PR0=HALF*(PRAD+PRADM)
         AB0=HALF*(ABROSD(ID)+ABROSD(ID-1))
         IF(ID.GE.ND-2.AND.ICONBE.EQ.0) GO TO 40
         CALL CONVEC(ID,T0,PT0,PG0,PR0,AB0,DELR,FLXCNV,VCON)
         IF(FLXCNV.EQ.0.) GO TO 40
         ICON0(ID)=1
         ICONBE=1
         if(id.eq.nd) then
            pip=(ptot+ptotm)/(ptot-ptotm)
            t=tm*(pip+delr)/(pip-delr)
            go to 40
         end if
         CALL CUBIC(DELTA0)
         FAC=DELTA0*(PTOT-PTOTM)/(PTOT+PTOTM)
         T=TM*(UN+FAC)/(UN-FAC)
         IF(T.LT.TM) T=TM
         IF(ABS(UN-T/TOLD).GT.ERRT.AND.J.LT.10) GO TO 30
C
C     Store the final quantitites
C
   40    IF(ID.GT.1.AND.ICON0(ID).EQ.0.AND.ICON0(ID-1).EQ.1)
     *      DELTC=DELT0
         if(id.eq.nd) then
            pip=(ptot+ptotm)/(ptot-ptotm)
            t=tm*(pip+delr)/(pip-delr)
         end if
         DELT0=TEMP(ID)-T
         PRADT(ID)=PRADT(ID)*(T/TEMP(ID))**4
         DENS(ID)=DENS(ID)*(TEMP(ID)/T)
         IF(TEMP(ID).NE.0.) CHANT0=ABS((T-TEMP(ID))/TEMP(ID))
         IF(CHANT0.GT.CHANTM) CHANTM=CHANT0
         TEMP(ID)=T
         IF(ICONIT.GT.1.AND.ICON0(ID).EQ.0.AND.ICONBE.EQ.1)
     *      TEMP(ID)=T-DELTC
         PRADM=PRADT(ID)
      END DO   
C
C     Diagnostic outprint
C
      IF(IPRING.EQ.2) THEN
         WRITE(6,600) ICONIT
         CALL CONOUT(1,IPRING)
      END IF
  600 FORMAT(1H1,' CONVECTIVE FLUX: AT CONTMD, ITER=',I2/)
C
C     2. New values of electron density and density
C
      CALL HESOL6
C
C     Evaluation of the Rosseland and Planck mean opacities
C
      DO ID=1,ND
         T=TEMP(ID)
         CALL WNSTOR(ID)
         CALL STEQEQ(ID,POP,1)
         CALL OPACF0(ID,NFREQ)
         CALL MEANOP(T,ABSO,SCAT,OPROS,OPPLA)
         ABROS=OPROS/DENS(ID)
         ABPLA=OPPLA/DENS(ID)
         ABROSD(ID)=ABROS
         ABPLAD(ID)=ABPLA
      END DO
      IF(CHANTM.GT.ERRT.AND.ICONIT.LT.NCONIT) GO TO 20
      RETURN
      END