1135 lines
33 KiB
Fortran
1135 lines
33 KiB
Fortran
SUBROUTINE ALIFR1(IJ)
|
|
C =====================
|
|
C
|
|
C hydrostatic and radiative equilibrium quantities -
|
|
C derivatives of the total heating and cooling rates in the
|
|
C ALI points with respect to the
|
|
C temperature, electron density, and populations
|
|
C a variant for consistent tridiagonal operator
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
PARAMETER (T23=TWO/3.D0, T43=4.D0/3.D0)
|
|
DIMENSION DSFP1(MLVEXP),DSFP1M(MLVEXP),DSFP1D(MLVEXP),
|
|
* DSFP1P(MLVEXP),DSFPMM(MLVEXP)
|
|
C
|
|
IF(IFALI.LE.1) RETURN
|
|
IF(IFALI.GT.5) THEN
|
|
CALL ALIFR3(IJ)
|
|
RETURN
|
|
END IF
|
|
|
|
WW=WC(IJ)
|
|
C
|
|
DSFT1M=0.
|
|
DSFN1M=0.
|
|
DSFM1M=0.
|
|
DSFT1D=0.
|
|
DSFN1D=0.
|
|
DSFM1D=0.
|
|
DSFX1=0.
|
|
DSFX1D=0.
|
|
DO II=1,NLVEXP
|
|
DSFP1M(II)=0.
|
|
DSFP1D(II)=0.
|
|
END DO
|
|
|
|
c Non-standard value of ILASCT & ILMCOR
|
|
|
|
IF(ILMCOR.NE.3) GO TO 199
|
|
c
|
|
c ****1. Special expressions for the first depth - id=1
|
|
c
|
|
ID=1
|
|
LNSKIP=.NOT.LSKIP(ID,IJ)
|
|
c
|
|
c Basic auxliliary quantities - derivatives of the source function
|
|
c
|
|
if(emis1(id).lt.1.e-35) emis1(id)=1.e-35
|
|
EMISIV=UN/EMIS1(ID)
|
|
ABST=UN/ABSO1(ID)
|
|
S0=EMIS1(ID)*ABST
|
|
c SC=ELEC(ID)*SIGEC(IJ)
|
|
c
|
|
c contribution from the improved boundary condition
|
|
c
|
|
dt=(abso1(id)/dens(id)+abso1(id+1)/dens(id+1))*
|
|
* (dm(id+1)-dm(id))*half
|
|
sa=s0*(un+4./dt*q0(ij))
|
|
s0=s0*(un+two/dt*q0(ij))
|
|
c sa=s0
|
|
sc=scat1(id)
|
|
SCT=SC*ABST
|
|
ST=S0+SCT*RAD1(ID)
|
|
CORR=UN/(UN-ALI1(ID)*SCT)
|
|
DSFT1=CORR*(S0*DEMT1(ID)*EMISIV-S0*DABT1(ID)*ABST)
|
|
DSFN1=CORR*(S0*DEMN1(ID)*EMISIV
|
|
c +SIGEC(IJ)*RAD1(ID)*ABST-
|
|
* -S0*DABN1(ID)*ABST)
|
|
DSFM1=CORR*(S0*DEMM1(ID)*EMISIV-S0*DABM1(ID)*ABST)
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=CORR*(S0*DEMP1(II,ID)*EMISIV-SA*DABP1(II,ID)*ABST)
|
|
END DO
|
|
if(emis1(id+1).lt.1.e-35) emis1(id+1)=1.e-35
|
|
EMISIP=UN/EMIS1(ID+1)
|
|
ABSTP=UN/ABSO1(ID+1)
|
|
S0P=EMIS1(ID+1)*ABSTP
|
|
scp=scat1(id+1)
|
|
SCTP=SCP*ABSTP
|
|
STP=S0P+SCTP*RAD1(ID+1)
|
|
CORRP=UN/(UN-ALI1(ID+1)*SCTP)
|
|
DSFT1P=CORRP*(S0P*DEMT1(ID+1)*EMISIP-S0P*DABT1(ID+1)*ABSTP)
|
|
DSFN1P=CORRP*(S0P*DEMN1(ID+1)*EMISIP
|
|
* -S0P*DABN1(ID+1)*ABSTP)
|
|
DSFM1P=CORRP*(S0P*DEMM1(ID+1)*EMISIP-S0P*DABM1(ID+1)*ABSTP)
|
|
DO II=1,NLVEXP
|
|
DSFP1P(II)=CORRP*(S0P*DEMP1(II,ID+1)*EMISIP-
|
|
* S0P*DABP1(II,ID+1)*ABSTP)
|
|
END DO
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDT(ID)=DSFT1*ALI1(ID)
|
|
DSFDN(ID)=DSFN1*ALI1(ID)
|
|
DSFDM(ID)=DSFM1*ALI1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Hydrostatic equilibrium quantities
|
|
c
|
|
WF=WW*FH(IJ)
|
|
IF(LNSKIP) THEN
|
|
FPRD(ID)=FPRD(ID)+WF*ABSO1(ID)*RAD1(ID)-
|
|
* WW*HEXTRD(IJ)*ABSO1(ID)
|
|
E0=WF*RAD1(ID)
|
|
D0=WF*ABSO1(ID)*ALI1(ID)
|
|
HEIT(ID)=HEIT(ID)+D0*DSFT1+E0*DABT1(ID)
|
|
HEIN(ID)=HEIN(ID)+D0*DSFN1+E0*DABN1(ID)
|
|
HEIM(ID)=HEIM(ID)+D0*DSFM1+E0*DABM1(ID)
|
|
DO II=1,NLVEXP
|
|
HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Differential equation part of radiative equilibrium
|
|
c
|
|
FLFIX(ID)=FLFIX(ID)+WF*RAD1(ID)-WW*HEXTRD(IJ)
|
|
FLRD(ID)=FLRD(ID)+W(IJ)*FH(IJ)*RAD1(ID)-W(IJ)*HALF*EXTRAD(IJ)
|
|
IF(REDIF(ID).GT.0.) THEN
|
|
WF=WF*ALI1(ID)
|
|
REDT(ID)=REDT(ID)+WF*DSFT1
|
|
REDN(ID)=REDN(ID)+WF*DSFN1
|
|
REDM(ID)=REDM(ID)+WF*DSFM1
|
|
REDX(ID)=REDX(ID)+WF*DSFX1
|
|
DO II=1,NLVEXP
|
|
REDP(II,ID)=REDP(II,ID)+WF*DSFP1(II)
|
|
END DO
|
|
END IF
|
|
c
|
|
C Integral equation part of the radiative equilibrium
|
|
C
|
|
IF(REINT(ID).GT.0) THEN
|
|
ABST=ABSO1(ID)-ELSCAT(ID)
|
|
D0=ABST*ALI1(ID)
|
|
FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABST*RAD1(ID))
|
|
REIT(ID)=REIT(ID)+WW*(D0*DSFT1+RAD1(ID)*DABT1(ID)-DEMT1(ID))
|
|
REIN(ID)=REIN(ID)+WW*(D0*DSFN1+
|
|
& RAD1(ID)*(DABN1(ID)-SIGEC(IJ))-DEMN1(ID))
|
|
REIM(ID)=REIM(ID)+WW*(D0*DSFM1+RAD1(ID)*DABM1(ID)-DEMM1(ID))
|
|
DO II=1,NLVEXP
|
|
REIP(II,ID)=REIP(II,ID)+WW*(D0*DSFP1(II)+
|
|
& RAD1(ID)*DABP1(II,ID)-DEMP1(II,ID))
|
|
END DO
|
|
END IF
|
|
|
|
C
|
|
c ****2. loop over depths
|
|
c
|
|
DO ID=2,ND-1
|
|
|
|
LNSKIP=.NOT.LSKIP(ID,IJ)
|
|
|
|
DSFTMM=DSFT1M
|
|
DSFNMM=DSFN1M
|
|
DSFMMM=DSFM1M
|
|
DO II=1,NLVEXP
|
|
DSFPMM(II)=DSFP1M(II)
|
|
END DO
|
|
DSFT1M=DSFT1
|
|
DSFN1M=DSFN1
|
|
DSFM1M=DSFM1
|
|
DO II=1,NLVEXP
|
|
DSFP1M(II)=DSFP1(II)
|
|
END DO
|
|
S0=S0P
|
|
DSFT1=DSFT1P
|
|
DSFN1=DSFN1P
|
|
DSFM1=DSFM1P
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=DSFP1P(II)
|
|
END DO
|
|
c
|
|
if(emis1(id+1).lt.1.e-35) emis1(id+1)=1.e-35
|
|
EMISIP=UN/EMIS1(ID+1)
|
|
ABSTP=UN/ABSO1(ID+1)
|
|
S0P=EMIS1(ID+1)*ABSTP
|
|
scp=scat1(id+1)
|
|
SCTP=SCP*ABSTP
|
|
STP=S0P+SCTP*RAD1(ID+1)
|
|
CORRP=UN/(UN-ALI1(ID+1)*SCTP)
|
|
DSFT1P=CORRP*(S0P*DEMT1(ID+1)*EMISIP-S0P*DABT1(ID+1)*ABSTP)
|
|
DSFN1P=CORRP*(S0P*DEMN1(ID+1)*EMISIP
|
|
* -S0P*DABN1(ID+1)*ABSTP)
|
|
DSFM1P=CORRP*(S0P*DEMM1(ID+1)*EMISIP-STP*DABM1(ID+1)*ABSTP)
|
|
DO II=1,NLVEXP
|
|
DSFP1P(II)=CORRP*(S0P*DEMP1(II,ID+1)*EMISIP-
|
|
* S0P*DABP1(II,ID+1)*ABSTP)
|
|
END DO
|
|
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDT(ID)=DSFT1*ALI1(ID)
|
|
DSFDN(ID)=DSFN1*ALI1(ID)
|
|
DSFDM(ID)=DSFM1*ALI1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Hydrostatic equilibrium equation
|
|
c
|
|
IF(LNSKIP) THEN
|
|
D0=WW*FAK1(ID)
|
|
A0=WW*FAK1(ID-1)
|
|
FPRD(ID)=FPRD(ID)+D0*RAD1(ID)-A0*RAD1(ID-1)
|
|
E0=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
HEIT(ID)=HEIT(ID)+D0*DSFT1
|
|
HEIN(ID)=HEIN(ID)+D0*DSFN1
|
|
HEIM(ID)=HEIM(ID)+D0*DSFM1
|
|
HEITM(ID)=HEITM(ID)+E0*DSFT1M
|
|
HEINM(ID)=HEINM(ID)+E0*DSFN1M
|
|
HEIMM(ID)=HEIMM(ID)+E0*DSFM1M
|
|
DO II=1,NLVEXP
|
|
HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)
|
|
HEIPM(II,ID)=HEIPM(II,ID)+E0*DSFP1M(II)
|
|
END DO
|
|
END IF
|
|
C
|
|
C Differential equation part of radiative equilibrium
|
|
C
|
|
DDT=UN/(ABSOT(ID)+ABSOT(ID-1))
|
|
DT=DDT/DELDMZ(ID-1)
|
|
FL=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))*DT
|
|
FLFIX(ID)=FLFIX(ID)+WW*FL
|
|
FLRD(ID)=FLRD(ID)+W(IJ)*FL
|
|
IF(REDIF(ID).GT.0) THEN
|
|
if(ifalih.eq.0) then
|
|
D0=WW*FAK1(ID)*DT
|
|
A0=WW*FAK1(ID-1)*DT
|
|
D0M=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0P=D0*ALIP1(ID)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
E0=WW*FL*DDT
|
|
REDX(ID)=REDX(ID)+E0*ABSO1(ID)
|
|
REDXM(ID)=REDXM(ID)+E0*ABSO1(ID-1)
|
|
E0M=E0*DENSI(ID-1)
|
|
E0=E0*DENSI(ID)
|
|
REDT(ID)=REDT(ID)+D0*DSFT1-E0*DABT1(ID)
|
|
REDTM(ID)=REDTM(ID)+D0M*DSFT1M-E0M*DABT1(ID-1)
|
|
REDN(ID)=REDN(ID)+D0*DSFN1-E0*DABN1(ID)
|
|
REDNM(ID)=REDNM(ID)+D0M*DSFN1M-E0M*DABN1(ID-1)
|
|
REDM(ID)=REDM(ID)+D0*DSFM1-E0*DABM1(ID)
|
|
REDMM(ID)=REDMM(ID)+D0M*DSFM1M-E0M*DABM1(ID-1)
|
|
DO II=1,NLVEXP
|
|
REDP(II,ID)=REDP(II,ID)+D0*DSFP1(II)-E0*DABP1(II,ID)
|
|
REDPM(II,ID)=REDPM(II,ID)+D0M*DSFP1M(II)-
|
|
* E0M*DABP1(II,ID-1)
|
|
END DO
|
|
else
|
|
d0=ww*alih1(id)
|
|
REDT(ID)=REDT(ID)+D0*DSFT1
|
|
REDN(ID)=REDN(ID)+D0*DSFN1
|
|
REDM(ID)=REDM(ID)+D0*DSFM1
|
|
DO II=1,NLVEXP
|
|
REDP(II,ID)=REDP(II,ID)+D0*DSFP1(II)
|
|
END DO
|
|
end if
|
|
END IF
|
|
c
|
|
C Integral equation part of the radiative equilibrium
|
|
C
|
|
IF(REINT(ID).GT.0) THEN
|
|
ABST=ABSO1(ID)-ELSCAT(ID)
|
|
D0=ABST*ALI1(ID)
|
|
FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABST*RAD1(ID))
|
|
REIN(ID)=REIN(ID)+WW*(D0*DSFN1+
|
|
* RAD1(ID)*(DABN1(ID)-SIGEC(IJ))-DEMN1(ID))
|
|
DO II=1,NLVEXP
|
|
REIP(II,ID)=REIP(II,ID)+WW*(D0*DSFP1(II)+
|
|
* RAD1(ID)*DABP1(II,ID)-DEMP1(II,ID))
|
|
END DO
|
|
REIT(ID)=REIT(ID)+WW*(D0*DSFT1+
|
|
* RAD1(ID)*DABT1(ID)-DEMT1(ID))
|
|
REIM(ID)=REIM(ID)+WW*(D0*DSFM1+
|
|
* RAD1(ID)*DABM1(ID)-DEMM1(ID))
|
|
END IF
|
|
|
|
END DO
|
|
|
|
C
|
|
c ****3. deepest point - ID=ND
|
|
c
|
|
ID=ND
|
|
|
|
LNSKIP=.NOT.LSKIP(ID,IJ)
|
|
|
|
DSFTMM=DSFT1M
|
|
DSFNMM=DSFN1M
|
|
DSFMMM=DSFM1M
|
|
DO II=1,NLVEXP
|
|
DSFPMM(II)=DSFP1M(II)
|
|
END DO
|
|
DSFT1M=DSFT1
|
|
DSFN1M=DSFN1
|
|
DSFM1M=DSFM1
|
|
DO II=1,NLVEXP
|
|
DSFP1M(II)=DSFP1(II)
|
|
END DO
|
|
S0=S0P
|
|
DSFT1=DSFT1P
|
|
DSFN1=DSFN1P
|
|
DSFM1=DSFM1P
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=DSFP1P(II)
|
|
END DO
|
|
DSFT1D=0.
|
|
DSFN1D=0.
|
|
DSFM1D=0.
|
|
DO II=1,NLVEXP
|
|
DSFP1D(II)=0.
|
|
END DO
|
|
C
|
|
C Improved lower boundary condition
|
|
C
|
|
IF(IBC.GT.0.AND.IDISK.EQ.0) THEN
|
|
DT=UN/(DELDMZ(ID-1)*(ABSOT(ID)+ABSOT(ID-1)))
|
|
PLAD=XKFB(ID)/XKF1(ID)
|
|
DBDT=PLAD/XKF1(ID)*HKT21(ID)*FREQ(IJ)*DT
|
|
IF(IBC.EQ.1) THEN
|
|
DSFT1=DSFT1+DBDT
|
|
ELSE IF(IBC.GE.2) THEN
|
|
PLAM=XKFB(ID-1)/XKF1(ID-1)
|
|
TAU23=T23*DT
|
|
TAU43=T43*DT
|
|
D0=(PLAD*(UN+TAU43)-T43*PLAM*DT)*DT*DT
|
|
RHD=DELDMZ(ID-1)*DENSI(ID)
|
|
E0=D0*RHD
|
|
DSFT1=DSFT1+DBDT*(UN+TAU23)-E0*DABT1(ID)
|
|
DSFN1=DSFN1-E0*(DABN1(ID)+ABSO1(ID)*DENSIM(ID))
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=DSFP1(II)-E0*DABP1(II,ID)
|
|
END DO
|
|
IF(IBC.GE.3) THEN
|
|
DBDTM=PLAM/XKF1(ID-1)*HKT21(ID-1)*FREQ(IJ)*DT
|
|
RHD=DELDMZ(ID-1)*DENSI(ID-1)
|
|
E0=D0*RHD
|
|
DSFT1D=-DBDTM*DT*T23-E0*DABT1(ID-1)
|
|
DSFN1D=-E0*(DABN1(ID-1)+ABSO1(ID-1)*DENSIM(ID-1))
|
|
DO II=1,NLVEXP
|
|
DSFP1D(II)=-E0*DABP1(II,ID-1)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
END IF
|
|
C
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDT(ID)=DSFT1*ALI1(ID)
|
|
DSFDN(ID)=DSFN1*ALI1(ID)
|
|
DSFDM(ID)=DSFM1*ALI1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Hydrostatic equilibrium equation
|
|
c
|
|
IF(LNSKIP) THEN
|
|
D0=WW*FAK1(ID)
|
|
A0=WW*FAK1(ID-1)
|
|
FPRD(ID)=FPRD(ID)+D0*RAD1(ID)-A0*RAD1(ID-1)
|
|
E0=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
HEIT(ID)=HEIT(ID)+D0*DSFT1
|
|
HEIN(ID)=HEIN(ID)+D0*DSFN1
|
|
HEIM(ID)=HEIM(ID)+D0*DSFM1
|
|
HEITM(ID)=HEITM(ID)+E0*DSFT1M
|
|
HEINM(ID)=HEINM(ID)+E0*DSFN1M
|
|
HEIMM(ID)=HEIMM(ID)+E0*DSFM1M
|
|
DO II=1,NLVEXP
|
|
HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)
|
|
HEIPM(II,ID)=HEIPM(II,ID)+E0*DSFP1M(II)
|
|
END DO
|
|
IF(IBC.GE.3) THEN
|
|
HEITM(ID)=HEITM(ID)-D0*DSFT1D
|
|
HEINM(ID)=HEINM(ID)-D0*DSFN1D
|
|
HEIMM(ID)=HEIMM(ID)-D0*DSFM1D
|
|
DO II=1,NLVEXP
|
|
HEIPM(II,ID)=HEIPM(II,ID)-D0*DSFP1D(II)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
C
|
|
C Differential equation part of radiative equilibrium
|
|
C
|
|
DDT=UN/(ABSOT(ID)+ABSOT(ID-1))
|
|
DT=DDT/DELDMZ(ID-1)
|
|
FL=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))*DT
|
|
FLFIX(ID)=FLFIX(ID)+WW*FL
|
|
FLRD(ID)=FLRD(ID)+W(IJ)*FL
|
|
IF(REDIF(ID).GT.0) THEN
|
|
D0=WW*FAK1(ID)*DT
|
|
A0=WW*FAK1(ID-1)*DT
|
|
D0M=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0P=D0*ALIP1(ID)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
E0=WW*FL*DDT
|
|
REDX(ID)=REDX(ID)+E0*ABSO1(ID)
|
|
REDXM(ID)=REDXM(ID)+E0*ABSO1(ID-1)
|
|
E0M=E0*DENSI(ID-1)
|
|
E0=E0*DENSI(ID)
|
|
REDT(ID)=REDT(ID)+D0*DSFT1-E0*DABT1(ID)
|
|
REDTM(ID)=REDTM(ID)+D0M*DSFT1M-E0M*DABT1(ID-1)
|
|
REDN(ID)=REDN(ID)+D0*DSFN1-E0*DABN1(ID)
|
|
REDNM(ID)=REDNM(ID)+D0M*DSFN1M-E0M*DABN1(ID-1)
|
|
REDM(ID)=REDM(ID)+D0*DSFM1-E0*DABM1(ID)
|
|
REDMM(ID)=REDMM(ID)+D0M*DSFM1M-E0M*DABM1(ID-1)
|
|
DO II=1,NLVEXP
|
|
REDP(II,ID)=REDP(II,ID)+D0*DSFP1(II)-E0*DABP1(II,ID)
|
|
REDPM(II,ID)=REDPM(II,ID)+D0M*DSFP1M(II)-
|
|
* E0M*DABP1(II,ID-1)
|
|
END DO
|
|
IF(IBC.GE.3) THEN
|
|
REDTM(ID)=REDTM(ID)+D0*DSFT1D
|
|
REDNM(ID)=REDNM(ID)+D0*DSFN1D
|
|
REDMM(ID)=REDMM(ID)+D0*DSFM1D
|
|
DO II=1,NLVEXP
|
|
REDPM(II,ID)=REDPM(II,ID)+D0*DSFP1D(II)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
c
|
|
C Integral equation part of the radiative equilibrium
|
|
C
|
|
IF(REINT(ID).GT.0) THEN
|
|
ABST=ABSO1(ID)-ELSCAT(ID)
|
|
D0=ABST*ALI1(ID)
|
|
FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABST*RAD1(ID))
|
|
REIN(ID)=REIN(ID)+WW*(D0*DSFN1+
|
|
* RAD1(ID)*(DABN1(ID)-SIGEC(IJ))-DEMN1(ID))
|
|
REIM(ID)=REIM(ID)+WW*(D0*DSFM1+RAD1(ID)*DABM1(ID)-DEMM1(ID))
|
|
DO II=1,NLVEXP
|
|
REIP(II,ID)=REIP(II,ID)+WW*(D0*DSFP1(II)+
|
|
* RAD1(ID)*DABP1(II,ID)-DEMP1(II,ID))
|
|
END DO
|
|
IF(IBC.EQ.0) THEN
|
|
REIT(ID)=REIT(ID)+WW*(D0*DSFT1+
|
|
* RAD1(ID)*DABT1(ID)-DEMT1(ID))
|
|
ELSE
|
|
REIT(ID)=REIT(ID)+WW*(D0*(DSFT1-DBDT)+E0*DABT1(ID)+
|
|
* RAD1(ID)*DABT1(ID)-DEMT1(ID)+
|
|
* ALI1(ID)/ABST*DBDT)
|
|
END IF
|
|
END IF
|
|
|
|
RETURN
|
|
|
|
c Non-standard value of ILASCT & ILMCOR
|
|
|
|
199 IF(ILASCT.NE.0) GO TO 299
|
|
c
|
|
c ****1. Special expressions for the first depth - id=1
|
|
c
|
|
ID=1
|
|
|
|
LNSKIP=.NOT.LSKIP(ID,IJ)
|
|
c
|
|
c Basic auxliliary quantities - derivatives of the source function
|
|
c
|
|
EMISIV=UN/EMIS1(ID)
|
|
ABST=UN/(ABSO1(ID)-ELSCAT(ID))
|
|
S0=EMIS1(ID)*ABST
|
|
DSFN1=S0*(DEMN1(ID)*EMISIV-(DABN1(ID)-SIGEC(IJ))*ABST)
|
|
DSFT1=S0*(DEMT1(ID)*EMISIV-DABT1(ID)*ABST)
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=S0*(DEMP1(II,ID)*EMISIV-DABP1(II,ID)*ABST)
|
|
END DO
|
|
EMISIP=UN/EMIS1(ID+1)
|
|
ABSTP=UN/(ABSO1(ID+1)-ELSCAT(ID+1))
|
|
S0P=EMIS1(ID+1)*ABSTP
|
|
DSFN1P=S0P*(DEMN1(ID+1)*EMISIP-(DABN1(ID+1)-SIGEC(IJ))*ABSTP)
|
|
DSFT1P=S0P*(DEMT1(ID+1)*EMISIP-DABT1(ID+1)*ABSTP)
|
|
DO II=1,NLVEXP
|
|
DSFP1P(II)=S0P*(DEMP1(II,ID+1)*EMISIP-DABP1(II,ID+1)*ABSTP)
|
|
END DO
|
|
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDT(ID)=DSFT1*ALI1(ID)
|
|
DSFDN(ID)=DSFN1*ALI1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Hydrostatic equilibrium quantities
|
|
c
|
|
WF=WW*FH(IJ)
|
|
IF(LNSKIP) THEN
|
|
FPRD(ID)=FPRD(ID)+WF*ABSO1(ID)*RAD1(ID)-
|
|
* WW*HEXTRD(IJ)*ABSO1(ID)
|
|
E0=WF*RAD1(ID)
|
|
D0=WF*ABSO1(ID)*ALI1(ID)
|
|
HEIT(ID)=HEIT(ID)+D0*DSFT1+E0*DABT1(ID)
|
|
HEIN(ID)=HEIN(ID)+D0*DSFN1+E0*DABN1(ID)
|
|
DO II=1,NLVEXP
|
|
HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Differential equation part of radiative equilibrium
|
|
c
|
|
FLFIX(ID)=FLFIX(ID)+WF*RAD1(ID)-WW*HEXTRD(IJ)
|
|
FLRD(ID)=FLRD(ID)+W(IJ)*FH(IJ)*RAD1(ID)-W(IJ)*HALF*EXTRAD(IJ)
|
|
IF(REDIF(ID).GT.0.) THEN
|
|
WF=WF*ALI1(ID)
|
|
REDT(ID)=REDT(ID)+WF*DSFT1
|
|
REDN(ID)=REDN(ID)+WF*DSFN1
|
|
DO II=1,NLVEXP
|
|
REDP(II,ID)=REDP(II,ID)+WF*DSFP1(II)
|
|
END DO
|
|
END IF
|
|
c
|
|
C Integral equation part of the radiative equilibrium
|
|
C
|
|
IF(REINT(ID).GT.0) THEN
|
|
ABST=ABSO1(ID)-ELSCAT(ID)
|
|
WWK=WW*ABST
|
|
FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABST*RAD1(ID))
|
|
D0=WW*(ALI1(ID)-UN)*ABST
|
|
E0=WW*(RAD1(ID)-S0)
|
|
REIN(ID)=REIN(ID)+D0*DSFN1+E0*(DABN1(ID)-SIGEC(IJ))
|
|
REIT(ID)=REIT(ID)+D0*DSFT1+E0*DABT1(ID)
|
|
DO II=1,NLVEXP
|
|
REIP(II,ID)=REIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
|
|
END DO
|
|
END IF
|
|
|
|
C
|
|
c ****2. loop over depths
|
|
c
|
|
DO ID=2,ND-1
|
|
|
|
LNSKIP=.NOT.LSKIP(ID,IJ)
|
|
|
|
DSFTMM=DSFT1M
|
|
DSFNMM=DSFN1M
|
|
DO II=1,NLVEXP
|
|
DSFPMM(II)=DSFP1M(II)
|
|
END DO
|
|
DSFT1M=DSFT1
|
|
DSFN1M=DSFN1
|
|
DO II=1,NLVEXP
|
|
DSFP1M(II)=DSFP1(II)
|
|
END DO
|
|
S0=S0P
|
|
DSFT1=DSFT1P
|
|
DSFN1=DSFN1P
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=DSFP1P(II)
|
|
END DO
|
|
|
|
EMISIP=UN/EMIS1(ID+1)
|
|
ABSTP=UN/(ABSO1(ID+1)-ELSCAT(ID+1))
|
|
S0P=EMIS1(ID+1)*ABSTP
|
|
DSFN1P=S0P*(DEMN1(ID+1)*EMISIP-(DABN1(ID+1)-SIGEC(IJ))*ABSTP)
|
|
DSFT1P=S0P*(DEMT1(ID+1)*EMISIP-DABT1(ID+1)*ABSTP)
|
|
DO II=1,NLVEXP
|
|
DSFP1P(II)=S0P*(DEMP1(II,ID+1)*EMISIP-DABP1(II,ID+1)*ABSTP)
|
|
END DO
|
|
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDT(ID)=DSFT1*ALI1(ID)
|
|
DSFDN(ID)=DSFN1*ALI1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Hydrostatic equilibrium equation
|
|
c
|
|
IF(LNSKIP) THEN
|
|
D0=WW*FAK1(ID)
|
|
A0=WW*FAK1(ID-1)
|
|
FPRD(ID)=FPRD(ID)+D0*RAD1(ID)-A0*RAD1(ID-1)
|
|
E0=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
HEIT(ID)=HEIT(ID)+D0*DSFT1
|
|
HEIN(ID)=HEIN(ID)+D0*DSFN1
|
|
HEITM(ID)=HEITM(ID)+E0*DSFT1M
|
|
HEINM(ID)=HEINM(ID)+E0*DSFN1M
|
|
DO II=1,NLVEXP
|
|
HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)
|
|
HEIPM(II,ID)=HEIPM(II,ID)+E0*DSFP1M(II)
|
|
END DO
|
|
END IF
|
|
C
|
|
C Differential equation part of radiative equilibrium
|
|
C
|
|
DDT=UN/(ABSOT(ID)+ABSOT(ID-1))
|
|
DT=DDT/DELDMZ(ID-1)
|
|
FL=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))*DT
|
|
FLFIX(ID)=FLFIX(ID)+WW*FL
|
|
IF(REDIF(ID).GT.0) THEN
|
|
D0=WW*FAK1(ID)*DT
|
|
A0=WW*FAK1(ID-1)*DT
|
|
D0M=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0P=D0*ALIP1(ID)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
E0=WW*FL*DDT
|
|
REDX(ID)=REDX(ID)+E0*ABSO1(ID)
|
|
REDXM(ID)=REDXM(ID)+E0*ABSO1(ID-1)
|
|
E0M=E0*DENSI(ID-1)
|
|
E0=E0*DENSI(ID)
|
|
REDT(ID)=REDT(ID)+D0*DSFT1-E0*DABT1(ID)
|
|
REDTM(ID)=REDTM(ID)+D0M*DSFT1M-E0M*DABT1(ID-1)
|
|
REDN(ID)=REDN(ID)+D0*DSFN1-E0*DABN1(ID)
|
|
REDNM(ID)=REDNM(ID)+D0M*DSFN1M-E0M*DABN1(ID-1)
|
|
DO II=1,NLVEXP
|
|
REDP(II,ID)=REDP(II,ID)+D0*DSFP1(II)-E0*DABP1(II,ID)
|
|
REDPM(II,ID)=REDPM(II,ID)+D0M*DSFP1M(II)-
|
|
* E0M*DABP1(II,ID-1)
|
|
END DO
|
|
END IF
|
|
c
|
|
C Integral equation part of the radiative equilibrium
|
|
C
|
|
IF(REINT(ID).GT.0) THEN
|
|
ABST=ABSO1(ID)-ELSCAT(ID)
|
|
WWK=WW*ABST
|
|
FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABST*RAD1(ID))
|
|
D0=WW*(ALI1(ID)-UN)*ABST
|
|
E0=WW*(RAD1(ID)-S0)
|
|
REIN(ID)=REIN(ID)+D0*DSFN1+E0*(DABN1(ID)-SIGEC(IJ))
|
|
REIT(ID)=REIT(ID)+D0*DSFT1+E0*DABT1(ID)
|
|
DO II=1,NLVEXP
|
|
REIP(II,ID)=REIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
|
|
END DO
|
|
END IF
|
|
|
|
END DO
|
|
|
|
C
|
|
c ****3. deepest point - ID=ND
|
|
c
|
|
ID=ND
|
|
|
|
LNSKIP=.NOT.LSKIP(ID,IJ)
|
|
|
|
DSFTMM=DSFT1M
|
|
DSFNMM=DSFN1M
|
|
DO II=1,NLVEXP
|
|
DSFPMM(II)=DSFP1M(II)
|
|
END DO
|
|
DSFT1M=DSFT1
|
|
DSFN1M=DSFN1
|
|
DO II=1,NLVEXP
|
|
DSFP1M(II)=DSFP1(II)
|
|
END DO
|
|
S0=S0P
|
|
DSFT1=DSFT1P
|
|
DSFN1=DSFN1P
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=DSFP1P(II)
|
|
END DO
|
|
C
|
|
C Improved lower boundary condition
|
|
C
|
|
IF(IBC.GT.0.AND.IDISK.EQ.0) THEN
|
|
DT=UN/(DELDMZ(ID-1)*(ABSOT(ID)+ABSOT(ID-1)))
|
|
PLAD=XKFB(ID)/XKF1(ID)
|
|
DBDT=PLAD/XKF1(ID)*HKT21(ID)*FREQ(IJ)*DT
|
|
IF(IBC.EQ.1) THEN
|
|
DSFT1=DSFT1+DBDT
|
|
ELSE IF(IBC.GE.2) THEN
|
|
PLAM=XKFB(ID-1)/XKF1(ID-1)
|
|
TAU23=T23*DT
|
|
TAU43=T43*DT
|
|
D0=(PLAD*(UN+TAU43)-T43*PLAM*DT)*DT*DT
|
|
RHD=DELDMZ(ID-1)*DENSI(ID)
|
|
E0=D0*RHD
|
|
DSFT1=DSFT1+DBDT*(UN+TAU23)-E0*DABT1(ID)
|
|
DSFN1=DSFN1-E0*(DABN1(ID)+ABSO1(ID)*DENSIM(ID))
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=DSFP1(II)-E0*DABP1(II,ID)
|
|
END DO
|
|
IF(IBC.GE.3) THEN
|
|
DBDTM=PLAM/XKF1(ID-1)*HKT21(ID-1)*FREQ(IJ)*DT
|
|
RHD=DELDMZ(ID-1)*DENSI(ID-1)
|
|
E0=D0*RHD
|
|
DSFT1D=-DBDTM*DT*T23-E0*DABT1(ID-1)
|
|
DSFN1D=-E0*(DABN1(ID-1)+ABSO1(ID-1)*DENSIM(ID-1))
|
|
DO II=1,NLVEXP
|
|
DSFP1D(II)=-E0*DABP1(II,ID-1)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
END IF
|
|
C
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDT(ID)=DSFT1*ALI1(ID)
|
|
DSFDN(ID)=DSFN1*ALI1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Hydrostatic equilibrium equation
|
|
c
|
|
IF(LNSKIP) THEN
|
|
D0=WW*FAK1(ID)
|
|
A0=WW*FAK1(ID-1)
|
|
FPRD(ID)=FPRD(ID)+D0*RAD1(ID)-A0*RAD1(ID-1)
|
|
E0=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
HEIT(ID)=HEIT(ID)+D0*DSFT1
|
|
HEIN(ID)=HEIN(ID)+D0*DSFN1
|
|
HEITM(ID)=HEITM(ID)+E0*DSFT1M
|
|
HEINM(ID)=HEINM(ID)+E0*DSFN1M
|
|
DO II=1,NLVEXP
|
|
HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)
|
|
HEIPM(II,ID)=HEIPM(II,ID)+E0*DSFP1M(II)
|
|
END DO
|
|
IF(IBC.GE.3) THEN
|
|
HEITM(ID)=HEITM(ID)-D0*DSFT1D
|
|
HEINM(ID)=HEINM(ID)-D0*DSFN1D
|
|
DO II=1,NLVEXP
|
|
HEIPM(II,ID)=HEIPM(II,ID)-D0*DSFP1D(II)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
C
|
|
C Differential equation part of radiative equilibrium
|
|
C
|
|
DDT=UN/(ABSOT(ID)+ABSOT(ID-1))
|
|
DT=DDT/DELDMZ(ID-1)
|
|
FL=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))*DT
|
|
FLFIX(ID)=FLFIX(ID)+WW*FL
|
|
IF(REDIF(ID).GT.0) THEN
|
|
D0=WW*FAK1(ID)*DT
|
|
A0=WW*FAK1(ID-1)*DT
|
|
D0M=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0P=D0*ALIP1(ID)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
E0=WW*FL*DDT
|
|
REDX(ID)=REDX(ID)+E0*ABSO1(ID)
|
|
REDXM(ID)=REDXM(ID)+E0*ABSO1(ID-1)
|
|
E0M=E0*DENSI(ID-1)
|
|
E0=E0*DENSI(ID)
|
|
REDT(ID)=REDT(ID)+D0*DSFT1-E0*DABT1(ID)
|
|
REDTM(ID)=REDTM(ID)+D0M*DSFT1M-E0M*DABT1(ID-1)
|
|
REDN(ID)=REDN(ID)+D0*DSFN1-E0*DABN1(ID)
|
|
REDNM(ID)=REDNM(ID)+D0M*DSFN1M-E0M*DABN1(ID-1)
|
|
DO II=1,NLVEXP
|
|
REDP(II,ID)=REDP(II,ID)+D0*DSFP1(II)-E0*DABP1(II,ID)
|
|
REDPM(II,ID)=REDPM(II,ID)+D0M*DSFP1M(II)-
|
|
* E0M*DABP1(II,ID-1)
|
|
END DO
|
|
IF(IBC.GE.3) THEN
|
|
REDTM(ID)=REDTM(ID)+D0*DSFT1D
|
|
REDNM(ID)=REDNM(ID)+D0*DSFN1D
|
|
DO II=1,NLVEXP
|
|
REDPM(II,ID)=REDPM(II,ID)+D0*DSFP1D(II)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
c
|
|
C Integral equation part of the radiative equilibrium
|
|
C
|
|
IF(REINT(ID).GT.0) THEN
|
|
ABST=ABSO1(ID)-ELSCAT(ID)
|
|
WWK=WW*ABST
|
|
FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABST*RAD1(ID))
|
|
D0=WW*(ALI1(ID)-UN)*ABST
|
|
E0=WW*(RAD1(ID)-S0)
|
|
REIN(ID)=REIN(ID)+D0*DSFN1+E0*(DABN1(ID)-SIGEC(IJ))
|
|
IF(IBC.EQ.0) THEN
|
|
REIT(ID)=REIT(ID)+D0*DSFT1+E0*DABT1(ID)
|
|
ELSE
|
|
REIT(ID)=REIT(ID)+D0*(DSFT1-DBDT)+E0*DABT1(ID)+
|
|
* ALI1(ID)/ABST*DBDT
|
|
END IF
|
|
DO II=1,NLVEXP
|
|
REIP(II,ID)=REIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
|
|
END DO
|
|
END IF
|
|
|
|
RETURN
|
|
|
|
C Non-standard value of ILASCT (ILASCT != 0)
|
|
c
|
|
c ****1. Special expressions for the first depth - id=1
|
|
c
|
|
299 ID=1
|
|
|
|
LNSKIP=.NOT.LSKIP(ID,IJ)
|
|
c
|
|
c Basic auxliliary quantities - derivatives of the source function
|
|
c
|
|
EMISIV=UN/EMIS1(ID)
|
|
ABST=UN/ABSO1(ID)
|
|
S0=EMIS1(ID)*ABST
|
|
DSFN1=S0*(DEMN1(ID)*EMISIV-DABN1(ID)*ABST)
|
|
DSFT1=S0*(DEMT1(ID)*EMISIV-DABT1(ID)*ABST)
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=S0*(DEMP1(II,ID)*EMISIV-DABP1(II,ID)*ABST)
|
|
END DO
|
|
EMISIP=UN/EMIS1(ID+1)
|
|
ABSTP=UN/ABSO1(ID+1)
|
|
S0P=EMIS1(ID+1)*ABSTP
|
|
DSFN1P=S0P*(DEMN1(ID+1)*EMISIP-DABN1(ID+1)*ABSTP)
|
|
DSFT1P=S0P*(DEMT1(ID+1)*EMISIP-DABT1(ID+1)*ABSTP)
|
|
DO II=1,NLVEXP
|
|
DSFP1P(II)=S0P*(DEMP1(II,ID+1)*EMISIP-DABP1(II,ID+1)*ABSTP)
|
|
END DO
|
|
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDT(ID)=DSFT1*ALI1(ID)
|
|
DSFDN(ID)=DSFN1*ALI1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
|
|
END DO
|
|
END IF
|
|
|
|
c
|
|
c Hydrostatic equilibrium quantities
|
|
c
|
|
WF=WW*FH(IJ)
|
|
IF(LNSKIP) THEN
|
|
FPRD(ID)=FPRD(ID)+WF*ABSO1(ID)*RAD1(ID)-
|
|
* WW*HEXTRD(IJ)*ABSO1(ID)
|
|
E0=WF*RAD1(ID)
|
|
D0=WF*ABSO1(ID)*ALI1(ID)
|
|
HEIT(ID)=HEIT(ID)+D0*DSFT1+E0*DABT1(ID)
|
|
HEIN(ID)=HEIN(ID)+D0*DSFN1+E0*DABN1(ID)
|
|
DO II=1,NLVEXP
|
|
HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Differential equation part of radiative equilibrium
|
|
c
|
|
FLFIX(ID)=FLFIX(ID)+WF*RAD1(ID)-WW*HEXTRD(IJ)
|
|
FLRD(ID)=FLRD(ID)+W(IJ)*FH(IJ)*RAD1(ID)-W(IJ)*HALF*EXTRAD(IJ)
|
|
IF(REDIF(ID).GT.0.) THEN
|
|
WF=WF*ALI1(ID)
|
|
REDT(ID)=REDT(ID)+WF*DSFT1
|
|
REDN(ID)=REDN(ID)+WF*DSFN1
|
|
DO II=1,NLVEXP
|
|
REDP(II,ID)=REDP(II,ID)+WF*DSFP1(II)
|
|
END DO
|
|
END IF
|
|
c
|
|
C Integral equation part of the radiative equilibrium
|
|
C
|
|
IF(REINT(ID).GT.0) THEN
|
|
SRH=SIGE*DENS1(ID)
|
|
ABST=ABSO1(ID)
|
|
ABSTE=ABST-ELSCAT(ID)
|
|
WWK=WW*ABST
|
|
FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABSTE*RAD1(ID))
|
|
D0=WW*(ABSTE*ALI1(ID)-ABST)
|
|
E0=WW*(RAD1(ID)-S0)
|
|
REIN(ID)=REIN(ID)+D0*DSFN1+E0*DABN1(ID)-WW*SIGEC(IJ)*RAD1(ID)
|
|
REIT(ID)=REIT(ID)+D0*DSFT1+E0*DABT1(ID)
|
|
DO II=1,NLVEXP
|
|
REIP(II,ID)=REIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
|
|
END DO
|
|
END IF
|
|
|
|
C
|
|
c ****2. loop over depths
|
|
c
|
|
DO ID=2,ND-1
|
|
|
|
LNSKIP=.NOT.LSKIP(ID,IJ)
|
|
|
|
DSFTMM=DSFT1M
|
|
DSFNMM=DSFN1M
|
|
DO II=1,NLVEXP
|
|
DSFPMM(II)=DSFP1M(II)
|
|
END DO
|
|
DSFT1M=DSFT1
|
|
DSFN1M=DSFN1
|
|
DO II=1,NLVEXP
|
|
DSFP1M(II)=DSFP1(II)
|
|
END DO
|
|
S0=S0P
|
|
DSFT1=DSFT1P
|
|
DSFN1=DSFN1P
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=DSFP1P(II)
|
|
END DO
|
|
|
|
EMISIP=UN/EMIS1(ID+1)
|
|
ABSTP=UN/ABSO1(ID+1)
|
|
S0P=EMIS1(ID+1)*ABSTP
|
|
DSFN1P=S0P*(DEMN1(ID+1)*EMISIP-DABN1(ID+1)*ABSTP)
|
|
DSFT1P=S0P*(DEMT1(ID+1)*EMISIP-DABT1(ID+1)*ABSTP)
|
|
DO II=1,NLVEXP
|
|
DSFP1P(II)=S0P*(DEMP1(II,ID+1)*EMISIP-DABP1(II,ID+1)*ABSTP)
|
|
END DO
|
|
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDT(ID)=DSFT1*ALI1(ID)
|
|
DSFDN(ID)=DSFN1*ALI1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Hydrostatic equilibrium equation
|
|
c
|
|
IF(LNSKIP) THEN
|
|
D0=WW*FAK1(ID)
|
|
A0=WW*FAK1(ID-1)
|
|
FPRD(ID)=FPRD(ID)+D0*RAD1(ID)-A0*RAD1(ID-1)
|
|
E0=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
HEIT(ID)=HEIT(ID)+D0*DSFT1
|
|
HEIN(ID)=HEIN(ID)+D0*DSFN1
|
|
HEITM(ID)=HEITM(ID)+E0*DSFT1M
|
|
HEINM(ID)=HEINM(ID)+E0*DSFN1M
|
|
DO II=1,NLVEXP
|
|
HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)
|
|
HEIPM(II,ID)=HEIPM(II,ID)+E0*DSFP1M(II)
|
|
END DO
|
|
END IF
|
|
C
|
|
C Differential equation part of radiative equilibrium
|
|
C
|
|
DDT=UN/(ABSOT(ID)+ABSOT(ID-1))
|
|
DT=DDT/DELDMZ(ID-1)
|
|
FL=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))*DT
|
|
FLFIX(ID)=FLFIX(ID)+WW*FL
|
|
IF(REDIF(ID).GT.0) THEN
|
|
D0=WW*FAK1(ID)*DT
|
|
A0=WW*FAK1(ID-1)*DT
|
|
D0M=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0P=D0*ALIP1(ID)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
E0=WW*FL*DDT
|
|
REDX(ID)=REDX(ID)+E0*ABSO1(ID)
|
|
REDXM(ID)=REDXM(ID)+E0*ABSO1(ID-1)
|
|
E0M=E0*DENSI(ID-1)
|
|
E0=E0*DENSI(ID)
|
|
REDT(ID)=REDT(ID)+D0*DSFT1-E0*DABT1(ID)
|
|
REDTM(ID)=REDTM(ID)+D0M*DSFT1M-E0M*DABT1(ID-1)
|
|
REDN(ID)=REDN(ID)+D0*DSFN1-E0*DABN1(ID)
|
|
REDNM(ID)=REDNM(ID)+D0M*DSFN1M-E0M*DABN1(ID-1)
|
|
DO II=1,NLVEXP
|
|
REDP(II,ID)=REDP(II,ID)+D0*DSFP1(II)-E0*DABP1(II,ID)
|
|
REDPM(II,ID)=REDPM(II,ID)+D0M*DSFP1M(II)-
|
|
* E0M*DABP1(II,ID-1)
|
|
END DO
|
|
END IF
|
|
c
|
|
C Integral equation part of the radiative equilibrium
|
|
C
|
|
IF(REINT(ID).GT.0) THEN
|
|
SRH=SIGE*DENS1(ID)
|
|
ABST=ABSO1(ID)
|
|
ABSTE=ABST-ELSCAT(ID)
|
|
WWK=WW*ABST
|
|
FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABSTE*RAD1(ID))
|
|
D0=WW*(ABSTE*ALI1(ID)-ABST)
|
|
E0=WW*(RAD1(ID)-S0)
|
|
REIN(ID)=REIN(ID)+D0*DSFN1+E0*DABN1(ID)-WW*SIGEC(IJ)*RAD1(ID)
|
|
REIT(ID)=REIT(ID)+D0*DSFT1+E0*DABT1(ID)
|
|
DO II=1,NLVEXP
|
|
REIP(II,ID)=REIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
|
|
END DO
|
|
END IF
|
|
|
|
END DO
|
|
|
|
C
|
|
c ****3. deepest point - ID=ND
|
|
c
|
|
ID=ND
|
|
|
|
LNSKIP=.NOT.LSKIP(ID,IJ)
|
|
|
|
DSFTMM=DSFT1M
|
|
DSFNMM=DSFN1M
|
|
DO II=1,NLVEXP
|
|
DSFPMM(II)=DSFP1M(II)
|
|
END DO
|
|
DSFT1M=DSFT1
|
|
DSFN1M=DSFN1
|
|
DO II=1,NLVEXP
|
|
DSFP1M(II)=DSFP1(II)
|
|
END DO
|
|
S0=S0P
|
|
DSFT1=DSFT1P
|
|
DSFN1=DSFN1P
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=DSFP1P(II)
|
|
END DO
|
|
C
|
|
C Improved lower boundary condition
|
|
C
|
|
IF(IBC.GT.0.AND.IDISK.EQ.0) THEN
|
|
DT=UN/(DELDMZ(ID-1)*(ABSOT(ID)+ABSOT(ID-1)))
|
|
PLAD=XKFB(ID)/XKF1(ID)
|
|
DBDT=PLAD/XKF1(ID)*HKT21(ID)*FREQ(IJ)*DT
|
|
IF(IBC.EQ.1) THEN
|
|
DSFT1=DSFT1+DBDT
|
|
ELSE IF(IBC.GE.2) THEN
|
|
PLAM=XKFB(ID-1)/XKF1(ID-1)
|
|
TAU23=T23*DT
|
|
TAU43=T43*DT
|
|
D0=(PLAD*(UN+TAU43)-T43*PLAM*DT)*DT*DT
|
|
RHD=DELDMZ(ID-1)*DENSI(ID)
|
|
E0=D0*RHD
|
|
DSFT1=DSFT1+DBDT*(UN+TAU23)-E0*DABT1(ID)
|
|
DSFN1=DSFN1-E0*(DABN1(ID)+ABSO1(ID)*DENSIM(ID))
|
|
DO II=1,NLVEXP
|
|
DSFP1(II)=DSFP1(II)-E0*DABP1(II,ID)
|
|
END DO
|
|
IF(IBC.GE.3) THEN
|
|
DBDTM=PLAM/XKF1(ID-1)*HKT21(ID-1)*FREQ(IJ)*DT
|
|
RHD=DELDMZ(ID-1)*DENSI(ID-1)
|
|
E0=D0*RHD
|
|
DSFT1D=-DBDTM*DT*T23-E0*DABT1(ID-1)
|
|
DSFN1D=-E0*(DABN1(ID-1)+ABSO1(ID-1)*DENSIM(ID-1))
|
|
DO II=1,NLVEXP
|
|
DSFP1D(II)=-E0*DABP1(II,ID-1)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
END IF
|
|
C
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDT(ID)=DSFT1*ALI1(ID)
|
|
DSFDN(ID)=DSFN1*ALI1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDP(II,ID)=DSFP1(II)*ALI1(ID)
|
|
END DO
|
|
END IF
|
|
c
|
|
c Hydrostatic equilibrium equation
|
|
c
|
|
IF(LNSKIP) THEN
|
|
D0=WW*FAK1(ID)
|
|
A0=WW*FAK1(ID-1)
|
|
FPRD(ID)=FPRD(ID)+D0*RAD1(ID)-A0*RAD1(ID-1)
|
|
c F0=D0*ALIP1(ID)
|
|
E0=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
HEIT(ID)=HEIT(ID)+D0*DSFT1
|
|
HEIN(ID)=HEIN(ID)+D0*DSFN1
|
|
HEITM(ID)=HEITM(ID)+E0*DSFT1M
|
|
HEINM(ID)=HEINM(ID)+E0*DSFN1M
|
|
DO II=1,NLVEXP
|
|
HEIP(II,ID)=HEIP(II,ID)+D0*DSFP1(II)
|
|
HEIPM(II,ID)=HEIPM(II,ID)+E0*DSFP1M(II)
|
|
END DO
|
|
IF(IBC.GE.3) THEN
|
|
HEITM(ID)=HEITM(ID)-D0*DSFT1D
|
|
HEINM(ID)=HEINM(ID)-D0*DSFN1D
|
|
DO II=1,NLVEXP
|
|
HEIPM(II,ID)=HEIPM(II,ID)-D0*DSFP1D(II)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
C
|
|
C Differential equation part of radiative equilibrium
|
|
C
|
|
DDT=UN/(ABSOT(ID)+ABSOT(ID-1))
|
|
DT=DDT/DELDMZ(ID-1)
|
|
FL=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))*DT
|
|
FLFIX(ID)=FLFIX(ID)+WW*FL
|
|
IF(REDIF(ID).GT.0) THEN
|
|
D0=WW*FAK1(ID)*DT
|
|
A0=WW*FAK1(ID-1)*DT
|
|
D0M=D0*ALIM1(ID)-A0*ALI1(ID-1)
|
|
D0P=D0*ALIP1(ID)
|
|
D0=D0*ALI1(ID)-A0*ALIP1(ID-1)
|
|
E0=WW*FL*DDT
|
|
REDX(ID)=REDX(ID)+E0*ABSO1(ID)
|
|
REDXM(ID)=REDXM(ID)+E0*ABSO1(ID-1)
|
|
E0M=E0*DENSI(ID-1)
|
|
E0=E0*DENSI(ID)
|
|
REDT(ID)=REDT(ID)+D0*DSFT1-E0*DABT1(ID)
|
|
REDTM(ID)=REDTM(ID)+D0M*DSFT1M-E0M*DABT1(ID-1)
|
|
REDN(ID)=REDN(ID)+D0*DSFN1-E0*DABN1(ID)
|
|
REDNM(ID)=REDNM(ID)+D0M*DSFN1M-E0M*DABN1(ID-1)
|
|
DO II=1,NLVEXP
|
|
REDP(II,ID)=REDP(II,ID)+D0*DSFP1(II)-E0*DABP1(II,ID)
|
|
REDPM(II,ID)=REDPM(II,ID)+D0M*DSFP1M(II)-
|
|
* E0M*DABP1(II,ID-1)
|
|
END DO
|
|
IF(IBC.GE.3) THEN
|
|
REDTM(ID)=REDTM(ID)+D0*DSFT1D
|
|
REDNM(ID)=REDNM(ID)+D0*DSFN1D
|
|
DO II=1,NLVEXP
|
|
REDPM(II,ID)=REDPM(II,ID)+D0*DSFP1D(II)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
c
|
|
C Integral equation part of the radiative equilibrium
|
|
C
|
|
IF(REINT(ID).GT.0) THEN
|
|
SRH=SIGE*DENS1(ID)
|
|
ABST=ABSO1(ID)
|
|
ABSTE=ABST-ELSCAT(ID)
|
|
WWK=WW*ABST
|
|
FCOOLI(ID)=FCOOLI(ID)+WW*(EMIS1(ID)-ABSTE*RAD1(ID))
|
|
D0=WW*(ABSTE*ALI1(ID)-ABST)
|
|
E0=WW*(RAD1(ID)-S0)
|
|
REIN(ID)=REIN(ID)+D0*DSFN1+E0*DABN1(ID)-WW*SIGEC(IJ)*RAD1(ID)
|
|
IF(IBC.EQ.0) THEN
|
|
REIT(ID)=REIT(ID)+D0*DSFT1+E0*DABT1(ID)
|
|
ELSE
|
|
REIT(ID)=REIT(ID)+D0*(DSFT1-DBDT)+E0*DABT1(ID)+
|
|
* ALI1(ID)/ABST*DBDT
|
|
END IF
|
|
DO II=1,NLVEXP
|
|
REIP(II,ID)=REIP(II,ID)+D0*DSFP1(II)+E0*DABP1(II,ID)
|
|
END DO
|
|
END IF
|
|
|
|
RETURN
|
|
END
|