169 lines
4.7 KiB
Fortran
169 lines
4.7 KiB
Fortran
SUBROUTINE RTEDFE
|
|
C =================
|
|
C
|
|
C Solution of the radiative transfer equation - frequency by
|
|
C frequency - for the known source function.
|
|
C
|
|
C The numerical method used:
|
|
c Discontinuous Finite Element (DFE) method
|
|
c Castor, Dykema, Klein, 1992, ApJ 387, 561.
|
|
C
|
|
C Input through blank COMMON block:
|
|
C CH - two-dimensional array absorption coefficient (frequency,
|
|
C depth)
|
|
C ET - emission coefficient (frequency, depth)
|
|
C
|
|
INCLUDE 'PARAMS.FOR'
|
|
INCLUDE 'MODELP.FOR'
|
|
INCLUDE 'SYNTHP.FOR'
|
|
PARAMETER (ONE=1.,TWO=2.,HALF=0.5)
|
|
PARAMETER (TAUREF = 0.6666666666667)
|
|
DIMENSION DT(MDEPTH),ST0(MDEPTH),AB0(MDEPTH),DELDM(MDEPTH),
|
|
* dtau(mdepth),rip(mdepth),rim(mdepth),riup(mdepth),
|
|
* AMU(3),WTMU(3),RINT1(MMU),
|
|
* AMUI(MMU),AMUW(MMU),TAU(MDEPTH),SS0(MDEPTH)
|
|
COMMON/RTEOPA/CH(MFREQ,MDEPTH),ET(MFREQ,MDEPTH),
|
|
* SC(MFREQ,MDEPTH)
|
|
COMMON/EMFLUX/FLUX(MFREQ),FLUXC(MFREQC)
|
|
COMMON/CONSCA/SCC1(mdepth),SCC2(MDEPTH)
|
|
COMMON/REFDEP/IREFD(MFREQ)
|
|
C
|
|
C angle points (AMU) and angular integration weights (WTMU)
|
|
C
|
|
DATA AMU/.887298334620742D0,.5D0,.112701665379258D0/,
|
|
* WTMU/.277777777777778D0,.444444444444444D0,.277777777777778D0/
|
|
C
|
|
DO I=1,ND-1
|
|
DELDM(I)=HALF*(DM(I+1)-DM(I))
|
|
END DO
|
|
C
|
|
c angle points
|
|
C
|
|
IF(IFLUX.EQ.0) THEN
|
|
NMUS=NMU
|
|
do i=1,nmu
|
|
amui(i)=amu(i)
|
|
amuw(i)=amu(i)*wtmu(i)
|
|
end do
|
|
ELSE IF(IFLUX.EQ.1) THEN
|
|
NMUS=NMU0
|
|
do i=1,nmus
|
|
amui(i)=angl(i)
|
|
amuw(i)=angl(i)*wangl(i)
|
|
end do
|
|
END IF
|
|
C
|
|
C overall loop over frequencies
|
|
C
|
|
DO IJ=1,NFREQ
|
|
FR=FREQ(IJ)
|
|
C
|
|
C total source function
|
|
C
|
|
DO ID=1,ND
|
|
AB0(ID)=CH(IJ,ID)
|
|
SCT=FRX1(IJ)*SCC2(ID)+FRX2(IJ)*SCC1(ID)
|
|
ST0(ID)=(ET(IJ,ID)+SCT)/AB0(ID)
|
|
SS0(ID)=-SCT/AB0(ID)
|
|
END DO
|
|
AH=0.
|
|
C
|
|
C optical depth scale
|
|
C
|
|
TAU(1)=0.
|
|
IREF=1
|
|
DO ID=1,ND-1
|
|
DT(ID)=DELDM(ID)*(AB0(ID+1)/DENS(ID+1)+AB0(ID)/DENS(ID))
|
|
TAU(ID+1)=TAU(ID)+DT(ID)
|
|
IF(TAU(ID).LE.TAUREF.AND.TAU(ID+1).GT.TAUREF) IREF=ID
|
|
END DO
|
|
IREFD(IJ)=IREF
|
|
C
|
|
C quantities for the lower boundary condition
|
|
C
|
|
FR15=FR*1.D-15
|
|
BNU=BN*FR15*FR15*FR15
|
|
PLAND=BNU/(EXP(HK*FR/TEMP(ND))-ONE)
|
|
DPLAN=BNU/(EXP(HK*FR/TEMP(ND-1))-ONE)
|
|
DPLAN=(PLAND-DPLAN)/DT(ND-1)
|
|
c
|
|
c loop over angle poits
|
|
c
|
|
DO I=1,NMUS
|
|
do id=1,nd-1
|
|
dtau(id)=dt(id)/amui(i)
|
|
enddo
|
|
C
|
|
c outgoing intensity
|
|
c
|
|
rip(nd)=PLAND+AMUI(I)*DPLAN
|
|
id=nd-1
|
|
dt0=dtau(id)
|
|
dtaup1=dt0+one
|
|
dtau2=dt0*dt0
|
|
bb=two*dtaup1
|
|
cc=dt0*dtaup1
|
|
aa=dtau2+bb
|
|
rim(id+1)=(aa*rip(id+1)-cc*st0(id+1)+dt0*st0(id))/bb
|
|
do id=nd-1,1,-1
|
|
dt0=dtau(id)
|
|
dtaup1=dt0+one
|
|
dtau2=dt0*dt0
|
|
bb=two*dtaup1
|
|
cc=dt0*dtaup1
|
|
aa=one/(dtau2+bb)
|
|
rim(id)=(two*rim(id+1)+dt0*st0(id+1)+cc*st0(id))*aa
|
|
rip(id+1)=(bb*rim(id+1)+cc*st0(id+1)-dt0*st0(id))*aa
|
|
enddo
|
|
do id=2,nd-1
|
|
riup(id)=(rim(id)*dtau(id-1)+rip(id)*dtau(id))/
|
|
* (dtau(id-1)+dtau(id))
|
|
enddo
|
|
riup(1)=rim(1)
|
|
riup(nd)=rip(nd)
|
|
c
|
|
AH=AH+AMUW(I)*RIUP(1)
|
|
RINT1(I)=RIUP(1)
|
|
rint1(i)=max(rint1(i),1.e-40)
|
|
c
|
|
c end of the loop over angle points
|
|
c
|
|
END DO
|
|
c
|
|
FLUX(IJ)=AH*HALF
|
|
if(iflux.ge.1) then
|
|
C
|
|
C output of emergent specific intensities to Unit 10 (line points)
|
|
C or 18 (two continuum points)
|
|
C
|
|
IF(IJ.GT.2) THEN
|
|
WRITE(10,618) WLAM(IJ),FLUX(IJ),(RINT1(IMU),IMU=1,NMUS)
|
|
ELSE
|
|
WRITE(18,618) WLAM(IJ),FLUX(IJ),(RINT1(IMU),IMU=1,NMUS)
|
|
END IF
|
|
end if
|
|
618 FORMAT(1H ,f10.3,1pe15.5/(1P5E15.5))
|
|
C
|
|
C if needed (if iprin.ge.3), output of interesting physical
|
|
C quantities at the monochromatic optical depth tau(nu)=2/3
|
|
C
|
|
IF(IPRIN.GE.3) THEN
|
|
T0=LOG(TAU(IREF+1)/TAU(IREF))
|
|
X0=LOG(TAU(IREF+1)/TAUREF)/T0
|
|
X1=LOG(TAUREF/TAU(IREF))/T0
|
|
DMREF=EXP(LOG(DM(IREF))*X0+LOG(DM(IREF+1))*X1)
|
|
TREF=EXP(LOG(TEMP(IREF))*X0+LOG(TEMP(IREF+1))*X1)
|
|
STREF=EXP(LOG(ST0(IREF))*X0+LOG(ST0(IREF+1))*X1)
|
|
SSREF=EXP(LOG(-SS0(IREF))*X0+LOG(-SS0(IREF+1))*X1)
|
|
SREF=STREF+SSREF
|
|
ALM=2.997925E18/FREQ(IJ)
|
|
WRITE(96,636) IJ,ALM,IREF,DMREF,TREF,STREF,SSREF,SREF
|
|
636 FORMAT(1H ,I3,F10.3,I4,1PE10.3,0PF10.1,1X,1P3E10.3)
|
|
END IF
|
|
C
|
|
C end of the loop over frequencies
|
|
C
|
|
END DO
|
|
RETURN
|
|
END
|