619 lines
17 KiB
Fortran
619 lines
17 KiB
Fortran
SUBROUTINE BRTE(ID)
|
|
C ===================
|
|
C
|
|
C The part of matrices A,B,C corresponding to the linearized
|
|
C radiative transfer equation
|
|
C i.e. the first NFREQE rows
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ARRAY1.FOR'
|
|
PARAMETER (XCON=8.0935D-21,YCON=1.68638E-10)
|
|
PARAMETER (SIXTH=UN/6.D0,
|
|
* THIRD=UN/3.D0)
|
|
C
|
|
IF(NFREQE.LE.0) RETURN
|
|
ispl=isplin
|
|
if(isplin.ge.5) isplin=isplin-5
|
|
NHE=NFREQE+INHE
|
|
NRE=NFREQE+INRE
|
|
NPC=NFREQE+INPC
|
|
NSE=NFREQE+INSE-1
|
|
NMP=NFREQE+INMP
|
|
C
|
|
GP=0.
|
|
GN=UN
|
|
IF(INMP.GT.0) THEN
|
|
GP=UN
|
|
GN=0.
|
|
END IF
|
|
c
|
|
c in the case of Compton scattering - boundary condition
|
|
c for the highest frequency
|
|
C
|
|
IJ1=1
|
|
if(icompt.gt.0.and.icombc.gt.0.and.ijex(1).gt.0) then
|
|
IJ1=2
|
|
ij=1
|
|
iji=nfreq
|
|
zj1=exp(-hk*freq(ij)/temp(id))
|
|
zj2=exp(-hk*freq(ij+1)/temp(id))
|
|
dlt=delj(iji-1,id)
|
|
if(ichcoo.eq.0) then
|
|
zj0=un/(hk*sqrt(freq(ij)*freq(ij+1))/temp(id))
|
|
zxx=un-3.*zj0+(un-dlt)*zj1+dlt*zj2
|
|
combid=zj0/dlnfr(iji-1)+(un-dlt)*zxx
|
|
comaid=-zj0/dlnfr(iji-1)+dlt*zxx
|
|
else
|
|
e2=ycon*temp(id)
|
|
zxx0=xcon*freq(ij)*(un+zj1)-3.*e2
|
|
zxxm=xcon*freq(ij+1)*(un+zj2)-3.*e2
|
|
zxx=(un-dlt)*zxx0+dlt*zxxm
|
|
combid=e2/dlnfr(iji-1)+(un-dlt)*zxx
|
|
comaid=-e2/dlnfr(iji-1)+dlt*zxx
|
|
end if
|
|
b(ij,ij)=combid
|
|
b(ij,ij+1)=comaid
|
|
vecl(ij)=-b(ij,ij)*rad(iji,id)-b(ij,ij+1)*rad(iji-1,id)
|
|
end if
|
|
C
|
|
C ----------------------------------------
|
|
C For ID = 1 - upper boundary condition
|
|
C ----------------------------------------
|
|
C
|
|
IF(ID.GT.1) GO TO 50
|
|
DDP=(DM(2)-DM(1))*HALF
|
|
DO IJ=IJ1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
OMEG0=ABSO0(IJ)/DENS(ID)
|
|
OMEGP=ABSOP(IJ)/DENS(ID+1)
|
|
DZP=OMEG0+OMEGP
|
|
DTAUP=DZP*DDP
|
|
ALF1=(FK0(IJ)*RAD0(IJ)-FKP(IJ)*RADP(IJ))/DTAUP
|
|
CHIEL0=SCAT0(IJ)
|
|
CHIELP=SCATP(IJ)
|
|
S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
|
|
BS=HALF*DTAUP
|
|
CS=0.
|
|
C2=0.
|
|
GAM2=0.
|
|
BET2=0.
|
|
SP=0.
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.0) then
|
|
call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
|
|
s0=s0+cms
|
|
end if
|
|
C
|
|
IF(MOD(ISPLIN,3).GT.0) THEN
|
|
C
|
|
C Spline collocation and/or Hermitian method (ISPLIN=1 or 2) -
|
|
C both give the same expression for the boundary conditions
|
|
C
|
|
BS=DTAUP*THIRD
|
|
CS=HALF *BS
|
|
SP=(EMISP(IJ)+CHIELP*RADP(IJ))/ABSOP(IJ)
|
|
C2=CS/ABSOP(IJ)
|
|
GAM2=CS*(RADP(IJ)-SP)
|
|
END IF
|
|
C
|
|
C auxiliary quantities
|
|
C
|
|
ALF2=BS*(RAD0(IJ)-S0)
|
|
BET2=ALF2+GAM2
|
|
X1=(ALF1-BET2)/DZP
|
|
B2=(BS+Q0(IJT))/ABSO0(IJ)
|
|
B1=X1/DENS(1)
|
|
B1=B1+UU0(IJT)*S0*DM(1)*HALF/DENS(1)
|
|
C1=X1/DENS(2)
|
|
C
|
|
C *** elements of the IJ-th row of matrices B and C
|
|
C
|
|
RTN=OMEG0*WMM(1)*B1
|
|
B(IJ,NHE)=-GN*RTN
|
|
B1=B1-B2*S0
|
|
C
|
|
RTNC=OMEGP*WMM(2)*C1
|
|
C(IJ,NHE)=-GN*RTNC
|
|
C1=C1-C2*SP
|
|
C
|
|
B(IJ,NRE)=B1*DABT0(IJ)+B2*(DEMT0(IJ)+DST*RAD0(IJ))
|
|
C(IJ,NRE)=C1*DABTP(IJ)+C2*(DEMTP(IJ)+DST*RADP(IJ))
|
|
B(IJ,NPC)=B1*DABN0(IJ)+
|
|
* B2*(DEMN0(IJ)+(DSN+SIGEC(IJT))*RAD0(IJ))+
|
|
* GN*RTN
|
|
C(IJ,NPC)=C1*DABNP(IJ)+
|
|
* C2*(DEMNP(IJ)+(DSN+SIGEC(IJT))*RADP(IJ))+
|
|
* GN*RTNC
|
|
B(IJ,NMP)=B1*DABM0(IJ)+B2*DEMM0(IJ)-GP*RTN
|
|
C(IJ,NMP)=C1*DABMP(IJ)+C2*DEMMP(IJ)-GP*RTNC
|
|
DO II=1,NLVEXP
|
|
B(IJ,NSE+II)=B(IJ,NSE+II)+
|
|
* B1*DRCH0(II,IJ)+B2*DRET0(II,IJ)
|
|
C(IJ,NSE+II)=C(IJ,NSE+II)+
|
|
* C1*DRCHP(II,IJ)+C2*DRETP(II,IJ)
|
|
END DO
|
|
B(IJ,NFREQE)=0.
|
|
B(IJ,IJ)=-FK0(IJ)/DTAUP-FH(IJT)-BS*(UN-CHIEL0/ABSO0(IJ))
|
|
* +Q0(IJT)*CHIEL0/ABSO0(IJ)
|
|
C(IJ,NFREQE)=0.
|
|
C(IJ,IJ)=FKP(IJ)/DTAUP-CS*(UN-CHIELP/ABSOP(IJ))
|
|
C
|
|
C *** the IJ-th element of the rhs vector
|
|
C
|
|
VECL(IJ)=ALF1+BET2+FH(IJT)*RAD0(IJ)
|
|
* -S0*Q0(IJT)
|
|
IF(IWINBL.LT.0) VECL(IJ)=VECL(IJ)-HEXTRD(IJT)
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.4) then
|
|
iji=nfreq-kij(ijt)+1
|
|
b(ij,ij)=b(ij,ij)+bs*(cmb+cme)
|
|
if(iji.gt.1) then
|
|
ijm=ijex(ijorig(iji-1))
|
|
if(ijm.gt.0) b(ij,ijm)=b(ij,ijm)+bs*cma
|
|
end if
|
|
if(iji.lt.nfreq) then
|
|
ijp=ijex(ijorig(iji+1))
|
|
if(ijp.gt.0) b(ij,ijp)=b(ij,ijp)+bs*cmc
|
|
end if
|
|
if(inre.gt.0) b(ij,nre)=b(ij,nre)+cmd*bs
|
|
if(inpc.gt.0) b(ij,npc)=b(ij,npc)+cms*bs/elec(id)
|
|
end if
|
|
c
|
|
END DO
|
|
isplin=ispl
|
|
go to 500
|
|
C
|
|
C ---------------------------------------
|
|
C For 1 < ID < ND - normal depth point
|
|
C ---------------------------------------
|
|
C
|
|
50 DDM=(DM(ID)-DM(ID-1))*HALF
|
|
IF(ID.EQ.ND) GO TO 150
|
|
DDP=(DM(ID+1)-DM(ID))*HALF
|
|
DO IJ=IJ1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
OMEG0=ABSO0(IJ)/DENS(ID)
|
|
OMEGP=ABSOP(IJ)/DENS(ID+1)
|
|
OMEGM=ABSOM(IJ)/DENS(ID-1)
|
|
DZP=OMEG0+OMEGP
|
|
DZM=OMEG0+OMEGM
|
|
DTAUP=DZP*DDP
|
|
DTAUM=DZM*DDM
|
|
DTAU0=HALF *(DTAUP+DTAUM)
|
|
FRD=FK0(IJ)*RAD0(IJ)
|
|
ALF1=(FRD-FKP(IJ)*RADP(IJ))/DTAUP/DTAU0
|
|
GAM1=(FRD-FKM(IJ)*RADM(IJ))/DTAUM/DTAU0
|
|
BET1=ALF1+GAM1
|
|
X1=HALF *BET1/DTAU0
|
|
A1=(GAM1+X1*DTAUM)/DZM
|
|
C1=(ALF1+X1*DTAUP)/DZP
|
|
B1=(A1+C1)/DENS(ID)
|
|
A1=A1/DENS(ID-1)
|
|
C1=C1/DENS(ID+1)
|
|
BS=UN
|
|
CHIELM=SCATM(IJ)
|
|
CHIEL0=SCAT0(IJ)
|
|
CHIELP=SCATP(IJ)
|
|
S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
|
|
AS=0.
|
|
CS=0.
|
|
A2=0.
|
|
C2=0.
|
|
A3=0.
|
|
C3=0.
|
|
BET2=0.
|
|
SM=0.
|
|
SP=0.
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.0) then
|
|
call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
|
|
s0=s0+cms
|
|
end if
|
|
C
|
|
IF(MOD(ISPLIN,3).EQ.0) GO TO 60
|
|
SM=(EMISM(IJ)+RADM(IJ)*CHIELM)/ABSOM(IJ)
|
|
SP=(EMISP(IJ)+RADP(IJ)*CHIELP)/ABSOP(IJ)
|
|
IF(ISPLIN.EQ.1) THEN
|
|
C
|
|
C spline collocation (ISPLIN=1)
|
|
C
|
|
AS=DTAUM/DTAU0*SIXTH
|
|
CS=DTAUP/DTAU0*SIXTH
|
|
BS=0.666666666666667D0
|
|
ALF2=AS*(RADM(IJ)-SM)
|
|
GAM2=CS*(RADP(IJ)-SP)
|
|
BET2=ALF2+GAM2
|
|
X =HALF *BET2/DTAU0
|
|
A2=(GAM2-X*DTAUM)/DZM
|
|
C2=(ALF2-X*DTAUP)/DZP
|
|
ELSE
|
|
C
|
|
C Hermitian method (ISPLIN=2)
|
|
C
|
|
AS=DTAUP*DTAUP/DTAUM/DTAU0
|
|
CS=DTAUM*DTAUM/DTAUP/DTAU0
|
|
AL3=(RADP(IJ)-SP-RAD0(IJ)+S0)*SIXTH
|
|
GA3=(RADM(IJ)-SM-RAD0(IJ)+S0)*SIXTH
|
|
AV=AL3*CS
|
|
CV=GA3*AS
|
|
AS=(UN-HALF *AS)*SIXTH
|
|
CS=(UN-HALF *CS)*SIXTH
|
|
BS=UN-AS-CS
|
|
X=(AV+CV)/DTAU0/4.D0
|
|
A2=(X*DTAUM+HALF *CV-AV)/DZM
|
|
C2=(X*DTAUP+HALF *AV-CV)/DZP
|
|
BET2=AS*(RADM(IJ)-SM)+CS*(RADP(IJ)-SP)
|
|
END IF
|
|
C
|
|
C auxiliary quantities
|
|
C
|
|
B1=B1-(A2+C2)/DENS(ID)
|
|
A1=A1-A2/DENS(ID-1)
|
|
C1=C1-C2/DENS(ID+1)
|
|
A2=AS/ABSOM(IJ)
|
|
C2=CS/ABSOP(IJ)
|
|
A3=A2*SM
|
|
C3=C2*SP
|
|
60 B2=BS/ABSO0(IJ)
|
|
B3=B2*S0
|
|
C
|
|
C *** elements of the IJ-th row of matrices A, B, and C
|
|
C
|
|
RTNA=OMEGM*WMM(ID-1)*A1
|
|
A(IJ,NHE)=-GN*RTNA
|
|
A1=A1-A3
|
|
C
|
|
RTN=OMEG0*WMM(ID)*B1
|
|
B(IJ,NHE)=-GN*RTN
|
|
B1=B1-B3
|
|
C
|
|
RTNC=OMEGP*WMM(ID+1)*C1
|
|
C(IJ,NHE)=-GN*RTNC
|
|
C1=C1-C3
|
|
C
|
|
A(IJ,NRE)= A1*DABTM(IJ)+A2*(DEMTM(IJ)+DST*RADM(IJ))
|
|
B(IJ,NRE)= B1*DABT0(IJ)+B2*(DEMT0(IJ)+DST*RAD0(IJ))
|
|
C(IJ,NRE)= C1*DABTP(IJ)+C2*(DEMTP(IJ)+DST*RADP(IJ))
|
|
A(IJ,NPC)= A1*DABNM(IJ)+
|
|
* A2*(DEMNM(IJ)+(DSN+SIGEC(IJT))*RADM(IJ))+
|
|
* GN*RTNA
|
|
B(IJ,NPC)= B1*DABN0(IJ)+
|
|
* B2*(DEMN0(IJ)+(DSN+SIGEC(IJT))*RAD0(IJ))+
|
|
* GN*RTN
|
|
C(IJ,NPC)= C1*DABNP(IJ)+
|
|
* C2*(DEMNP(IJ)+(DSN+SIGEC(IJT))*RADP(IJ))+
|
|
* GN*RTNC
|
|
A(IJ,NMP)= A1*DABMM(IJ)+A2*DEMMM(IJ)-GP*RTNA
|
|
B(IJ,NMP)= B1*DABM0(IJ)+B2*DEMM0(IJ)-GP*RTN
|
|
C(IJ,NMP)= C1*DABMP(IJ)+C2*DEMMP(IJ)-GP*RTNC
|
|
DO II=1,NLVEXP
|
|
A(IJ,NSE+II)=A(IJ,NSE+II)+
|
|
* A1*DRCHM(II,IJ)+A2*DRETM(II,IJ)
|
|
B(IJ,NSE+II)=B(IJ,NSE+II)+
|
|
* B1*DRCH0(II,IJ)+B2*DRET0(II,IJ)
|
|
C(IJ,NSE+II)=C(IJ,NSE+II)+
|
|
* C1*DRCHP(II,IJ)+C2*DRETP(II,IJ)
|
|
END DO
|
|
A(IJ,NFREQE)=0.
|
|
A(IJ,IJ)=FKM(IJ)/DTAUM/DTAU0-AS*(UN-CHIELM/ABSOM(IJ))
|
|
B(IJ,NFREQE)=0.
|
|
B(IJ,IJ)=-FK0(IJ)/DTAU0*(UN/DTAUP+UN/DTAUM)-
|
|
* BS*(UN-CHIEL0/ABSO0(IJ))
|
|
C(IJ,NFREQE)=0.
|
|
C(IJ,IJ)=FKP(IJ)/DTAUP/DTAU0-CS*(UN-CHIELP/ABSOP(IJ))
|
|
C
|
|
C *** the IJ-th element of the rhs vector
|
|
C
|
|
VECL(IJ)=BET1+BET2+BS*(RAD0(IJ)-S0)
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.4) then
|
|
iji=nfreq-kij(ijt)+1
|
|
b(ij,ij)=b(ij,ij)+bs*(cmb+cme)
|
|
if(iji.gt.1) then
|
|
ijm=ijex(ijorig(iji-1))
|
|
if(ijm.gt.0) b(ij,ijm)=b(ij,ijm)+bs*cma
|
|
end if
|
|
if(iji.lt.nfreq) then
|
|
ijp=ijex(ijorig(iji+1))
|
|
if(ijp.gt.0) b(ij,ijp)=b(ij,ijp)+bs*cmc
|
|
end if
|
|
if(inre.gt.0) b(ij,nre)=b(ij,nre)+cmd*bs
|
|
if(inpc.gt.0) b(ij,npc)=b(ij,npc)+cms*bs/elec(id)
|
|
end if
|
|
c
|
|
END DO
|
|
isplin=ispl
|
|
go to 500
|
|
C
|
|
C --------------------------------------
|
|
C For ID=ND - lower boundary condition
|
|
C --------------------------------------
|
|
C
|
|
150 CONTINUE
|
|
IF(IDISK.EQ.0.OR.IFZ0.LT.0) THEN
|
|
T=TEMP(ID)
|
|
TM=TEMP(ID-1)
|
|
IF(TEMPBD.NE.0.) THEN
|
|
T=TEMPBD
|
|
TM=T
|
|
END IF
|
|
HKT=HK/T
|
|
HKTM=HK/TM
|
|
C
|
|
C auxiliary quantites
|
|
C
|
|
DO IJ=IJ1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
CHIELM=SCATM(IJ)
|
|
CHIEL0=SCAT0(IJ)
|
|
OMEGM=ABSOM(IJ)/DENS(ID-1)
|
|
OMEG0=ABSO0(IJ)/DENS(ID)
|
|
DZM=OMEG0+OMEGM
|
|
DTAUM=DZM*DDM
|
|
FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ)
|
|
GAM1=FRD/DTAUM
|
|
A1=GAM1/DZM
|
|
AS=0.
|
|
BS=0.
|
|
A2=0.
|
|
B2=0.
|
|
A3=0.
|
|
B3=0.
|
|
ALF2=0.
|
|
BET2=0.
|
|
GAM2=0.
|
|
C
|
|
C second-order boundary condition
|
|
C
|
|
IF(IBC.GT.0.AND.IBC.LT.4) THEN
|
|
BS=DTAUM*HALF
|
|
S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.0) then
|
|
call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
|
|
s0=s0+cms
|
|
end if
|
|
C
|
|
GAM2=BS*(RAD0(IJ)-S0)
|
|
BET2=GAM2
|
|
X1=BET2/DZM
|
|
A1=A1-X1
|
|
B2=BS/ABSO0(IJ)
|
|
B3=B2*S0
|
|
END IF
|
|
C
|
|
C auxiliary parameters
|
|
C
|
|
FR=FREQ(IJT)
|
|
FR15=FR*1.D-15
|
|
X=HKT*FR
|
|
EX=EXP(X)
|
|
XM=HKTM*FR
|
|
EXM=EXP(XM)
|
|
PLAN=BN*FR15*FR15*FR15/(EX-UN)*RRDIL
|
|
IF(INRE.EQ.0.OR.ID.GE.NDRE) THEN
|
|
PLANM=BN*FR15*FR15*FR15/(EXM-UN)*RRDIL
|
|
GAM3=(PLAN-PLANM)/DTAUM*THIRD
|
|
A1=A1-GAM3/DZM
|
|
GAM1=GAM1-GAM3
|
|
END IF
|
|
C1=A1
|
|
A1=C1/DENS(ID-1)
|
|
B1=C1/DENS(ID)
|
|
C
|
|
C *** elements of the IJ-th row of matrices A and B
|
|
C
|
|
RTNA=OMEGM*WMM(ID-1)*A1
|
|
A(IJ,NHE)=-GN*RTNA
|
|
A1=A1-A3
|
|
C
|
|
RTN=OMEG0*WMM(ID)*B1
|
|
B(IJ,NHE)=-GN*RTN
|
|
B1=B1-B3
|
|
C
|
|
DPLANM=PLANM*XM/TM/(UN-UN/EXM)
|
|
A(IJ,NRE)=A1*DABTM(IJ)+A2*(DEMTM(IJ)+DST*RADM(IJ))-
|
|
* DPLANM/DTAUM*THIRD
|
|
A(IJ,NPC)=A1*DABNM(IJ)+
|
|
* A2*(DEMNM(IJ)+(DSN+SIGEC(IJT))*RADM(IJ))+
|
|
* GN*RTNA
|
|
BB=HALF+THIRD/DTAUM
|
|
DPLAN=PLAN*X/T/(UN-UN/EX)
|
|
B(IJ,NRE)= B1*DABT0(IJ)+B2*(DEMT0(IJ)+DST*RAD0(IJ))+
|
|
* BB*DPLAN
|
|
B(IJ,NPC)= B1*DABN0(IJ)+
|
|
* B2*(DEMN0(IJ)+(DSN+SIGEC(IJT))*RAD0(IJ))+
|
|
* GN*RTN
|
|
A(IJ,NMP)= A1*DABMM(IJ)+A2*DEMMM(IJ)-GP*RTNA
|
|
B(IJ,NMP)= B1*DABM0(IJ)+B2*DEMM0(IJ)-GP*RTN
|
|
DO II=1,NLVEXP
|
|
A(IJ,NSE+II)=A(IJ,NSE+II)+
|
|
* A1*DRCHM(II,IJ)+A2*DRETM(II,IJ)
|
|
B(IJ,NSE+II)=B(IJ,NSE+II)+
|
|
* B1*DRCH0(II,IJ)+B2*DRET0(II,IJ)
|
|
END DO
|
|
A(IJ,NFREQE)=0.
|
|
A(IJ,IJ)=FKM(IJ)/DTAUM-AS*(UN-CHIELM/ABSOM(IJ))
|
|
B(IJ,NFREQE)=0.
|
|
C
|
|
C *** the IJ-th element of the rhs vector
|
|
C
|
|
IF(IBC.EQ.0.OR.IBC.EQ.4) THEN
|
|
B(IJ,IJ)=B(IJ,IJ)-FK0(IJ)/DTAUM-
|
|
* BS*(UN-CHIEL0/ABSO0(IJ))-HALF
|
|
VECL(IJ)=GAM1+BET2-HALF*(PLAN-RAD0(IJ))
|
|
ELSE
|
|
B(IJ,IJ)=B(IJ,IJ)-FK0(IJ)/DTAUM-
|
|
* BS*(UN-CHIEL0/ABSO0(IJ))-FHD(IJT)
|
|
VECL(IJ)=GAM1+BET2-HALF*PLAN+FHD(IJT)*RAD0(IJ)
|
|
END IF
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.4) then
|
|
iji=nfreq-kij(ijt)+1
|
|
b(ij,ij)=b(ij,ij)+bs*(cmb+cme)
|
|
if(iji.gt.1) then
|
|
ijm=ijex(ijorig(iji-1))
|
|
if(ijm.gt.0) b(ij,ijm)=b(ij,ijm)+bs*cma
|
|
end if
|
|
if(iji.lt.nfreq) then
|
|
ijp=ijex(ijorig(iji+1))
|
|
if(ijp.gt.0) b(ij,ijp)=b(ij,ijp)+bs*cmc
|
|
end if
|
|
if(inre.gt.0) b(ij,nre)=b(ij,nre)+cmd*bs
|
|
if(inpc.gt.0) b(ij,npc)=b(ij,npc)+cms*bs/elec(id)
|
|
end if
|
|
c
|
|
END DO
|
|
C
|
|
ELSE
|
|
C
|
|
C --------------------------------------
|
|
C For ID=ND - lower boundary condition
|
|
C --------------------------------------
|
|
C
|
|
C for disks -
|
|
C lower b.c. expresses just I(taumax,-mu,nu)=I(taumax,+mu,nu)
|
|
C
|
|
DO IJ=IJ1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
CHIELM=SCATM(IJ)
|
|
CHIEL0=SCAT0(IJ)
|
|
OMEGM=ABSOM(IJ)/DENS(ID-1)
|
|
OMEG0=ABSO0(IJ)/DENS(ID)
|
|
DZM=OMEG0+OMEGM
|
|
DTAUM=DZM*DDM
|
|
FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ)
|
|
GAM1=FRD/DTAUM
|
|
A1=GAM1/DZM
|
|
AS=0.
|
|
A2=0.
|
|
A3=0.
|
|
ALF2=0.
|
|
GAM2=0.
|
|
BS=DTAUM*HALF
|
|
S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.0) then
|
|
call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
|
|
s0=s0+cms
|
|
end if
|
|
C
|
|
GAM2=BS*(RAD0(IJ)-S0)
|
|
BET2=ALF2+GAM2
|
|
X1=BET2/DZM
|
|
A1=A1-X1
|
|
B2=BS/ABSO0(IJ)
|
|
B3=B2*S0
|
|
C1=A1
|
|
A1=C1/DENS(ID-1)
|
|
B1=C1/DENS(ID)
|
|
C
|
|
C *** elements of the IJ-th row of matrix A
|
|
C
|
|
RTN=OMEGM*WMM(ID)*A1
|
|
A(IJ,NHE)=-GN*RTN
|
|
A(IJ,NMP)=-GP*RTN
|
|
A1=A1-A3
|
|
A(IJ,NRE)=A1*DABTM(IJ)+A2*(DEMTM(IJ)+DST*RADM(IJ))
|
|
A(IJ,NPC)=A1*DABNM(IJ)+
|
|
* A2*(DEMNM(IJ)+(DSN+SIGEC(IJT))*RADM(IJ))+
|
|
* GN*RTN
|
|
A(IJ,NMP)= A1*DABMM(IJ)+A2*DEMMM(IJ)-GP*RTNA
|
|
DO I=1,NLVEXP
|
|
A(IJ,NSE+I)=A1*DRCHM(I,IJ)+A2*DRETM(I,IJ)
|
|
END DO
|
|
A(IJ,NFREQE)=0.
|
|
A(IJ,IJ)=FKM(IJ)/DTAUM-AS*(UN-CHIELM/ABSOM(IJ))
|
|
C
|
|
C *** elements of the IJ-th row of matrix B
|
|
C
|
|
RTN=OMEG0*WMM(ID)*B1
|
|
B(IJ,NHE)=-GN*RTN
|
|
B(IJ,NMP)=-GP*RTN
|
|
B1=B1-B3
|
|
B(IJ,NRE)=B1*DABT0(IJ)+B2*(DEMT0(IJ)+DST*RAD0(IJ))
|
|
B(IJ,NPC)=B1*DABN0(IJ)+
|
|
* B2*(DEMN0(IJ)+(DSN+SIGEC(IJT))*RAD0(IJ))+
|
|
* GN*RTN
|
|
B(IJ,NMP)= B1*DABM0(IJ)+B2*DEMM0(IJ)-GP*RTN
|
|
DO I=1,NLVEXP
|
|
B(IJ,NSE+I)=B1*DRCH0(I,IJ)+B2*DRET0(I,IJ)
|
|
END DO
|
|
B(IJ,NFREQE)=0.
|
|
B(IJ,IJ)=-FK0(IJ)/DTAUM-BS*(UN-CHIEL0/ABSO0(IJ))
|
|
C
|
|
C *** the IJ-th element of the rhs vector
|
|
C
|
|
VECL(IJ)=GAM1+BET2
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.4) then
|
|
iji=nfreq-kij(ijt)+1
|
|
b(ij,ij)=b(ij,ij)+bs*(cmb+cme)
|
|
if(iji.gt.1) then
|
|
ijm=ijex(ijorig(iji-1))
|
|
if(ijm.gt.0) b(ij,ijm)=b(ij,ijm)+bs*cma
|
|
end if
|
|
if(iji.lt.nfreq) then
|
|
ijp=ijex(ijorig(iji+1))
|
|
if(ijp.gt.0) b(ij,ijp)=b(ij,ijp)+bs*cmc
|
|
end if
|
|
if(inre.gt.0) b(ij,nre)=b(ij,nre)+cmd*bs
|
|
if(inpc.gt.0) b(ij,npc)=b(ij,npc)+cms*bs/elec(id)
|
|
end if
|
|
c
|
|
END DO
|
|
END IF
|
|
isplin=ispl
|
|
500 CONTINUE
|
|
c
|
|
c zeroing radiation field for very low intensities (if required)
|
|
c
|
|
if(radzer.gt.0) then
|
|
c
|
|
C find the peak in nu*rad_nu:
|
|
c
|
|
radsum=0.
|
|
DO IJ=IJ1,NFREQE
|
|
radsum=max(freq(ij)*radex(ij,id),radsum)
|
|
END DO
|
|
C
|
|
C if much smaller than peak in nu*rad_nu, then set to zero:
|
|
C
|
|
DO IJ=IJ1,NFREQE
|
|
if(freq(ij)*radex(ij,id).lt.radzer*radsum) then
|
|
do ii=1,nn0
|
|
a(ij,ii)=0.
|
|
b(ij,ii)=0.
|
|
c(ij,ii)=0.
|
|
end do
|
|
vecl(ij)=0.
|
|
b(ij,ij)=un
|
|
end if
|
|
end do
|
|
end if
|
|
isplin=ispl
|
|
isplin=ispl
|
|
c
|
|
RETURN
|
|
END
|