subroutine hedif
c     ================
c
c subroutine to calculate the depth dependent abundance profile for
c a layered H+He atmosphere.
c
      INCLUDE 'IMPLIC.FOR'                                                      
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ATOMIC.FOR'
      common/hediff/ hcmass,radstr
c     real depth(mdepth+1),qs(mdepth+1),
      dimension depth(mdepth+1),qs(mdepth+1),
     *     ps(mdepth+1),gams(mdepth+1),
     *     abunds(mdepth+1),hms(mdepth+1)
c
      data smas,srad /1.9891e33,6.9599e10/
      data z1,z2,a1,a2 /1.,2.,1.,4./
      data bigg,pi  / 6.6732e-8,3.141592654/
c
c    Set up starting values
c
      do id=1,nd
         depth(id+1)=dm(id)
      end do
      if(radstr.lt.1.e3) radius=radstr*srad
      if(hcmass.gt.1.e-10) hcmass=hcmass*1.e-13
c
      gam=1.e-30
      gams(1)=1.e-30
c
   10 continue
      depth(1)=1.e-10
      q1=(depth(1)*4*pi*radius**2/smas)    
      p1=(q1*grav**2/(4*pi*bigg))
      ps(1)=p1
      qs(1)=q1
      hms(1)=0
      abunds(1)=0.0
      dpsl=-6
      hm=0.0
      do i=2,nd+1,1
         q2=(depth(i)*4*pi*radius**2/smas)    
         p2=(q2*grav**2/(4*pi*bigg))
         dp=p2-p1
         dlp=log(p2)-log(p1)
         gam=gam+raph(gam,z1,z2,a1,a2)*dlp
         abun0=gam
         hm=hm+(q2-q1)/((1+gam*a2/a1))
         p1=p2
         ps(i)=p2
         qs(i)=q2
         gams(i)=gam
         abunds(i)=abun0
         hms(i)=hm
      end do
c
      dh1=(log10(hcmass)-log10(hms(nd+1)))
      dh=hcmass/hms(nd+1)
      if(dh.ge.0.99) go to 20
      gam=gams(1)*1.1
      gams(1)=gam
      hm=0.0
      go to 10
c
c     Now work backwards to get the full profiles
c
   20 continue
      q1=(depth(nd+1)*4*pi*radius**2/smas)    
      p1=(q1*grav**2/(4*pi*bigg))
c
c     store new helium abundance and corresponding new YTOT, MMY, WMM
c
      write(6,600)
      do id=1,nd
         aheold=0.
         ahenew=abunds(id+1)
         if(iathe.gt.0) then
            aheold=abund(iathe,id)
            abund(iathe,id)=ahenew
         end if
         ytot(id)=ytot(id)-aheold+ahenew
         wmy(id)=wmy(id)+(ahenew-aheold)*4.003
         wmm(id)=wmy(id)*hmass/ytot(id)
         write(6,601) id,aheold,ahenew,ytot(id),wmy(id),wmm(id)
      end do
  600 format(' stratified helium abundance'/
     *'  id     He(old)     He(new)    ytot     wmm     wmy'/)
  601 format(i4,1p5e11.3)
c
      return
      end