SUBROUTINE TEMCOR
C     =================
C
C     Auxiliary procedure for INILAM.
C     Tests wheather the convective flux corresponding to newly
C     determined temperature and logarithmic gradient DELTA
C     is larger than the total flux. If so, the temperature is modified
C     by an iterative procedure for determining new temperature
C     that yields convective flux less than SIG4P*TEFF**4
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ARRAY1.FOR'
      INCLUDE 'ALIPAR.FOR'
      COMMON/CUBCON/ACNV,BCNV,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
C
      IF(ICONV.LE.0.AND.INDL.EQ.0) RETURN
      FLXTO0=SIG4P*TEFF**4
      ANEREL=ELEC(1)/(DENS(1)/WMM(1)+ELEC(1))
      DLTND=DELTA(ND)
      IFNDM1=0
C
      DO 100 ID=2,ND
         flxtot=flxto0
         if(idisk.eq.1) then
            flxtot=flxto0*(un-thetav(id))
            gradv=zd(id)*qgrav
         end if
         T=TEMP(ID)
         P=PTOTAL(ID)
         PG=PGS(ID)
         PRAD=P-PG-HALF*DENS(ID)*VTURB(ID)**2
         TM=TEMP(ID-1)
         PM=PTOTAL(ID-1)
         PGM=PGS(ID-1)
         PRADM=PM-PGM-HALF*DENS(ID-1)*VTURB(ID-1)**2
         IF(ID.EQ.ND.AND.IFNDM1.EQ.1) THEN
            FAC=DLTND*(P-PM)/(P+PM)
            T=TM*(UN+FAC)/(UN-FAC)
         END IF
         KKK=0
   10    KKK=KKK+1
         T0=HALF*(T+TM)
         P0=HALF*(P+PM)
         PG0=HALF*(PG+PGM)
         PR0=HALF*(PRAD+PRADM) 
         AB0=HALF*(ABROSD(ID)+ABROSD(ID-1))
         DLT=(T-TM)/(P-PM)*P0/T0
         DELTA(ID)=DLT
C
C        convective flux and its derivatives
C
         CALL CONVEC(ID,T0,P0,PG0,PR0,AB0,DLT,FLXCNV,VCON)
         FLXC(ID)=FLXCNV
         IF(FLXCNV.LT.0.999999*FLXTOT) GO TO 100
C
C        in case that convective flux it too large:
C        iterative procedure for determining new temperature such
C        that yields convective flux less than SIG4P*TEFF**4
C
C        Basically by the Newton-Raphson method;
C        derivative of the convective flux wrt. T calculated
C        numerically; derivative wrt. DELTA calculated analytically
C
         DHCDD=0.
         IF(DELMDE.GT.0.) DHCDD=1.5D0/DELMDE*FLXCNV
         T1=1.001D0*T0
         CALL CONVEC(ID,T1,P0,PG0,PR0,AB0,DLT,FLXC1,VCON)
         DHCDT=(FLXC1-FLXCNV)*1.D3/T0*HALF
         TT=T*T-TM*TM
         DDT0= DLT/HALF*TM/TT
         DFLCDT=DHCDT+DHCDD*DDT0
         DELTEM=(FLXTOT-FLXCNV)/DFLCDT
         T1=T+DELTEM
         IF(T1.LT.TM+HALF*(T-TM)) T1=TM+HALF*(T-TM)
         T=T1
         write(94,601) iter,id,kkk,delta(id),temp(id),t,deltem,flxcnv
  601    format(3i4,f10.4,2f12.3,f12.5,1pe15.6)
         TEMP(ID)=T
         IF(ID.EQ.ND-1) IFNDM1=1
         IF(KKK.LE.10.and.abs(deltem/t).gt.1.e-5) GO TO 10
C
C        Determination of electron density from the total pressure
C
         AN=PG/T/BOLK
         CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1)
         RHO=WMM(ID)*(AN-ANE)
         DENS(ID)=RHO
         ELEC(ID)=ANE
         CALL WNSTOR(ID)
         CALL STEQEQ(ID,POP,1)
         CALL OPACF0(ID,NFREQ)
         CALL MEANOP(T,ABSO,SCAT,OPROS,OPPLA)
         ABROS=OPROS/DENS(ID)
         ABROSD(ID)=ABROS
  100 CONTINUE
      RETURN
      END