391 lines
12 KiB
Fortran
391 lines
12 KiB
Fortran
SUBROUTINE RYBMAT(IJ)
|
|
C =====================
|
|
C
|
|
C evaluation of the complete-linearization matrices in the Rybicki
|
|
c formalism
|
|
c
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ARRAY1.FOR'
|
|
COMMON/RYBMTX/RA(MDEPTH),RB(MDEPTH),RC(MDEPTH),VR(MDEPTH),
|
|
* UA(MDEPTH),UB(MDEPTH),UC(MDEPTH),
|
|
* VA(MDEPTH),VB(MDEPTH),VC(MDEPTH),WR(MDEPTH),
|
|
* WM(MDEPTH,MDEPTH)
|
|
common/dsctva/dsct1(mdepth),dscn1(mdepth)
|
|
PARAMETER (SIXTH=UN/6.D0,
|
|
* THIRD=UN/3.D0,
|
|
* TWOTHR=TWO/3.D0)
|
|
C
|
|
IJT=IJFR(IJ)
|
|
C
|
|
C ==============================================================
|
|
C 1. components corresponding to the radiative transfer equation
|
|
C ==============================================================
|
|
C
|
|
C -----------------------------------
|
|
C ID = 1 - upper boundary condition
|
|
C
|
|
ID=1
|
|
DDM=(DM(ID+1)-DM(ID))*HALF
|
|
DTM=UN/((ABSO1(ID)+ABSO1(ID+1))*DDM)
|
|
DTM2=DTM*DTM
|
|
ALF=DTM*DDM
|
|
FD=TWO*FH(IJ)
|
|
EXTI=EXTRAD(IJ)
|
|
BET=(EXTI-FD*RAD1(ID))*DTM
|
|
GAM=(FAK1(ID)*RAD1(ID)-FAK1(ID+1)*RAD1(ID+1))*TWO*DTM2
|
|
S0=(EMIS1(ID)+SCAT1(ID)*RAD1(ID))/ABSO1(ID)
|
|
C1=ALF*(TWO*GAM-BET)
|
|
B1=C1-S0/ABSO1(ID)
|
|
UNQ=UN+TWO*DTM*Q0(IJ)
|
|
c unq=un
|
|
RB(ID)=-(UN+DTM*(FD+TWO*FAK1(ID)*DTM))+
|
|
* SCAT1(ID)/ABSO1(ID)*UNQ
|
|
RC(ID)= TWO*FAK1(ID+1)*DTM2
|
|
VR(ID)= GAM-BET+RAD1(ID)-S0*UNQ
|
|
UB(ID)= B1*DABT1(ID)+(DEMT1(ID)+
|
|
* DSCT1(ID)*RAD1(ID))/ABSO1(ID)*UNQ-
|
|
* emis1(id)/abso1(id)**2*dabt1(id)*two*dtm*q0(ij)+
|
|
* dtm*s0*(uu0(ij)*dm(1)*dabt1(id)-
|
|
* two*q0(ij)*dtm*ddm*dabt1(id))
|
|
UC(ID)= C1*DABT1(ID+1)-two*dtm*q0(ij)*s0*dtm*ddm*dabt1(id+1)
|
|
c
|
|
if(iubc.gt.0) then
|
|
corf=half/dtm
|
|
rb(id)=rb(id)*corf
|
|
rc(id)=rc(id)*corf
|
|
vr(id)=vr(id)*corf
|
|
c1=(gam-rad1(id)+s0)*corf*alf
|
|
b1=c1-corf*s0/abso1(id)
|
|
UB(ID)= B1*DABT1(ID)+(DEMT1(ID)+
|
|
* DSCT1(ID)*RAD1(ID))/ABSO1(ID)*corf
|
|
UC(ID)= C1*DABT1(ID+1)
|
|
end if
|
|
C
|
|
C
|
|
C ----------------------------------
|
|
C 1 < ID < ND - normal depth point
|
|
C
|
|
DO ID=2,ND-1
|
|
DDM=(DM(ID)-DM(ID-1))*HALF
|
|
DDP=(DM(ID+1)-DM(ID))*HALF
|
|
DZM=ABSO1(ID)+ABSO1(ID-1)
|
|
DZP=ABSO1(ID)+ABSO1(ID+1)
|
|
DTAUP=DZP*DDP
|
|
DTAUM=DZM*DDM
|
|
DTAU0=HALF *(DTAUP+DTAUM)
|
|
FRD=FAK1(ID)*RAD1(ID)
|
|
ALF1=(FRD-FAK1(ID+1)*RAD1(ID+1))/DTAUP/DTAU0
|
|
GAM1=(FRD-FAK1(ID-1)*RAD1(ID-1))/DTAUM/DTAU0
|
|
BET1=ALF1+GAM1
|
|
X1=HALF *BET1/DTAU0
|
|
A1=(GAM1+X1*DTAUM)/DZM
|
|
C1=(ALF1+X1*DTAUP)/DZP
|
|
B1=A1+C1
|
|
BS=UN
|
|
CHIELM=SCAT1(ID-1)
|
|
CHIEL0=SCAT1(ID)
|
|
CHIELP=SCAT1(ID+1)
|
|
S0=(EMIS1(ID)+CHIEL0*RAD1(ID))/ABSO1(ID)
|
|
AS=0.
|
|
CS=0.
|
|
A2=0.
|
|
C2=0.
|
|
A3=0.
|
|
C3=0.
|
|
BET2=0.
|
|
SM=0.
|
|
SP=0.
|
|
IF(MOD(ISPLIN,3).EQ.0) GO TO 60
|
|
SM=(EMIS1(ID-1)+RAD1(ID-1)*CHIELM)/ABSO1(ID-1)
|
|
SP=(EMIS1(ID+1)+RAD1(ID+1)*CHIELP)/ABSO1(ID+1)
|
|
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*(RAD1(ID-1)-SM)
|
|
GAM2=CS*(RAD1(ID+1)-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=(RAD1(ID+1)-SP-RAD1(ID)+S0)*SIXTH
|
|
GA3=(RAD1(ID-1)-SM-RAD1(ID)+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*(RAD1(ID-1)-SM)+CS*(RAD1(ID+1)-SP)
|
|
END IF
|
|
C
|
|
C auxiliary quantities
|
|
C
|
|
B1=B1-(A2+C2)
|
|
A1=A1-A2
|
|
C1=C1-C2
|
|
A2=AS/ABSO1(ID-1)
|
|
C2=CS/ABSO1(ID+1)
|
|
A3=A2*SM
|
|
C3=C2*SP
|
|
C
|
|
60 CONTINUE
|
|
B2=BS/ABSO1(ID)
|
|
B3=B2*S0
|
|
A1=A1-A3
|
|
B1=B1-B3
|
|
C1=C1-C3
|
|
C
|
|
C *** elements of the matrices
|
|
C
|
|
RA(ID)= FAK1(ID-1)/DTAUM/DTAU0-AS*(UN-CHIELM/ABSO1(ID-1))
|
|
RB(ID)=-FAK1(ID)/DTAU0*(UN/DTAUP+UN/DTAUM)-
|
|
* BS*(UN-CHIEL0/ABSO1(ID))
|
|
RC(ID)= FAK1(ID+1)/DTAUP/DTAU0-CS*(UN-CHIELP/ABSO1(ID+1))
|
|
VR(ID)= BET1+BET2+BS*(RAD1(ID)-S0)
|
|
UA(ID)= A1*DABT1(ID-1)+A2*(DEMT1(ID-1)+DSCT1(ID-1)*RAD1(ID-1))
|
|
UB(ID)= B1*DABT1(ID)+B2*(DEMT1(ID)+DSCT1(ID)*RAD1(ID))
|
|
UC(ID)= C1*DABT1(ID+1)+C2*(DEMT1(ID+1)+DSCT1(ID+1)*RAD1(ID+1))
|
|
END DO
|
|
C
|
|
C ----------------------------------
|
|
C ID=ND - lower boundary condition
|
|
C
|
|
ID=ND
|
|
DDM=HALF*(DM(ID)-DM(ID-1))
|
|
T0=TEMP(ID)
|
|
TM=TEMP(ID-1)
|
|
IF(IBC.GT.0.AND.IBC.LT.4.AND.IDISK.EQ.0) THEN
|
|
DTM=UN/((ABSO1(ID-1)+ABSO1(ID))*DDM)
|
|
DTM2=DTM*DTM
|
|
FD=TWO*FHD(IJT)
|
|
FR=FREQ(IJT)
|
|
FR15=FR*1.D-15
|
|
BNU=BN*FR15*FR15*FR15
|
|
X0=HK*FR/T0
|
|
XM=HK*FR/TM
|
|
PLAND=BNU/(EXP(X0)-UN)
|
|
PLANM=BNU/(EXP(XM)-UN)
|
|
DPLDT0=PLAND/(UN-EXP(-X0))*X0/T0
|
|
DPLDTM=PLANM/(UN-EXP(-XM))*XM/TM
|
|
DPLAN=(PLAND-PLANM)*DTM
|
|
ALF=DTM*DDM
|
|
BET=(PLAND-FD*RAD1(ID))*DTM
|
|
GAM=(FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1)-
|
|
* THIRD*(PLAND-PLANM))*TWO*DTM2
|
|
S0=(EMIS1(ID)+SCAT1(ID)*RAD1(ID))/ABSO1(ID)
|
|
A1=ALF*(TWO*GAM-BET)
|
|
B1=A1-S0/ABSO1(ID)
|
|
RA(ID)= TWO*FAK1(ID-1)*DTM2
|
|
RB(ID)=-(UN+DTM*(FD+TWO*FAK1(ID)*DTM))+
|
|
* SCAT1(ID)/ABSO1(ID)
|
|
VR(ID)= GAM-BET+RAD1(ID)-S0
|
|
UA(ID)= B1*DABT1(ID-1)
|
|
* +(DEMT1(ID-1)+
|
|
* DSCT1(ID-1)*RAD1(ID-1))/ABSO1(ID-1)
|
|
* -DPLDTM*DTM2*TWOTHR
|
|
UB(ID)= B1*DABT1(ID)+(DEMT1(ID)+
|
|
* DSCT1(ID)*RAD1(ID))/ABSO1(ID)+
|
|
* DPLDT0*DTM*(UN+TWOTHR*DTM)
|
|
C
|
|
if(ifryb.gt.0) then
|
|
DTM=UN/((ABSO1(ID-1)+ABSO1(ID))*DDM)
|
|
FR=FREQ(IJT)
|
|
FR15=FR*1.D-15
|
|
BNU=BN*FR15*FR15*FR15
|
|
X0=HK*FR/T0
|
|
XM=HK*FR/TM
|
|
PLAND=BNU/(EXP(X0)-UN)
|
|
PLANM=BNU/(EXP(XM)-UN)
|
|
DPLDT0=PLAND/(UN-EXP(-X0))*X0/T0
|
|
DPLDTM=PLANM/(UN-EXP(-XM))*XM/TM
|
|
GAM=(FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1)-
|
|
* THIRD*(PLAND-PLANM))*DTM
|
|
S0=(EMIS1(ID)+SCAT1(ID)*RAD1(ID))/ABSO1(ID)
|
|
BS=HALF/DTM
|
|
BET=BS*(RAD1(ID)-S0)
|
|
A1=(GAM-BET)*DTM*DDM
|
|
B1=A1-BS*S0/ABSO1(ID)
|
|
RA(ID)=FAK1(ID-1)*DTM
|
|
RB(ID)=-FAK1(ID)*DTM-BS*(UN-SCAT1(ID)/ABSO1(ID))-FHD(IJT)
|
|
VR(ID)=GAM+BET-HALF*PLAND+FHD(IJT)*RAD1(ID)
|
|
UA(ID)=A1*DABT1(ID-1)-DPLDTM*DTM*THIRD
|
|
UB(ID)=B1*DABT1(ID)+
|
|
* BS*(DEMT1(ID)+DSCT1(ID)*RAD1(ID))/ABSO1(ID)+
|
|
* (HALF+THIRD*DTM)*DPLDT0
|
|
end if
|
|
C
|
|
ELSE
|
|
C
|
|
C for disks -
|
|
C lower b.c. expresses just I(taumax,-mu,nu)=I(taumax,+mu,nu)
|
|
C
|
|
DZM=ABSO1(ID)+ABSO1(ID-1)
|
|
FRD=FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1)
|
|
DTAUM=DZM*DDM
|
|
GAM1=FRD/DTAUM
|
|
A1=GAM1/DZM
|
|
AS=0.
|
|
BS=DTAUM*HALF
|
|
S0=(EMIS1(ID)+SCAT1(ID)*RAD1(ID))/ABSO1(ID)
|
|
GAM2=BS*(RAD1(ID)-S0)
|
|
X1=GAM2/DZM
|
|
A1=A1-X1
|
|
B2=BS/ABSO1(ID)
|
|
B1=A1-B2*S0
|
|
RA(ID)=FAK1(ID-1)/DTAUM-AS*(UN-SCAT1(ID-1)/ABSO1(ID-1))
|
|
RB(ID)=-FAK1(ID)/DTAUM-BS*(UN-SCAT1(ID)/ABSO1(ID))
|
|
UA(ID)=A1*DABT1(ID-1)+A2*(DEMT1(ID-1)+DSCT1(ID-1)*RAD1(ID-1))
|
|
UB(ID)=B1*DABT1(ID)+B2*(DEMT1(ID)+DSCT1(ID)*RAD1(ID))
|
|
VR(ID)=GAM1+GAM2
|
|
END IF
|
|
C
|
|
C
|
|
C =====================================================
|
|
C components corresponding to the radiative equilibrium
|
|
C =====================================================
|
|
C
|
|
DO ID=1,ND
|
|
C
|
|
C ********* integral equation part
|
|
C
|
|
IF(REINT(ID).GT.0.) THEN
|
|
HEAT = ABSO1(ID)-SCAT1(ID)
|
|
DHEAT = (DABT1(ID)-DSCT1(ID))*RAD1(ID)
|
|
WDR = W(IJ)*dens(id)*reint(id)
|
|
VB(ID) = HEAT*WDR
|
|
WM(ID,ID)= WM(ID,ID)+(DHEAT-DEMT1(ID))*WDR
|
|
WR(ID) = WR(ID)-(HEAT*RAD1(ID)-EMIS1(ID))*WDR
|
|
END IF
|
|
END DO
|
|
C
|
|
C ********* differential equation part
|
|
C
|
|
ID=1
|
|
IF(REDIF(ID).GT.0.) THEN
|
|
WF=W(IJ)*FH(IJT)*REDIF(ID)
|
|
VB(ID)=VB(ID)+WF
|
|
WR(ID)=WR(ID)-WF*RAD1(ID)
|
|
END IF
|
|
c
|
|
c original variant with 1st-order mid-point differences
|
|
c
|
|
if(icentr.eq.0) then
|
|
nd1=nd
|
|
if(ilbc.ne.0) nd1=nd-1
|
|
DO ID=2,nd1
|
|
IF(REDIF(ID).GT.0.) THEN
|
|
DDM=(DM(ID)-DM(ID-1))*HALF
|
|
OMEG0=ABSO1(ID)
|
|
OMEGM=ABSO1(ID-1)
|
|
DTAUM=(OMEG0+OMEGM)*DDM
|
|
FRD=FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1)
|
|
GAMR=FRD/DTAUM*redif(id)
|
|
A1=GAMR/(OMEG0+OMEGM)
|
|
VA(ID) =-W(IJ)*FAK1(ID-1)/DTAUM*redif(id)
|
|
VB(ID) = VB(ID)+W(IJ)*FAK1(ID)/DTAUM*redif(id)
|
|
WM(ID,ID-1)= WM(ID,ID-1)-A1*W(IJ)*DABT1(ID-1)
|
|
WM(ID,ID) = WM(ID,ID)-A1*W(IJ)*DABT1(ID)
|
|
WR(ID) = WR(ID)-W(IJ)*GAMR
|
|
END IF
|
|
END DO
|
|
c
|
|
if(ilbc.gt.0) then
|
|
id=nd
|
|
IF(REDIF(ID).GT.0.) THEN
|
|
DDM=(DM(ID)-DM(ID-1))*HALF
|
|
DTAUM=(ABSO1(ID)+ABSO1(ID-1))*DDM*3.D0
|
|
FR=FREQ(IJT)
|
|
FR15=FR*1.D-15
|
|
BNU=BN*FR15*FR15*FR15
|
|
X0=HK*FR/TEMP(ID)
|
|
XM=HK*FR/TEMP(ID-1)
|
|
PLAND=BNU/(EXP(X0)-UN)
|
|
PLANM=BNU/(EXP(XM)-UN)
|
|
DPLDT0=PLAND/(UN-EXP(-X0))*X0/TEMP(ID)
|
|
DPLDTM=PLANM/(UN-EXP(-XM))*XM/TEMP(ID-1)
|
|
FLX=(PLAND-PLANM)/DTAUM*redif(id)
|
|
A1=FLX/(ABSO1(ID)+ABSO1(ID-1))
|
|
WM(ID,ID-1)= WM(ID,ID-1)+W(IJ)*(DPLDTM-A1*DABT1(ID-1))
|
|
WM(ID,ID) = WM(ID,ID)+W(IJ)*(DPLDT0-A1*DABT1(ID))
|
|
WR(ID) = WR(ID)-W(IJ)*FLX
|
|
END IF
|
|
end if
|
|
c
|
|
c centered difference variant
|
|
c
|
|
else
|
|
do id=2,nd-1
|
|
if(redif(id).gt.0) then
|
|
wwr=w(ij)*redif(id)
|
|
ddm=half*(dm(id)-dm(id-1))
|
|
ddp=half*(dm(id+1)-dm(id))
|
|
dtm=(abso1(id-1)+abso1(id))*ddm
|
|
dtp=(abso1(id+1)+abso1(id))*ddp
|
|
dt0=dtm+dtp
|
|
frm=fak1(id)*rad1(id)-fak1(id-1)*rad1(id-1)
|
|
frp=fak1(id+1)*rad1(id+1)-fak1(id)*rad1(id)
|
|
alp=dtp/dtm/dt0
|
|
gam=dtm/dtp/dt0
|
|
flx=alp*frm+gam*frp
|
|
c
|
|
c matrix elements
|
|
c
|
|
va(id)=-wwr*fak1(id-1)*alp
|
|
vb(id)=vb(id)+wwr*fak1(id)*(alp-gam)
|
|
vc(id)=wwr*fak1(id+1)*gam
|
|
c
|
|
dmtm=ddm/(dtm*dtm)
|
|
dmtp=ddp/(dtp*dtp)
|
|
rm=dtm/dt0
|
|
rp=dtp/dt0
|
|
wm(id,id-1) = wm(id,id-1)+wwr*dabt1(id-1)*dmtm*
|
|
* (rm*rm*frp-(un+rm)*frm)
|
|
wm(id,id) = wm(id,id)+wwr*dabt1(id)*
|
|
* ((dmtp*rp*rp-dmtm*rp*(un+rm))*frm +
|
|
* (dmtm*rm*rm-dmtp*rm*(un+rp)*frp))
|
|
wr(id)=wr(id)-wwr*flx
|
|
end if
|
|
end do
|
|
c
|
|
id=nd
|
|
IF(REDIF(ID).GT.0.) THEN
|
|
DDM=(DM(ID)-DM(ID-1))*HALF
|
|
DTAUM=(ABSO1(ID)+ABSO1(ID-1))*DDM
|
|
FRD=FAK1(ID)*RAD1(ID)-FAK1(ID-1)*RAD1(ID-1)
|
|
FLX=FRD/DTAUM*redif(id)
|
|
A1=FLX/(ABSO1(ID)+ABSO1(ID-1))
|
|
VA(ID) =-W(IJ)*FAK1(ID-1)/DTAUM*redif(id)
|
|
VB(ID) = VB(ID)+W(IJ)*FAK1(ID)/DTAUM*redif(id)
|
|
WM(ID,ID-1)= WM(ID,ID-1)-A1*W(IJ)*DABT1(ID-1)
|
|
WM(ID,ID) = WM(ID,ID)-A1*W(IJ)*DABT1(ID)
|
|
WR(ID) = WR(ID)-W(IJ)*FLX
|
|
END IF
|
|
end if
|
|
c
|
|
if(nretc.gt.0) then
|
|
do id=1,nretc
|
|
wr(id)=0.
|
|
vb(id)=0.
|
|
va(id)=0.
|
|
wm(id,id)=1.
|
|
if(id.gt.1) wm(id,id-1)=0.
|
|
end do
|
|
end if
|
|
RETURN
|
|
END
|