SUBROUTINE CONTMP C ================= C C Auxiliary procedure for LTEGR C Determination of temperature in convectively unstable layers 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/ichndm/ichanm 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. C PRAD0=DPRAD/1.732D0 DO ID=1,ND TEMPR(ID)=TEMP(ID) IF(ID.EQ.1) THEN DELTR(ID)=0. ELSE if(ilgder.eq.0) then DELTR(ID)= * (TEMP(ID)-TEMP(ID-1))/(PTOTAL(ID)-PTOTAL(ID-1))* * (PTOTAL(ID)+PTOTAL(ID-1))/(TEMP(ID)+TEMP(ID-1)) else DELTR(ID)= * LOG(TEMP(ID)/TEMP(ID-1))/LOG(PTOTAL(ID)/PTOTAL(ID-1)) end if END IF END DO ICONIT=0 C C ------------------------------------------------------ C Global iteration loop for calculating convective model C ------------------------------------------------------ C 20 ICONIT=ICONIT+1 ICONBE=0 CHANTM=0. DO ID=1,ND T=TEMP(ID) PTOT=PTOTAL(ID) PGAS=PGS(ID) PTURB=HALF*DENS(ID)*VTURB(ID)*VTURB(ID) PRAD=PTOT-PGAS-PTURB PRADT(ID)=PRAD FLXTOT=FLXTO0 IF(IDISK.EQ.1) THEN FLXTOT=FLXTO0*(UN-THETA(ID)) GRAVD=ZD(ID)*QGRAV END IF ICON0(ID)=0 C IF(ID.EQ.1) GO TO 45 J=0 IF(ICONIT.EQ.1) T=T-TEMPR(ID-1)+TEMP(ID-1) TM=TEMP(ID-1) IF(T.LT.0.) T=TM PTOTM=PTOTAL(ID-1) if(ilgder.eq.0) then PT0=HALF*(PTOT+PTOTM) else pt0=sqrt(ptot*ptotm) end if DELR=DELTR(ID) C C Inner iteration loop for determining temperature in the C conectively unstable layers C 30 J=J+1 TOLD=T if(ilgder.eq.0) then T0=HALF*(T+TM) PG0=HALF*(PGAS+PGM) PR0=HALF*(PRAD+PRADM) AB0=HALF*(ABROSD(ID)+ABROSD(ID-1)) else t0=sqrt(t*tm) pg0=sqrt(pgas*pgm) pr0=sqrt(prad*pradm) ab0=sqrt(abrosd(id)*abrosd(id-1)) end if 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..or.id.lt.idconz) GO TO 40 ICON0(ID)=1 ICONBE=1 CALL CUBIC(DELTA0) REFF=DELTA0/DELR PRAD=PRADM+(TAUROS(ID)-TAUROS(ID-1))*DPRAD*REFF PRADT(ID)=PRAD PGAS=PTOT-PRAD-PTURB if(ilgder.eq.0) then IF(REFF.GT.UN) REFF=UN IF(REFF.LT.0.) REFF=0. FAC=DELTA0*(PTOT-PTOTM)/(PTOT+PTOTM) T=TM*(UN+FAC)/(UN-FAC) IF(T.LT.TM) T=TM else T=TM*(PTOT/PTOTM)**DELTA0 IF(T.LT.TM) T=TM*1.0001 end if 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 45 DELT0=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 PGM=PGAS PRADM=PRAD PGS(ID)=PGAS 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 CONTMP, ITER=',I2/) C C 2. New values of electron density, density, sound spped, C and mean opacities and optical depths C c ANEREL=ELEC(1)/(DENS(1)/WMM(1)+ELEC(1)) DO 70 ID=1,ND T=TEMP(ID) P=PTOTAL(ID) ITINT=0 60 ITINT=ITINT+1 if(ioptab.ge.-1) then AN=PGS(ID)/T/BOLK CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1) ELEC(ID)=ANE DENS(ID)=WMM(ID)*(AN-ANE) PHMOL(ID)=AHMOL C C Corresponding LTE populations C if(ioptab.ge.0) then CALL WNSTOR(ID) CALL STEQEQ(ID,POP,1) C C Evaluation of the Rosseland and Planck mean opacities C CALL OPACF0(ID,NFREQ) CALL MEANOP(T,ABSO,SCAT,OPROS,OPPLA) ABROS=OPROS/DENS(ID) ABPLA=OPPLA/DENS(ID) else rho=dens(id) call meanopt(t,id,rho,abros,abpla) abrosd(id)=abros abplad(id)=abpla end if else rho=rhoeos(t,p) dens(id)=rho call meanopt(t,id,rho,abros,abpla) abrosd(id)=abros abplad(id)=abpla end if C C New values of the the column mass C PTOLD=PTOTAL(ID) if(idisk.eq.0.and.ichanm.gt.0) then IF(ID.EQ.1) THEN DM(ID)=TAUROS(ID)/ABROS PTOTAL(ID)=DM(ID)*GRAV+PRAD0 ELSE DM(ID)=DM(ID-1)+(TAUROS(ID)-TAUROS(ID-1))/ * HALF/(ABROSD(ID-1)+ABROS) PTOTAL(ID)=DM(ID)*GRAV+PRAD0 END IF C C Store the final quantitites C PTURB=HALF*DENS(ID)*VTURB(ID)*VTURB(ID) PGS(ID)=PTOTAL(ID)-PRADT(ID)-PTURB end if ABROSD(ID)=ABROS ABPLAD(ID)=ABPLA IF((PTOTAL(ID)-PTOLD)/PTOLD.LT.1.D-3) GO TO 70 IF(ITINT.GT.5) THEN WRITE(6,601) ID,PTOLD,PTOTAL(ID) GO TO 70 ELSE GO TO 60 END IF 70 CONTINUE C *** TEMPORARY C IF(ICONIT.LT.NCONIT) GO TO 20 601 FORMAT(1H0,'SLOW CONVERGENCE OF INTERNAL ITERATIONS IN', * ' CONTMP: ID, PTOT(OLD), PTOT(NEW) ='/I3,1P2D10.2/) RETURN END