76 lines
1.9 KiB
Fortran
76 lines
1.9 KiB
Fortran
subroutine rtefe2(dtau,s,rup,rdown,ri)
|
|
C ======================================
|
|
C
|
|
C formal solver of the radiative transfer equation
|
|
C for one frequency, angle, and for completely known source function;
|
|
c original Feautrier (second-order) scheme
|
|
c
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
parameter (one=1.d0)
|
|
dimension dtau(mdepth),s(mdepth),ri(mdepth)
|
|
dimension a(mdepth),b(mdepth),c(mdepth),d(mdepth),
|
|
* f(mdepth),v(mdepth),z(mdepth)
|
|
c
|
|
c set up the global tridiagonal matrix
|
|
c
|
|
c upper boundary condition
|
|
c
|
|
id=1
|
|
cc=two/dtau(id)
|
|
c(id)=cc/dtau(id)
|
|
b(id)=one+cc+c(id)
|
|
a(id)=0.
|
|
v(id)=s(id)+cc*rup
|
|
c
|
|
c normal depth points
|
|
c
|
|
do id=2,nd-1
|
|
dtinv=two/(dtau(id-1)+dtau(id))
|
|
a(id)=dtinv/dtau(id-1)
|
|
c(id)=dtinv/dtau(id)
|
|
b(id)=one+a(id)+c(id)
|
|
v(id)=s(id)
|
|
enddo
|
|
c
|
|
c lower boundary condition
|
|
c
|
|
id=nd
|
|
aa=two/dtau(id-1)
|
|
a(id)=aa/dtau(id-1)
|
|
b(id)=one+aa+a(id)
|
|
if(rdown.eq.0.) b(id)=one+a(id)
|
|
c(id)=0.
|
|
v(id)=s(id)+aa*rdown
|
|
c
|
|
c ---------------------------------------------------
|
|
c solution by elimination
|
|
c 1. forward sweep
|
|
c
|
|
f(1)=(b(1)-c(1))/c(1)
|
|
d(1)=one/(one+f(1))
|
|
z(1)=v(1)/b(1)
|
|
c
|
|
c ii) normal depth points
|
|
c
|
|
do id=2,nd-1
|
|
f(id)=(b(id)-a(id)-c(id)+a(id)*f(id-1)*d(id-1))/c(id)
|
|
d(id)=one/(one+f(id))
|
|
z(id)=(v(id)+a(id)*z(id-1))*d(id)/c(id)
|
|
enddo
|
|
c
|
|
c iii) upper boundary
|
|
c
|
|
id=nd
|
|
z(id)=(v(id)+a(id)*z(id-1))/(b(id)-a(id)*d(id-1))
|
|
c
|
|
c 2. backward elimination
|
|
c
|
|
ri(nd)=z(nd)
|
|
do id=nd-1,1,-1
|
|
ri(id)=ri(id+1)*d(id)+z(id)
|
|
enddo
|
|
c
|
|
return
|
|
end
|