595 lines
16 KiB
Fortran
595 lines
16 KiB
Fortran
SUBROUTINE RTE
|
|
C
|
|
C solution of the radiative transfer equation by Feautrier method
|
|
C
|
|
INCLUDE 'PARAMS.FOR'
|
|
INCLUDE 'MODELP.FOR'
|
|
INCLUDE 'SYNTHP.FOR'
|
|
INCLUDE 'LINDAT.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)
|
|
CHARACTER*4 TYPION(9)
|
|
COMMON/RTEOPA/CH(MFREQ,MDEPTH),ET(MFREQ,MDEPTH),
|
|
* SC(MFREQ,MDEPTH)
|
|
COMMON/EMFLUX/FLUX(MFREQ),FLUXC(MFREQC)
|
|
COMMON/BLAPAR/RELOP,SPACE0,CUTOF0,TSTD,DSTD,ALAMC
|
|
COMMON/CTRFUN/CINT1(MDEPTH),CINT2(MDEPTH),
|
|
* CTRI(MDEPTH),CTRR(MDEPTH),XKAR(MDEPTH),
|
|
* ABXLI(MFREQ),EMXLI(MFREQ),IJCTR(MFREQ)
|
|
COMMON/REFDEP/IREFD(MFREQ)
|
|
COMMON/CENTRL/ZND,IFZ0
|
|
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 /
|
|
DATA TYPION /' I ',' II ',' III',' IV ',' V ',
|
|
* ' VI ',' VII','VIII',' IX '/
|
|
C
|
|
NMU=3
|
|
ND1=ND-1
|
|
C
|
|
C Overall loop over frequencies
|
|
C
|
|
DO IJ=1,NFREQ
|
|
TAUMIN=CH(IJ,1)/DENS(1)*DM(1)*HALF
|
|
TAU(1)=TAUMIN
|
|
IREF=1
|
|
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
|
|
IREFD(IJ)=IREF
|
|
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
|
|
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
|
|
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.EQ.0) THEN
|
|
B=DTP1*HALF
|
|
A=0.
|
|
DO I=1,NMU
|
|
BI=B/AMU(I)
|
|
AI=A/AMU(I)
|
|
VL(I)=ST0(ID)*BI+ST0(ID-1)*AI
|
|
DO J=1,NMU
|
|
AA(I,J)=-AI*SS0(ID-1)*WTMU(J)
|
|
BB(I,J)=BI*SS0(ID)*WTMU(J)
|
|
END DO
|
|
AA(I,I)=AA(I,I)+AMU(I)/DTP1-AI
|
|
BB(I,I)=BB(I,I)+AMU(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
|
|
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
|
|
END IF
|
|
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
|
|
ID=ND
|
|
FKK(ND)=THIRD
|
|
AJ=0.
|
|
AK=0.
|
|
DO I=1,NMU
|
|
RMU=WTMU(I)*ANU(I,ID)
|
|
AJ=AJ+RMU
|
|
AK=AK+RMU*AMU(I)*AMU(I)
|
|
END DO
|
|
RDD(ID)=AJ
|
|
FKK(ND)=AK/AJ
|
|
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
|
|
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
|
|
C
|
|
C Lower boundary condition
|
|
C 1.option - different from stellar atmospheres
|
|
C
|
|
IF(IFZ0.EQ.0) THEN
|
|
B=DTP1*HALF
|
|
BBB=FKK(ND)/DTP1+B*(UN+SS0(ND))
|
|
AAA=FKK(ND-1)/DTP1
|
|
VLL=B*ST0(ND)
|
|
ELSE
|
|
C
|
|
C Lower boundary condition
|
|
C 2.option - stellar atmospheric
|
|
C
|
|
BBB=FKK(ND)/DTP1+HALF
|
|
AAA=FKK(ND1)/DTP1
|
|
VLL=HALF*PLAND+DPLAN*THIRD
|
|
END IF
|
|
BBB=BBB-AAA*DDD(ND1)
|
|
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
|
|
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 THIRD PART - DETERMINATION OF THE SPECIFIC INTENSITIES
|
|
C RECALCULATION OF THE TRANSFER EQUATION WITH GIVEN SOURCE FUNCTION
|
|
C
|
|
if(iflux.eq.0) return
|
|
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 1.option - different from stellar atmospheres
|
|
C
|
|
IF(IFZ0.EQ.0) THEN
|
|
B=DTP1*HALF/ANX
|
|
BBB=ANX/DTP1+B*(UN+SS0(ND))
|
|
AAA=ANX/DTP1
|
|
VLL=B*ST0(ND)
|
|
ELSE
|
|
C
|
|
C Lower boundary condition
|
|
C 2.option - stellar atmospheric
|
|
C
|
|
AAA=ANX/DTP1
|
|
BBB=AAA+UN
|
|
VLL=PLAND+ANX*DPLAN
|
|
END IF
|
|
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 FLUX(IJ)=FLX
|
|
C
|
|
C output of emergent specific intensities to Unit 10
|
|
C and 18 (continuum)
|
|
C
|
|
IF(IJ.GT.2) THEN
|
|
WRITE(10,641) WLAM(IJ),FLX,(RINT(1,IMU),IMU=1,NMU0)
|
|
ELSE
|
|
WRITE(18,641) WLAM(IJ),FLX,(RINT(1,IMU),IMU=1,NMU0)
|
|
END IF
|
|
c
|
|
if(iprin.eq.4) then
|
|
c
|
|
c compute contribution function C_i (ctri) and C_r (ctrr)
|
|
c following Magain (1986, A&A 163, 135)
|
|
c
|
|
if(ijctr(ij).gt.0) then
|
|
xfr0=(freq(ij)-freq(2))/(freq(1)-freq(2))
|
|
tauc=ch(1,1)/dens(1)*dm(1)*half
|
|
do id=1,nd
|
|
chc1=ch(1,id)
|
|
chc2=ch(2,id)
|
|
chcc=chc2+xfr0*(chc1-chc2)
|
|
etc1=et(1,id)
|
|
etc2=et(2,id)
|
|
etcc=etc2+xfr0*(etc1-etc2)
|
|
stcc=etcc/chcc
|
|
cint=cint2(id)+xfr0*(cint1(id)-cint2(id))
|
|
avx=(chc1+chc2)*0.5*relop
|
|
call linop(id,abxli,emxli,avx)
|
|
sli0=emxli(ij)/abxli(ij)
|
|
abt0=ch(ij,id)
|
|
emt0=et(ij,id)
|
|
stt0=emt0/abt0
|
|
Xkar(id)=abxli(ij)+chcc*stcc/cint
|
|
ctri(id)=tauc*abt0/chc1*stt0*exp(-tau(id))
|
|
if(tau(id).gt.70.) ctri(id)=0.
|
|
ctrr(id)=tauc/chc1*abxli(ij)*(un-sli0/cint)
|
|
if(id.lt.nd) then
|
|
dtc=(ch(1,id+1)/dens(id+1)+ch(1,id)/dens(id))
|
|
tauc=tauc+half*dtc*(dm(id+1)-dm(id))
|
|
endif
|
|
end do
|
|
taurs=Xkar(1)/dens(1)*dm(1)*half
|
|
xcti=ctri(1)*half*(dm(2)-dm(1))
|
|
xctr=ctrr(1)*half*(dm(2)-dm(1))
|
|
do i=1,nd-1
|
|
ctrr(i)=ctrr(i)*exp(-taurs)
|
|
if(i.eq.1) xctr=xctr*exp(-taurs)
|
|
if(i.gt.1) then
|
|
xcti=xcti+ctri(i)*half*(dm(i+1)-dm(i-1))
|
|
xctr=xctr+ctrr(i)*half*(dm(i+1)-dm(i-1))
|
|
endif
|
|
if(taurs.gt.70.) ctrr(i)=0.
|
|
dtrs=(dm(i+1)-dm(i))*(Xkar(i+1)/dens(i+1)+Xkar(i)/dens(i))
|
|
taurs=taurs+half*dtrs
|
|
end do
|
|
ctrr(nd)=0.
|
|
alam=2.997925d18/freq(ij)
|
|
il0=ijctr(ij)
|
|
il=indlin(il0)
|
|
iat=indat(il)/100
|
|
ion=mod(indat(il),100)
|
|
write(97,376) il,alam,typat(iat),typion(ion),iref,dmref,tref
|
|
376 format(i5,f11.4,2x,2a4,i8,1pe12.4,0pf10.1)
|
|
do id=1,nd
|
|
ctrip=ctri(id)/xcti
|
|
ctrrp=ctrr(id)/xctr
|
|
write(97,377) id,dm(id),tau(id),ctrip,ctrrp
|
|
377 format(i4,1p4e12.4)
|
|
end do
|
|
else if(ij.eq.1) then
|
|
do id=1,nd
|
|
cint1(id)=rint(id,nmu0)
|
|
end do
|
|
else if(ij.eq.2) then
|
|
do id=1,nd
|
|
cint2(id)=rint(id,nmu0)
|
|
end do
|
|
endif
|
|
endif
|
|
641 FORMAT(1H ,f10.3,1pe15.5/(1P5E15.5))
|
|
c
|
|
c end of the global loop over frequencies
|
|
c
|
|
END DO
|
|
RETURN
|
|
END
|