654 lines
20 KiB
Fortran
654 lines
20 KiB
Fortran
SUBROUTINE ALIFR6(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
|
|
WW=WC(IJ)
|
|
C
|
|
DSFT1M=0.
|
|
DSFN1M=0.
|
|
DO II=1,NLVEXP
|
|
DSFP1M(II)=0.
|
|
END DO
|
|
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)
|
|
if(ilmcor.ne.3) then
|
|
if(ilasct.eq.0) then
|
|
ABST=UN/(ABSO1(ID)-ELSCAT(ID))
|
|
S0=EMIS1(ID)*ABST
|
|
DSFN1=S0*(DEMN1(ID)*EMISIV-(DABN1(ID)-SIGEC(IJ))*ABST)
|
|
else
|
|
ABST=UN/ABSO1(ID)
|
|
S0=EMIS1(ID)*ABST
|
|
DSFN1=S0*(DEMN1(ID)*EMISIV-DABN1(ID)*ABST)
|
|
end if
|
|
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
|
|
C
|
|
c correction for electron scattering source function contribution
|
|
c
|
|
else
|
|
abst=un/abso1(id)
|
|
s0=emis1(id)*abst
|
|
sc=elec(id)*sigeC(IJ)
|
|
sct=sc*abst
|
|
st=s0+sct*rad1(id)
|
|
corr=un/(un-ali1(id)*sct)
|
|
dsft1=corr*(s0*demt1(id)*emisiv-st*dabt1(id)*abst)
|
|
dsfn1=corr*(s0*demn1(id)*emisiv+sigeC(IJ)*rad1(id)*abst-
|
|
* st*dabn1(id)*abst)
|
|
do ii=1,nlvexp
|
|
dsfp1(ii)=corr*(s0*demp1(ii,id)*emisiv-
|
|
* st*dabp1(ii,id)*abst)
|
|
end do
|
|
end if
|
|
|
|
EMISIP=UN/EMIS1(ID+1)
|
|
if(ilmcor.ne.3) then
|
|
if(ilasct.eq.0) then
|
|
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)
|
|
else
|
|
ABSTP=UN/ABSO1(ID+1)
|
|
S0P=EMIS1(ID+1)*ABSTP
|
|
DSFN1P=S0P*(DEMN1(ID+1)*EMISIP-DABN1(ID+1)*ABSTP)
|
|
end if
|
|
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
|
|
else
|
|
abstp=un/abso1(id+1)
|
|
s0p=emis1(id+1)*abstp
|
|
scp=elec(id+1)*sigeC(IJ)
|
|
sctp=scp*abstp
|
|
stp=s0p+sctp*rad1(id+1)
|
|
corrp=un/(un-ali1(id+1)*sctp)
|
|
dsft1p=corrp*(s0p*demt1(id+1)*emisip-stp*dabt1(id+1)*abstp)
|
|
dsfn1p=corrp*(s0p*demn1(id+1)*emisip+sigeC(IJ)*rad1(id+1)*abstp-
|
|
* stp*dabn1(id+1)*abstp)
|
|
do ii=1,nlvexp
|
|
dsfp1p(ii)=corrp*(s0p*demp1(ii,id+1)*emisip-
|
|
* stp*dabp1(ii,id+1)*abstp)
|
|
end do
|
|
end if
|
|
|
|
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
|
|
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDTM(ID)=DSFT1M*ALIM1(ID)
|
|
DSFDNM(ID)=DSFN1M*ALIM1(ID)
|
|
DSFDTP(ID)=DSFT1P*ALIP1(ID)
|
|
DSFDNP(ID)=DSFN1P*ALIP1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDPM(II,ID)=DSFP1M(II)*ALIM1(ID)
|
|
DSFDPP(II,ID)=DSFP1P(II)*ALIP1(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
|
|
IF(IFALI.GE.7) THEN
|
|
D0P=WF*ABSO1(ID)*ALIP1(ID)
|
|
HEITP(ID)=HEITP(ID)+D0P*DSFT1P
|
|
HEINP(ID)=HEINP(ID)+D0P*DSFN1P
|
|
DO II=1,NLVEXP
|
|
HEIPP(II,ID)=HEIPP(II,ID)+D0P*DSFP1P(II)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
c
|
|
c Differential equation part of radiative equilibrium
|
|
c
|
|
FLFIX(ID)=FLFIX(ID)+WF*RAD1(ID)-WW*HEXTRD(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
|
|
if(ilmcor.ne.3) then
|
|
if(ilasct.eq.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))
|
|
else
|
|
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)
|
|
end if
|
|
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
|
|
else
|
|
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))
|
|
end if
|
|
c
|
|
c the following are needed only for tridiagonal Lambda^star
|
|
C
|
|
C Upper sub-diagonal band
|
|
C
|
|
WWKC=WWK*ALIP1(ID)
|
|
CREIT(ID)=CREIT(ID)+WWKC*DSFT1P
|
|
CREIN(ID)=CREIN(ID)+WWKC*DSFN1P
|
|
DO II=1,NLVEXP
|
|
CREIP(II,ID)=CREIP(II,ID)+WWKC*DSFP1P(II)
|
|
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)
|
|
if(ilmcor.ne.3) then
|
|
if(ilasct.eq.0) then
|
|
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)
|
|
else
|
|
ABSTP=UN/ABSO1(ID+1)
|
|
S0P=EMIS1(ID+1)*ABSTP
|
|
DSFN1P=S0P*(DEMN1(ID+1)*EMISIP-DABN1(ID+1)*ABSTP)
|
|
end if
|
|
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
|
|
else
|
|
abstp=un/abso1(id+1)
|
|
s0p=emis1(id+1)*abstp
|
|
scp=elec(id+1)*sigeC(IJ)
|
|
sctp=scp*abstp
|
|
stp=s0p+sctp*rad1(id+1)
|
|
corrp=un/(un-ali1(id+1)*sctp)
|
|
dsft1p=corrp*(s0p*demt1(id+1)*emisip-stp*dabt1(id+1)*abstp)
|
|
dsfn1p=corrp*(s0p*demn1(id+1)*emisip+sigeC(IJ)*rad1(id+1)*abstp-
|
|
* stp*dabn1(id+1)*abstp)
|
|
do ii=1,nlvexp
|
|
dsfp1p(ii)=corrp*(s0p*demp1(ii,id+1)*emisip-
|
|
* stp*dabp1(ii,id+1)*abstp)
|
|
end do
|
|
end if
|
|
|
|
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
|
|
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDTM(ID)=DSFT1M*ALIM1(ID)
|
|
DSFDNM(ID)=DSFN1M*ALIM1(ID)
|
|
DSFDTP(ID)=DSFT1P*ALIP1(ID)
|
|
DSFDNP(ID)=DSFN1P*ALIP1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDPM(II,ID)=DSFP1M(II)*ALIM1(ID)
|
|
DSFDPP(II,ID)=DSFP1P(II)*ALIP1(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)
|
|
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(IFALI.GE.7) THEN
|
|
HEITP(ID)=HEITP(ID)+F0*DSFT1P
|
|
HEINP(ID)=HEINP(ID)+F0*DSFN1P
|
|
DO II=1,NLVEXP
|
|
HEIPP(II,ID)=HEIPP(II,ID)+F0*DSFP1P(II)
|
|
END DO
|
|
A0M=A0*ALIM1(ID-1)
|
|
EHET(ID)=EHET(ID)-A0M*DSFTMM
|
|
EHEN(ID)=EHEN(ID)-A0M*DSFNMM
|
|
DO II=1,NLVEXP
|
|
EHEP(II,ID)=EHEP(II,ID)-A0M*DSFPMM(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(IFALI.GE.7) THEN
|
|
REDTP(ID)=REDTP(ID)+D0P*DSFT1M
|
|
REDNP(ID)=REDNP(ID)+D0P*DSFN1M
|
|
DO II=1,NLVEXP
|
|
REDPP(II,ID)=REDPP(II,ID)+D0P*DSFP1P(II)
|
|
END DO
|
|
A0M=A0*ALIM1(ID-1)
|
|
ERET(ID)=ERET(ID)-A0M*DSFTMM
|
|
EREN(ID)=EREN(ID)-A0M*DSFNMM
|
|
DO II=1,NLVEXP
|
|
EREP(II,ID)=EREP(II,ID)-A0M*DSFPMM(II)
|
|
END DO
|
|
END IF
|
|
END IF
|
|
c
|
|
C Integral equation part of the radiative equilibrium
|
|
C
|
|
IF(REINT(ID).GT.0) THEN
|
|
if(ilmcor.ne.3) then
|
|
if(ilasct.eq.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))
|
|
else
|
|
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*SIGE*RAD1(ID)
|
|
end if
|
|
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
|
|
else
|
|
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)-sige)-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))
|
|
end if
|
|
c
|
|
c the following are needed only for tridiagonal Lambda^star
|
|
C
|
|
C Lower sub-diagonal band
|
|
C
|
|
WWKA=WWK*ALIM1(ID)
|
|
AREIT(ID)=AREIT(ID)+WWKA*DSFT1M
|
|
AREIN(ID)=AREIN(ID)+WWKA*DSFN1M
|
|
DO II=1,NLVEXP
|
|
AREIP(II,ID)=AREIP(II,ID)+WWKA*DSFP1M(II)
|
|
END DO
|
|
C
|
|
C Upper sub-diagonal band
|
|
C
|
|
WWKC=WWK*ALIP1(ID)
|
|
CREIT(ID)=CREIT(ID)+WWKC*DSFT1P
|
|
CREIN(ID)=CREIN(ID)+WWKC*DSFN1P
|
|
DO II=1,NLVEXP
|
|
CREIP(II,ID)=CREIP(II,ID)+WWKC*DSFP1P(II)
|
|
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
|
|
|
|
IF(IRDER.EQ.1.OR.IRDER.EQ.3) THEN
|
|
DSFDTM(ID)=DSFT1M*ALIM1(ID)
|
|
DSFDNM(ID)=DSFN1M*ALIM1(ID)
|
|
DSFDTP(ID)=DSFT1P*ALIP1(ID)
|
|
DSFDNP(ID)=DSFN1P*ALIP1(ID)
|
|
END IF
|
|
IF(IRDER.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
DSFDPM(II,ID)=DSFP1M(II)*ALIM1(ID)
|
|
DSFDPP(II,ID)=DSFP1P(II)*ALIP1(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)
|
|
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(IFALI.GE.7) THEN
|
|
HEITP(ID)=HEITP(ID)+F0*DSFT1P
|
|
HEINP(ID)=HEINP(ID)+F0*DSFN1P
|
|
DO II=1,NLVEXP
|
|
HEIPP(II,ID)=HEIPP(II,ID)+F0*DSFP1P(II)
|
|
END DO
|
|
A0M=A0*ALIM1(ID-1)
|
|
EHET(ID)=EHET(ID)-A0M*DSFTMM
|
|
EHEN(ID)=EHEN(ID)-A0M*DSFNMM
|
|
DO II=1,NLVEXP
|
|
EHEP(II,ID)=EHEP(II,ID)-A0M*DSFPMM(II)
|
|
END DO
|
|
END IF
|
|
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(IFALI.GE.7) THEN
|
|
REDTP(ID)=REDTP(ID)+D0P*DSFT1M
|
|
REDNP(ID)=REDNP(ID)+D0P*DSFN1M
|
|
DO II=1,NLVEXP
|
|
REDPP(II,ID)=REDPP(II,ID)+D0P*DSFP1P(II)
|
|
END DO
|
|
A0M=A0*ALIM1(ID-1)
|
|
ERET(ID)=ERET(ID)-A0M*DSFTMM
|
|
EREN(ID)=EREN(ID)-A0M*DSFNMM
|
|
DO II=1,NLVEXP
|
|
EREP(II,ID)=EREP(II,ID)-A0M*DSFPMM(II)
|
|
END DO
|
|
END IF
|
|
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
|
|
if(ilmcor.ne.3) then
|
|
if(ilasct.eq.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))
|
|
else
|
|
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)
|
|
end if
|
|
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
|
|
else
|
|
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
|
|
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
|
|
c
|
|
c the following are needed only for tridiagonal Lambda^star
|
|
C
|
|
C Lower sub-diagonal band
|
|
C
|
|
WWKA=WWK*ALIM1(ID)
|
|
AREIT(ID)=AREIT(ID)+WWKA*DSFT1M
|
|
AREIN(ID)=AREIN(ID)+WWKA*DSFN1M
|
|
DO II=1,NLVEXP
|
|
AREIP(II,ID)=AREIP(II,ID)+WWKA*DSFP1M(II)
|
|
END DO
|
|
END IF
|
|
|
|
RETURN
|
|
END
|