145 lines
4.2 KiB
Fortran
145 lines
4.2 KiB
Fortran
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
|