SUBROUTINE CONVC1(ID,T,PTOT,PG,PRAD,ABROS,DELTA,FLXCNV,FC0)
C     ===========================================================
C
C     Determination of the mixing-lengths convective flux
C
C     Input:  T     - temperature
C             PTOT  - total pressure
C             PG    - gas pressure
C             PRAD  - radiation pressure
C             ABROS - Rosseland opacity (per gram)
C             DELTA - corresponding temperature gradient
C     Output: FLXCNV - convective flux (expressed as H, ie F/4/pi)
C             VCONV  - convective velocity
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      COMMON/CUBCON/A,B,DDEL,GRDADB,DLT,RHO,FLXTOT,GRAVD
C
      VCONV=0.
      FLXCNV=0.
      DLT=0.
      IF(HMIX0.LT.0.) RETURN
C 
C     Thermodynamic derivatives
C
      if(ioptab.ge.-1) then
         CALL TRMDER(ID,T,PG,PRAD,TAURS(ID),HEATCP,DLRDLT,GRDADB,RHO)
       else
         call trmdrt(id,t,ptot,heatcp,dlrdlt,grdadb,rho)
      end if
      DDEL=DELTA-GRDADB
C
C     Convective instability criterion
C
      if(idisk.eq.0) then
         HSCALE=PTOT/RHO/GRAV
       else
         if(gravd.eq.0.) return
         hscale=ptot/rho/gravd
      end if
      HMIX=HMIX0
      if(hmix0.eq.0.) hmix=1.
      VCO=HMIX*SQRT(ABS(aconml*PTOT/RHO*DLRDLT))
      FLCO=bconml*RHO*HEATCP*T*HMIX/12.5664
      FC0=FLCO*VCO
      IF(DDEL.LT.0.) RETURN
c
      TAUE=HMIX*ABROS*RHO*HSCALE
      FAC=TAUE/(UN+HALF *TAUE*TAUE)
C
C     Set up parameters A and B (see Mihalas, Eq. 7-76, 7-79, etc)
C
      B=5.67d-5*T**3/(rho*heatcp*VCO)*FAC*cconml*half
      IF(FLXTOT.GT.0.) A=FLCO*VCO/FLXTOT*DELTA
C
C     Determination of   Delta - Delta(E)
C
      D=B*B/2.D0
      DLT=D+DDEL-B*SQRT(D/2.D0+DDEL)
      IF(DLT.LT.0.) DLT=0.
C
C     Resulting convective velocity VCONV and flux FLXCNV
C
      VCONV=VCO*SQRT(DLT)
      FLXCNV=FLCO*VCONV*DLT
      RETURN
      END