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