223 lines
6.4 KiB
Fortran
223 lines
6.4 KiB
Fortran
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
|