subroutine env * * calculates white dwarf envelopes * * the variables are r the radius, m the log of the mass within a * sphere of radius r, p the log of the pressure, t the log of the * temperature, el1 the degeneracy parameter, ol1 the log of the * total opacity, rl1 the log of the density, atg the adiabatic * gradient, rtg the radiative gradient, phs the mixing length, * drdtp the derivative of the density(log) with respect to the * temperature(log) at constant pressure, cp the specific heat per * gram at constant pressure, ttg the actual temperature gradient. * * the control parameters are: * nmax, maximum number of shells * nsuff, indicates if max number of shells is exceeded * nlim1, shell number at the bottom of the convection zone * kind, shell number at the top of the convection zone * nlim, total number of shells for a model * 1/nhp, fraction of appropriate scale height used as integration step * iter, total number of iterations * nlum, number of luminosities * lmax, number of isotherms in eos table * ml, number of density points per isotherm (eos) * lmam, number of isotherms for which optical depths are computed * * the units are cgs throughout * implicit double precision(a-h,o-z) double precision m logical writeit common/shells/r(800),m(800),p(800),t(800),el1(800),ol1(800), 1 rl1(800),atg(800),rtg(800),phs(800) common/tenv/rhok(3,60,35),eta(3,60,35),pp(3,60,35),datg(3,60,35), 1 od(3,18,35),tk(3,60),uu(3,60,35),xtt(3,60,35),xrt(3,60,35) common/thresh/thresh common/je/je common/comp/amhyhe,amheca,x(4) common/xcomp/xcomp(400,3) common/slr/sl(11) common/c/az(3),z(3),taux,ts,ps,sm,rs,err,te,gs,rho,ol1x,rl1x,el1x 1 ,atgx,rtgx,tl,pl,xmass,acu,am,bk,an0,pi common/d/ ml(3,60),nmax,kind,lmax,nsuff,nlim,nlim1,j,iter,lmam, 1 nlum,ko,nhp common/grad/ drdtp(400),cp(400),ttg(400),drdtpx,cpx,pecx,ttgx, 1 deptx,dpdr,dtdr,dmdr,dr,xr,xt,cv,drdptx common/surf/u(2,3),v(2,3),w(2,3),rm,bm,is,ks,ls,ms,kstart,npass common/contrl/ds,g,smass,wc,idummy(7) common/terp/stpms,rat1 common/wprop/ bn2(400),ac2(400) common/writeit/writeit common/oldval/bold,teold common/fin/ifin common/dfin/rfin,blfin common/wprep/iprep common/dprep/aa(20,500), ecv(400), ext(400), exr(400) common /mixl/ aml,mls data dlummax, dtemax /0.1d0, 0.5d0/ * * fundamental constants (Pi,Avagadro's number,Boltzmann const.) * bk=1.380622e-16 ano=6.022169e+23 pi=3.1415927e0 stop=stpms+smass-33.298635 sm=10.**(smass-33.298635) gslast=1e8 xmass = sm if ( ifin .eq. 0 ) then rs=10.**(w(2,is)-10.8425968) sl(1)=10.**w(1,is) writeit = .false. elseif ( ifin .eq. 1 ) then rs = 10.**(rfin-10.8425968) sl(1) = 10.**blfin endif * * compute the surface gravity gs and scale the optical depth tables * accordingly * gs=27398.*sm/(rs*rs) glog=dlog10(gslast/gs) do 360 jm=1,3 do 260 l=1,lmam kmax=ml(jm,l) do 160 k=1,kmax od(jm,l,k)=od(jm,l,k)+glog 160 continue 260 continue 360 continue * * set chem comp for gray atm integration * specify a luminosity and compute the effective temperature te * call ccomp j = 1 nsuff=1 nlim1=0 am=(x(1)*az(1) + x(2)*az(2) * x(3)*az(3)) * 1.660531e-24 te=5800.*(sl(j)/(rs*rs))**0.25 * * if stepped far enough in h-r diagram, then write out envelope * to tape6 * if ( ifin .eq. 1 ) then dlum = dabs( bold - bm ) dteff = dabs( teold - dlog10(te) ) if ( dlum .gt. dlummax .or. dteff .gt. dtemax ) then bold = bm teold = dlog10( te ) writeit = .true. endif endif 103 format(1x,' amhyhe=',1pe14.4,' amheca=',1pe14.4) * * integrate gray atmosphere * call gray it=2 nhp=8 iter=1 * * start the inward integration with values provided by the gray * subroutine, if there is convection take the geometrical depth * as zero * p(1)=ps t(1)=ts r(1)=rs*6.96e+10 m(1)=dlog10(sm) dept=1e-20 deptx=dept kind1=0 do 3 n=1,nmax ko=n tl=t(n) pl=p(n) xmass=10.**m(n) call ccomp am=(x(1)*az(1) + x(2)*az(2) * x(3)*az(3)) * 1.660531e-24 do 105 iii = 1,3 xcomp(n,iii) = x(iii) 105 continue * * check that our (p,t) point is on the eeos table * (inserted 6/13/88) * skip over check if teff > thresh (30,000) * if (te .gt. thresh) goto 987 itchk = (tl - tk(je,1)) / (tk(je,2) - tk(je,1)) + 1 if ( tl .lt. tk(je,1) .or. tl .gt. tk(je,lmax) ) then if ( tl .lt. tk(je,1) ) itchk = 1 if ( tl .gt. tk(je,lmax) ) itchk = lmax endif mltmp = ml(je,itchk) * * end eeos check * 987 continue * * given p and t at r call the subroutine state and store the values * of el1,rl1,ol1,atg,rtg,drdtp,cp * call estatem ol1x = opac3(rl1x,tl) rtgx = 7.732749e+9*10.**(ol1x+pl-4.*tl)*sl(j)/xmass ol1(n)=ol1x rtg(n)=rtgx el1(n)=el1x rl1(n)=rl1x atg(n)=atgx drdtp(n)=drdtpx cp(n)=cpx gam1=((xt**2)*10.e0**(p(n)-rl1x-t(n))/cv)+xr bn2(n)=-drdtpx*(rtgx-atgx)*((6.67e-8*1.989e+33* 1 10.e0**m(n)/r(n)**2)**2)*10.e0**(rl1x-p(n)) ac2(n)=gam1*10.e0**(p(n)-rl1x)/r(n)**2 dpdr=-27398.*(10.**m(n))*rho*(6.96e+10/r(n))**2 dmdr=4.*pi*rho*r(n)**2 if ( ifin .eq. 1 ) then ecv(n) = cv ext(n) = xt exr(n) = xr endif * * define the integration step as a fraction of the pressure scale * height * dr=10.**p(n)/nhp/dpdr * * Schwarzschild's test, if radiative equilibrium take rtg * as the actual gradient, if convective define the mixing length * as the smaller of either the local pressure scale height or the * geometrical depth from the top of the convection zone and call * gradt to evaluate the true gradient * if ( atg(n) .le. rtg(n)) then nlim1=n kind=n-kind1 kind1=kind1+1 phs(n)=dabs(dr*nhp) dept=dept+dabs(dr) deptx=dept pecx=9./8./(1.601743e+07*cpx*10.**(ol1x+2.5*rl1x+m(n) 1 -3.*t(n)-0.5*p(n))*(dabs(drdtpx))**0.5 2 *(phs(n)**2)*(6.96e+10/r(n))**2) pecx = pecx/aml/aml if(mls.eq.0)then * * call Bohm-Vitense version * call gradt else pecx = pecx * (8./9.)/dsqrt(2.d0) call gradt2 endif ttg(n)=ttgx else ttg(n)=rtg(n) kind1=0 dept=0. endif dtdr=ttg(n)*dpdr*10.**(t(n)-p(n)) * * evaluate r,p,t,m at the next shell and quit if the envelope * mass fraction exceeds -stpms or the temperature exceeds tk(lmax) * call intg if ( t(n+1) .gt. tk(je,lmax) ) then goto 13 elseif ( stop .gt. m(n) ) then goto 13 elseif ( n+2 .gt. nmax ) then goto 11 endif nlim=n+1 3 continue * * end of large loop (envelope integration) at statement 3 * if we go to 11, then we ran out of memory, and didn't integrate * all the way to the fitting point. * 11 nsuff=-1 13 plim2=10.**p(nlim) plim1=plim2 acu=0. n = n-1 call write4 if ( writeit ) writeit = .false. gslast=gs return end ************************************************************************