        subroutine ccomp
*
*  this routine returns the h, he, and c mass fractions given as
*  input:  xmass (position in the envelope), amhyhe, amheca.
*
*  this is the first diffusion version.  the profiles returned are
*  derived under several assumptions (some questionable).  in
*  short, we can do better, but this is going to be it for now.
*  we're doing a pulsation study, so it's mostly important that
*  the profiles be smooth.  this version will not mix h into he
*  by convection
*
*  Modified by PAB to remove discontinuity in He abundance when
* helium layer is only 100x thicker than hydrogen layer. 09 05 91
* Here, I let both H and C be trace elements between zone(1,2)
* and zone(2,1). 09 09 91 (This one works in giving a smooth profile.)
*
        implicit double precision(a-h,o-z)

        common/je/je
        common/comp/amhyhe,amheca,x(4)
        common/c/az(3),z(3),taux,ts,ps,sm,rs,err,te,gs,rho, 
     1      ol1x,rl1x,el1x,atgx,rtgx,tl,pl,xmass,acu,am,bk,an0,pi
        common/alpha/alph(4),qmix1,qmix2,zone(2,2)
        common/idiffus/idiffus

      ame = 1.e0-(xmass/sm)
*
* determine surface composition
*
      if(amhyhe.lt.1.0e-15.and.amheca.lt.1.0e-15) je=3
      if(amhyhe.lt.1.0e-15.and.amheca.gt.1.0e-15) je=2
      if(amhyhe.gt.1.0e-15) je=1
*
*  if it's a diffusion equilibrium, then proceed, otherwise, jump
*  down to the discontinuity code
*
*  determine local composition and load comp vector x accordingly

*** hydrogen surface layer ***
      if ( idiffus .eq. 1 ) then
        if ( je .eq. 1 ) then 
          if(zone(1,2) .lt. zone(2,1))then
*
*  Normal case. We compute profiles without a care in the world.
*  pure h 
*
           if ( ame .le. zone(1,1) ) then
              x(1) = 1.0
              x(2) = 0.0
              x(3) = 0.0
              x(4) = 0.0
*
*  lower h layer (includes he tail)
*
           elseif ( ame .gt. zone(1,1) .and. ame .le. qmix1 ) then
              xx = ame/qmix1
              yy = func1(xx)
              x(1) = yy
              x(2) = 1.0  - yy
              x(3) = 0.0
              x(4) = 0.0
*
*  upper he layer
*
           elseif ( ame .gt. qmix1  .and. ame .le. zone(1,2) ) then
              xx = ame/qmix1
              yy = func2(xx)
              x(1) = yy
              x(2) = 1.0  - yy
              x(3) = 0.0
              x(4) = 0.0
*
*  pure he
*
           elseif ( ame .gt. zone(1,2) .and. ame .le. zone(2,1) ) then
              x(1) = 0.
              x(2) = 1.0
              x(3) = 0.0
              x(4) = 0.0
*
*  lower he layer
*
           elseif ( ame .gt. zone(2,1) .and. ame .le. qmix2) then
              xx = ame/qmix2
              yy = func3(xx)
              x(1) = 0.
              x(2) = yy
              x(3) = 1.0  - yy
              x(4) = 0.0
*
*  upper c layer
*
           elseif ( ame .gt. qmix2 .and. ame .le. zone(2,2) ) then
              xx = ame/qmix2
              yy = func4(xx)
              x(1) = 0.
              x(2) = yy
              x(3) = 1.0  - yy
              x(4) = 0.0
*
*  pure c 
*
           elseif ( ame .gt. zone(2,2) ) then
              x(1) = 0.0
              x(2) = 0.0
              x(3) = 1.0
              x(4) = 0.0
           endif
          else
*
* Special case where M_H is only 100x less than M_He. 
* We compute profiles very carefully to avoid discontinuous slopes.
*
* compute changeover point from H as trace to C as trace
*
           zond1 = dlog10(zone(1,2))
           zond2 = dlog10(zone(2,1))
           zond3 = (zond1 + zond2) / 2.
           change = 10.**(zond3)
*
*  pure h  (same as before)
*
           if ( ame .le. zone(1,1) ) then
              x(1) = 1.0
              x(2) = 0.0
              x(3) = 0.0
              x(4) = 0.0
*
*  lower h layer (includes he tail) (again, same as before)
*
           elseif ( ame .gt. zone(1,1) .and. ame .le. qmix1 ) then
              xx = ame/qmix1
              yy = func1(xx)
              x(1) = yy
              x(2) = 1.0  - yy
              x(3) = 0.0
              x(4) = 0.0
*
*  upper he layer (up to point where C should become a trace element)
*
           elseif ( ame .gt. qmix1  .and. ame .le. zone(2,1) ) then
              xx = ame/qmix1
              yy = func2(xx)
              x(1) = yy
              x(2) = 1.0  - yy
              x(3) = 0.0
              x(4) = 0.0
*
*  Middle of He layer (between zone(1,2) and Zone (2,1))
*
           elseif ( ame .gt. zone(2,1)  .and. ame .le. zone(1,2) ) then
              xx = ame/qmix1
              yy = func2(xx)
              x(1) = yy
              yhe = 1.0 - yy
              zz = ame/qmix2
              ww = func3(zz)
	      if (zz.le.1.0) then
                ww = func3(zz)
	      else
                ww = func4(zz)
	      endif
*
* changing following lines so that there is no negative helium abundance
*
              x(2) = ww - yy
              x(3) = 1.0 - ww
              x(4) = 0.0
*
*  lower He layer. (Below where H is a trace element)
*
           elseif ( ame .gt. zone(1,2) .and. ame .le. qmix2) then
              xx = ame/qmix2
              yy = func3(xx)
              x(1) = 0.
              x(2) = yy
              x(3) = 1.0  - yy
*
*  upper c layer (same as before)
*
           elseif ( ame .gt. qmix2 .and. ame .le. zone(2,2) ) then
              xx = ame/qmix2
              yy = func4(xx)
              x(1) = 0.
              x(2) = yy
              x(3) = 1.0  - yy
              x(4) = 0.0
*
*  pure c  (same as before):
*
           elseif ( ame .gt. zone(2,2) ) then
              x(1) = 0.0
              x(2) = 0.0
              x(3) = 1.0
              x(4) = 0.0
           endif
          endif

*** helium surface layer ***

        elseif ( je .eq. 2 ) then
*
*  pure he
*
           if ( ame .le. zone(2,1) ) then
              x(1) = 0.
              x(2) = 1.0
              x(3) = 0.0
              x(4) = 0.0
*
*  lower he layer
*
           elseif ( ame .gt. zone(2,1) .and. ame .le. qmix2) then
              xx = ame/qmix2
              yy = func3(xx)
              x(1) = 0.
              x(2) = yy
              x(3) = 1.0  - yy
              x(4) = 0.0
*
*  upper c layer
*
           elseif ( ame .gt. qmix2 .and. ame .le. zone(2,2) ) then
              xx = ame/qmix2
              yy = func4(xx)
              x(1) = 0.
              x(2) = yy
              x(3) = 1.0  - yy
              x(4) = 0.0
*
*  pure c 
*
           elseif ( ame .gt. zone(2,2) ) then
              x(1) = 0.0
              x(2) = 0.0
              x(3) = 1.0
              x(4) = 0.0
           endif

*** carbon surface layer ***

        elseif ( je .eq. 3 ) then
           x(1) = 0.0
           x(2) = 0.0
           x(3) = 1.0
           x(4) = 0.0
        endif
      else
*
*  come here for discontinuity profiles
*
* determine local composition
*
      if(je.ne.1) goto 100
      if(ame.gt.amhyhe.and.amheca.gt.1.0e-15) je=2
      if(ame.gt.amhyhe.and.amheca.lt.1.0e-15) je=3
100   continue
      if(je.ne.2) goto 101
      if(ame.gt.amheca) je=3
101   continue
*
* load comp vector X accordingly
*
      if ( je .eq. 1 ) then
        x(1) = 1.0
        x(2) = 0.0
        x(3) = 0.0
        x(4) = 0.0
      else if ( je .eq. 2 ) then
        x(1) = 0.0
        x(2) = 1.0
        x(3) = 0.0
        x(4) = 0.0
      else if ( je .eq. 3 ) then
        x(1) = 0.0
        x(2) = 0.0
        x(3) = 1.0
        x(4) = 0.0
      endif
      endif

        return
        end

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

