subroutine gray
*
*     gray atmosphere stratification,   the pressure(log)  is patm, the
*     temperature(log) is tatm and the optical depth(log) is tau
*
      implicit double precision(a-h,o-z)

      logical writeit

      common/vca/modnr
      common/dvca/sgtmp
      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/je/je
      common/comp/amhyhe,amheca,x(4)
      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/writeit/writeit
      common/oldval/bold,teold 
      common/thresh/thresh
*
*     interpolate linearly in tk-od table to find the pressure
*     corresponding to initial tau and tatm
*
      tau=-4.
      tatm=dlog10(te*(0.75*0.5773)**0.25)
      patm = dlog10( gs ) + tau + 1.
98    continue
      patm1 = patm
*
*  integrate equation of hydrostatic equilibrium using log pressure
*  as independant variable,  quit when convection sets in or when tau=1
*  if radiative equilibrium
*
*  integration loops back to statement 34 after incrementing pressure 
*
      xmass=sm
34    tl=tatm
      pl=patm
*
* trying out stop at tau = 0.1 instead of tau = 10.0
*
      if ( tau .ge. -1. ) then 
         tauz=10.**tau
         ps=patm
         ts=tatm
         return
      endif
*
*  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
*
*  ccomp called before gray, no need to call here 
*
      call estatem
      ol1x = opac3(rl1x,tl)
      rtgx = 7.732749e+9*10.**(ol1x+pl-4.*tl)*sl(j)/xmass
 
      diff=rtgx-atgx
      if ( atgx .le. 0.0 ) diff = rtgx - 0.40
*
* return if convective
*
      if(diff .gt. 0.) then
         if ( patm .eq. patm1 ) then
	    patm = patm - .3
            goto 98
         endif
         taux=10.**t1+diff1*(10.**tau-10.**t1)/(diff1-diff) 
         ps=patm
         ts=tatm
         return
      endif
*
* otherwise continue with the atm integration
*
      diff1=diff
      t1=tau
      dtaudp=(2.30258*10.**(pl+ol1x))/gs
      tau1=10.**tau+0.05*dtaudp
      tau=dlog10(tau1)
      qtau=0.7104-0.1331*dexp(-3.4488d0*tau1)
      tatm=dlog10(te*((tau1+qtau)*0.75)**0.25)
      patm=patm+0.05
      go to 34
*
*  goto 34 ==> loop back to top and continue integration
*
      end

************************************************************************
*
*     subroutine gshell is entry gshell in subroutine cshell
*
************************************************************************