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 

************************************************************************