358 lines
9.1 KiB
Fortran
358 lines
9.1 KiB
Fortran
SUBROUTINE RTEINT
|
|
C =================
|
|
C
|
|
C Solution of the radiative transfer equation - for one frequency
|
|
C - for the known source function.
|
|
C
|
|
C Solution for specific intensities
|
|
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
|
|
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)
|
|
PARAMETER (MMA=20)
|
|
DIMENSION ST0(MDEPTH),SS0(MDEPTH),AB0(MDEPTH),
|
|
* AA(MMA,MMA),BB(MMA,MMA),CC(MMA,MMA),VL(MMA),
|
|
* FFD(MMA,MMA),FF0D(MMA,MMA),FFPD(MMA,MMA),
|
|
* D(MMA,MMA,MDEPTH),ANU(MMA,MDEPTH)
|
|
DIMENSION tau(mdepth),angl(mma),WANG(MMA)
|
|
COMMON/OPTDPT/DT(MDEPTH)
|
|
C
|
|
nmuf=nmu
|
|
nmu=intens
|
|
do imu=1,nmu
|
|
angl(imu)=0.1+float(imu-1)*0.9/float(nmu-1)
|
|
wang(imu)=0.9/float(nmu-1)
|
|
end do
|
|
wang(1)=wang(1)*0.5
|
|
wang(nmu)=WANG(NMU)*0.5
|
|
C
|
|
DO 500 IJO=1,NFREQ
|
|
IJ=IJO
|
|
IF(ispodf.eq.0) IJ=JIK(IJO)
|
|
IF(IJX(IJ).EQ.-1) GO TO 500
|
|
call opacf1(IJ)
|
|
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)
|
|
SS0(ID)=-SCAT1(ID)/AB0(ID)
|
|
RAD1(ID)=0.
|
|
END DO
|
|
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.
|
|
taumin=absot(1)*dm(1)/2.
|
|
C
|
|
ALB1=0.
|
|
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/ANGL(I)
|
|
EX=EXP(-TAMM)
|
|
P0=UN-EX
|
|
QQ0=QQ0+P0*ANGL(I)*WANG(I)
|
|
U0=U0+EX*WANG(I)
|
|
US0=US0+P0/TAMM*WANG(I)
|
|
END IF
|
|
C
|
|
BI=B/ANGL(I)
|
|
CI=C/ANGL(I)
|
|
VL(I)=(BI+P0)*ST0(ID)+CI*ST0(ID+1)
|
|
IF(IWINBL.LT.0) VL(I)=VL(I)+EXTRAD(IJ)
|
|
DO J=1,NMU
|
|
BB(I,J)=SS0(ID)*WANG(J)*(BI+P0)-ALB1*WANG(J)
|
|
CC(I,J)=-CI*SS0(ID+1)*WANG(J)
|
|
END DO
|
|
BB(I,I)=BB(I,I)+ANGL(I)/DTP1+UN+BI
|
|
CC(I,I)=CC(I,I)+ANGL(I)/DTP1-CI
|
|
ANU(I,ID)=0.
|
|
END DO
|
|
C
|
|
IF(ISPLIN.LE.2) THEN
|
|
CALL MATINV(BB,NMU,MMA)
|
|
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
|
|
CALL MATINV(BB,NMU,MMA)
|
|
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
|
|
ENDIF
|
|
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)*WANG(J)
|
|
CC(I,J)=-C*SS0(ID+1)*WANG(J)
|
|
BB(I,J)=B*SS0(ID)*WANG(J)
|
|
END DO
|
|
END DO
|
|
DO I=1,NMU
|
|
VL(I)=VL0
|
|
DIV=ANGL(I)*ANGL(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,MMA)
|
|
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,MMA)
|
|
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/ANGL(I)
|
|
AI=A/ANGL(I)
|
|
VL(I)=ST0(ID)*BI+ST0(ID-1)*AI
|
|
DO J=1,NMU
|
|
AA(I,J)=-AI*SS0(ID-1)*WANG(J)
|
|
BB(I,J)=BI*SS0(ID)*WANG(J)
|
|
END DO
|
|
AA(I,I)=AA(I,I)+ANGL(I)/DTP1-AI
|
|
BB(I,I)=BB(I,I)+ANGL(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)=ANGL(I)/DTP1
|
|
VL(I)=PLAND+ANGL(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=ANGL(I)/DTP1
|
|
B=HALF/A
|
|
AA(I,I)=A
|
|
VL(I)=B*ST0(ID)+PLAND+ANGL(I)*DPLAN+AA(I,I)*ANU(I,ID-1)
|
|
DO J=1,NMU
|
|
BB(I,J)=B*SS0(ID)*WANG(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,MMA)
|
|
C
|
|
DO 230 I=1,NMU
|
|
ANU(I,ID)=0.
|
|
DO 230 J=1,NMU
|
|
D(I,J,ID)=0.
|
|
ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
|
|
230 CONTINUE
|
|
C
|
|
C ***************** Backsolution
|
|
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
|
|
END DO
|
|
c
|
|
sum=0.
|
|
sua=0.
|
|
do imu=1,nmu
|
|
sum=sum+anu(imu,1)*angl(imu)*wang(imu)
|
|
sua=sua+angl(imu)*wang(imu)
|
|
end do
|
|
C
|
|
wlam=2.997925e18/freq(ij)
|
|
WRITE(18,641) WLAM,flux(ij),sum,sua,(2.*ANU(IMU,1),IMU=1,NMU)
|
|
641 format(f11.3,(1p13e11.3))
|
|
c
|
|
500 continue
|
|
nmu=nmuf
|
|
c
|
|
return
|
|
end
|