295 lines
9.8 KiB
Fortran
295 lines
9.8 KiB
Fortran
SUBROUTINE LTEGR
|
|
C ================
|
|
C
|
|
C Driving procedure for computing the initial LTE-grey model
|
|
C atmosphere
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
C
|
|
COMMON ESEMAT(MLEVEL,MLEVEL),BESE(MLEVEL),
|
|
* DEPTH(MDEPTH),DEPTH0(MDEPTH),TAU(MDEPTH),TAU0(MDEPTH),
|
|
* TEMP0(MDEPTH),ELEC0(MDEPTH),DENS0(MDEPTH),DM0(MDEPTH)
|
|
C
|
|
C ----------------
|
|
C Input parameters
|
|
C ----------------
|
|
C
|
|
C NDEPTH - number of depth points for evaluating LTE-grey model
|
|
C = 0 - NDEPTH is taken to be ND-1
|
|
C TAUFIR - Rosseland optical depth in the first depth point
|
|
C TAULAS - Rosseland optical depth in the last depth point
|
|
C ABROS0 - an estimate of the Rosseland opacity (per gram) at the
|
|
C first depth point
|
|
C TSURF = 0 - surface temperature and the Hopf function are
|
|
C evaluated exactly
|
|
C > 0 - value of surface temperature is set to TSURF, and
|
|
C the Hopf function is assumed to be constant,
|
|
C corresponding to TSURF
|
|
C ALBAVE = 0 - wind blanketing is not considered
|
|
C > 0 - wind blanketing is considered; the averaged value
|
|
C of albedo [precisely the quantity (1+rho)/(1-rho)
|
|
C in the notation of Hummer (Ap.J. 257, 724, 1982)
|
|
C see his Eq. (3.1)] is ALBAVE
|
|
C DION0 - initial estimate of the degree of ionization at
|
|
C the first depth point (=1 for completely ionized;
|
|
C =1/2 for completely neutral)
|
|
C
|
|
C --------------------------------------------------------------------
|
|
C IDEPTH - mode of determining the mass-depth scale to be used
|
|
C in linearization
|
|
C = 0 - depth scale DM (in g*cm**-2) is evaluated as mass
|
|
C corresponding to Rosseland optical depths which
|
|
C are equidistantly spaced in logarithms between
|
|
C the first point TAUFIR and the last point TAULAS
|
|
C the last-but-one point is, however, set to
|
|
C TAULAS-1.
|
|
C = 2 - depth scale DM is evaluated as that corresponding
|
|
C to input values of Rosseland optical depth -
|
|
C array TAU0(ID), ID=1,ND
|
|
C = 1 - similar, but now DM is evaluated as mass
|
|
C corresponding to Rosseland optical depths which
|
|
C are equidistantly spaced in logarithms between
|
|
C the first point TAU1 and the last-but-one point
|
|
C TAU2; the last point is TAUL
|
|
C (i.e. similar to option 0, now TAU1 and TAUL may
|
|
C be different from TAUFIR nad TAULAS)
|
|
C = 3 - depth scale DM has already been read in START
|
|
C NCONIT - number of internal iterations for calculating the
|
|
C gray model with convection
|
|
C = 0 - and HMIX0>0, then NCONIT is set to 10
|
|
C IPRING - controls diagnostic output of the LTE-gray model
|
|
C calculations
|
|
C = 0 - no output
|
|
C = 1 - only final LTE-gray model
|
|
C = 2 - results of all internal iterations
|
|
C IHM > 0 - negative hydrogen ion considered in particle and
|
|
C charge conservation in ELDENS
|
|
C IH2 > 0 - hydrogen molecule considered in particle
|
|
C conservation in ELDENS
|
|
C IH2P > 0 - ionized hydrogen molecule considered in particle
|
|
C and charge conservation in ELDENS
|
|
C
|
|
IF(NDGREY.EQ.0) THEN
|
|
NDEPTH=ND
|
|
ELSE
|
|
NDEPTH=NDGREY
|
|
END IF
|
|
IF(NDEPTH.GT.MDEPTH)
|
|
* CALL QUIT('ndepth.gt.mdepth in LTEGR',ndepth,mdepth)
|
|
IDEPTH=IDGREY
|
|
IF(ALBAVE.GT.0.AND.TSURF.EQ.0.) TSURF=(0.433*ALBAVE)**0.25
|
|
HOPF0=0.
|
|
IF(TSURF.NE.0.) HOPF0=4.D0*TSURF**4/3.D0
|
|
T4=TEFF**4
|
|
ANEREL=(DION0-HALF )/DION0
|
|
IF(NCONIT.EQ.0.AND.HMIX0.GT.0.) NCONIT=10
|
|
IF(ANEREL.LT.1.D-3) ANEREL=1.D-3
|
|
LTE0=LTE
|
|
LTE=.TRUE.
|
|
IF(NDEPTH.EQ.0) NDEPTH=ND-1
|
|
ND0=ND
|
|
ND=NDEPTH
|
|
DO I=1,ND0
|
|
DM0(I)=DM(I)
|
|
END DO
|
|
C
|
|
C tau(ross) scale - logarithmically equidistant points between
|
|
C input TAUFIR and TAULAS
|
|
C
|
|
DML0=LOG(TAUFIR)
|
|
DLGM=(LOG(TAULAS)-DML0)/(NDEPTH-1)
|
|
DO I=1,NDEPTH
|
|
TAU0(I)=DML0+(I-1)*DLGM
|
|
TAU(I)=EXP(TAU0(I))
|
|
TAUROS(I)=TAU(I)
|
|
END DO
|
|
C
|
|
DPRAD=1.891204931D-15*T4
|
|
if(ifprad.eq.0) dprad=0.
|
|
C
|
|
PRAD0=DPRAD/1.732D0
|
|
ABROS=ABROS0
|
|
PLOG1=0.
|
|
PLOG2=0.
|
|
PLOG3=0.
|
|
DPLOG1=0.
|
|
DPLOG2=0.
|
|
IF(IPRING.GT.0) WRITE(6,601)
|
|
C
|
|
C -------------------------------------------------------------------
|
|
C
|
|
C 1.part
|
|
C Integration of the hydrostatic equilibrium equation on the
|
|
C tau(ross) scale;
|
|
C basically by a predictor-corrector method
|
|
C (similar to Kurucz's ATLAS code)
|
|
C
|
|
DO I=1,NDEPTH
|
|
J=0
|
|
TAUR=TAU(I)
|
|
C
|
|
C predictor step
|
|
C
|
|
IF(I.EQ.1) PLOG=LOG(GRAV/ABROS*TAUR+prad0)
|
|
IF(I.GT.1.AND.I.LE.4) PLOG=PLOG1+DPLOG1
|
|
IF(I.GT.4) PLOG=(3.*PLOG4+8.*DPLOG1-4.*DPLOG2+8.*DPLOG3)/3.
|
|
ERROR=1.
|
|
GO TO 40
|
|
C
|
|
C corrector step
|
|
C ------ iterate between hydrostatic equilibrium (which determines an
|
|
C increment of the total pressure) and state equations (which
|
|
C determine relevant number densities and then the Rosseland
|
|
C opacity)
|
|
C
|
|
30 IF(I.EQ.1) PNEW=LOG(GRAV/ABROS*TAUR+prad0)
|
|
IF(I.GT.1.AND.I.LE.4) PNEW=(PLOG+2.*PLOG1+DPLOG+DPLOG1)/3.
|
|
IF(I.GT.4) PNEW=(126.*PLOG1-14.*PLOG3+9.*PLOG4+42.*DPLOG+
|
|
* 108.*DPLOG1-54.*DPLOG2+24.*DPLOG3)/121.
|
|
ERROR=ABS(PNEW-PLOG)
|
|
PLOG=PNEW
|
|
40 PTOT=EXP(PLOG)
|
|
P=PTOT-TAUR*DPRAD-prad0
|
|
J=J+1
|
|
CALL ROSSOP(I,P,TAUR,HOPF0,T4,T,ANE,ABROS)
|
|
DPLOG=GRAV/ABROS*TAUR/PTOT*DLGM
|
|
IF(ERROR.GT.1.D-4.AND.J.LT.10) GO TO 30
|
|
C
|
|
C ------ end of the iteration loop;
|
|
C set up necessary quantities for the next depth step of the
|
|
C hydrostatic equilibrium
|
|
C
|
|
PLOG4=PLOG3
|
|
PLOG3=PLOG2
|
|
PLOG2=PLOG1
|
|
PLOG1=PLOG
|
|
DPLOG3=DPLOG2
|
|
DPLOG2=DPLOG1
|
|
DPLOG1=DPLOG
|
|
TEMP0(I)=T
|
|
ELEC0(I)=ANE
|
|
AN=P/T/BOLK
|
|
DEPTH(I)=(PTOT-prad0)/GRAV
|
|
DM(I)=DEPTH(I)
|
|
DENS0(I)=WMM(I)*(AN-ANE)
|
|
IF(IPRING.GT.0) WRITE(6,602) I,TAU(I),DEPTH(I),T,
|
|
* AN,ANE,P,ABROS
|
|
PTOTAL(I)=PTOT
|
|
PGS(I)=P
|
|
END DO
|
|
C
|
|
C -------------------------------------------------------------------
|
|
C
|
|
C 2. Second part - taking into account convection
|
|
C
|
|
IF(HMIX0.GT.0.) THEN
|
|
CALL CONTMP
|
|
GO TO 110
|
|
END IF
|
|
C
|
|
C -------------------------------------------------------------------
|
|
C
|
|
C 3. Third part
|
|
C Interpolation of the computed model to the depth scale which is
|
|
C going to be used in the subsequent - complete-linearization -
|
|
C part of the model atmosphere construction
|
|
C
|
|
ND=ND0
|
|
C
|
|
C First option - logarithmically equidistant Rosseland opt.depths
|
|
C the same first and last depth as in the first part
|
|
C
|
|
IF(IDEPTH.EQ.0) THEN
|
|
TAU1=TAUFIR
|
|
TAUL=TAULAS
|
|
TAU2=TAULAS-1.
|
|
C
|
|
C Second option - logarithmically equidistant Rosseland opt.depths
|
|
C the first, last-but-one, and last depths are read
|
|
C
|
|
ELSE IF(IDEPTH.EQ.1) THEN
|
|
READ(IBUFF,*) TAU1,TAU2,TAUL
|
|
END IF
|
|
C
|
|
IF(IDEPTH.LE.1) THEN
|
|
DML0=LOG(TAU1)
|
|
IF(TAUL.GT.0.) THEN
|
|
DLGM=(LOG(TAU2)-DML0)/(ND-2)
|
|
DO I=1,ND-1
|
|
TAU0(I)=DML0+(I-1)*DLGM
|
|
END DO
|
|
TAU0(ND)=LOG(TAUL)
|
|
ELSE
|
|
DLGM=(LOG(TAU2)-DML0)/(ND-1)
|
|
DO I=1,ND
|
|
TAU0(I)=DML0+(I-1)*DLGM
|
|
END DO
|
|
END IF
|
|
ELSE IF(IDEPTH.EQ.2) THEN
|
|
C
|
|
C Third option - prescribed set of Rosseland optical depths
|
|
C
|
|
READ(IBUFF,*) (TAU0(I),I=1,ND)
|
|
DO I=1,ND
|
|
TAU0(I)=LOG(TAU0(I))
|
|
END DO
|
|
ELSE if(idepth.eq.3) then
|
|
C
|
|
C Fourth option - interpolation to the prescribed mass scale DM
|
|
C
|
|
DO I=1,ND
|
|
DM(I)=DM0(I)
|
|
DM0(I)=LOG(DM(I))
|
|
END DO
|
|
CALL INTERP(DEPTH0,TEMP0,DM0,TEMP,NDEPTH,ND,2,0,0)
|
|
CALL INTERP(DEPTH0,ELEC0,DM0,ELEC,NDEPTH,ND,2,0,1)
|
|
CALL INTERP(DEPTH0,DENS0,DM0,DENS,NDEPTH,ND,2,0,1)
|
|
END IF
|
|
C
|
|
C in the first three options - interpolation from the previous
|
|
C Rosseland opacity scale to the new scale and from the previous
|
|
C mass depth scale to the new one
|
|
C
|
|
IF(IDEPTH.LE.2) THEN
|
|
DO I=1,NDEPTH
|
|
TEMP0(I)=TEMP(I)
|
|
ELEC0(I)=ELEC(I)
|
|
DENS0(I)=DENS(I)
|
|
TAU(I)=LOG(TAUROS(I))
|
|
DEPTH0(I)=LOG(DM(I))
|
|
END DO
|
|
CALL INTERP(TAU,DEPTH0,TAU0,DM0,NDEPTH,ND,3,0,0)
|
|
CALL INTERP(TAU,TEMP0,TAU0,TEMP,NDEPTH,ND,3,0,0)
|
|
CALL INTERP(TAU,ELEC0,TAU0,ELEC,NDEPTH,ND,3,0,1)
|
|
CALL INTERP(TAU,DENS0,TAU0,DENS,NDEPTH,ND,3,0,1)
|
|
DO ID=1,ND
|
|
DM(ID) = EXP(DM0(ID))
|
|
TOTN(ID) = DENS(ID)/WMM(ID)+ELEC(ID)
|
|
PTOTAL(ID) = DM(ID)*GRAV+PRAD0
|
|
PGS(ID) = TOTN(ID)*BOLK*TEMP(ID)
|
|
END DO
|
|
END IF
|
|
C
|
|
C Recalculation of the populations
|
|
C
|
|
DO ID=1,ND
|
|
CALL WNSTOR(ID)
|
|
CALL STEQEQ(ID,POP,1)
|
|
END DO
|
|
C
|
|
110 CONTINUE
|
|
IF(HMIX0.GE.0.) THEN
|
|
CALL CONOUT(2,IPRING)
|
|
END IF
|
|
LCHC=LCHC0
|
|
LTE=LTE0
|
|
IRSPLT=IRSPL0
|
|
601 FORMAT(1H1, 'COMPUTED LTE-GREY MODEL'//' ID TAU',7X,
|
|
* 'MASS',5X,'TEMP',7X,'N',10X,'NE',9X,'P',9X,'ROSS.OP'/)
|
|
602 FORMAT(1H ,I4,1P2D11.3,0PF8.0,1P4D11.3)
|
|
RETURN
|
|
END
|