452 lines
12 KiB
Fortran
452 lines
12 KiB
Fortran
SUBROUTINE RTECF1(IJ)
|
|
C =====================
|
|
C
|
|
C Solution of the radiative transfer equation with Compton scattering
|
|
C for one frequency (assuming the radiation intensity in i
|
|
C other frequencies is given
|
|
C solution is done for individual angles, and new Eddington factors
|
|
C are determined
|
|
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)
|
|
COMMON/OPTDPT/DT(MDEPTH)
|
|
COMMON/SURFEX/EXTJ(MFREQ),EXTH(MFREQ)
|
|
COMMON/EXTINT/WANGLE,EXTIN(MFREQ)
|
|
COMMON/AUXRTE/
|
|
* COMA(MDEPTH),COMB(MDEPTH),COMC(MDEPTH),VL(MDEPTH),
|
|
* COME(MDEPTH),U(MDEPTH),V(MDEPTH),BS(MDEPTH),
|
|
* AL(MDEPTH),BE(MDEPTH),GA(MDEPTH)
|
|
common/comgfs/gfm(mfreq,mdeptc),gfp(mfreq,mdeptc)
|
|
DIMENSION RI(MDEPTH),RDH(MDEPTH),RDK(MDEPTH),RDN(MDEPTH),
|
|
* DTAU(MDEPTH),ST0(MDEPTH),RDWN(MMUC),
|
|
* ali(mdepth)
|
|
DIMENSION AANU(MDEPTH),DDD(MDEPTH),FKK(MDEPTH),ali0(mdepth),
|
|
* SS0(MDEPTH),
|
|
* AAA(MDEPTH),BBB(MDEPTH),CCC(MDEPTH),EEE(MDEPTH),
|
|
* ZZZ(MDEPTH),ALRH(MDEPTH),ALRM(MDEPTH),ALRP(MDEPTH),
|
|
* ss0c(mdepth)
|
|
C
|
|
IF(IJ.EQ.1) THEN
|
|
if(icompt.gt.0.and.icombc.gt.0) then
|
|
IJE=IJEX(IJ)
|
|
DO ID=1,ND
|
|
rad1(id)=rad(nfreq,id)
|
|
fak1(id)=0.333333
|
|
ali1(id)=0.
|
|
if(ije.gt.0) then
|
|
RADEX(IJE,ID)=rad1(id)
|
|
FAKEX(IJE,ID)=fak1(id)
|
|
END IF
|
|
END DO
|
|
return
|
|
end if
|
|
END IF
|
|
C
|
|
WW=W(IJ)
|
|
IJI=NFREQ-KIJ(IJ)+1
|
|
FR=FREQ(IJ)
|
|
CALL RTECF0(IJ)
|
|
c
|
|
do id=1,nd
|
|
rad1(id)=0.
|
|
ali1(id)=0.
|
|
rdh(id)=0.
|
|
rdk(id)=0.
|
|
rdn(id)=0.
|
|
st0(id)=vl(id)+(comb(id)+bs(id))*rad(iji,id)
|
|
ss0(id)=0.
|
|
end do
|
|
rdh1=0.
|
|
rdhd=0.
|
|
c
|
|
if(iji.gt.1) then
|
|
do id=1,nd
|
|
st0(id)=st0(id)+coma(id)*rad(iji-1,id)
|
|
end do
|
|
end if
|
|
if(iji.lt.nfreq) then
|
|
do id=1,nd
|
|
st0(id)=st0(id)+comc(id)*rad(iji+1,id)
|
|
end do
|
|
end if
|
|
c
|
|
if(idisk.eq.0.or.ifz0.lt.0) then
|
|
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
|
|
ENDIF
|
|
DPLAN=(PLAND-DPLAN)/DT(ND-1)
|
|
end if
|
|
c
|
|
if(icomrt.eq.0) then
|
|
c
|
|
c ========================================================
|
|
c Formal angle-dependent solution done by Feautrier scheme
|
|
c ========================================================
|
|
c
|
|
c loop over angles points
|
|
c
|
|
do i=1,nmu
|
|
do id=1,nd-1
|
|
dtau(id)=dt(id)/amu(i)
|
|
end do
|
|
c
|
|
c boundary conditions
|
|
c
|
|
rup=0.
|
|
rdown=0.
|
|
rup=extint(ij,i)
|
|
if(idisk.eq.0.or.ifz0.lt.0) rdown=pland+amu(i)*dplan
|
|
c
|
|
c solution of the transfer equation
|
|
c
|
|
call rtefe2(dtau,st0,rup,rdown,ri)
|
|
ttau=0.
|
|
do id=1,nd
|
|
riid=wtmu(i)*ri(id)
|
|
rad1(id)=rad1(id)+riid
|
|
rdk(id)=rdk(id)+amu(i)*amu(i)*riid
|
|
end do
|
|
rdh1=rdh1+amu(i)*wtmu(i)*ri(1)
|
|
rdhd=rdhd+amu(i)*wtmu(i)*ri(nd)
|
|
end do
|
|
rdh1=rdh1-half*hextrd(ij)
|
|
c
|
|
c ----------------------
|
|
c end of the loop over angle points
|
|
c
|
|
c ===========================================
|
|
c Formal angle-dependent solution done by DFE
|
|
c ===========================================
|
|
c
|
|
else
|
|
c
|
|
c loop over angle points
|
|
c ----------------------
|
|
c
|
|
do i=1,nmuc
|
|
do id=1,nd-1
|
|
dtau(id)=dt(id)/abs(amuc(i))
|
|
end do
|
|
c
|
|
c boundary conditions
|
|
c
|
|
rup=0.
|
|
rdown=0.
|
|
if(amuc(i).lt.0.) rup=extint(ij,i)
|
|
C
|
|
C diffusion approximation for semi-infinite atmospheres
|
|
C
|
|
if(idisk.eq.0.or.ifz0.lt.0) rdown=pland+amuc(i)*dplan
|
|
c
|
|
c the case of finite slab - irradiation of the back side
|
|
c
|
|
if(amuc(i).gt.0.) rdown=rdwn(nmuc-i+1)
|
|
c
|
|
c solution of the transfer equation
|
|
c
|
|
call rtesol(dtau,st0,rup,rdown,amuc(i),ri,ali)
|
|
ttau=0.
|
|
do id=1,nd
|
|
riid=ri(id)*half
|
|
rad1(id)=rad1(id)+wtmuc(i)*riid
|
|
ali1(id)=ali1(id)+wtmuc(i)*ali(id)
|
|
rdh(id)=rdh(id)+amuc1(i)*riid
|
|
rdk(id)=rdk(id)+amuc2(i)*riid
|
|
rdn(id)=rdn(id)+amuc3(i)*riid
|
|
end do
|
|
rdwn(i)=ri(nd)
|
|
if(amuc(i).gt.0.) rdh1=rdh1+amuc1(i)*ri(1)*half
|
|
rdhd=rdhd+abs(amuc1(i))*ri(nd)*half
|
|
end do
|
|
c
|
|
c ----------------------
|
|
c end of the loop over angle points
|
|
c
|
|
end if
|
|
c
|
|
do id=1,nd
|
|
fak1(id)=fak(ij,id)
|
|
radk(ij,id)=rdk(id)
|
|
if(icomve .gt. 0) then
|
|
fkk(id)=rdk(id)/rad1(id)
|
|
else
|
|
fkk(id)=fak(ij,id)
|
|
endif
|
|
ss0(id)=0.
|
|
end do
|
|
if(icomve.gt.0) then
|
|
do id=1,nd
|
|
fak(ij,id)=rdk(id)/rad1(id)
|
|
fak1(id)=fak(ij,id)
|
|
fkk(id)=fak(ij,id)
|
|
end do
|
|
end if
|
|
if(rad1(1).gt.0.) then
|
|
flux(ij)=rdh1
|
|
fhd(ij)=rdhd/rad1(nd)
|
|
end if
|
|
c
|
|
ah=rdh1
|
|
if(iwinbl.lt.0) ah=ah+half*hextrd(ij)
|
|
aj=rad1(1)
|
|
fh(ij)=ah/aj
|
|
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
|
|
C Upper boundary condition
|
|
C
|
|
U0=0.
|
|
QQ0=0.
|
|
US0=0.
|
|
TAUMIN=ABSO1(1)*DEDM1
|
|
NMU=3
|
|
DO I=1,NMU
|
|
IF(IWINBL.EQ.0.AND.WANGLE.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)
|
|
if(tamm.gt.0.) US0=US0+P0/TAMM*WTMU(I)
|
|
END IF
|
|
END DO
|
|
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+FH(IJ)+B)*BQ+SS0(ID)
|
|
CCC(ID)=(FKK(ID+1)/DTP1)*BQ-CQ*(UN+SS0(ID+1))
|
|
ZZZ(ID)=UN/BBB(ID)
|
|
VLL=ST0(ID)+CQ*ST0(ID+1)
|
|
c IF(IWINBL.LT.0) VLL=VLL+TWO*HEXTRD(IJ)/DTP1
|
|
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)
|
|
ALP=UN/DTM1*DT0
|
|
GAM=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*ALP)*SIXTH
|
|
C=(UN-HALF*DTM1*DTM1*GAM)*SIXTH
|
|
END IF
|
|
AAA(ID)=ALP*FKK(ID-1)-A*(UN+SS0(ID-1))
|
|
CCC(ID)=GAM*FKK(ID+1)-C*(UN+SS0(ID+1))
|
|
BBB(ID)=(ALP+GAM)*FKK(ID)+(UN-A-C)*(UN+SS0(ID))
|
|
VLL=A*ST0(ID-1)+C*ST0(ID+1)+(UN-A-C)*ST0(ID)
|
|
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)
|
|
ENDIF
|
|
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
|
|
C
|
|
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)
|
|
frd=bbb(nd)*rad1(nd)-aaa(nd)*rad1(nd-1)
|
|
frd1=(bbb(nd)-un)*rad1(nd)-aaa(nd)*rad1(nd-1)
|
|
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)
|
|
C write(42,642),ij,id,rad1(id),st0(id),fak1(id)
|
|
ALRH(ID)=ZZZ(ID)/(UN-DDD(ID)*EEE(ID+1))
|
|
ALRM(ID)=0
|
|
ALRP(ID)=0
|
|
END DO
|
|
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
|
|
ELSE IF(ILMCOR.EQ.3) THEN
|
|
DO ID=1,ND
|
|
ALI1(ID)=ALI1(ID)/(UN+SS0C(ID)*ALI1(ID))
|
|
ALIM1(ID)=ALIM1(ID)/(UN+SS0C(ID)*ALIM1(ID))
|
|
ALIP1(ID)=ALIP1(ID)/(UN+SS0C(ID)*ALIP1(ID))
|
|
END DO
|
|
END IF
|
|
C
|
|
DO ID=1,ND
|
|
radcm(iji,id)=rad1(id)
|
|
END DO
|
|
C
|
|
C radiation pressure
|
|
C
|
|
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
|
|
c
|
|
if(chmax.ge.1.91e-3.and.chmax.le.2.03e-3) then
|
|
tauij=taumin
|
|
do id=1,nd
|
|
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
|
|
do id=1,nd
|
|
fak(ij,id)=fak1(id)
|
|
end do
|
|
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
|
|
RETURN
|
|
END
|