274 lines
7.7 KiB
Fortran
274 lines
7.7 KiB
Fortran
SUBROUTINE BREZ(ID)
|
|
C ===================
|
|
C
|
|
C The part of matrices A and B corresponding to the radiative
|
|
C equilibrium equation
|
|
C i.e. the (NFREQE+INRE)-th row
|
|
C
|
|
C Input: ID - depth index
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ARRAY1.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
DIMENSION REXB(MLEVEL)
|
|
EQUIVALENCE (REX(1),REXB(1))
|
|
C
|
|
NRE=NFREQE+INRE
|
|
NHE=NFREQE+INHE
|
|
NPC=NFREQE+INPC
|
|
NMP=NFREQE+INMP
|
|
NSE=NFREQE+INSE-1
|
|
IJ1=1
|
|
if(icompt.gt.0.and.icombc.gt.0.and.ijex(1).gt.0) IJ1=2
|
|
C
|
|
ittc=abs(nretc)/100
|
|
if(iter.gt.ittc) then
|
|
if(id.le.mod(abs(nretc),100)) then
|
|
b(nre,nre)=1.
|
|
if(nretc.lt.0) then
|
|
c(nre,nre)=-1.
|
|
vecl(nre)=temp(id+1)-temp(id)
|
|
end if
|
|
return
|
|
end if
|
|
end if
|
|
C
|
|
C the rhs vector accounts for total net cooling in ALI
|
|
C transitions (FCOOL)
|
|
C
|
|
VECL(NRE)=FCOOL(ID)-reint(id)*TVISC(ID)
|
|
if(reint(id).le.0) go to 100
|
|
C
|
|
C ********* integral equation part of the radiative
|
|
C equilibrium equation
|
|
C
|
|
BREPC=0.
|
|
BREMP=0.
|
|
DO I=1,NLVEXP
|
|
REXB(I)=0.
|
|
END DO
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=IJ1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
BREPC=BREPC+((DABN0(IJ)-SIGEC(IJT))*RAD0(IJ)-
|
|
* DEMN0(IJ))*WDEP0(IJ)
|
|
BREMP=BREMP+(DABM0(IJ)*RAD0(IJ)-DEMM0(IJ))*WDEP0(IJ)
|
|
DO I=1,NLVEXP
|
|
REXB(I)=REXB(I)+(DRCH0(I,IJ)*RAD0(IJ)-
|
|
* DRET0(I,IJ))*WDEP0(IJ)
|
|
END DO
|
|
B(NRE,NRE)=B(NRE,NRE)+(DABT0(IJ)*RAD0(IJ)-
|
|
* DEMT0(IJ))*WDEP0(IJ)*reint(id)
|
|
HEAT=ABSO0(IJ)-SCAT0(IJ)
|
|
B(NRE,IJ)=WDEP0(IJ)*HEAT*reint(id)
|
|
VECL(NRE)=VECL(NRE)-(HEAT*RAD0(IJ)-EMIS0(IJ))*WDEP0(IJ)*
|
|
* reint(id)
|
|
c
|
|
c additional terms for Compton scattering
|
|
c
|
|
if(icompt.gt.5) then
|
|
ijt=ijfr(ij)
|
|
call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
|
|
vecl(nre)=vecl(nre)+abso0(ij)*cms*wdep0(ij)*reint(id)
|
|
if(icompt.gt.6) then
|
|
if(icmdra.gt.0) then
|
|
b(nre,ij)=b(nre,ij)-abso0(ij)*(cmb+cme)*wdep0(ij)*reint(id)
|
|
else
|
|
b(nre,ij)=b(nre,ij)-abso0(ij)*(cmb+cme)*reint(id)
|
|
end if
|
|
iji=nfreq-kij(ijt)+1
|
|
if(iji.gt.1) then
|
|
ijm=ijex(ijorig(iji-1))
|
|
if(ijm.gt.0) then
|
|
if(icmdra.gt.0) then
|
|
b(nre,ijm)=b(nre,ijm)-abso0(ij)*cma*wdep0(ij)*reint(id)
|
|
else
|
|
b(nre,ijm)=b(nre,ijm)-abso0(ij)*cma*reint(id)
|
|
end if
|
|
end if
|
|
end if
|
|
if(iji.lt.nfreq) then
|
|
ijp=ijex(ijorig(iji+1))
|
|
if(ijp.gt.0) then
|
|
if(icmdra.gt.0) then
|
|
b(nre,ijp)=b(nre,ijp)-abso0(ij)*cmc*wdep0(ij)*reint(id)
|
|
else
|
|
b(nre,ijp)=b(nre,ijp)-abso0(ij)*cmc*reint(id)
|
|
end if
|
|
end if
|
|
end if
|
|
b(nre,nre)=b(nre,nre)-cmd*abso0(ij)*wdep0(ij)*reint(id)
|
|
b(nre,npc)=b(nre,npc)-cms*abso0(ij)/elec(id)*wdep0(ij)*
|
|
* reint(id)
|
|
end if
|
|
end if
|
|
C
|
|
END DO
|
|
END IF
|
|
C
|
|
C corrections for ALI frequency points
|
|
C
|
|
B(NRE,NRE)=B(NRE,NRE)+REIT(ID)*reint(id)
|
|
IF(INPC.GT.0) B(NRE,NPC)=B(NRE,NPC)+(BREPC+REIN(ID))*reint(id)
|
|
IF(INPC.GT.0) B(NRE,NPC)=B(NRE,NPC)+(BREPC+REIN(ID))*reint(id)
|
|
IF(INMP.GT.0) B(NRE,NMP)=B(NRE,NMP)+(BREMP+REIM(ID))*reint(id)
|
|
IF(INHE.GT.0) B(NRE,NHE)=REIX(ID)*reint(id)
|
|
A(NRE,NRE)=AREIT(ID)*reint(id)
|
|
IF(INPC.GT.0) A(NRE,NPC)=AREIN(ID)*reint(id)
|
|
C(NRE,NRE)=CREIT(ID)*reint(id)
|
|
IF(INPC.GT.0) C(NRE,NPC)=CREIN(ID)*reint(id)
|
|
IF(INMP.GT.0) C(NRE,NMP)=CREIM(ID)*reint(id)
|
|
IF(INHE.GT.0) C(NRE,NHE)=CREIX(ID)*reint(id)
|
|
c END IF
|
|
C
|
|
C terms arising because of viscosity
|
|
C
|
|
B(NRE,NRE)=B(NRE,NRE)+DTVIST(ID)*reint(id)
|
|
IF(INPC.GT.0) B(NRE,NPC)=B(NRE,NPC)-DTVISR(ID)*reint(id)
|
|
IF(INHE.GT.0) B(NRE,NFREQE+INHE)=
|
|
* (DTVISR(ID)+DTVISN(ID))*reint(id)
|
|
IF(INMP.GT.0) B(NRE,NFREQE+INMP)=
|
|
* DTVISR(ID)*HMASS/WMM(ID)*reint(id)
|
|
C
|
|
DO II=1,NLVEXP
|
|
B(NRE,NSE+II)=B(NRE,NSE+II)+(REXB(II)+REIP(II,ID))*reint(id)
|
|
END DO
|
|
IF(IFALI.GT.5.AND.ID.GT.1) THEN
|
|
DO II=1,NLVEXP
|
|
A(NRE,NSE+II)=A(NRE,NSE+II)+AREIP(II,ID)*reint(id)
|
|
END DO
|
|
END IF
|
|
IF(IFALI.GT.5.AND.ID.LT.ND) THEN
|
|
DO II=1,NLVEXP
|
|
C(NRE,NSE+II)=C(NRE,NSE+II)+CREIP(II,ID)*reint(id)
|
|
END DO
|
|
END IF
|
|
C
|
|
C ********* differential equation part of the
|
|
C radiative equilibrium equation
|
|
C
|
|
100 CONTINUE
|
|
if(redif(id).eq.0) return
|
|
C
|
|
TEFFD=TEFF**4*(UN-THETAV(ID))
|
|
VECL(NRE)=VECL(NRE)+SIG4P*TEFFD*redif(id)
|
|
if(id.eq.1) go to 200
|
|
C
|
|
DDM=(ZD(ID-1)-ZD(ID))*HALF
|
|
AREN=0.
|
|
BREN=0.
|
|
AREPC=0.
|
|
BREPC=0.
|
|
C
|
|
GP=0.
|
|
GN=UN
|
|
IF(INMP.GT.0) THEN
|
|
GP=UN
|
|
GN=0.
|
|
END IF
|
|
C
|
|
DO I=1,NLVEXP
|
|
REXB(I)=0.
|
|
REXA(I)=0.
|
|
END DO
|
|
C
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
OMEG0=ABSO0(IJ)
|
|
OMEGM=ABSOM(IJ)
|
|
DTAUM=(OMEG0+OMEGM)*DDM
|
|
FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ)
|
|
GAMR=FRD/DTAUM
|
|
A1=GAMR/(OMEG0+OMEGM)*WDEP0(IJ)
|
|
C
|
|
C Corresponding elements of matrix A
|
|
C
|
|
A(NRE,IJ)=-WDEP0(IJ)*FKM(IJ)/DTAUM*redif(id)
|
|
AREPC=AREPC-A1*DABNM(IJ)
|
|
A(NRE,NRE)=A(NRE,NRE)-A1*DABTM(IJ)*redif(id)
|
|
C
|
|
C Corresponding elements of matrix B
|
|
C Columns corresponding to mean intensities
|
|
C
|
|
B(NRE,IJ)=B(NRE,IJ)+WDEP0(IJ)*FK0(IJ)/DTAUM*redif(id)
|
|
BREPC=BREPC-A1*DABN0(IJ)
|
|
C
|
|
C Column corresponding to temperature
|
|
C
|
|
B(NRE,NRE)=B(NRE,NRE)-A1*DABT0(IJ)*redif(id)
|
|
C
|
|
C auxiliary vectors for columns corresponding to populations
|
|
C
|
|
DO I=1,NLVEXP
|
|
REXA(I)=REXA(I)-A1*DRCHM(I,IJ)
|
|
REXB(I)=REXB(I)-A1*DRCH0(I,IJ)
|
|
END DO
|
|
C
|
|
C The rhs vector
|
|
C
|
|
VECL(NRE)=VECL(NRE)-WDEP0(IJ)*GAMR*redif(id)
|
|
END DO
|
|
END IF
|
|
C
|
|
C Column corresponding to temperature
|
|
C
|
|
A(NRE,NRE)=A(NRE,NRE)+REDTM(ID)*REDIF(ID)
|
|
B(NRE,NRE)=B(NRE,NRE)+REDT(ID)*REDIF(ID)
|
|
C(NRE,NRE)=C(NRE,NRE)+REDTP(ID)*REDIF(ID)
|
|
C
|
|
C Column corresponding to electron density (for matrices A and B)
|
|
C
|
|
IF(INPC.NE.0) THEN
|
|
A(NRE,NPC)=A(NRE,NPC)+(AREPC+REDNM(ID))*redif(id)
|
|
B(NRE,NPC)=B(NRE,NPC)+(BREPC+REDN(ID))*redif(id)
|
|
C(NRE,NPC)=C(NRE,NPC)+REDNP(ID)*redif(id)
|
|
END IF
|
|
C
|
|
C Columns corresponding to populations (for matrices A and B)
|
|
C Note that auxiliary arrays REXA and REX, which contain derivatives
|
|
C of the absorption and emission coefficients wrt populations, have
|
|
C been generated by BRTE
|
|
C
|
|
DO II=1,NLVEXP
|
|
A(NRE,NSE+II)=A(NRE,NSE+II)+(REXA(II)+REDPM(II,ID))*redif(id)
|
|
B(NRE,NSE+II)=B(NRE,NSE+II)+(REXB(II)+REDP(II,ID))*redif(id)
|
|
C(NRE,NSE+II)=C(NRE,NSE+II)+REDPP(II,ID)*redif(id)
|
|
END DO
|
|
RETURN
|
|
c
|
|
C upper boundary condition for the differential form
|
|
c
|
|
200 CONTINUE
|
|
C
|
|
C Columns corresponding to mean intensities; rhs vector
|
|
C
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
WF=WDEP0(IJ)*FH(IJT)*REDIF(ID)
|
|
B(NRE,IJ)=B(NRE,IJ)+WF
|
|
VECL(NRE)=VECL(NRE)-WF*RAD0(IJ)
|
|
END DO
|
|
END IF
|
|
C
|
|
C Column corresponding to temperature
|
|
C
|
|
B(NRE,NRE)=B(NRE,NRE)+REDT(ID)*REDIF(ID)
|
|
C
|
|
C Column corresponding to electron density
|
|
C
|
|
IF(INPC.NE.0) THEN
|
|
B(NRE,NPC)=B(NRE,NPC)+REDN(ID)*redif(id)
|
|
END IF
|
|
C
|
|
C Columns corresponding to populations
|
|
C
|
|
DO II=1,NLVEXP
|
|
B(NRE,NSE+II)=B(NRE,NSE+II)+REDP(II,ID)*redif(id)
|
|
END DO
|
|
RETURN
|
|
END
|