SUBROUTINE LTEGRD C ================= C C Driving procedure for computing the initial LTE-grey disk model C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'MODELQ.FOR' C PARAMETER (ERRT=1.D-3, THIRD=UN/3.D0, FOUR=4.D0) DIMENSION TAU0(MDEPTH),TEMP0(MDEPTH),ELEC0(MDEPTH), * DENS0(MDEPTH),ZD0(MDEPTH),DM0(MDEPTH) COMMON/PRSAUX/VSND2(MDEPTH),HG1,HR1,RR1 COMMON/FACTRS/GAMJ(MDEPTH),GAMH,FAK0 COMMON/TOTJHK/TOTJ(MDEPTH),TOTH(MDEPTH),TOTK(MDEPTH), * RDOPAC(MDEPTH),FLOPAC(MDEPTH) COMMON/FLXAUX/T4,PGAS,PRAD,PGM,PRADM,ITGMAX,ITGMX0 COMMON/CUBCON/A,B,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD C C ---------------- C Input parameters C ---------------- C C NDEPTH - number of depth point for evaluating LTE-grey model C < 0 - the program computes an isothermal structure; C the temperature is specified by an extra record C (see below) C DM1 - mass at the first depth point C ABROS0 - initial estimate of the Rosseland opacity (per gram) C at the first depth point C ABPLA0 - initial estimate of the Planck mean opacity (per C gram) at the first depth point C DION0 - initial estimate of the degree of ionization at the C first depth point (=1 for completely ionized; =1/2 for C completely neutral) C ITGMAX - maximum number of global iterations of the procedure C for determining the Eddington and opacity factors C - number of iterations in recalculating new depth scale C in order to get a better coverage of optical depths C > 0 - subroutine NEWDM (blind addition of points) C < 0 - subroutine NEWDMT - new depths determined as the C equal-segment endpoints of the curve y=tauros(m) C C -------------------------------------------------------------------- C C IDEPTH - mode of determining the mass-depth scale to be used C in linearization C = 0 - depth scale DM is set up by the program C = 1 - interpolation to the depth scale DM (in g*cm**-2), C which has been read in START C = 2 - DM is evaluated as mass corresponding to Rosseland C optical depths which are equidistantly spaced in C logarithms between the first point TAU1 and the C last-but-one point TAU2=0.99*TAUMAX; the last point C is TAUMAX, where TAUMAX is the Ross. optical depth C corresponding to the last depth point (midplane), C set up by the program C NCONIT - number of internal iterations for calculating the C LTE-gray model with convection C = 0 - if HMIX0>0, program sets NCONIT=10. C IPRING - switch for determining an amount of output from C the calculation of LTE-gray model C = 0 - only final LTE-gray model is printed C = 1 - more tables are printed C = 2 - complete output 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 C -------------------------------------------------------------------- C IF(NDGREY.EQ.0) THEN NDEPTH=ND ELSE NDEPTH=NDGREY END IF IF(NDEPTH.GT.MDEPTH) call quit(' NDEPTH too large in LTEGR', * ndept,mdepth) IDEPTH=IDGREY ITGMAX=ITGMX0 IF(HMIX0.GT.0..AND.NCONIT.EQ.0) NCONIT=10 if(dion0.lt.0) then abpmin=-dion0 dion0=1. end if T4=TEFF**4 TOTF=SIG4P*T4 ABFL0=SIGE/WMM(1) if(idmfix.eq.1) then t0=teff DMTOT=TOTF/EDISC else call greyd t0=temp(1) EDISC=TOTF/DMTOT end if c ZND=0. VSND20=2.76D-16*T0/WMM(1)*DION0+VTB*VTB HSCALG=SQRT(TWO*VSND20/QGRAV) HSCALR=4.19168946D-10*TOTF*ABFL0/QGRAV R=HSCALR/HSCALG WRITE(6,615) HSCALG,HSCALR,R 615 FORMAT(/' GAS PRESSURE SCALE HEIGHT = ',1PD10.3/ * ' RAD.PRESSURE SCALE HEIGHT = ',1PD10.3/ * ' RATIO = ',1PD10.3/) GAMH=UN FAK0=THIRD ANEREL=(DION0-HALF)/DION0 IF(ANEREL.LT.ERRT) ANEREL=ERRT IF(NDEPTH.EQ.0) NDEPTH=ND LCHC0=LCHC LCHC=.TRUE. ND0=ND ND=NDEPTH DO ID=1,ND0 DM0(ID)=DM(ID) END DO C C mass-depth scale C Initial estimate of the density, geometrical distance z, and C pressure C CALL ZMRHO(R,HSCALG) C if(ipring.eq.2) then xdm=dm(1) DO ID=1,ND if(id.gt.1) xdm=xdm-half*(dens(id)+dens(id-1))* * (zd(id)-zd(id-1)) WRITE(6,602) ID,DM(ID),TAUROS(ID), * TEMP(ID),ELEC(ID),PTOTAL(ID), * ZD(ID),ABROSD(ID),ABPLAD(ID) * ,dens(id),xdm end do end if c ITGREY=-1 AMUV0=DMVISC**(ZETA0+UN) AMUV1=UN-AMUV0 DO ID=1,ND PGS(ID)=DENS(ID)*VSND20 IF(DM(ID).LE.DMVISC*DM(ND)) THEN VISCD(ID)=(UN-FRACTV)*(ZETA1+UN)/ * DMVISC**(ZETA1+UN)*(DM(ID)/DM(ND))**ZETA1 THETA(ID)=(UN-FRACTV)*(DM(ID)/DMVISC/DM(ND))**(ZETA1+UN) ELSE VISCD(ID)=FRACTV*(ZETA0+UN)/AMUV1* * (DM(ID)/DM(ND))**ZETA0 THETA(ID)=(UN-FRACTV)+FRACTV*((DM(ID)/DM(ND))**(ZETA0+UN)- * AMUV0)/AMUV1 END IF GAMJ(ID)=UN C C First estimates of the values of the Rosseland opacity C and function tauthe C IF(ID.EQ.1) THEN TAUR=DM(ID)*ABROS0 TAUTHE(ID)=TAUR*THETA(ID)/(ZETA1+TWO) ABROSD(ID)=ABROS0 ABPLAD(ID)=ABPLA0 ELSE DDM=DM(ID)-DM(ID-1) TAUR=TAUROS(ID-1)+DDM*ABROSD(ID-1) TAUTHE(ID)=TAUTHE(ID-1)+DDM*ABROSD(ID-1)*THETA(ID) ABROSD(ID)=ABROSD(ID-1) ABPLAD(ID)=ABPLAD(ID-1) END IF C do ii=1,nlevel wop(ii,id)=un end do C C Determination of temperature and mean opacities C CALL TEMPER(ID,TAUR,ITGREY) END DO C IF(IPRING.GE.2) THEN WRITE(6,601) xdm=dm(1) DO ID=1,ND if(id.gt.1) xdm=xdm-half*(dens(id)+dens(id-1))* * (zd(id)-zd(id-1)) WRITE(6,602) ID,DM(ID),TAUROS(ID), * TEMP(ID),ELEC(ID),PTOTAL(ID), * ZD(ID),ABROSD(ID),ABPLAD(ID) * ,dens(id),xdm END DO END IF C C C Simultaneous solution of the hydrostatic equilibrium and the C z-m relation, assuming sound speed fixed C if(nconit.ge.0) CALL HESOLV C if(nconit.lt.-2) then do id=2,nd dm(id)=dm(id-1)-half*(dens(id)+dens(id-1))* * (zd(id)-zd(id-1)) end do end if IF(IPRING.GE.2) THEN xdm=dm(1) WRITE(6,601) DO ID=1,ND if(id.gt.1) xdm=xdm-half*(dens(id)+dens(id-1))* * (zd(id)-zd(id-1)) WRITE(6,602) ID,DM(ID),TAUROS(ID), * TEMP(ID),ELEC(ID),PTOTAL(ID), * ZD(ID),ABROSD(ID),ABPLAD(ID) * ,dens(id),xdm END DO END IF C C ------------------------------------------------------------------- C C Outer iteration loop for the pseudo-grey model; C basically generalized Unsold-Lucy procedure C C 1.iteration - assumes that C Rosseland opacity = flux mean opacity C Planck mean opac = absorption mean opacity C Eddington factors = 1/3 and 1/sgrt(3) C C next iterations C improvement of mean opacities and Eddington factors; C corrections of the temperature (generalized Unsold-Lucy) C C C 1.part C 100 ITGREY=ITGREY+1 C ANEREL=ELEC(1)/(DENS(1)/WMM(1)+ELEC(1)) DO ID=1,ND TAUR=TAUROS(ID) IF(ITGREY.GT.1) TAUR=TAUFLX(ID) CALL TEMPER(ID,TAUR,ITGREY) END DO C C Again simultaneous solution of the hydrostatic equilibrium C and the z-m relation, assuming sound speed fixed C if(nconit.ge.0) CALL HESOLV C IF(IPRING.GE.1) THEN WRITE(6,601) xdm=dm(1) DO ID=1,ND WRITE(6,602) ID,DM(ID),TAUROS(ID), * TEMP(ID),ELEC(ID),PTOTAL(ID), * ZD(ID),ABROSD(ID),ABPLAD(ID) * ,dens(id),xdm END DO END IF C 601 FORMAT(1H1,' ID DM TAUROSS TEMP NE P', * 8X,'ZD ROSS.MEAN PLANCK',' dens '/) 602 FORMAT(1H ,I3,1P2D9.2,0PF11.0,1P3D9.2,2X,2D9.2,2x,2d11.4,d9.2) C C If required, modification of the depth scale (logarithmically C equidistant not in m, but in Tau(ross) C C ***************************************************** C IF(NNEWD.GT.0.AND.ITGREY.EQ.0.AND.TAUROS(ND).GT.10.) THEN DO III=1,NNEWD CALL NEWDM END DO END IF C C another modification of the depth scale C IF(NNEWD.LT.0) THEN DO III=1,-NNEWD CALL NEWDMT END DO END IF C IF(HMIX0.GT.0.) THEN CALL CONTMD GO TO 200 END IF C C ***************************************************** C C If ITGMAX = 0 - no iterative improvement of the pseudo-grey C model is required C IF(ITGMAX.EQ.0) GO TO 200 IF(ITGREY.EQ.0) ITGREY=1 C C Opacities and mean intensities in all frequency points ; C evaluation of appropriate integrals over frequency C CALL RADTOT C C Interpolation of TOTH and FLOPAC, which are determined by RTE C at the intermediate depth ponts DM(ID+1/2) to the grid DM C DO ID=2,ND-1 A1=DM(ID+1)-DM(ID-1) A0=(DM(ID)-DM(ID-1))/A1 A1=(DM(ID+1)-DM(ID))/A1 TOTH(ID)=A0*TOTH(ID+1)+A1*TOTH(ID) FLOPAC(ID)=A0*FLOPAC(ID+1)+A1*FLOPAC(ID) END DO TOTH(ND)=0. FLOPAC(ND)=FLOPAC(ND-1) C C Determination of new temperature C IF(IPRING.GE.1) WRITE(6,613) ITGREY DO ID=1,ND HMECH=TOTF*(UN-THETA(ID)) DFLUX=TOTH(ID)-HMECH FAKK=TOTK(ID)/TOTJ(ID) ABRAD=RDOPAC(ID)/DENS(ID)/TOTJ(ID) GAMJ(ID)=ABRAD/ABPLAD(ID)/FAKK*THIRD IF(ID.NE.ND) THEN ABFLX=FLOPAC(ID)/TOTH(ID) ELSE ABFLX=ABFLXM END IF abflx=abrosd(id) IF(ID.EQ.1) THEN FHH=TOTH(ID)/TOTJ(ID) GAMH=FAKK/FHH/5.7753D-1 TAUFLX(ID)=ABFLX*DM(ID) TAUTHE(ID)=TAUFLX(ID)*THETA(ID)/(ZETA1+TWO) DFINT=TAUFLX(ID)*DFLUX DB0=FAKK/FHH*DFLUX ELSE IF(DM(ID).LE.DMVISC*DM(ND)) THEN ZETAD=ZETA1 ELSE ZETAD=ZETA0 END IF DDM=DM(ID)-DM(ID-1) A0=(ABFLXM*DM(ID)-ABFLX*DM(ID-1))/DDM/(ZETAD+TWO) A1=(ABFLX-ABFLXM)/DDM/(ZETAD+3.D0) TAUFLX(ID)=TAUFLX(ID-1)+DDM*HALF*(ABFLXM+ABFLX) TAUTHE(ID)=TAUTHE(ID-1)+ * A0*(THETA(ID)*DM(ID)-THETA(ID-1)*DM(ID-1))+ * A1*(THETA(ID)*DM(ID)**2-THETA(ID-1)*DM(ID-1)**2) DFINT=DFINT+DDM*HALF*(ABFLXM*DFLUXM+ABFLX*DFLUX) END IF ABFLXM=ABFLX DFLUXM=DFLUX IF(ITGMAX.GE.0) THEN C C generalized Unsold-Lucy procedure C B0=FOUR*SIG4P*TEMP(ID)**4 DIS=TOTF*VISCD(ID)/ABPLAD(ID)/DM(ND) DB1=ABRAD/ABPLAD(ID)*TOTJ(ID)-B0+DIS DB=DB1-3.D0*GAMJ(ID)*(DB0+DFINT) BNEW=FOUR*SIG4P*TEMP(ID)**4+DB TEMP(ID)=SQRT(SQRT(BNEW/4.D0/SIG4P)) BREL=DB/B0 END IF C C diagnostic output of iterative improvement C R2=ABFLX/ABROSD(ID) R3=ABRAD/ABPLAD(ID) IF(IPRING.GE.1) WRITE(6,614) ID,FAKK,TAUROS(ID), * ABROSD(ID),ABFLX,R2, * ABRAD,ABPLAD(ID),R3,TOTH(ID),HMECH,BREL END DO C 613 FORMAT(1H1,'ITERATIVE IMPROVEMENT, ITGREY =',I2// * ' ID FK TAUROS ABROS',4X, * 'ABFLUX RATIO ABRAD ABPLA RATIO',10X,'FLUX', * 7X,'MECH',4X,'DELTA(B)/B'/) 614 FORMAT(1H ,I3,1P2D9.2,1X,3D9.2,1X,3D9.2,3X,2D13.5,3X,D10.2) C IF(ITGREY.LE.IABS(ITGMAX)) GO TO 100 C C End of iteration loop for the pseudo-grey model C ------------------------------------------------------------------- C C C 2. The final part C Interpolation of the computed model to the depth scale which is C going to be used in the subsequent - complete-linearization - C step of the model construction C C C First option - no interpolation C C Second option - interpolation to the prescribed mass scale DM C 200 CONTINUE IF(IDEPTH.EQ.1) THEN CALL INTERP(DM,TEMP,DM0,TEMP0,ND,ND0,2,1,0) CALL INTERP(DM,ELEC,DM0,ELEC0,ND,ND0,2,1,1) CALL INTERP(DM,DENS,DM0,DENS0,ND,ND0,2,1,1) CALL INTERP(DM,ZD,DM0,ZD0,ND,ND0,2,1,0) END IF C C Third option - logarithmically equidistant Rosseland opt.depths C IF(IDEPTH.EQ.2) THEN READ(IBUFF,*) TAU1 TAU0(ND0)=TAUROS(ND) TAU2=TAU0(ND0)*0.99 DML0=LOG(TAU1) DLGM=(LOG(TAU2)-DML0)/(ND0-2) DO I=1,ND0-1 TAU0(I)=EXP(DML0+(I-1)*DLGM) END DO CALL INTERP(TAUROS,DM,TAU0,DM0,ND,ND0,2,1,0) CALL INTERP(TAUROS,TEMP,TAU0,TEMP0,ND,ND0,2,1,0) CALL INTERP(TAUROS,ELEC,TAU0,ELEC0,ND,ND0,2,1,1) CALL INTERP(TAUROS,DENS,TAU0,DENS0,ND,ND0,2,1,1) CALL INTERP(TAUROS,ZD,TAU0,ZD0,ND,ND0,2,1,0) END IF C C Fourth option - truncation of the disk and computing only C a disk atmosphere C IF(IDEPTH.GE.3) THEN TAU1=TAUROS(1) IF(IDEPTH.EQ.3) THEN READ(IBUFF,*) TDIV ELSE IF(IDEPTH.EQ.4) THEN READ(IBUFF,*) TAU0(ND) ELSE IF(IDEPTH.EQ.5) THEN READ(IBUFF,*) TDIV ELSE IF(IDEPTH.EQ.6) THEN READ(IBUFF,*) TAU0(1),TAU0(ND) END IF IF(IDEPTH.EQ.3.OR.IDEPTH.EQ.5) THEN DO ID=1,ND IF(TAUROS(ID).LE.TDIV.AND.TAUROS(ID+1).GT.TDIV) * ID1=ID END DO END IF if(tauros(nd).le.tdiv) ID1=ND IF(IDEPTH.EQ.3) THEN ND0=ID1 DO ID=1,ND0 DM0(ID)=DM(ID) TEMP0(ID)=TEMP(ID) ELEC0(ID)=ELEC(ID) DENS0(ID)=DENS(ID) ZD0(ID)=ZD(ID) END DO nd=nd0 ELSE IF(IDEPTH.EQ.5) THEN TAU0(ND0)=TAUROS(ID1) END IF IF(IDEPTH.GE.4) THEN TAU2=TAU0(ND0)*0.99 DML0=LOG(TAU1) DLGM=(LOG(TAU2)-DML0)/(ND0-2) DO I=1,ND0-1 TAU0(I)=EXP(DML0+(I-1)*DLGM) END DO CALL INTERP(TAUROS,DM,TAU0,DM0,ND,ND0,2,1,0) CALL INTERP(TAUROS,TEMP,TAU0,TEMP0,ND,ND0,2,1,0) CALL INTERP(TAUROS,ELEC,TAU0,ELEC0,ND,ND0,2,1,1) CALL INTERP(TAUROS,DENS,TAU0,DENS0,ND,ND0,2,1,1) CALL INTERP(TAUROS,ZD,TAU0,ZD0,ND,ND0,2,1,0) END IF c IFZ0=-1 ZND=ZD0(ND0) IF(INZD.GT.0) THEN INZD=0 IF(INSE.GT.0) INSE=INSE-1 END IF END IF C C in the two last 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.GT.0) THEN ND=ND0 DO I=1,ND DM(I)=DM0(I) TEMP(I)=TEMP0(I) ELEC(I)=ELEC0(I) DENS(I)=DENS0(I) ZD(I)=ZD0(I) END DO END IF C C Recalculation of the populations C DO ID=1,ND AN=DENS(ID)/WMM(ID)+ELEC(ID) ANEREL=ELEC(ID)/AN CALL ELDENS(ID,TEMP(ID),AN,ANE,ENRG,ENTT,WM,1) ELEC(ID)=ANE DENS(ID)=WMM(ID)*(AN-ANE) PGS(ID)=AN*BOLK*TEMP(ID) PHMOL(ID)=AHMOL CALL WNSTOR(ID) CALL STEQEQ(ID,POP,1) END DO if(nconit.lt.0) CALL PSOLVE IF(HMIX0.GE.0.AND.IPRING.GT.0) CALL CONOUT(2,IPRING) LCHC=LCHC0 RETURN END