SUBROUTINE CONOUT(IMOD,IPRIN)
C     =============================
C
C     Diagnostic outprint of temperature gradients, convective flux,
C     and their derivatives
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      COMMON/CUBCON/A,B,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
C
      IF(IPRIN.GT.0) WRITE(6,600)
      ANEREL=ELEC(1)/(DENS(1)/WMM(1)+ELEC(1))
      ICBEG=0
      FLXTO0=SIG4P*TEFF**4
      DO ID=1,ND
         T=TEMP(ID)
         PTOT=PTOTAL(ID)
         PG=PGS(ID)
         PRAD=PTOT-PG-HALF*DENS(ID)*VTURB(ID)**2
         if(prad.lt.0.) prad=0.
         IF(IMOD.EQ.2) THEN
            if(ioptab.ge.0) then
            CALL OPACF0(ID,NFREQ)
            CALL MEANOP(T,ABSO,SCAT,OPROS,OPPLA)
            ABROSD(ID)=OPROS/DENS(ID)
            else
            call meanopt(t,id,dens(id),opros,oppla)
            if(hmix0.gt.0.) abrosd(id)=opros
            end if
         END IF
         FLXTOT=flxto0
         if(idisk.eq.1) then
            flxtot=flxto0*(un-thetav(id))
            gravd=zd(id)*qgrav
            prad=pradt(id)
         end if
         IF(ID.EQ.1) THEN
            TAU=DM(ID)*ABROSD(ID)
            FLXCR=0.
            GRDADB=0.
            DELTA(ID)=0.
            FLXC(ID)=0.
            FLXCNV=0.
          ELSE
            TM=TEMP(ID-1)
            TAU=TAUM+HALF*(DM(ID)-DM(ID-1))*(ABROSD(ID)+ABROSD(ID-1))
            PTOTM=PTOTAL(ID-1)
            PGM=PGS(ID-1)
            PRADM=PTOTM-PGM-HALF*DENS(ID-1)*VTURB(ID-1)**2
            if(idisk.eq.1) pradm=pradt(id-1)
            if(pradm.lt.0.) pradm=0.
            if(ilgder.eq.0) then
               T0=HALF*(T+TM)
               PT0=HALF*(PTOT+PTOTM)
               PG0=HALF*(PG+PGM)
               PR0=HALF*(PRAD+PRADM)
               AB0=HALF*(ABROSD(ID)+ABROSD(ID-1))
               DLT=(T-TM)/(PTOT-PTOTM)*PT0/T0
             else
               T0=SQRT(T*TM)
               PT0=SQRT(PTOT*PTOTM)
               PG0=SQRT(PG*PGM)
               PR0=SQRT(PRAD*PRADM)
               AB0=SQRT(ABROSD(ID)*ABROSD(ID-1))
               DLT=LOG(T/TM)/LOG(PTOT/PTOTM)  
            end if       
            DELTA(ID)=DLT
            flxcnv=0.
            if(idisk.ne.1.or.id.lt.nd)
     *      CALL CONVEC(ID,T0,PT0,PG0,PR0,AB0,DLT,FLXCNV,VCON)
            if(hmix0.gt.0.) FLXC(ID)=FLXCNV
            IF(ICBEG.EQ.0.AND.FLXC(ID).GT.0..AND.FLXC(ID-1).EQ.0..
     *         AND.ID.GT.25) ICBEG=ID
            if(icbeg.gt.0.and.flxc(id).gt.0.) icend=id
         END IF
         PRADR=PRAD/PTOT
         conrel=0.
         radrel=1.
         if(flxtot.gt.0.) then
           conrel=FLXCNV/FLXTOT
           radrel=flrd(id)/flxtot
         end if
         IF(IPRIN.GT.0) WRITE(6,601) ID,TAU,T,DELTA(ID),
     *                  GRDADB,conrel,radrel,conrel+radrel

         TAUM=TAU
      END DO
      if(iprin.gt.0) write(6,603) icbeg,icend
  603 format(/' convective zone between depths (inclusive)  ',2i4/)
C
c     if ICONV=3, then:
c     in the convective zone the radiative (+ convective) equilibrium
c     is taken obligatorily in the differential form, i.e. 
c     NDRE is modified to have the value just below the beginning of
c     convection zone, and
c     rediff(id) has to be set to unity, and reint(id) to 0 for id => ndre
c
      IF(ICBEG.GT.3.AND.ICONV.EQ.3) THEN 
         NDRE=ICBEG-1
         DO ID=1,ND
            IF(ID.GE.NDRE) THEN
               REINT(ID)=0.
               REDIF(ID)=1.
             else
               reint(id)=1.
               redif(id)=0.
            END IF
         END DO
         WRITE(6,602) NDRE
      END IF
c
      IF(ICBEG.GT.3.AND.ICONV.EQ.2) THEN 
         NDRE=ICBEG-1
         DO ID=1,ND
            IF(ID.GE.NDRE) THEN
               REDIF(ID)=1.
            END IF
         END DO
         WRITE(6,602) NDRE
      END IF
c
  600 FORMAT(//'  ID',3X,'TAUR',7X,'TEMP',6X,
     *       'DELTA',2X,'DELTA(AD)',2X,'CON/TOT  RAD/TOT  (C+R)/TOT'//)
  601 FORMAT(1H ,I4,1PD9.2,0PF9.1,1P5D10.2)
  602 FORMAT(//' NDRE IS RESET IN CONOUT DUE TO THE EXISTENCE OF'
     *       ,' CONVECTIVE ZONE'/'  NDRE= ',I3/)
      RETURN
      END