c
      subroutine atmint(fl,Teff,alfa,x,y,z,c12,o16,fmtot,fmenv,Ro,peff)
      implicit double precision (a-h,o-z)
c
c  Integrate an atmosphere with given L/Lo (fl) and Teff.
c  from tau=0 to tau=2/3 to obtain pressure at photosphere
c  and mass (in solar masses) above the photosphere, and the
c  radius (in solar radii) at the top of the atmosphere
c
c    INPUT:     fl = L/Lsun at surface
c             Teff = Effective temperature
c             alfa = mixing length to pressure scale height ratio
c                x = hydrogen mass fraction
c                y = helium4 mass fraction
c              c12 = C12 mass fraction
c              o16 = O16 mass fraction
c            fmtot = total mass of star, solar units
c     
c    OUTPUT:
c            fmenv = mass between tau=0 and photosphere, solar units
c               Ro = estimate of radius at tau=0, solar units
c             Peff = pressure at tau=2/3
c
c Uses Numerical Recipes integrator, and a subroutine (DIFATM)
c to provide r.h.s.'s to continuity, hse, radiative diffusion,
c and optical depth equation.  Allows for convective transport, too.
c
c Assumes that the atmosphere has a constant gravity with
c R=Reff.  The T(tau) relation is set by the equation of radiative
c diffusion or convection, and turns out to be equivalent to the Eddington
c approximation when radiative.
c
c Common block CONST contains fundamental physical constants
      COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR
c Common block ATMIN contains things needed to evaluate atmosphere
c structure equations in DIFATM, but not in calling arguments.  
c Implemented to preserve structure of Numerical Recipes routines.
      common/atmin/xc,yc,zc,cc,oc,flc,fmc,reff,grav,alfac
c Array struc(4) contains r/Rsun, lnP, lnT, and S (=1-Mr/Msun)
      dimension struc(4),dydx(4)
      external deratm
c
c load common block atmin for transferring variable to r.h.s. subroutine
      xc=x
      yc=y
      zc=z
      cc=c12
      oc=o16
      flc=fl
      fmc=fmtot
      alfac=alfa
c
c  Conditions at tau=0 (Eddington approximation)
c
c  To=Te/2^1/4
c  Po= radiation pressure at that T, times 1.1
c  So=0.d0 by definition
c  Ro is about equal to Reff; it doesn't matter at this level of approximation
c    we keep track of it though to estimate the `true' Ro.
c
      struc(4)=0.d0
      To=Teff/(2.d0**0.25d0)
      struc(3)=dlog(To)
      Po=1.1d0*4.d0*csb/3.d0/ccl*To**4
      struc(2)=dlog(po)
      reff=dsqrt(cls*fl/cp4/csb/Teff**4)/crs
      struc(1)=reff
      tau=1.d-10
      dtau0=1.d-10
      taup=2.d0/3.d0
      grav=cg*fmtot*cms/reff**2/crs**2
c
c Okay, now integrate downwards to the photosphere
c
      iastep=0
      do while (tau.lt.taup)
        call DERATM(tau,struc,dydx)
c Integrate down a step, keeping all changes to less than 1.d-3
        eps=1.d-4
        call BSSTEP(struc,dydx,4,tau,dtau0,eps,struc,dtau,
     1              dtaunxt,deratm)
        dtau0=dtaunxt
c Check to see if we'll pass the photosphere at the next step. If so
c shorten up on the step to ensure that we don't overshoot.
        if (tau+dtau0.gt.taup) dtau0=taup-tau
        istep=istep+1
      enddo
c
c Evaluate the radius at the top of the atmosphere
c
      Ro=Reff/struc(1)*Reff
      peff=dexp(struc(2))
      fmenv=struc(4)
      return
      end
c
      subroutine DERATM(tau,struc,dydx)
      implicit double precision (a-h,o-z)
      dimension struc(4),dydx(4)
c Common block CTRL contains most of the control parameters for ISUEVO
      common/CTRL/dllenv,dltenv,tcut,dtfac,sstep(7),tstep(7),
     1            nmod,nsav,nprint,nzp,itstep
c Common block ATMIN contains things needed to evaluate atmosphere
c structure equations but not in calling arguments.  Implemented to
c preserve structure of Numerical Recipes routines.
      common/ATMIN/x,y,z,c12,o16,fl,fmtot,reff,grav,alfa
c Common block CONST contains fundamental physical constants
      COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR
c Common block MLT contains switch for ML parameters
      COMMON/MLT/infmlt
c
c Differential equations used to integrate
c the atmosphere down from tau=0 to tau=2/3 to find the
c photospheric pressure.   For use in the Numerical Recipes routines
c such as RK4, RKQC, or BSSTEP
c
      r=struc(1)*crs
      p=dexp(struc(2))
      t=dexp(struc(3))
      tcut=5.5d0
      call EOS(p,t,x,y,c12,tcut,q,cp,ga,d,xd,xt,u,cv,eta,fmu)
      call OPACITY(dlog10(d),dlog10(t),x,y,Z,c12,o16,fk,0
     1          ,DLKLT,DLKLD)
c Definition of optical depth gives dr/dtau
      dydx(1)=-1.d0/fk/d/crs
c H.S.E. gives  d(lnP)/dtau
      dydx(2)=grav/fk/p
c Radiation diffusion or convection gives d(lnT)/dtau
      condl=3.d0*cls*fl/(16.d0*cp4*crs**2*csb)
c Radiative temperature gradient:
      gr=condl*fK*P/(grav*reff**2*T**4)
c check for convective stability
      IF(gr.LE.GA) THEN
c   stable against convection; set gradient to radiative gradient
        GRAD=GR
      ELSE
c   Convectively unstable; use mixing length theory to evaluate
c   the temperature gradient...
        call tgrad(fmtot,reff,t,p,d,alfa,infmlt,q,cp,fk,ga,gr,grad)
      endif
c
      dydx(3)=grav*grad/p/fk
      dydx(4)=cp4*Reff**2/fk*crs**2/cms
c      print *,t,p,dydx(1),dydx(2),dydx(3),dydx(4),fk
      return
      end