753 lines
19 KiB
Fortran
753 lines
19 KiB
Fortran
SUBROUTINE RTEFR1(IJ)
|
|
C =====================
|
|
C
|
|
C Solution of the radiative transfer equation - for one frequency
|
|
C - for the known source function.
|
|
C Determination of the radiation field and variable Eddington
|
|
C factors.
|
|
C
|
|
C The numerical method used:
|
|
C for ISPLIN = 0 - the ordinary Feautrier scheme
|
|
C = 1 - the spline collocation method
|
|
C = 2 - Hermitian fourth-order method
|
|
C = 3 - improved Feautrier scheme
|
|
C (Rybicki & Hummer 1991, A&A 245, 171.)
|
|
C
|
|
C In all cases, the overall matrix system is solved by the standard
|
|
C Gaussian elimination, analogous to that described in SOLVE
|
|
C (auxiliary matrix D is called ALF in SOLVE; auxiliary vector ANU
|
|
C is called BET in SOLVE)
|
|
C
|
|
C U0 - derivative of Q0 wrt taumin;
|
|
C for "fixed" frequencies, U0 has the meaning of
|
|
C absorption coefficient * second moment H
|
|
C ( a quantity needed for lower boundary condition of the
|
|
C hydrostatic equilibrium equation, specifically for
|
|
C accounting for an effect of fixed-option transitions)
|
|
C
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ITERAT.FOR'
|
|
PARAMETER (SIXTH=UN/6.D0,
|
|
* THIRD=UN/3.D0,
|
|
* TWOTHR=TWO/3.D0)
|
|
DIMENSION AANU(MDEPTH),DDD(MDEPTH),FKK(MDEPTH),
|
|
* RDD(MDEPTH),ST0(MDEPTH),SS0(MDEPTH),AB0(MDEPTH),
|
|
* AA(MMU,MMU),BB(MMU,MMU),CC(MMU,MMU),VL(MMU),
|
|
* FFD(MMU,MMU),FF0D(MMU,MMU),
|
|
* FFPD(MMU,MMU),ali0(mdepth),ss0c(mdepth),
|
|
* AAA(MDEPTH),BBB(MDEPTH),CCC(MDEPTH),EEE(MDEPTH),
|
|
* ZZZ(MDEPTH),ALRH(MDEPTH),ALRM(MDEPTH),ALRP(MDEPTH),
|
|
* D(MMU,MMU,MDEPTH),ANU(MMU,MDEPTH),scor(mdepth)
|
|
DIMENSION rmmu(2*MMU),wmmu(2*MMU),rwmu(2*MMU),
|
|
* dtau(mdepth),ri(mdepth),ali(mdepth),alij1(mdepth)
|
|
dimension tau(mdepth)
|
|
COMMON/OPTDPT/DT(MDEPTH)
|
|
C
|
|
WW=W(IJ)
|
|
ISPL=ISPLIN
|
|
IF(ISPLIN.GE.5) THEN
|
|
ISPLIN=ISPL-5
|
|
IF(IJALI(IJ).GT.0) THEN
|
|
IF(IRTE.EQ.0) THEN
|
|
CALL RTEDF1(IJ)
|
|
ELSE
|
|
CALL RTEDF2(IJ)
|
|
END IF
|
|
ISPLIN=ISPL
|
|
if(ifprad.eq.0) return
|
|
DO ID=1,ND
|
|
if(.not.lskip(ID,IJ))
|
|
* PRADT(ID)=PRADT(ID)+RAD1(ID)*FAK1(ID)*W(ij)
|
|
END DO
|
|
if(.not.lskip(1,IJ))
|
|
* PRD0=PRD0+ABSO1(1)*W(IJ)*(RAD1(1)*FH(IJ)-HEXTRD(IJ))
|
|
DO ID=1,ND
|
|
PRADA(ID)=PRADA(ID)+RAD1(ID)*FAK1(ID)*WW
|
|
END DO
|
|
RETURN
|
|
END IF
|
|
END IF
|
|
C
|
|
if(icompt.gt.0.and.(iter.gt.1.or.ilam.gt.0)) then
|
|
call rtecf1(ij)
|
|
return
|
|
end if
|
|
c
|
|
C
|
|
FR=FREQ(IJ)
|
|
C
|
|
C total source function
|
|
C
|
|
AH=0.
|
|
DO ID=1,ND
|
|
AB0(ID)=ABSO1(ID)
|
|
ST0(ID)=EMIS1(ID)/AB0(ID)
|
|
RAD1(ID)=0.
|
|
ALI1(ID)=0.
|
|
END DO
|
|
C
|
|
C non-coherent electron scattering by lambda iteration
|
|
C
|
|
IF(NELSC.LE.0) THEN
|
|
DO ID=1,ND
|
|
SS0(ID)=-SCAT1(ID)/AB0(ID)
|
|
END DO
|
|
ELSE
|
|
DO ID=1,ND
|
|
ST0(ID)=ST0(ID)+SCAT1(ID)*EMEL1(ID)*RAD1(ID)/AB0(ID)
|
|
SS0(ID)=0.
|
|
END DO
|
|
END IF
|
|
C
|
|
C optical depth scale
|
|
C
|
|
tau(1)=absot(1)*dm(1)
|
|
DO ID=1,ND-1
|
|
DT(ID)=DELDMZ(ID)*(ABSOT(ID+1)+ABSOT(ID))
|
|
tau(id+1)=tau(id)+dt(id)
|
|
END DO
|
|
c
|
|
U0=0.
|
|
QQ0=0.
|
|
US0=0.
|
|
c TAUMIN=ABSO1(1)*DEDM1
|
|
taumin=absot(1)*dm(1)/2.
|
|
ALB1=0.
|
|
C
|
|
C Allowance for wind blanketing
|
|
C
|
|
IF(IWINBL.GT.0) ALB1=TWO*ALBE(IJ)/(UN+ALBE(IJ))
|
|
C
|
|
C ************** Forward elimination
|
|
C
|
|
C Upper boundary condition
|
|
C
|
|
ID=1
|
|
DTP1=DT(1)
|
|
P0=0.
|
|
EX=UN
|
|
IF(MOD(ISPLIN,3).EQ.0) THEN
|
|
B=DTP1*HALF
|
|
C=0.
|
|
ELSE
|
|
B=DTP1*THIRD
|
|
C=B*HALF
|
|
END IF
|
|
QQ0=0.
|
|
US0=0.
|
|
DO I=1,NMU
|
|
IF(IDISK.EQ.0) THEN
|
|
C
|
|
C allowance for non-zero optical depth at the first depth point
|
|
C
|
|
TAMM=TAUMIN/AMU(I)
|
|
EX=EXP(-TAMM)
|
|
P0=UN-EX
|
|
QQ0=QQ0+P0*AMU(I)*WTMU(I)
|
|
U0=U0+EX*WTMU(I)
|
|
US0=US0+P0/TAMM*WTMU(I)
|
|
END IF
|
|
C
|
|
BI=B/AMU(I)
|
|
CI=C/AMU(I)
|
|
VL(I)=(BI+P0)*ST0(ID)+CI*ST0(ID+1)
|
|
IF(IWINBL.LT.0) VL(I)=VL(I)+EXTINT(IJ,I)
|
|
DO J=1,NMU
|
|
BB(I,J)=SS0(ID)*WTMU(J)*(BI+P0)-ALB1*WTMU(J)
|
|
CC(I,J)=-CI*SS0(ID+1)*WTMU(J)
|
|
END DO
|
|
BB(I,I)=BB(I,I)+AMU(I)/DTP1+UN+BI
|
|
CC(I,I)=CC(I,I)+AMU(I)/DTP1-CI
|
|
ANU(I,ID)=0.
|
|
END DO
|
|
C
|
|
IF(ISPLIN.LE.2) THEN
|
|
CALL MATINV(BB,NMU,MMU)
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
D(I,J,ID)=0.
|
|
DO K=1,NMU
|
|
D(I,J,ID)=D(I,J,ID)+BB(I,K)*CC(K,J)
|
|
END DO
|
|
ANU(I,1)=ANU(I,1)+BB(I,J)*VL(J)
|
|
END DO
|
|
END DO
|
|
ELSE
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
FF0D(I,J)=BB(I,J)/CC(I,I)
|
|
END DO
|
|
FF0D(I,I)=FF0D(I,I)-UN
|
|
END DO
|
|
C
|
|
c CALL MINV3(BB)
|
|
CALL MATINV(BB,NMU,MMU)
|
|
DO I=1,NMU
|
|
ANU(I,ID)=0.
|
|
DO J=1,NMU
|
|
D(I,J,ID)=BB(I,J)*CC(J,J)
|
|
ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
|
|
END DO
|
|
END DO
|
|
END IF
|
|
C
|
|
C Normal depth points 1 < ID < ND
|
|
C
|
|
DO ID=2,ND-1
|
|
DTM1=DTP1
|
|
DTP1=DT(ID)
|
|
DT0=TWO/(DTM1+DTP1)
|
|
AL=UN/DTM1*DT0
|
|
GA=UN/DTP1*DT0
|
|
BE=AL+GA
|
|
IF(MOD(ISPLIN,3).EQ.0) THEN
|
|
A=0.
|
|
C=0.
|
|
ELSE IF(ISPLIN.EQ.1) THEN
|
|
A=DTM1*DT0*SIXTH
|
|
C=DTP1*DT0*SIXTH
|
|
ELSE
|
|
A=(UN-HALF*AL*DTP1*DTP1)*SIXTH
|
|
C=(UN-HALF*GA*DTM1*DTM1)*SIXTH
|
|
END IF
|
|
B=UN-A-C
|
|
VL0=A*ST0(ID-1)+B*ST0(ID)+C*ST0(ID+1)
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
AA(I,J)=-A*SS0(ID-1)*WTMU(J)
|
|
CC(I,J)=-C*SS0(ID+1)*WTMU(J)
|
|
BB(I,J)=B*SS0(ID)*WTMU(J)
|
|
END DO
|
|
END DO
|
|
DO I=1,NMU
|
|
VL(I)=VL0
|
|
DIV=AMU(I)*AMU(I)
|
|
AA(I,I)=AA(I,I)+DIV*AL-A
|
|
CC(I,I)=CC(I,I)+DIV*GA-C
|
|
BB(I,I)=BB(I,I)+DIV*BE+B
|
|
END DO
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
VL(I)=VL(I)+AA(I,J)*ANU(J,ID-1)
|
|
END DO
|
|
END DO
|
|
IF(ISPLIN.LE.2) THEN
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
S=0.
|
|
DO K=1,NMU
|
|
S=S+AA(I,K)*D(K,J,ID-1)
|
|
END DO
|
|
BB(I,J)=BB(I,J)-S
|
|
END DO
|
|
END DO
|
|
C
|
|
CALL MATINV(BB,NMU,MMU)
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
D(I,J,ID)=0.
|
|
DO K=1,NMU
|
|
D(I,J,ID)=D(I,J,ID)+BB(I,K)*CC(K,J)
|
|
END DO
|
|
END DO
|
|
END DO
|
|
ELSE
|
|
DO I=1,NMU
|
|
BB(I,I)=-AA(I,I)+BB(I,I)-CC(I,I)
|
|
DO J=1,NMU
|
|
FFPD(I,J)=AA(I,I)*FF0D(I,J)
|
|
END DO
|
|
END DO
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
S=0.
|
|
DO K=1,NMU
|
|
S=S+FFPD(I,K)*D(K,J,ID-1)
|
|
END DO
|
|
FFD(I,J)=(BB(I,J)+S)/CC(I,I)
|
|
END DO
|
|
END DO
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
FF0D(I,J)=FFD(I,J)
|
|
END DO
|
|
FFD(I,I)=FFD(I,I)+UN
|
|
END DO
|
|
C
|
|
CALL MATINV(FFD,NMU,MMU)
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
D(I,J,ID)=FFD(I,J)
|
|
BB(I,J)=FFD(I,J)/CC(J,J)
|
|
END DO
|
|
END DO
|
|
END IF
|
|
DO I=1,NMU
|
|
ANU(I,ID)=0.
|
|
DO J=1,NMU
|
|
ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
|
|
END DO
|
|
END DO
|
|
END DO
|
|
C
|
|
C Lower boundary condition
|
|
C
|
|
ID=ND
|
|
C
|
|
C First option:
|
|
C b.c. is different from stellar atmospheres; expresses symmetry
|
|
C at the central plane I(taumax,-mu,nu)=I(taumax,+mu,nu)
|
|
C
|
|
IF(IFZ0.GE.0.AND.IDISK.EQ.1) THEN
|
|
B=DTP1*HALF
|
|
A=0.
|
|
DO I=1,NMU
|
|
BI=B/AMU(I)
|
|
AI=A/AMU(I)
|
|
VL(I)=ST0(ID)*BI+ST0(ID-1)*AI
|
|
DO J=1,NMU
|
|
AA(I,J)=-AI*SS0(ID-1)*WTMU(J)
|
|
BB(I,J)=BI*SS0(ID)*WTMU(J)
|
|
END DO
|
|
AA(I,I)=AA(I,I)+AMU(I)/DTP1-AI
|
|
BB(I,I)=BB(I,I)+AMU(I)/DTP1+BI
|
|
END DO
|
|
DO I=1,NMU
|
|
S1=0.
|
|
DO J=1,NMU
|
|
S=0.
|
|
S1=S1+AA(I,J)*ANU(J,ID-1)
|
|
DO K=1,NMU
|
|
S=S+AA(I,K)*D(K,J,ID-1)
|
|
END DO
|
|
BB(I,J)=BB(I,J)-S
|
|
END DO
|
|
VL(I)=VL(I)+S1
|
|
END DO
|
|
C
|
|
C Second option:
|
|
C b.c. is the same as in stellar atmospheres - the last depth point
|
|
C is not at the central plane
|
|
C
|
|
ELSE
|
|
FR15=FR*1.D-15
|
|
BNU=BN*FR15*FR15*FR15
|
|
PLAND=BNU/(EXP(HK*FR/TEMP(ND ))-UN)*RRDIL
|
|
DPLAN=BNU/(EXP(HK*FR/TEMP(ND-1))-UN)*RRDIL
|
|
IF(TEMPBD.GT.0.) THEN
|
|
PLAND=BNU/(EXP(HK*FR/TEMPBD)-UN)*RRDIL
|
|
DPLAN=BNU/(EXP(HK*FR/TEMPBD)-UN)*RRDIL
|
|
END IF
|
|
DPLAN=(PLAND-DPLAN)/DT(ND-1)
|
|
IF(IBC.EQ.0.OR.IBC.EQ.4) THEN
|
|
DO I=1,NMU
|
|
AA(I,I)=AMU(I)/DTP1
|
|
VL(I)=PLAND+AMU(I)*DPLAN+AA(I,I)*ANU(I,ID-1)
|
|
DO J=1,NMU
|
|
BB(I,J)=-AA(I,I)*D(I,J,ID-1)
|
|
END DO
|
|
BB(I,I)=BB(I,I)+AA(I,I)+UN
|
|
END DO
|
|
ELSE
|
|
DO I=1,NMU
|
|
A=AMU(I)/DTP1
|
|
B=HALF/A
|
|
AA(I,I)=A
|
|
VL(I)=B*ST0(ID)+PLAND+AMU(I)*DPLAN+AA(I,I)*ANU(I,ID-1)
|
|
DO J=1,NMU
|
|
BB(I,J)=B*SS0(ID)*WTMU(J)-AA(I,I)*D(I,J,ID-1)
|
|
END DO
|
|
BB(I,I)=BB(I,I)+A+B+UN
|
|
END DO
|
|
END IF
|
|
END IF
|
|
C
|
|
CALL MATINV(BB,NMU,MMU)
|
|
C
|
|
DO I=1,NMU
|
|
ANU(I,ID)=0.
|
|
DO J=1,NMU
|
|
D(I,J,ID)=0.
|
|
ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
|
|
END DO
|
|
END DO
|
|
C
|
|
C ***************** Backsolution
|
|
C
|
|
FKK(ND)=THIRD
|
|
AJ=0.
|
|
AH=0.
|
|
AK=0.
|
|
DO I=1,NMU
|
|
RMU=WTMU(I)*ANU(I,ID)
|
|
AJ=AJ+RMU
|
|
AH=AH+RMU*AMU(I)
|
|
AK=AK+RMU*AMU(I)*AMU(I)
|
|
END DO
|
|
RDD(ID)=AJ
|
|
IF(IBC.EQ.0) THEN
|
|
FKK(ND)=THIRD
|
|
ELSE
|
|
FKK(ID)=AK/AJ
|
|
FHD(IJ)=AH/AJ
|
|
END IF
|
|
C
|
|
DO ID=ND-1,1,-1
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
ANU(I,ID)=ANU(I,ID)+D(I,J,ID)*ANU(J,ID+1)
|
|
END DO
|
|
END DO
|
|
AJ=0.
|
|
AK=0.
|
|
DO I=1,NMU
|
|
RMU=WTMU(I)*ANU(I,ID)
|
|
AJ=AJ+RMU
|
|
AK=AK+RMU*AMU(I)*AMU(I)
|
|
END DO
|
|
C
|
|
C solution of the transfer equation
|
|
C Variables:
|
|
C ANU - Feautrier intensity
|
|
C RDD - mean intensity
|
|
C FKK - Eddington factor f(K) = K/J
|
|
C
|
|
FKK(ID)=AK/AJ
|
|
RDD(ID)=AJ
|
|
END DO
|
|
C
|
|
if(idisk.eq.1) then
|
|
do id=1,nd
|
|
fak(ij,id)=fkk(id)
|
|
end do
|
|
end if
|
|
C
|
|
C the "surface" Eddington factor fH
|
|
C
|
|
AH=0.
|
|
DO I=1,NMU
|
|
AH=AH+WTMU(I)*AMU(I)*ANU(I,1)
|
|
END DO
|
|
FH0=AH/AJ-HALF*ALB1
|
|
FH(IJ)=FH0
|
|
q0(ij)=qq0
|
|
uu0(ij)=u0
|
|
|
|
c q0(ij)=0.
|
|
c uu0(ij)=0.
|
|
|
|
C
|
|
C ********************
|
|
C
|
|
C Again solution of the transfer equation, now with Eddington
|
|
C FKK and FH determined above, to insure strict consistency of the
|
|
C radiation field and Eddington factors
|
|
C
|
|
if(ilmcor.eq.2) then
|
|
do id=1,nd
|
|
scor(id)=un/(un+ss0(id))
|
|
end do
|
|
else if(ilmcor.eq.3) then
|
|
do id=1,nd
|
|
ss0c(id)=ss0(id)
|
|
st0(id)=st0(id)-ss0(id)*rdd(id)
|
|
ss0(id)=0.
|
|
end do
|
|
end if
|
|
C
|
|
C Upper boundary condition
|
|
C
|
|
ID=1
|
|
DTP1=DT(ID)
|
|
IF(MOD(ISPLIN,3).EQ.0) THEN
|
|
B=DTP1*HALF
|
|
C=0.
|
|
ELSE
|
|
B=DTP1*THIRD
|
|
C=B*HALF
|
|
END IF
|
|
BQ=UN/(B+QQ0)
|
|
CQ=C*BQ
|
|
BBB(ID)=(FKK(ID)/DTP1+FH0+B)*BQ+SS0(ID)
|
|
CCC(ID)=(FKK(ID+1)/DTP1)*BQ-CQ*(UN+SS0(ID+1))
|
|
VLL=ST0(ID)+CQ*ST0(ID+1)
|
|
IF(IWINBL.LT.0) VLL=VLL+HEXTRD(IJ)*BQ
|
|
if(ilmcor.eq.2) then
|
|
bbb(id)=bbb(id)*scor(id)
|
|
ccc(id)=ccc(id)*scor(id)
|
|
vll=vll*scor(id)
|
|
end if
|
|
ZZZ(ID)=UN/BBB(ID)
|
|
AANU(ID)=VLL*ZZZ(ID)
|
|
DDD(ID)=CCC(ID)*ZZZ(ID)
|
|
IF(ISPLIN.GT.2) FFF=BBB(ID)/CCC(ID)-UN
|
|
C
|
|
C Normal depth point
|
|
C
|
|
DO ID=2,ND-1
|
|
DTM1=DTP1
|
|
DTP1=DT(ID)
|
|
DT0=TWO/(DTP1+DTM1)
|
|
AL=UN/DTM1*DT0
|
|
GA=UN/DTP1*DT0
|
|
IF(MOD(ISPLIN,3).EQ.0) THEN
|
|
A=0.
|
|
C=0.
|
|
ELSE IF(ISPLIN.EQ.1) THEN
|
|
A=DTM1*DT0*SIXTH
|
|
C=DTP1*DT0*SIXTH
|
|
ELSE
|
|
A=(UN-HALF*DTP1*DTP1*AL)*SIXTH
|
|
C=(UN-HALF*DTM1*DTM1*GA)*SIXTH
|
|
END IF
|
|
AAA(ID)=AL*FKK(ID-1)-A*(UN+SS0(ID-1))
|
|
CCC(ID)=GA*FKK(ID+1)-C*(UN+SS0(ID+1))
|
|
BBB(ID)=(AL+GA)*FKK(ID)+(UN-A-C)*(UN+SS0(ID))
|
|
VLL=A*ST0(ID-1)+C*ST0(ID+1)+(UN-A-C)*ST0(ID)
|
|
if(ilmcor.eq.2) then
|
|
aaa(id)=aaa(id)*scor(id)
|
|
bbb(id)=bbb(id)*scor(id)
|
|
ccc(id)=ccc(id)*scor(id)
|
|
vll=vll*scor(id)
|
|
end if
|
|
AANU(ID)=VLL+AAA(ID)*AANU(ID-1)
|
|
IF(ISPLIN.LE.2) THEN
|
|
ZZZ(ID)=UN/(BBB(ID)-AAA(ID)*DDD(ID-1))
|
|
DDD(ID)=CCC(ID)*ZZZ(ID)
|
|
AANU(ID)=AANU(ID)*ZZZ(ID)
|
|
ELSE
|
|
SUM=-AAA(ID)+BBB(ID)-CCC(ID)
|
|
FFF=(SUM+AAA(ID)*FFF*DDD(ID-1))/CCC(ID)
|
|
DDD(ID)=UN/(UN+FFF)
|
|
AANU(ID)=AANU(ID)*DDD(ID)/CCC(ID)
|
|
END IF
|
|
END DO
|
|
C
|
|
C Lower boundary condition
|
|
C
|
|
ID=ND
|
|
c
|
|
c stellar atmospheric
|
|
c
|
|
IF(IDISK.EQ.0.OR.IFZ0.LT.0) then
|
|
IF(IBC.EQ.0) THEN
|
|
BBB(ID)=FKK(ID)/DTP1+HALF
|
|
AAA(ID)=FKK(ID-1)/DTP1
|
|
VLL=HALF*PLAND+THIRD*DPLAN
|
|
ELSE IF(IBC.LT.4) THEN
|
|
B=UN/DTP1
|
|
A=TWO*B*B
|
|
BBB(ID)=UN+SS0(ID)+B*TWO*FHD(IJ)+A*FKK(ID)
|
|
AAA(ID)=A*FKK(ID-1)
|
|
VLL=ST0(ID)+B*(PLAND+TWOTHR*DPLAN)
|
|
ELSE
|
|
B=UN/DTP1
|
|
A=TWO*B*B
|
|
BBB(ID)=B+A*FKK(ID)
|
|
AAA(ID)=A*FKK(ID-1)
|
|
VLL=B*(PLAND+TWOTHR*DPLAN)
|
|
END IF
|
|
c
|
|
c accretion disk - symmetric boundary
|
|
c
|
|
ELSE
|
|
B=TWO/DTP1
|
|
BBB(ID)=FKK(ID)/DTP1*B+UN+SS0(ND)
|
|
AAA(ID)=FKK(ID-1)/DTP1*B
|
|
VLL=ST0(ID)
|
|
END IF
|
|
if(ilmcor.eq.2) then
|
|
aaa(id)=aaa(id)*scor(id)
|
|
bbb(id)=bbb(id)*scor(id)
|
|
vll=vll*scor(id)
|
|
end if
|
|
EEE(ND)=AAA(ID)/BBB(ID)
|
|
ZZZ(ID)=UN/(BBB(ID)-AAA(ID)*DDD(ID-1))
|
|
RAD1(ID)=(VLL+AAA(ID)*AANU(ID-1))*ZZZ(ID)
|
|
FAK1(ID)=FKK(ND)
|
|
ALRH(ID)=ZZZ(ID)
|
|
c
|
|
C Backsolution
|
|
C
|
|
DO ID=ND-1,1,-1
|
|
EEE(ID)=AAA(ID)/(BBB(ID)-CCC(ID)*EEE(ID+1))
|
|
RAD1(ID)=AANU(ID)+DDD(ID)*RAD1(ID+1)
|
|
FAK1(ID)=FKK(ID)
|
|
ALRH(ID)=ZZZ(ID)/(UN-DDD(ID)*EEE(ID+1))
|
|
ALRM(ID)=0
|
|
ALRP(ID)=0
|
|
END DO
|
|
FLUX(IJ)=FH(IJ)*RAD1(1)-half*hextrd(ij)
|
|
* -(st0(1)-ss0(1)*rad1(1))*q0(ij)
|
|
|
|
C
|
|
C evaluate approximate Lambda operator
|
|
C
|
|
C a) Rybicki-Hummer Lambda^star operator (diagonal)
|
|
C (for JALI = 1)
|
|
C
|
|
DO ID=1,ND
|
|
ALIM1(ID)=0.
|
|
ALIP1(ID)=0.
|
|
END DO
|
|
IF(JALI.EQ.1) THEN
|
|
DO ID=1,ND
|
|
ALI1(ID)=ALRH(ID)
|
|
END DO
|
|
c
|
|
IF(IBC.EQ.0) THEN
|
|
ali1(nd-1)=rad1(nd-1)/st0(nd-1)
|
|
ali1(nd)=rad1(nd)/st0(nd)
|
|
END IF
|
|
C
|
|
C for IFALI>5:
|
|
C tridiagonal Rybicki-Hummer operator (off-diagonal terms)
|
|
C
|
|
IF(IFALI.GE.6) THEN
|
|
ALIP1(1)=ALRH(2)*DDD(1)
|
|
DO ID=2,ND-1
|
|
ALIM1(ID)=ALRH(ID-1)*EEE(ID)
|
|
ALIP1(ID)=ALRH(ID+1)*DDD(ID)
|
|
END DO
|
|
ALIM1(ND)=ALRH(ND-1)*EEE(ND)
|
|
IF(IBC.EQ.0) THEN
|
|
ALIM1(nd)=0.
|
|
ALIM1(nd-1)=0.
|
|
ALIP1(nd)=0.
|
|
ALIP1(nd-1)=0.
|
|
END IF
|
|
END IF
|
|
c
|
|
C b) diagonal Olson-Kunasz Lambda^star operator,
|
|
C (for JALI = 2)
|
|
C
|
|
ELSE IF(JALI.EQ.2) THEN
|
|
DO ID=1,ND-1
|
|
ALI0(ID)=0.
|
|
DO I=1,NMU
|
|
DIV=DT(ID)/AMU(I)
|
|
ALI0(ID)=ALI0(ID)+(UN-EXP(-DIV))/DIV*WTMU(I)
|
|
END DO
|
|
END DO
|
|
DO ID=2,ND-1
|
|
ALI1(ID)=UN-HALF*(ALI0(ID)+ALI0(ID-1))
|
|
END DO
|
|
ALI1(1)=UN-HALF*(ALI0(1)+US0)
|
|
ALI1(ND)=UN-ALI0(ND-1)
|
|
ali1(nd-1)=rad1(nd-1)/st0(nd-1)
|
|
ali1(nd)=rad1(nd)/st0(nd)
|
|
END IF
|
|
C
|
|
C correction of Lambda^star for scattering
|
|
C
|
|
IF(ILMCOR.EQ.1) THEN
|
|
DO ID=1,ND
|
|
ALI1(ID)=ALI1(ID)*(UN+SS0(ID))
|
|
ALIM1(ID)=ALIM1(ID)*(UN+SS0(ID))
|
|
ALIP1(ID)=ALIP1(ID)*(UN+SS0(ID))
|
|
END DO
|
|
IF(IBC.EQ.4) THEN
|
|
ALI1(ND)= ALI1(ND)/(UN+SS0(ND))
|
|
ALIM1(ND)= ALIM1(ND)/(UN+SS0(ND))
|
|
ALIP1(ND)= ALIP1(ND)/(UN+SS0(ND))
|
|
END IF
|
|
END IF
|
|
C
|
|
if(ifalih.gt.0) then
|
|
c
|
|
c solution for the individual angles - to get Lambda^star_H
|
|
C
|
|
do id=1,nd
|
|
alih1(id)=0.
|
|
alij1(id)=0.
|
|
end do
|
|
nw=nmu
|
|
do i=1,nw
|
|
rmmu(i)=-amu(nw-i+1)
|
|
rmmu(i+nw)=amu(i)
|
|
wmmu(i)=wtmu(nw-i+1)
|
|
wmmu(i+nw)=wtmu(i)
|
|
end do
|
|
do i=1,2*nw
|
|
rwmu(i)=rmmu(i)*wmmu(i)*half
|
|
end do
|
|
C
|
|
c --------------------- loop over angles
|
|
c
|
|
do i=1,2*nw
|
|
do id=1,nd-1
|
|
dtau(id)=dt(id)/abs(rmmu(i))
|
|
end do
|
|
c
|
|
c boundary conditions
|
|
c
|
|
c rup=extrad(ij)
|
|
rup=extint(ij,i)
|
|
C
|
|
C diffusion approximation for semi-infinite atmospheres
|
|
C
|
|
rdown=pland+rmmu(i)*dplan
|
|
c
|
|
c solution of the transfer equation
|
|
c
|
|
call rtesol(dtau,st0,rup,rdown,rmmu(i),ri,ali)
|
|
c
|
|
DO ID=1,ND
|
|
alih1(id)=alih1(id)+rwmu(i)*ali(id)
|
|
alij1(id)=alij1(id)+wmmu(i)*ali(id)*half
|
|
END DO
|
|
end do
|
|
end if
|
|
C
|
|
c --------------------- end of loop over angles
|
|
c
|
|
ISPLIN=ISPL
|
|
c
|
|
if(idisk.ne.0) then
|
|
iji=nfreq-kij(ij)+1
|
|
DO ID=1,ND
|
|
rad(iji,id)=rad1(id)
|
|
END DO
|
|
end if
|
|
C
|
|
C radiation pressure
|
|
C
|
|
if(ifprad.gt.0) then
|
|
if(.not.lskip(1,IJ))
|
|
* PRD0=PRD0+ABSO1(1)*WW*(RAD1(1)*FH(IJ)-HEXTRD(IJ))
|
|
DO ID=1,ND
|
|
if(.not.lskip(ID,IJ))
|
|
* PRADT(ID)=PRADT(ID)+RAD1(ID)*FAK1(ID)*WW
|
|
PRADA(ID)=PRADA(ID)+RAD1(ID)*FAK1(ID)*WW
|
|
END DO
|
|
end if
|
|
c
|
|
if(chmax.ge.1.91e-3.and.chmax.le.2.03e-3) then
|
|
tauij=taumin
|
|
do id=1,nd
|
|
if(ilmcor.eq.3) ss0(id)=ss0c(id)
|
|
if(id.gt.1) tauij=tauij+dt(id-1)
|
|
write(97,697) ij,id,tauij,rad1(id),st0(id)/(un+ss0(id)),st0(id),
|
|
* un+ss0(id),ali1(id)
|
|
end do
|
|
697 format(2i4,1p6e12.4)
|
|
end if
|
|
|
|
C
|
|
C store quantities for explicit (linearized) frequencies
|
|
C
|
|
IF(IJEX(IJ).LE.0) RETURN
|
|
IJE=IJEX(IJ)
|
|
DO ID=1,ND
|
|
RADEX(IJE,ID)=RAD1(ID)
|
|
FAKEX(IJE,ID)=FAK1(ID)
|
|
END DO
|
|
c Q0(IJE)=QQ0
|
|
c UU0(IJE)=U0
|
|
RETURN
|
|
END
|