284 lines
7.3 KiB
Fortran
284 lines
7.3 KiB
Fortran
SUBROUTINE BHEZ(ID)
|
|
C ==================
|
|
C
|
|
C The part of matrices A and B corresponding to the hydrostatic
|
|
C equilibrium equation,
|
|
C i.e. the (NFREQE+INHE)-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'
|
|
COMMON/SURFEX/EXTJ(MFREQ),EXTH(MFREQ)
|
|
C
|
|
NHE=NFREQE+INHE
|
|
NRE=NFREQE+INRE
|
|
NPC=NFREQE+INPC
|
|
NZD=NFREQE+INZD
|
|
NSE=NFREQE+INSE-1
|
|
c
|
|
if(inhe.le.0) return
|
|
IJ1=1
|
|
C
|
|
C *********** Linearized equation for the fictitious massive particle
|
|
C density
|
|
C
|
|
IF(INMP.GT.0) THEN
|
|
NMP=NFREQE+INMP
|
|
B(NMP,NMP)=-UN
|
|
B(NMP,NHE)=UN
|
|
IF(INPC.GT.0) B(NMP,NPC)=-UN
|
|
END IF
|
|
C
|
|
C *********** Linearized hydrostatic equilibrium
|
|
C
|
|
HEXT=0.
|
|
HEXN=0.
|
|
GRD=0.
|
|
FLUXW=0.
|
|
DO I=1,NLVEXP
|
|
HEX(I)=0.
|
|
END DO
|
|
C
|
|
IF(ID.GT.1) GO TO 50
|
|
C
|
|
C *** Upper boundary condition (ID=1)
|
|
C
|
|
C 1. possibility - the same as in stellar atmospheres
|
|
C Basically, linearized eq. (7-10) of Mihalas (1978)
|
|
C
|
|
IF(IBCHE.EQ.0) THEN
|
|
X1=PCK/DENS(ID)
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=IJ1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
IF(.NOT.LSKIP(ID,IJT)) THEN
|
|
FLUXW=W(IJT)*(FH(IJT)*RAD0(IJ)-HEXTRD(IJT))
|
|
GRD=GRD+FLUXW*ABSO0(IJ)
|
|
HEXN=HEXN+FLUXW*DABN0(IJ)
|
|
HEXT=HEXT+FLUXW*DABT0(IJ)
|
|
DO I=1,NLVEXP
|
|
HEX(I)=HEX(I)+FLUXW*DRCH0(I,IJ)
|
|
END DO
|
|
C
|
|
C Columns corresponding to mean intensities
|
|
C
|
|
B(NHE,IJ)=X1*WDEP0(IJ)*FH(IJT)*ABSO0(IJ)
|
|
END IF
|
|
END DO
|
|
END IF
|
|
C
|
|
RTN=X1*WMM(ID)/DENS(ID)*(GRD+FPRD(ID))
|
|
VT0=HALF*VTURB(ID)*VTURB(ID)/DM(ID)*WMM(ID)
|
|
C
|
|
C columns corresponding to total particle density, fictitious
|
|
C massive particle density, temperature, and electron density,
|
|
C respectively
|
|
C
|
|
B(NHE,NHE)=BOLK*TEMP(ID)/DM(ID)-GN*(RTN-VT0)
|
|
IF(INMP.GT.0) B(NHE,NFREQE+INMP)=GP*(VT0-RTN)
|
|
IF(INRE.GT.0) THEN
|
|
B(NHE,NRE)=BOLK*PSI0(NHE)/DM(1)+X1*(HEXT+HEIT(ID))
|
|
C(NHE,NRE)=X1*HEITP(ID)
|
|
END IF
|
|
IF(INPC.GT.0) THEN
|
|
B(NHE,NPC)=X1*(HEXN+HEIN(ID))+GN*(RTN-VT0)
|
|
C(NHE,NPC)=X1*HEINP(ID)
|
|
END IF
|
|
C
|
|
C Columns corresponding to populations
|
|
C
|
|
DO II=1,NLVEXP
|
|
B(NHE,NSE+II)=B(NHE,NSE+II)+X1*(HEX(II)+HEIP(II,ID))
|
|
C(NHE,NSE+II)=C(NHE,NSE+II)+X1*HEIPP(II,ID)
|
|
END DO
|
|
C
|
|
C The rhs vector also accounts for the total radiation pressure in
|
|
C the fixed-option transitions (array FPRD)
|
|
C
|
|
GRAV=QGRAV*ZD(1)
|
|
VECL(NHE)=GRAV-BOLK*TEMP(ID)*PSI0(NHE)/DM(ID)-
|
|
* X1*(GRD+FPRD(ID))-VT0/WMM(ID)*DENS(ID)
|
|
C
|
|
RETURN
|
|
ELSE IF(IBCHE.EQ.1) THEN
|
|
C
|
|
C 2. possibility - specifically disk - Hubeny (1990), Eq. (4.19)
|
|
C newer variant
|
|
C
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=IJ1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
IF(.NOT.LSKIP(ID,IJT)) THEN
|
|
FLUXW=W(IJT)*(FH(IJT)*RAD0(IJ)-HEXTRD(IJT))
|
|
GRD=GRD+FLUXW*ABSO0(IJ)
|
|
HEXN=HEXN+FLUXW*DABN0(IJ)
|
|
HEXT=HEXT+FLUXW*DABT0(IJ)
|
|
DO I=1,NLVEXP
|
|
HEX(I)=HEX(I)+FLUXW*DRCH0(I,IJ)
|
|
END DO
|
|
END IF
|
|
END DO
|
|
END IF
|
|
c
|
|
CCC=PCK/QGRAV
|
|
HR1=CCC*(GRD+FPRD(1))/DENS(1)
|
|
PG1=BOLK*PSI0(NHE)*TEMP(1)
|
|
HG1=SQRT(TWO*PG1/DENS(1)/QGRAV)
|
|
X=(ZD(1)-HR1)/HG1
|
|
IF(X.LT.3.) THEN
|
|
IF(X.LT.0.) X=0.
|
|
F1=8.86226925D-1*EXP(X*X)*ERFCX(X)
|
|
ELSE
|
|
F1=HALF*(UN-HALF/X/X)/X
|
|
END IF
|
|
X1=X*1.01
|
|
F1D=0.
|
|
IF(X1.LT.3.) THEN
|
|
F1D=8.86226925D-1*EXP(X1*X1)*ERFCX(X1)
|
|
ELSE
|
|
F1D=HALF*(UN-HALF/X1/X1)/X1
|
|
END IF
|
|
IF(X.GT.0.) F1D=(F1D-F1)*100./X
|
|
GGG=DENS(1)*HG1*F1
|
|
RF1=DENS(1)*F1D
|
|
CCD=CCC*F1D
|
|
C
|
|
DO IJ=1,NFREQE
|
|
B(NHE,IJ)=-CCD*WDEP0(IJ)*FH(IJFR(IJ))*ABSO0(IJ)
|
|
END DO
|
|
C
|
|
C columns corresponding to total particle density and temperature
|
|
C
|
|
B(NHE,NHE)=B(NHE,NHE)+(GGG+HR1*RF1)/PSI0(NHE)
|
|
IF(INRE.GT.0) B(NHE,NRE)=
|
|
* (GGG-RF1*ZD(1)+RF1*HR1)*HALF/TEMP(1)-CCD*(HEXT+HEIT(ID))
|
|
IF(INZD.GT.0) B(NHE,NZD)=RF1
|
|
IF(INPC.GT.0) B(NHE,NPC)=-CCD*(HEXN+HEIN(ID))
|
|
DO II=1,NLVEXP
|
|
B(NHE,NSE+II)=-CCD*(HEX(II)+HEIP(II,ID))
|
|
END DO
|
|
C
|
|
C The rhs vector
|
|
C
|
|
VECL(NHE)=DM(1)-GGG
|
|
RETURN
|
|
ELSE IF(IBCHE.EQ.2) THEN
|
|
C
|
|
C 3. possibility - specifically disk - Hubeny (1990), Eq. (4.19)
|
|
C older variant
|
|
C
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=IJ1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
IF(.NOT.LSKIP(ID,IJT)) THEN
|
|
FLUXW=W(IJT)*(FH(IJT)*RAD0(IJ)-HEXTRD(IJT))
|
|
GRD=GRD+FLUXW*ABSO0(IJ)
|
|
END IF
|
|
END DO
|
|
END IF
|
|
CCC=PCK/QGRAV
|
|
PR1=CCC*(GRD+FPRD(1))/DENS(1)
|
|
PG1=BOLK*PSI0(NHE)*TEMP(1)
|
|
HG1=SQRT(TWO*PG1/DENS(1)/QGRAV)
|
|
X=(ZD(1)-PR1)/HG1
|
|
IF(X.LT.3.) THEN
|
|
IF(X.LT.0.) X=0.
|
|
F1=8.86226925D-1*EXP(X*X)*ERFCX(X)
|
|
ELSE
|
|
F1=HALF*(UN-HALF/X/X)/X
|
|
END IF
|
|
GGG=HG1*QGRAV*HALF/F1
|
|
C
|
|
C columns corresponding to total particle density and temperature
|
|
C
|
|
B(NHE,NHE)=BOLK*TEMP(1)
|
|
IF(INRE.GT.0) B(NHE,NFREQE+INRE)=PG1/TEMP(1)
|
|
C
|
|
C The rhs vector
|
|
C
|
|
VECL(NHE)=DM(1)*GGG-PG1
|
|
RETURN
|
|
ELSE
|
|
C
|
|
C 4. a simple form of the bounary condition P_gas(ID=1)=PGAS0,
|
|
C where PGAS0 is an input parameter
|
|
C
|
|
B(NHE,NHE)=BOLK*TEMP(1)
|
|
IF(INRE.GT.0) B(NHE,NRE)=BOLK*PSI0(NHE)
|
|
VECL(NHE)=PGAS0-BOLK*TEMP(1)*PSI0(NHE)
|
|
RETURN
|
|
END IF
|
|
C
|
|
C *** Normal depth point (ID > 1)
|
|
C
|
|
C Columns (for matrices A and B) corresponding to mean intensities
|
|
C
|
|
50 CONTINUE
|
|
GRAV=QGRAV*(ZD(ID)+ZD(ID-1))*HALF
|
|
GRAVZ=GRAV*(ZD(ID)-ZD(ID-1))
|
|
DGRV=GRAVZ*HALF*WMM(ID)
|
|
GRD=0.
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=IJ1,NFREQE
|
|
IF(.NOT.LSKIP(ID,IJFR(IJ))) THEN
|
|
GRD=GRD+(FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ))*W(IJFR(IJ))
|
|
A(NHE,IJ)=-PCK*W(IJFR(IJ))*FKM(IJ)
|
|
B(NHE,IJ)=PCK*W(IJFR(IJ))*FK0(IJ)
|
|
END IF
|
|
END DO
|
|
END IF
|
|
C
|
|
VT0=HALF*VTURB(ID)*VTURB(ID)*WMM(ID)
|
|
VTM=HALF*VTURB(ID-1)*VTURB(ID-1)*WMM(ID)
|
|
C
|
|
C columns corresponding to total particle density
|
|
C
|
|
A(NHE,NHE)=-BOLK*TEMP(ID-1)-GN*(VTM+DGRV)
|
|
B(NHE,NHE)=BOLK*TEMP(ID)+GN*(VT0+DGRV)
|
|
C
|
|
C columns corresponding to temperature
|
|
C
|
|
IF(INRE.GT.0) THEN
|
|
A(NHE,NRE)=-BOLK*PSIM(NHE)+PCK*HEITM(ID)
|
|
B(NHE,NRE)=BOLK*PSI0(NHE)+PCK*HEIT(ID)
|
|
C(NHE,NRE)=PCK*HEITP(ID)
|
|
END IF
|
|
C
|
|
C columns corresponding to electron density
|
|
C
|
|
IF(INPC.GT.0) THEN
|
|
A(NHE,NPC)=GN*(VTM+DGRV)+PCK*HEINM(ID)
|
|
B(NHE,NPC)=-GN*(VT0+DGRV)+PCK*HEIN(ID)
|
|
C(NHE,NPC)=PCK*HEINP(ID)
|
|
END IF
|
|
C
|
|
C columns corresponding to NMP
|
|
C
|
|
IF(INMP.GT.0) THEN
|
|
A(NHE,NFREQE+INMP)=-GP*(VTM+DGRV)
|
|
B(NHE,NFREQE+INMP)=GP*(VT0+DGRV)
|
|
END IF
|
|
C
|
|
C columns corresponding to populations
|
|
C
|
|
DO II=1,NLVEXP
|
|
A(NHE,NSE+II)=A(NHE,NSE+II)+PCK*HEIPM(II,ID)
|
|
B(NHE,NSE+II)=B(NHE,NSE+II)+PCK*HEIP(II,ID)
|
|
C(NHE,NSE+II)=C(NHE,NSE+II)+PCK*HEIPP(II,ID)
|
|
END DO
|
|
C
|
|
C the rhs vector
|
|
C
|
|
VECL(NHE)=-GRAVZ*(DENS(ID)+DENS(ID-1))*HALF-
|
|
* BOLK*(TEMP(ID)*PSI0(NHE)-TEMP(ID-1)*PSIM(NHE))-
|
|
* PCK*(GRD+FPRD(ID))-
|
|
* VT0/WMM(ID)*DENS(ID)+VTM/WMM(ID)*DENS(ID-1)
|
|
C
|
|
RETURN
|
|
END
|