100 lines
3.1 KiB
Fortran
100 lines
3.1 KiB
Fortran
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
|