453 lines
12 KiB
Fortran
453 lines
12 KiB
Fortran
SUBROUTINE RTECD
|
|
C ================
|
|
C
|
|
C solution of the radiative transfer equation by Feautrier method
|
|
C for two continuum points
|
|
C used when one employs RTEDFE, ie. the DFE method for the
|
|
C transfer equation for the inner frequency points
|
|
C
|
|
INCLUDE 'PARAMS.FOR'
|
|
INCLUDE 'MODELP.FOR'
|
|
INCLUDE 'SYNTHP.FOR'
|
|
DIMENSION D(3,3,MDEPTH),ANU(3,MDEPTH),AANU(MDEPTH),DDD(MDEPTH),
|
|
* AA(3,3),BB(3,3),CC(3,3),VL(3),AMU(3),WTMU(3),
|
|
* DT(MDEPTH),TAU(MDEPTH),
|
|
* RDD(MDEPTH),FKK(MDEPTH),ST0(MDEPTH),SS0(MDEPTH),
|
|
* RINT(MDEPTH,MMU)
|
|
COMMON/RTEOPA/CH(MFREQ,MDEPTH),ET(MFREQ,MDEPTH),
|
|
* SC(MFREQ,MDEPTH)
|
|
COMMON/EMFLUX/FLUX(MFREQ),FLUXC(MFREQC)
|
|
COMMON/CONSCA/SCC1(mdepth),SCC2(MDEPTH)
|
|
PARAMETER (UN=1.D0, HALF=0.5D0)
|
|
PARAMETER (THIRD=UN/3., QUART=UN/4., SIXTH=UN/6.D0)
|
|
PARAMETER (TAUREF = 0.6666666666667)
|
|
DATA AMU/.887298334620742D0,.5D0,.112701665379258D0/,
|
|
1 WTMU/.277777777777778D0,.444444444444444D0,.277777777777778D0
|
|
1 /
|
|
C
|
|
NMU=3
|
|
ND1=ND-1
|
|
C
|
|
C loop over two continuum frequencies
|
|
C
|
|
DO 100 IJ=1,2
|
|
TAUMIN=CH(IJ,1)/DENS(1)*DM(1)*HALF
|
|
TAU(1)=TAUMIN
|
|
DO I=1,ND1
|
|
DT(I)=(DM(I+1)-DM(I))*(CH(IJ,I+1)/DENS(I+1)+CH(IJ,I)/DENS(I))*
|
|
* HALF
|
|
ST0(I)=ET(IJ,I)/CH(IJ,I)
|
|
SS0(I)=-SC(IJ,I)/CH(IJ,I)
|
|
TAU(I+1)=TAU(I)+DT(I)
|
|
IF(TAU(I).LE.TAUREF.AND.TAU(I+1).GT.TAUREF) IREF=I
|
|
END DO
|
|
ST0(ND)=ET(IJ,ND)/CH(IJ,ND)
|
|
SS0(ND)=-SC(IJ,ND)/CH(IJ,ND)
|
|
FR=FREQ(IJ)
|
|
BNU=BN*(FR*1.E-15)**3
|
|
PLAND=BNU/(EXP(HK*FR/TEMP(ND ))-UN)
|
|
DPLAN=BNU/(EXP(HK*FR/TEMP(ND-1))-UN)
|
|
DPLAN=(PLAND-DPLAN)/DT(ND1)
|
|
C
|
|
C +++++++++++++++++++++++++++++++++++++++++
|
|
C FIRST PART - VARIABLE EDDINGTON FACTORS
|
|
C +++++++++++++++++++++++++++++++++++++++++
|
|
C
|
|
C Allowance for wind blanketing
|
|
C
|
|
ALB1=0.
|
|
DO I=1,NMU
|
|
C
|
|
C ************************
|
|
C UPPER BOUNDARY CONDITION
|
|
C ************************
|
|
C
|
|
ID=1
|
|
DTP1=DT(1)
|
|
Q0=0.
|
|
P0=0.
|
|
C
|
|
C allowance for non-zero optical depth at the first depth point
|
|
C
|
|
TAMM=TAUMIN/AMU(I)
|
|
IF(TAMM.GT.0.01) THEN
|
|
P0=UN-EXP(-TAMM)
|
|
ELSE
|
|
P0=TAMM*(UN-HALF*TAMM*(UN-TAMM*THIRD*(UN-QUART*TAMM)))
|
|
END IF
|
|
EX=UN-P0
|
|
Q0=Q0+P0*AMU(I)*WTMU(I)
|
|
C
|
|
DIV=DTP1/AMU(I)*THIRD
|
|
VL(I)=DIV*(ST0(ID)+HALF*ST0(ID+1))+ST0(ID)*P0
|
|
DO J=1,NMU
|
|
BB(I,J)=SS0(ID)*WTMU(J)*(DIV+P0)-ALB1*WTMU(J)
|
|
CC(I,J)=-HALF*DIV*SS0(ID+1)*WTMU(J)
|
|
END DO
|
|
BB(I,I)=BB(I,I)+AMU(I)/DTP1+UN+DIV
|
|
CC(I,I)=CC(I,I)+AMU(I)/DTP1-HALF*DIV
|
|
ANU(I,ID)=0.
|
|
END DO
|
|
C
|
|
C Matrix inversion: instead of calling MATINV, a very fast inlined
|
|
C routine MINV3 for a specific 3 x 3 matrix inversion
|
|
C
|
|
C CALL MATINV(BB,NMU,3)
|
|
C
|
|
C ******************************
|
|
BB(2,1)=BB(2,1)/BB(1,1)
|
|
BB(2,2)=BB(2,2)-BB(2,1)*BB(1,2)
|
|
BB(2,3)=BB(2,3)-BB(2,1)*BB(1,3)
|
|
BB(3,1)=BB(3,1)/BB(1,1)
|
|
BB(3,2)=(BB(3,2)-BB(3,1)*BB(1,2))/BB(2,2)
|
|
BB(3,3)=BB(3,3)-BB(3,1)*BB(1,3)-BB(3,2)*BB(2,3)
|
|
C
|
|
BB(3,2)=-BB(3,2)
|
|
BB(3,1)=-BB(3,1)-BB(3,2)*BB(2,1)
|
|
BB(2,1)=-BB(2,1)
|
|
C
|
|
BB(3,3)=UN/BB(3,3)
|
|
BB(2,3)=-BB(2,3)*BB(3,3)/BB(2,2)
|
|
BB(2,2)=UN/BB(2,2)
|
|
BB(1,3)=-(BB(1,2)*BB(2,3)+BB(1,3)*BB(3,3))/BB(1,1)
|
|
BB(1,2)=-BB(1,2)*BB(2,2)/BB(1,1)
|
|
BB(1,1)=UN/BB(1,1)
|
|
C
|
|
BB(1,1)=BB(1,1)+BB(1,2)*BB(2,1)+BB(1,3)*BB(3,1)
|
|
BB(1,2)=BB(1,2)+BB(1,3)*BB(3,2)
|
|
BB(2,1)=BB(2,2)*BB(2,1)+BB(2,3)*BB(3,1)
|
|
BB(2,2)=BB(2,2)+BB(2,3)*BB(3,2)
|
|
BB(3,1)=BB(3,3)*BB(3,1)
|
|
BB(3,2)=BB(3,3)*BB(3,2)
|
|
C ******************************
|
|
C
|
|
DO I=1,NMU
|
|
DO J=1,NMU
|
|
S=0.
|
|
DO K=1,NMU
|
|
S=S+BB(I,K)*CC(K,J)
|
|
END DO
|
|
D(I,J,ID)=S
|
|
ANU(I,1)=ANU(I,1)+BB(I,J)*VL(J)
|
|
END DO
|
|
END DO
|
|
C
|
|
C *******************
|
|
C NORMAL DEPTH POINTS
|
|
C *******************
|
|
C
|
|
DO ID=2,ND1
|
|
DTM1=DTP1
|
|
DTP1=DT(ID)
|
|
DT0=HALF*(DTM1+DTP1)
|
|
AL=UN/DTM1/DT0
|
|
GA=UN/DTP1/DT0
|
|
BE=AL+GA
|
|
A=(UN-HALF*AL*DTP1*DTP1)*SIXTH
|
|
C=(UN-HALF*GA*DTM1*DTM1)*SIXTH
|
|
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)*WTMU(J)
|
|
CC(I,J)=-C*SS0(ID+1)*WTMU(J)
|
|
BB(I,J)=B*SS0(ID)*WTMU(J)
|
|
END DO
|
|
END DO
|
|
DO I=1,NMU
|
|
DIV=AMU(I)**2
|
|
VL(I)=VL0
|
|
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
|
|
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 Matrix inversion: instead of calling MATINV, a very fast inlined
|
|
C routine MINV3 for a specific 3 x 3 matrix inversion
|
|
C
|
|
C CALL MATINV(BB,NMU,3)
|
|
C
|
|
C ******************************
|
|
BB(2,1)=BB(2,1)/BB(1,1)
|
|
BB(2,2)=BB(2,2)-BB(2,1)*BB(1,2)
|
|
BB(2,3)=BB(2,3)-BB(2,1)*BB(1,3)
|
|
BB(3,1)=BB(3,1)/BB(1,1)
|
|
BB(3,2)=(BB(3,2)-BB(3,1)*BB(1,2))/BB(2,2)
|
|
BB(3,3)=BB(3,3)-BB(3,1)*BB(1,3)-BB(3,2)*BB(2,3)
|
|
C
|
|
BB(3,2)=-BB(3,2)
|
|
BB(3,1)=-BB(3,1)-BB(3,2)*BB(2,1)
|
|
BB(2,1)=-BB(2,1)
|
|
C
|
|
BB(3,3)=UN/BB(3,3)
|
|
BB(2,3)=-BB(2,3)*BB(3,3)/BB(2,2)
|
|
BB(2,2)=UN/BB(2,2)
|
|
BB(1,3)=-(BB(1,2)*BB(2,3)+BB(1,3)*BB(3,3))/BB(1,1)
|
|
BB(1,2)=-BB(1,2)*BB(2,2)/BB(1,1)
|
|
BB(1,1)=UN/BB(1,1)
|
|
C
|
|
BB(1,1)=BB(1,1)+BB(1,2)*BB(2,1)+BB(1,3)*BB(3,1)
|
|
BB(1,2)=BB(1,2)+BB(1,3)*BB(3,2)
|
|
BB(2,1)=BB(2,2)*BB(2,1)+BB(2,3)*BB(3,1)
|
|
BB(2,2)=BB(2,2)+BB(2,3)*BB(3,2)
|
|
BB(3,1)=BB(3,3)*BB(3,1)
|
|
BB(3,2)=BB(3,3)*BB(3,2)
|
|
C ******************************
|
|
C
|
|
DO I=1,NMU
|
|
ANU(I,ID)=0.
|
|
DO J=1,NMU
|
|
S=0.
|
|
DO K=1,NMU
|
|
S=S+BB(I,K)*CC(K,J)
|
|
END DO
|
|
D(I,J,ID)=S
|
|
ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
|
|
END DO
|
|
END DO
|
|
END DO
|
|
C
|
|
C ************
|
|
C LOWER BOUNDARY CONDITION
|
|
C ************
|
|
C
|
|
ID=ND
|
|
DO I=1,NMU
|
|
AA(I,I)=AMU(I)/DTP1
|
|
VL(I)=PLAND+AMU(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
|
|
C
|
|
C Matrix inversion: instead of calling MATINV, a very fast inlined
|
|
C routine MINV3 for a specific 3 x 3 matrix inversion
|
|
C
|
|
C CALL MATINV(BB,NMU,3)
|
|
C
|
|
C ******************************
|
|
BB(2,1)=BB(2,1)/BB(1,1)
|
|
BB(2,2)=BB(2,2)-BB(2,1)*BB(1,2)
|
|
BB(2,3)=BB(2,3)-BB(2,1)*BB(1,3)
|
|
BB(3,1)=BB(3,1)/BB(1,1)
|
|
BB(3,2)=(BB(3,2)-BB(3,1)*BB(1,2))/BB(2,2)
|
|
BB(3,3)=BB(3,3)-BB(3,1)*BB(1,3)-BB(3,2)*BB(2,3)
|
|
C
|
|
BB(3,2)=-BB(3,2)
|
|
BB(3,1)=-BB(3,1)-BB(3,2)*BB(2,1)
|
|
BB(2,1)=-BB(2,1)
|
|
C
|
|
BB(3,3)=UN/BB(3,3)
|
|
BB(2,3)=-BB(2,3)*BB(3,3)/BB(2,2)
|
|
BB(2,2)=UN/BB(2,2)
|
|
BB(1,3)=-(BB(1,2)*BB(2,3)+BB(1,3)*BB(3,3))/BB(1,1)
|
|
BB(1,2)=-BB(1,2)*BB(2,2)/BB(1,1)
|
|
BB(1,1)=UN/BB(1,1)
|
|
C
|
|
BB(1,1)=BB(1,1)+BB(1,2)*BB(2,1)+BB(1,3)*BB(3,1)
|
|
BB(1,2)=BB(1,2)+BB(1,3)*BB(3,2)
|
|
BB(2,1)=BB(2,2)*BB(2,1)+BB(2,3)*BB(3,1)
|
|
BB(2,2)=BB(2,2)+BB(2,3)*BB(3,2)
|
|
BB(3,1)=BB(3,3)*BB(3,1)
|
|
BB(3,2)=BB(3,3)*BB(3,2)
|
|
C ******************************
|
|
C
|
|
DO I=1,NMU
|
|
ANU(I,ID)=0.
|
|
DO J=1,NMU
|
|
D(I,J,ID)=0.
|
|
ANU(I,ID)=ANU(I,ID)+BB(I,J)*VL(J)
|
|
END DO
|
|
END DO
|
|
C
|
|
C ************
|
|
C BACKSOLUTION
|
|
C ************
|
|
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
|
|
AJ=0.
|
|
AK=0.
|
|
DO I=1,NMU
|
|
DIV=WTMU(I)*ANU(I,ID)
|
|
AJ=AJ+DIV
|
|
AK=AK+DIV*AMU(I)**2
|
|
END DO
|
|
FKK(ID)=AK/AJ
|
|
END DO
|
|
C
|
|
C surface Eddington actor
|
|
C
|
|
AH=0.
|
|
DO I=1,NMU
|
|
AH=AH+WTMU(I)*AMU(I)*ANU(I,1)
|
|
END DO
|
|
FH=AH/AJ-HALF*ALB1
|
|
C
|
|
FKK(ND)=THIRD
|
|
C
|
|
C
|
|
C +++++++++++++++++++++++++++++++++++++++++
|
|
C SECOND PART - DETERMINATION OF THE MEAN INTENSITIES
|
|
C RECALCULATION OF THE TRANSFER EQUATION WITH GIVEN EDDINGTON FACTORS
|
|
C +++++++++++++++++++++++++++++++++++++++++
|
|
C
|
|
DTP1=DT(1)
|
|
DIV=DTP1*THIRD
|
|
BBB=FKK(1)/DTP1+FH+DIV+SS0(1)*(DIV+Q0)
|
|
CCC=FKK(2)/DTP1-HALF*DIV*(UN+SS0(2))
|
|
VLL=DIV*(ST0(1)+HALF*ST0(2))+ST0(1)*Q0
|
|
AANU(1)=VLL/BBB
|
|
DDD(1)=CCC/BBB
|
|
DO ID=2,ND1
|
|
DTM1=DTP1
|
|
DTP1=DT(ID)
|
|
DT0=HALF*(DTP1+DTM1)
|
|
AL=UN/DTM1/DT0
|
|
GA=UN/DTP1/DT0
|
|
A=(UN-HALF*DTP1*DTP1*AL)*SIXTH
|
|
C=(UN-HALF*DTM1*DTM1*GA)*SIXTH
|
|
AAA=AL*FKK(ID-1)-A*(UN+SS0(ID-1))
|
|
CCC=GA*FKK(ID+1)-C*(UN+SS0(ID+1))
|
|
BBB=(AL+GA)*FKK(ID)+(UN-A-C)*(UN+SS0(ID))
|
|
VLL=A*ST0(ID-1)+C*ST0(ID+1)+(UN-A-C)*ST0(ID)
|
|
BBB=BBB-AAA*DDD(ID-1)
|
|
DDD(ID)=CCC/BBB
|
|
AANU(ID)=(VLL+AAA*AANU(ID-1))/BBB
|
|
END DO
|
|
BBB=FKK(ND)/DTP1+HALF
|
|
AAA=FKK(ND1)/DTP1
|
|
BBB=BBB-AAA*DDD(ND1)
|
|
VLL=HALF*PLAND+DPLAN*THIRD
|
|
RDD(ND)=(VLL+AAA*AANU(ND1))/BBB
|
|
DO IID=1,ND1
|
|
ID=ND-IID
|
|
RDD(ID)=AANU(ID)+DDD(ID)*RDD(ID+1)
|
|
END DO
|
|
FLUX(IJ)=FH*RDD(1)
|
|
C
|
|
if(ij.eq.1) then
|
|
do id=1,nd
|
|
scc1(id)=-rdd(id)*ss0(id)*ch(1,id)
|
|
end do
|
|
else
|
|
do id=1,nd
|
|
scc2(id)=-rdd(id)*ss0(id)*ch(2,id)
|
|
end do
|
|
end if
|
|
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)
|
|
SCREF=EXP(LOG(-SS0(IREF))*X0+LOG(-SS0(IREF+1))*X1)
|
|
SSREF=EXP(LOG(-SS0(IREF)*RDD(IREF))*X0+
|
|
* LOG(-SS0(IREF+1)*RDD(IREF+1))*X1)
|
|
SREF=STREF+SSREF
|
|
ALM=2.997925E18/FREQ(IJ)
|
|
WRITE(96,636) IJ,ALM,IREF,DMREF,TREF,SCREF,STREF,SSREF,SREF
|
|
636 FORMAT(1H ,I3,F10.3,I4,1PE10.3,0PF10.1,1X,1P3E10.3,E11.3)
|
|
END IF
|
|
c
|
|
C ********************************************************************
|
|
C
|
|
C THIRD PART - DETERMINATION OF THE SPECIFIC INTENSITIES
|
|
C RECALCULATION OF THE TRANSFER EQUATION WITH GIVEN SOURCE FUNCTION
|
|
C
|
|
if(iflux.eq.0) go to 100
|
|
DO IMU=1,NMU0
|
|
ANX=ANGL(IMU)
|
|
DTP1=DT(1)
|
|
DIV=DTP1*THIRD/ANX
|
|
C
|
|
TAMM=TAUMIN/ANX
|
|
IF(TAMM.LT.0.01) THEN
|
|
P0=TAMM*(UN-HALF*TAMM*(UN-TAMM*THIRD*(UN-QUART*TAMM)))
|
|
ELSE
|
|
P0=UN-EXP(-TAMM)
|
|
END IF
|
|
C
|
|
BBB=ANX/DTP1+UN+DIV
|
|
CCC=ANX/DTP1-HALF*DIV
|
|
VLL=(DIV+P0)*(ST0(1)-SS0(1)*RDD(1))
|
|
* +HALF*DIV*(ST0(2)-SS0(2)*RDD(2))
|
|
AANU(1)=VLL/BBB
|
|
DDD(1)=CCC/BBB
|
|
DIV=ANX*ANX
|
|
DO ID=2,ND1
|
|
DTM1=DT(ID-1)
|
|
DTP1=DT(ID)
|
|
DT0=HALF*(DTP1+DTM1)
|
|
AL=UN/DTM1/DT0
|
|
GA=UN/DTP1/DT0
|
|
A=(UN-HALF*DTP1*DTP1*AL)*SIXTH
|
|
C=(UN-HALF*DTM1*DTM1*GA)*SIXTH
|
|
AAA=DIV*AL-A
|
|
CCC=DIV*GA-C
|
|
BBB=DIV*(AL+GA)+UN-A-C
|
|
VLL=A*(ST0(ID-1)-SS0(ID-1)*RDD(ID-1))+
|
|
* C*(ST0(ID+1)-SS0(ID+1)*RDD(ID+1))+
|
|
* (UN-A-C)*(ST0(ID)-SS0(ID)*RDD(ID))
|
|
BBB=BBB-AAA*DDD(ID-1)
|
|
DDD(ID)=CCC/BBB
|
|
AANU(ID)=(VLL+AAA*AANU(ID-1))/BBB
|
|
END DO
|
|
C
|
|
C Lower boundary condition
|
|
C
|
|
AAA=ANX/DTP1
|
|
BBB=AAA+UN
|
|
VLL=PLAND+ANX*DPLAN
|
|
C
|
|
RINT(ND,IMU)=(VLL+AAA*AANU(ND1))/(BBB-AAA*DDD(ND1))
|
|
DO IID=1,ND1
|
|
ID=ND-IID
|
|
RINT(ID,IMU)=AANU(ID)+DDD(ID)*RINT(ID+1,IMU)
|
|
END DO
|
|
END DO
|
|
c
|
|
FLX=0.
|
|
DO IMU=1,NMU0
|
|
RINT(1,IMU)=RINT(1,IMU)/HALF
|
|
FLX=FLX+ANGL(IMU)*WANGL(IMU)*RINT(1,IMU)
|
|
END DO
|
|
FLX=FLX*HALF
|
|
C
|
|
C output of emergent specific intensities in continuum to Unit 18
|
|
C
|
|
if(iflux.ge.1) then
|
|
WRITE(18,641) WLAM(IJ),FLX,(RINT(1,IMU),IMU=1,NMU0)
|
|
end if
|
|
100 CONTINUE
|
|
641 FORMAT(1H ,f10.3,1pe15.5/(1P5E15.5))
|
|
c
|
|
c call rtedfe for the internal points
|
|
c
|
|
CALL RTEDFE
|
|
c
|
|
RETURN
|
|
END
|