172 lines
4.9 KiB
Fortran
172 lines
4.9 KiB
Fortran
SUBROUTINE RTECOM
|
|
C =================
|
|
C
|
|
C Solution of the radiative transfer equation with Compton scattering
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ITERAT.FOR'
|
|
COMMON/OPTDPT/DT(MDEPTH)
|
|
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 aa(mdepth),bb(mdepth),cc(mdepth),
|
|
* d(mdepth),f(mdepth),z(mdepth),rd(mdepth)
|
|
c
|
|
PRD0=0.
|
|
DO ID=1,ND
|
|
PRADT(ID)=0.
|
|
END DO
|
|
c
|
|
c ---------------------------------------------
|
|
c "1st formal solution" to update Eddington factors
|
|
c ---------------------------------------------
|
|
c
|
|
if(ncfor1.gt.0) then
|
|
do iform=1,ncfor1
|
|
ij0=1
|
|
if(icombc.gt.0) ij0=2
|
|
DO IJ=ij0,NFREQ
|
|
CALL OPACF1(IJ)
|
|
CALL RTECF1(IJ)
|
|
END DO
|
|
if(icombc.gt.0) then
|
|
ij=1
|
|
iji=nfreq
|
|
call rtecf0(ij)
|
|
do id=1,nd
|
|
rad(iji,id)=-rad(iji-1,id)*coma(id)/(comb(id)+bs(id))
|
|
end do
|
|
end if
|
|
end do
|
|
end if
|
|
c
|
|
c -----------------------------------------------
|
|
c coupled solution for the frequency derivatives terms
|
|
c -----------------------------------------------
|
|
c
|
|
c fully coupled treatment - traditional formulation
|
|
c
|
|
if(ncfull.gt.0) then
|
|
do icfull=1,ncfull
|
|
call rtecmc
|
|
c
|
|
c iterative treatment of the derivative terms
|
|
c
|
|
if(ncitot.gt.0) then
|
|
do ictot=1,ncitot
|
|
c
|
|
if(nccoup.gt.0) then
|
|
do iccoup=1,nccoup
|
|
do ij=1,nfreq
|
|
ijo=ijorig(ij)
|
|
fr=freq(ijo)
|
|
call opacf1(ijo)
|
|
call rtecf0(ijo)
|
|
do id=1,nd
|
|
comb(id)=comb(id)+bs(id)
|
|
bb(id)=be(id)+un-comb(id)
|
|
aa(id)=al(id)
|
|
cc(id)=ga(id)
|
|
vl(id)=vl(id)+
|
|
* (coma(id)*gfm(ij,id)+
|
|
* comc(id)*gfp(ij,id))*rad(ij,id)
|
|
end do
|
|
c
|
|
c ----------------
|
|
c forward sweep
|
|
c ----------------
|
|
c
|
|
c i) upper boundary
|
|
c
|
|
f(1)=(bb(1)-cc(1))/cc(1)
|
|
d(1)=un/(un+f(1))
|
|
z(1)=vl(1)/bb(1)
|
|
c
|
|
c ii) normal depth points
|
|
c
|
|
do id=2,nd-1
|
|
f(id)=(bb(id)-aa(id)-cc(id)+
|
|
* aa(id)*f(id-1)*d(id-1))/cc(id)
|
|
d(id)=un/(un+f(id))
|
|
z(id)=(vl(id)+aa(id)*z(id-1))*d(id)/cc(id)
|
|
end do
|
|
c
|
|
c iii) lower boundary
|
|
c
|
|
id=nd
|
|
z(id)=(vl(id)+aa(id)*z(id-1))/(bb(id)-aa(id)*d(id-1))
|
|
c
|
|
c --------------------
|
|
c backward elimination
|
|
c --------------------
|
|
c
|
|
rd(nd)=z(nd)
|
|
do id=nd-1,1,-1
|
|
rd(id)=rd(id+1)*d(id)+z(id)
|
|
end do
|
|
c
|
|
do id=1,nd
|
|
rad(ij,id)=rd(id)
|
|
end do
|
|
end do
|
|
end do
|
|
end if
|
|
c
|
|
c ---------------------------------------------
|
|
c "2nd formal solution" to update Eddington factors
|
|
c ---------------------------------------------
|
|
c
|
|
if(ncfor2.gt.0) then
|
|
do iform=1,ncfor2
|
|
ij0=nfreq
|
|
if(icombc.gt.0) ij0=nfreq-1
|
|
PRD0=0.
|
|
DO ID=1,ND
|
|
PRADT(ID)=0.
|
|
END DO
|
|
DO IJ=1,ij0
|
|
ijo=ijorig(ij)
|
|
CALL OPACF1(IJo)
|
|
CALL RTECF1(IJo)
|
|
END DO
|
|
PRD0=PRD0*PCK
|
|
DO ID=1,ND
|
|
PRADT(ID)=PRADT(ID)*PCK
|
|
END DO
|
|
if(icombc.gt.0) then
|
|
ij=1
|
|
iji=nfreq
|
|
call rtecf0(ij)
|
|
do id=1,nd
|
|
radcm(iji,id)=-radcm(iji-1,id)*coma(id)/
|
|
* (comb(id)+bs(id))
|
|
end do
|
|
flux(1)=radcm(iji,1)*fh(2)
|
|
end if
|
|
c
|
|
do id=1,nd
|
|
DO IJ=1,NFREQ
|
|
rad(ij,id)=radcm(ij,id)
|
|
END DO
|
|
end do
|
|
c
|
|
end do
|
|
end if
|
|
c
|
|
end do
|
|
end if
|
|
c
|
|
c ---------------------------------------------
|
|
c end of formal solutions
|
|
c ---------------------------------------------
|
|
c
|
|
end do
|
|
end if
|
|
return
|
|
end
|