723 lines
20 KiB
Fortran
723 lines
20 KiB
Fortran
SUBROUTINE RHSGEN(ID)
|
|
C =====================
|
|
C
|
|
C Auxiliary procedure for SOLVE
|
|
C controls evaluation of rhs vector when
|
|
C matrices A,B, and C are kept fixed.
|
|
C
|
|
C Input: ID - depth index
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ARRAY1.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
PARAMETER (XCON=8.0935D-21,YCON=1.68638E-10)
|
|
PARAMETER (SIXTH=UN/6.D0,
|
|
* THIRD=UN/3.D0)
|
|
C DIMENSION AJ(MLEVEL)
|
|
DIMENSION POPP(MLEVEL)
|
|
COMMON/CUBCON/ACNV,BCNV,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
|
|
C
|
|
ispl=isplin
|
|
if(isplin.ge.5) isplin=isplin-5
|
|
IJ1=1
|
|
C
|
|
C in the case of Compton scattering - boundary condition
|
|
C for the highest frequency
|
|
C
|
|
if(icompt.gt.0.and.icombc.gt.0.and.ijex(1).gt.0) then
|
|
IJ1=2
|
|
ij=1
|
|
iji=nfreq
|
|
zj1=exp(-hk*freq(ij)/temp(id))
|
|
zj2=exp(-hk*freq(ij+1)/temp(id))
|
|
dlt=delj(iji-1,id)
|
|
if(ichcoo.eq.0) then
|
|
zj0=un/(hk*sqrt(freq(ij)*freq(ij+1))/temp(id))
|
|
zxx=un-3.*zj0+(un-dlt)*zj1+dlt*zj2
|
|
combid=zj0/dlnfr(iji-1)+(un-dlt)*zxx
|
|
comaid=-zj0/dlnfr(iji-1)+dlt*zxx
|
|
else
|
|
e2=ycon*temp(id)
|
|
zxx0=xcon*freq(ij)*(un+zj1)-3.*e2
|
|
zxxm=xcon*freq(ij+1)*(un+zj2)-3.*e2
|
|
zxx=(un-dlt)*zxx0+dlt*zxxm
|
|
combid=e2/dlnfr(iji-1)+(un-dlt)*zxx
|
|
comaid=-e2/dlnfr(iji-1)+dlt*zxx
|
|
end if
|
|
vecl(ij)=-combid*rad(iji,id)-comaid*rad(iji-1,id)
|
|
end if
|
|
C
|
|
C evaluation of the opacity, emissivity, scattering, and
|
|
C their derivatives, at the current depth point ID;
|
|
C
|
|
DO I=1,NN
|
|
PSI0(I)=PSY0(I,ID)
|
|
END DO
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
WDEP0(IJ)=W(IJT)
|
|
RAD0(IJ)=RADEX(IJ,ID)
|
|
FK0(IJ)=FAKEX(IJ,ID)
|
|
ABSO0(IJ)=ABSOEX(IJ,ID)
|
|
EMIS0(IJ)=EMISEX(IJ,ID)
|
|
SCAT0(IJ)=SCATEX(IJ,ID)
|
|
END DO
|
|
END IF
|
|
C
|
|
IF(ID.GT.1) THEN
|
|
DO I=1,NN
|
|
PSIM(I)=PSY0(I,ID-1)
|
|
END DO
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
RADM(IJ)=RADEX(IJ,ID-1)
|
|
FKM(IJ)=FAKEX(IJ,ID-1)
|
|
ABSOM(IJ)=ABSOEX(IJ,ID-1)
|
|
EMISM(IJ)=EMISEX(IJ,ID-1)
|
|
SCATM(IJ)=SCATEX(IJ,ID-1)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
C
|
|
IF(ID.LT.ND) THEN
|
|
DO I=1,NN
|
|
PSIP(I)=PSY0(I,ID+1)
|
|
END DO
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
RADP(IJ)=RADEX(IJ,ID+1)
|
|
FKP(IJ)=FAKEX(IJ,ID+1)
|
|
ABSOP(IJ)=ABSOEX(IJ,ID+1)
|
|
EMISP(IJ)=EMISEX(IJ,ID+1)
|
|
SCATP(IJ)=SCATEX(IJ,ID+1)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
c
|
|
C
|
|
C ------------------------------------------------------------
|
|
C Actual evaluation of the rhs vector VECL
|
|
C ------------------------------------------------------------
|
|
C
|
|
DO I=1,NN
|
|
VECL(I)=0.
|
|
END DO
|
|
C
|
|
IF(NFREQE.LE.0) GO TO 100
|
|
C
|
|
C -----------------------------------------------------------
|
|
C 1. Radiative transfer components
|
|
C -----------------------------------------------------------
|
|
C
|
|
C For ID = 1 - upper boundary condition
|
|
C
|
|
IF(ID.GT.1) GO TO 50
|
|
DDP=DELDMZ(1)
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
IF(IZSCAL.EQ.0) THEN
|
|
OMEG0=ABSO0(IJ)/DENS(ID)
|
|
OMEGP=ABSOP(IJ)/DENS(ID+1)
|
|
ELSE
|
|
OMEG0=ABSO0(IJ)
|
|
OMEGP=ABSOP(IJ)
|
|
END IF
|
|
DZP=OMEG0+OMEGP
|
|
DTAUP=DZP*DDP
|
|
ALF1=(FK0(IJ)*RAD0(IJ)-FKP(IJ)*RADP(IJ))/DTAUP
|
|
CHIEL0=SCAT0(IJ)
|
|
CHIELP=SCATP(IJ)
|
|
S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
|
|
BS=HALF*DTAUP
|
|
CS=0.
|
|
C2=0.
|
|
SP=0.
|
|
BET2=0.
|
|
GAM2=0.
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.0) then
|
|
call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
|
|
s0=s0+cms
|
|
end if
|
|
C
|
|
IF(ISPLIN.GT.0) THEN
|
|
C
|
|
C Spline collocation and/or Hermitian method (ISPLIN=1 or 2) -
|
|
C both give the same expression for the boundary conditions
|
|
C
|
|
BS=DTAUP*THIRD
|
|
CS=HALF*BS
|
|
SP=(EMISP(IJ)+CHIELP*RADP(IJ))/ABSOP(IJ)
|
|
C2=CS/ABSOP(IJ)
|
|
GAM2=CS*(RADP(IJ)-SP)
|
|
END IF
|
|
C
|
|
C auxiliary quantities
|
|
C
|
|
ALF2=BS*(RAD0(IJ)-S0)
|
|
BET2=ALF2+GAM2
|
|
C
|
|
C *** the IJ-th element of the rhs vector
|
|
C
|
|
VECL(IJ)=ALF1+BET2+FH(IJT)*RAD0(IJ)-S0*Q0(IJT)
|
|
IF(IWINBL.LT.0) VECL(IJ)=VECL(IJ)-HEXTRD(IJT)
|
|
C
|
|
END DO
|
|
GO TO 100
|
|
C
|
|
C For 1 < ID < ND - normal depth point
|
|
C
|
|
50 DDM=DELDMZ(ID-1)
|
|
IF(ID.EQ.ND) GO TO 80
|
|
DDP=DELDMZ(ID)
|
|
DO IJ=1,NFREQE
|
|
IF(IZSCAL.EQ.0) THEN
|
|
OMEG0=ABSO0(IJ)/DENS(ID)
|
|
OMEGP=ABSOP(IJ)/DENS(ID+1)
|
|
OMEGM=ABSOM(IJ)/DENS(ID-1)
|
|
ELSE
|
|
OMEG0=ABSO0(IJ)
|
|
OMEGP=ABSOP(IJ)
|
|
OMEGM=ABSOM(IJ)
|
|
END IF
|
|
DZP=OMEG0+OMEGP
|
|
DZM=OMEG0+OMEGM
|
|
DTAUP=DZP*DDP
|
|
DTAUM=DZM*DDM
|
|
DTAU0=HALF*(DTAUP+DTAUM)
|
|
FRD=FK0(IJ)*RAD0(IJ)
|
|
ALF1=(FRD-FKP(IJ)*RADP(IJ))/DTAUP/DTAU0
|
|
GAM1=(FRD-FKM(IJ)*RADM(IJ))/DTAUM/DTAU0
|
|
BET1=ALF1+GAM1
|
|
BS=UN
|
|
CHIELM=SCATM(IJ)
|
|
CHIEL0=SCAT0(IJ)
|
|
CHIELP=SCATP(IJ)
|
|
S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
|
|
BET2=0.
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.0) then
|
|
call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
|
|
s0=s0+cms
|
|
end if
|
|
C
|
|
IF(ISPLIN.EQ.1) THEN
|
|
C
|
|
C spline collocation (ISPLIN=1)
|
|
C
|
|
AS=DTAUM/DTAU0*SIXTH
|
|
CS=DTAUP/DTAU0*SIXTH
|
|
BS=0.666666666666667D0
|
|
SM=(EMISM(IJ)+RADM(IJ)*CHIELM)/ABSOM(IJ)
|
|
SP=(EMISP(IJ)+RADP(IJ)*CHIELP)/ABSOP(IJ)
|
|
ALF2=AS*(RADM(IJ)-SM)
|
|
GAM2=CS*(RADP(IJ)-SP)
|
|
BET2=ALF2+GAM2
|
|
ELSE IF(ISPLIN.EQ.2) THEN
|
|
C
|
|
C Hermitian method (ISPLIN=2)
|
|
C
|
|
AS=DTAUP*DTAUP/DTAUM/DTAU0
|
|
CS=DTAUM*DTAUM/DTAUP/DTAU0
|
|
SM=(EMISM(IJ)+RADM(IJ)*CHIELM)/ABSOM(IJ)
|
|
SP=(EMISP(IJ)+RADP(IJ)*CHIELP)/ABSOP(IJ)
|
|
AS=(UN-HALF*AS)*SIXTH
|
|
CS=(UN-HALF*CS)*SIXTH
|
|
BS=UN-AS-CS
|
|
BET2=AS*(RADM(IJ)-SM)+CS*(RADP(IJ)-SP)
|
|
END IF
|
|
C
|
|
C *** the IJ-th element of the rhs vector
|
|
C
|
|
VECL(IJ)=BET1+BET2+BS*(RAD0(IJ)-S0)
|
|
C
|
|
END DO
|
|
GO TO 100
|
|
C
|
|
C For ID=ND - lower boundary condition
|
|
C
|
|
80 T=TEMP(ID)
|
|
IF(TEMPBD.GT.0.) T=TEMPBD
|
|
HKT=HK/T
|
|
C
|
|
c the case of stellar atmospheres
|
|
C
|
|
IF(IDISK.EQ.0.OR.IFZ0.LT.0) THEN
|
|
C
|
|
C 1. the case NDRE.le.ND, ie. radiative equilibrium equation at this
|
|
C depth is treated as differential equation.
|
|
C
|
|
C 2. the case NDRE > ND, ie. radiative equilibrium equation at this
|
|
C depth is treated as integral equation.
|
|
C Boundary condition of the transfer equation is then complex,
|
|
C because it has to incorporate the constraint of the total flux
|
|
C equal to SIG4P*Teff**4 (integral form of the radiative
|
|
C equation does not fix the total flux)
|
|
C
|
|
C Details of this approach are presented in Mihalas, Heasley, Auer
|
|
C NCAR-TN/STR-104 (1975).
|
|
C
|
|
C This option is presented here for the sake of continuity with
|
|
C previous work; otherwise the first option (NDRE.le.ND) is
|
|
C recommended due to its better numerical properties
|
|
C
|
|
C auxiliary quantities for the second option
|
|
C
|
|
IF(ID.LT.NDRE.AND.INRE.NE.0) THEN
|
|
SUMB=0.
|
|
SUMF=0.
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
FR=FREQ(IJT)
|
|
FR15=FR*1.D-15
|
|
WD=WDEP0(IJ)
|
|
X=HKT*FR
|
|
EX=EXP(X)
|
|
PLAN=BN*FR15*FR15*FR15/(EX-UN)*RRDIL
|
|
DPLAN=PLAN*X/T/(UN-UN/EX)*WD
|
|
SUMB=SUMB+(PLAN-RAD0(IJ))*WD
|
|
FI=DPLAN/ABSO0(IJ)
|
|
SUMF=SUMF+FI
|
|
END DO
|
|
FL=SIG4P*TEFF*TEFF*TEFF*TEFF
|
|
ZZ=(FL-HALF*SUMB)/SUMF
|
|
END IF
|
|
C
|
|
C auxiliary quantities for both options
|
|
C
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
CHIELM=SCATM(IJ)
|
|
CHIEL0=SCAT0(IJ)
|
|
IF(IZSCAL.EQ.0) THEN
|
|
OMEG0=ABSO0(IJ)/DENS(ID)
|
|
OMEGM=ABSOM(IJ)/DENS(ID-1)
|
|
ELSE
|
|
OMEG0=ABSO0(IJ)
|
|
OMEGM=ABSOM(IJ)
|
|
END IF
|
|
DZM=OMEG0+OMEGM
|
|
DTAUM=DZM*DDM
|
|
GAM1=(FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ))/DTAUM
|
|
BET2=0.
|
|
C
|
|
C second-order boundary condition
|
|
C
|
|
IF(IBC.GT.0.AND.IBC.LT.4) THEN
|
|
BS=DTAUM*HALF
|
|
S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.0) then
|
|
call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
|
|
s0=s0+cms
|
|
end if
|
|
C
|
|
BET2=BS*(RAD0(IJ)-S0)
|
|
END IF
|
|
C
|
|
FR=FREQ(IJT)
|
|
X=HKT*FR
|
|
EX=EXP(X)
|
|
PLAN=BN*(FR*1.D-15)**3/(EX-UN)*RRDIL
|
|
IF(INRE.EQ.0.OR.ID.GE.NDRE) THEN
|
|
DPLAN=BN*(FR*1.D-15)**3/(EXP(HK*FR/TEMP(ID-1))-UN)
|
|
GAM1=GAM1-(PLAN-DPLAN)/DTAUM*THIRD
|
|
ELSE
|
|
DPLAN=PLAN*X/T/(UN-UN/EX)
|
|
FI=DPLAN/ABSO0(IJ)
|
|
X1=FI*ZZ
|
|
GAM1=GAM1-X1
|
|
END IF
|
|
C
|
|
C *** the IJ-th element of the rhs vector
|
|
C
|
|
IF(IBC.EQ.0.OR.IBC.EQ.4) THEN
|
|
VECL(IJ)=GAM1+BET2-HALF*(PLAN-RAD0(IJ))
|
|
ELSE
|
|
VECL(IJ)=GAM1+BET2-HALF*PLAN+FHD(IJT)*RAD0(IJ)
|
|
END IF
|
|
END DO
|
|
C
|
|
C the case of accretion disks
|
|
C
|
|
ELSE
|
|
C
|
|
DO IJ=IJ1,NFREQE
|
|
CHIELM=SCATM(IJ)
|
|
CHIEL0=SCAT0(IJ)
|
|
IF(IZSCAL.EQ.0) THEN
|
|
OMEG0=ABSO0(IJ)/DENS(ID)
|
|
OMEGM=ABSOM(IJ)/DENS(ID-1)
|
|
ELSE
|
|
OMEG0=ABSO0(IJ)
|
|
OMEGM=ABSOM(IJ)
|
|
END IF
|
|
DZM=OMEG0+OMEGM
|
|
DTAUM=DZM*DDM
|
|
FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ)
|
|
GAM1=FRD/DTAUM
|
|
BS=DTAUM*HALF
|
|
S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
|
|
GAM2=BS*(RAD0(IJ)-S0)
|
|
C
|
|
C *** the IJ-th element of the rhs vector
|
|
C
|
|
VECL(IJ)=GAM1+GAM2
|
|
END DO
|
|
END IF
|
|
C
|
|
C -----------------------------------------------------------
|
|
C 2. Hydrostatic equilibrium
|
|
C -----------------------------------------------------------
|
|
C
|
|
100 isplin=ispl
|
|
IF(INHE.EQ.0) GO TO 170
|
|
NHE=NFREQE+INHE
|
|
GRD=0.
|
|
C
|
|
IF(ID.EQ.1) THEN
|
|
C
|
|
C Upper boundary condition (ID=1)
|
|
C Basically, linearized eq. (7-10) of Mihalas (1978)
|
|
C
|
|
IF(IDISK.EQ.0.OR.IBCHE.EQ.0) THEN
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
IF(.NOT.LSKIP(ID,IJT)) THEN
|
|
FLUXW=W(IJT)*FH(IJT)*RAD0(IJ)
|
|
GRD=GRD+FLUXW*ABSO0(IJ)
|
|
END IF
|
|
END DO
|
|
END IF
|
|
C
|
|
X1=PCK/DENS(ID)
|
|
VT0=HALF*VTURB(ID)*VTURB(ID)/DM(ID)*WMM(ID)
|
|
C
|
|
C The rhs vector also accounts for the total radiation pressure in
|
|
C the fixed-option transitions (array FPRD, generated by FIXLIN)
|
|
C
|
|
VECL(NHE)=GRAV-BOLK*TEMP(ID)*PSI0(NHE)/DM(ID)-
|
|
* X1*(GRD+FPRD(ID))-VT0/WMM(ID)*DENS(ID)
|
|
C
|
|
ELSE
|
|
IF(IBCHE.EQ.1) THEN
|
|
C
|
|
C specifically disk - newer variant
|
|
C
|
|
CCC=PCK/QGRAV
|
|
HR1=CCC*SIG4P*TEFF**4*ABROSD(1)
|
|
PG1=BOLK*PSI0(NHE)*TEMP(1)
|
|
HG1=SQRT(TWO*PG1/DENS(1)/QGRAV)
|
|
X=(ZD(1)-PR1)/HG1
|
|
IF(X.LT.3.) THEN
|
|
IF(X.LT.0.) X=0.
|
|
F1=8.86226925D-1*EXP(X*X)*ERFCX(X)
|
|
ELSE
|
|
F1=HALF*(UN-HALF/X/X)/X
|
|
END IF
|
|
GGG=DENS(1)*HG1*F1
|
|
VECL(NHE)=DM(1)-GGG
|
|
C
|
|
C specifically disk - older variant
|
|
C
|
|
ELSE IF(IBCHE.EQ.2) THEN
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=IJ1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
IF(.NOT.LSKIP(ID,IJT)) THEN
|
|
FLUXW=W(IJT)*FH(IJT)*RAD0(IJ)
|
|
GRD=GRD+FLUXW*ABSO0(IJ)
|
|
END IF
|
|
END DO
|
|
END IF
|
|
C
|
|
CCC=PCK/QGRAV
|
|
PR1=CCC*(GRD+FPRD(1))/DENS(1)
|
|
PG1=BOLK*PSI0(NHE)*TEMP(1)
|
|
HG1=SQRT(TWO*PG1/DENS(1)/QGRAV)
|
|
X=(ZD(1)-PR1)/HG1
|
|
IF(X.LT.3.) THEN
|
|
IF(X.LT.0.) X=0.
|
|
F1=8.86226925D-1*EXP(X*X)*ERFCX(X)
|
|
ELSE
|
|
F1=HALF*(UN-HALF/X/X)/X
|
|
END IF
|
|
GGG=HG1*QGRAV*HALF/F1
|
|
C
|
|
C The rhs vector
|
|
C
|
|
VECL(NHE)=DM(1)*GGG-PG1
|
|
ELSE
|
|
VECL(NHE)=PGAS0-BOLK*TEMP(1)*PSI0(NHE)
|
|
END IF
|
|
END IF
|
|
C
|
|
ELSE
|
|
C
|
|
C Normal depth point (ID > 1)
|
|
C
|
|
C
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
IF(.NOT.LSKIP(ID,IJFR(IJ)))
|
|
* GRD=GRD+(FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ))*W(IJFR(IJ))
|
|
END DO
|
|
END IF
|
|
C
|
|
VT0=HALF*VTURB(ID)*VTURB(ID)*WMM(ID)
|
|
VTM=HALF*VTURB(ID-1)*VTURB(ID-1)*WMM(ID-1)
|
|
C
|
|
C the rhs vector
|
|
C again, which accounts for the total radiation pressure in the
|
|
C fixed-option transitions (array FPRD, generated by FIXLIN)
|
|
C
|
|
IF(IDISK.EQ.1) GRAV=QGRAV*(ZD(ID)+ZD(ID-1))*HALF
|
|
IF(IZSCAL.EQ.0) THEN
|
|
VECL(NHE)=GRAV*(DM(ID)-DM(ID-1))-
|
|
* BOLK*(TEMP(ID)*PSI0(NHE)-TEMP(ID-1)*PSIM(NHE))-
|
|
* PCK*(GRD+FPRD(ID))-
|
|
* VT0/WMM(ID)*DENS(ID)+VTM/WMM(ID-1)*DENS(ID-1)
|
|
ELSE
|
|
GRAVZ=GRAV*(ZD(ID)-ZD(ID-1))
|
|
VECL(NHE)=-GRAVZ*(DENS(ID)+DENS(ID-1))*HALF-
|
|
* BOLK*(TEMP(ID)*PSI0(NHE)-TEMP(ID-1)*PSIM(NHE))-
|
|
* PCK*(GRD+FPRD(ID))-
|
|
* VT0/WMM(ID)*DENS(ID)+VTM/WMM(ID-1)*DENS(ID-1)
|
|
END IF
|
|
END IF
|
|
C
|
|
C -----------------------------------------------------------
|
|
C 2a. z-m relation
|
|
C -----------------------------------------------------------
|
|
C
|
|
170 IF(INZD.LE.0.OR.ID.EQ.ND.OR.IDISK.EQ.0) GO TO 200
|
|
NZD=NFREQE+INZD
|
|
DDP=(DM(ID+1)-DM(ID))*HALF
|
|
VECL(NZD)=ZD(ID+1)-ZD(ID)+DDP/DENS(ID)+DDP/DENS(ID+1)
|
|
C
|
|
C -----------------------------------------------------------
|
|
C 3. Radiative equilibrium
|
|
C -----------------------------------------------------------
|
|
C
|
|
200 IF(INRE.EQ.0) GO TO 250
|
|
NRE=NFREQE+INRE
|
|
c
|
|
ittc=abs(nretc)/100
|
|
if(iter.gt.ittc) then
|
|
if(id.le.mod(abs(nretc),100)) then
|
|
if(nretc.lt.0) vecl(nre)=temp(id+1)-temp(id)
|
|
go to 250
|
|
end if
|
|
end if
|
|
C
|
|
C integral equation part of the radiative
|
|
C equilibrium equation
|
|
C
|
|
C the rhs vector accounts for total net cooling in ALI transitions
|
|
C
|
|
VECL(NRE)=FCOOL(ID)
|
|
IF(IDISK.EQ.1) VECL(NRE)=FCOOL(ID)-TVISC(ID)*reint(id)
|
|
C
|
|
IF(REINT(ID).GT.0.AND.NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
HEAT=ABSO0(IJ)-SCAT0(IJ)
|
|
VECL(NRE)=VECL(NRE)-
|
|
* (HEAT*RAD0(IJ)-EMIS0(IJ))*WDEP0(IJ)*reint(id)
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.5) then
|
|
ijt=ijfr(ij)
|
|
call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
|
|
vecl(nre)=vecl(nre)+abso0(ij)*cms*wdep0(ij)*reint(id)
|
|
end if
|
|
c
|
|
END DO
|
|
END IF
|
|
C
|
|
C
|
|
C differential equation part of the
|
|
C radiative equilibrium equation
|
|
C
|
|
if(redif(id).gt.0) then
|
|
TEFFD=TEFF**4
|
|
IF(IDISK.EQ.1) TEFFD=TEFFD*(UN-THETAV(ID))
|
|
VECL(NRE)=vecl(nre)+SIG4P*TEFFD*redif(id)
|
|
IF(ID.GT.1) THEN
|
|
DDM=(DM(ID)-DM(ID-1))*HALF
|
|
DDM=DELDMZ(ID-1)
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
OMEG0=ABSO0(IJ)*DENSI(ID)
|
|
OMEGM=ABSOM(IJ)*DENSI(ID-1)
|
|
DTAUM=(OMEG0+OMEGM)*DDM
|
|
FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ)
|
|
GAMR=FRD/DTAUM
|
|
VECL(NRE)=VECL(NRE)-WDEP0(IJ)*GAMR*redif(id)
|
|
END DO
|
|
END IF
|
|
ELSE
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
VECL(NRE)=VECL(NRE)-
|
|
* WDEP0(IJ)*FH(IJT)*RAD0(IJ)*REDIF(ID)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
end if
|
|
C
|
|
C -----------------------------------------------------------
|
|
C 4. Statistical equilibrium
|
|
C -----------------------------------------------------------
|
|
C
|
|
250 IF(IABS(IFPOPR).GE.3.and.ifpopr.le.5) THEN
|
|
if(inse.gt.0) then
|
|
NSE=NFREQE+INSE-1
|
|
CALL SABOLF(ID)
|
|
CALL LEVGRP(ID,IIEXP,0,POPP)
|
|
CALL RATMAT(ID,IIEXP,0,ESEMAT,BESE)
|
|
IF(IFPOPR.LE.3) THEN
|
|
CALL MATINV(ESEMAT,NLVEXP,MLEVEL)
|
|
DO I=1,NLVEXP
|
|
SUM=0.
|
|
DO J=1,NLVEXP
|
|
SUM=SUM+ESEMAT(I,J)*BESE(J)
|
|
END DO
|
|
VECL(NSE+I)=POPGRP(I)-SUM
|
|
IF(IGZERO(I,ID).GT.0) VECL(NFREQE+INSE-1+I)=0.
|
|
END DO
|
|
ELSE
|
|
DO I=1,NLVEXP
|
|
SUM=0.
|
|
DO J=1,NLVEXP
|
|
SUM=SUM+ESEMAT(I,J)*POPGRP(J)
|
|
END DO
|
|
VECL(NSE+I)=BESE(I)-SUM
|
|
IF(IGZERO(I,ID).GT.0) VECL(NFREQE+INSE-1+I)=0.
|
|
END DO
|
|
END IF
|
|
end if
|
|
END IF
|
|
C
|
|
C -----------------------------------------------------------
|
|
C 5. charge conservation
|
|
C -----------------------------------------------------------
|
|
C
|
|
IF(INPC.EQ.0) GO TO 400
|
|
NPC=NFREQE+INPC
|
|
C
|
|
C This part is very similar to procedure ELCOR (obviously);
|
|
C array AJ has the meaning of coefficients of the charge conserv.
|
|
C ie. charge conservation is written
|
|
C AJ * (vector of populations) = electron density
|
|
C
|
|
T=TEMP(ID)
|
|
ANE=ELEC(ID)
|
|
CALL STATE(2,ID,T,ANE)
|
|
C
|
|
VPC=QFIX(ID)+Q*DENS(ID)/WMM(ID)/YTOT(ID)
|
|
IF(IOPTAB.EQ.0) VPC=VPC*ABUND(IATREF,ID)
|
|
DO IAT=1,NATOM
|
|
IF(IIFIX(IAT).NE.1) THEN
|
|
DO I=N0A(IAT),NKA(IAT)
|
|
IL=ILK(I)
|
|
CH=IZ(IEL(I))-1
|
|
IF(IL.GT.0) CH=IZ(IL)+(IZ(IL)-1)*ANE*USUMS(IL,ID)
|
|
IF(IMODL(I).GE.0) VPC=VPC+CH*POPUL(I,ID)
|
|
END DO
|
|
END IF
|
|
END DO
|
|
C
|
|
VECL(NPC)=ANE-VPC
|
|
C
|
|
C -----------------------------------------------------------
|
|
C 6. Convection
|
|
C -----------------------------------------------------------
|
|
C
|
|
400 IF(HMIX0.gt.0.) then
|
|
NRE=NFREQE+INRE
|
|
NDEL=NFREQE+INDL
|
|
C
|
|
C Upper boundary condition (ID=1)
|
|
C
|
|
ANEREL=ELEC(1)/(DENS(1)/WMM(1)+ELEC(1))
|
|
IF(ID.EQ.1) THEN
|
|
DELTA(ID)=0.
|
|
FLXC(ID)=0.
|
|
ELSE
|
|
C
|
|
C Normal depth point 1 < ID < ND
|
|
C
|
|
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
|
|
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
|
|
IF(INDL.GT.0) VECL(NDEL)=DELTA(ID)-DLT
|
|
IF(IDISK.EQ.1) GRAVD=ZD(ID)*QGRAV
|
|
C
|
|
C convective flux
|
|
C
|
|
CALL CONVEC(ID,T0,P0,PG0,PR0,AB0,DLT,FLXCNV,VCON)
|
|
FLXC(ID)=FLXCNV
|
|
if(redif(id).gt.0) then
|
|
VECL(NRE)=VECL(NRE)-FLXC(ID)*redif(id)
|
|
end if
|
|
C
|
|
if(reint(id).gt.0.AND.ICONV.LE.2) then
|
|
TP=TEMP(ID+1)
|
|
PM=PTOTAL(ID+1)
|
|
PGP=PGS(ID+1)
|
|
PRADP=PM-PGP-HALF*DENS(ID+1)*VTURB(ID+1)**2
|
|
T0=HALF*(T+TP)
|
|
P0=HALF*(P+PM)
|
|
PG0=HALF*(PG+PGP)
|
|
PR0=HALF*(PRAD+PRADP)
|
|
AB0=HALF*(ABROSD(ID)+ABROSD(ID+1))
|
|
DLT=(TP-T)/(PM-P)*(PM+P)/(TP+T)
|
|
CALL CONVEC(ID+1,T0,P0,PG0,PR0,AB0,DLT,FLXCP,VCON)
|
|
DELM=(DM(ID+1)-DM(ID))*HALF
|
|
RDELM=DENS(ID)/DELM
|
|
VECL(NRE)=VECL(NRE)-RDELM*(FLXCP-FLXCNV)*reint(id)
|
|
end if
|
|
END IF
|
|
END IF
|
|
C
|
|
C skip rows corresponding to fully-zeroed populations
|
|
C
|
|
NSE=NFREQE+INSE
|
|
INONZ=NSE
|
|
DO II=NSE,NN0
|
|
IF(IGZERT(II-NSE+1).EQ.0) THEN
|
|
IF(INONZ.NE.II) VECL(INONZ)=VECL(II)
|
|
INONZ=INONZ+1
|
|
END IF
|
|
END DO
|
|
C
|
|
RETURN
|
|
END
|