80 lines
2.2 KiB
Fortran
80 lines
2.2 KiB
Fortran
SUBROUTINE TRMDER(ID,T,PG,PRAD,TAU,HEATCP,DLRDLT,GRDADB,RHO)
|
|
C ============================================================
|
|
C
|
|
C Thermodynamic derivatives
|
|
C Evaluation similar as in Kurucz's ATLAS
|
|
C
|
|
C Input: T - temperature
|
|
C PG - gas pressure
|
|
C PRAD - radiation pressure
|
|
C Output: DEDT - d(energy)/d(T)
|
|
C DRDT - d(rho)/d(T)
|
|
C DEDPG - d(energy)/d(PG)
|
|
C DRDPG - d(rho)/d(PG)
|
|
C RHO - density
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
common/derdif/dift,difp
|
|
common/terden/rhoter,anta,entrp
|
|
common/adiaba/grdad0,itgrad
|
|
DIMENSION TT(5),PP(5),RHON(5),ENER(5),entr(5)
|
|
C
|
|
C NUMERICAL EVALUATION OF THERMODYNAMIC DERIVATIVES
|
|
C
|
|
c prad=0.
|
|
P=PG
|
|
TT(1)=T*(UN+DIFT)
|
|
TT(2)=T*(UN-DIFT)
|
|
TT(3)=T
|
|
TT(4)=T
|
|
TT(5)=T
|
|
PP(1)=P
|
|
PP(2)=P
|
|
PP(3)=P*(UN+DIFP)
|
|
PP(4)=P*(UN-DIFT)
|
|
PP(5)=P
|
|
tmoli0=tmolim
|
|
if(t.lt.tmolim) tmolim=t*(un+dift+0.001)
|
|
if(t.ge.tmolim) tmolim=t*(in-dift-0.001)
|
|
DO I=1,5
|
|
TE=TT(I)
|
|
TKN=TE*1.38054D-16
|
|
ANT=PP(I)/TKN
|
|
CALL ELDENS(ID,TE,ANT,ANE,ENRG,ENT,WM,0)
|
|
RHON(I)=rhoter
|
|
ENER(I)=(1.5D0*PP(I)+ENRG+3.D0*PRAD*(TE/T)**4)/RHON(I)
|
|
entr(i)=ent/RHON(I)
|
|
END DO
|
|
tmolim=tmoli0
|
|
entrp=entr(5)
|
|
c
|
|
DRDT=(RHON(1)-RHON(2))/(2.*T*DIFT)
|
|
DRDPG=(RHON(3)-RHON(4))/(2.*P*DIFP)
|
|
RHO=RHON(5)
|
|
DPDPG=1.
|
|
dpdt=0.
|
|
if(tau.lt.50.) DPDT=4.*PRAD/T*(un-exp(-tau))
|
|
DLRDLT=T/RHO*(DRDT-DRDPG*DPDT/DPDPG)
|
|
ptot=pg+prad
|
|
c
|
|
if(ifentr.le.0) then
|
|
DEDT= (ENER(1)-ENER(2))/(2.*T*DIFT)
|
|
DEDPG=(ENER(3)-ENER(4))/(2.*P*DIFP)
|
|
HEATCV=DEDT-DEDPG*DRDT/DRDPG
|
|
HEATCP=DEDT-DEDPG*DPDT/DPDPG-
|
|
* PTOT/RHO/RHO*(DRDT-DRDPG*DPDT/DPDPG)
|
|
GRDADB=-PTOT/RHO/T*DLRDLT/HEATCP
|
|
else
|
|
DSDT= (ENTR(1)-ENTR(2))/(2.*T*DIFT)
|
|
DSDP= (ENTR(3)-ENTR(4))/(2.*P*DIFP)
|
|
grdadb=-dsdp/dsdt*pg/t
|
|
c heatcp=-PTOT/RHO/T*DLRDLT/grdadb
|
|
heatcp=t*dsdt
|
|
end if
|
|
c
|
|
if(iter.le.itgrad.and.grdad0.gt.0)
|
|
* grdadb=grdad0*0.4+(un-grdad0)*grdadb
|
|
RETURN
|
|
END
|