      subroutine ledoux
*
* compute modified Ledoux B term to include composition gradient
* effects in the Brunt-Vaisala frequency
* Pressure is the dependent variable throught and I sum partial pressures.
* In this version I change only the He abundance at the H/He interface.
* while I keep the sum of the abundances unity at the He/C interface.
* This duplicates the numerical differencing results, but the physical 
* justification eludes me at present.
* PAB Sept 10, 1991
*
      implicit double precision(a-h,o-z)
*
* EOS data
*
      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)
*
* composition info
*
      common/je/je
      common/comp/amhyhe,amheca,x(4)
*
* need for lmax
*
      common/d/ ml(3,60),nmax,kind,lmax,nsuff,nlim,nlim1,j,iter,lmam, 
     1 nlum,ko,nhp
*
* has log T=tl, log Rho=rho and log P=pl          
* also has stellar mass=sm and current mass=xmass 
*
      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
*
* diffusion profile info
*
      common/alpha/alph(4),qmix1,qmix2,zone(2,2)
*
* ledoux term value
*
      common/bldx/bled
      real*8 rhot1(35),rhot2(35),rhot3(35),rhot4(35)
      real*8 ptmp1(35),ptmp2(35),ptmp3(35),ptmp4(35)
      real*8 y21(35),y22(35),y23(35),y24(35)
      real*8 temp(60)
*
*  initialize summation variables
*
      bled    = 0.0
*
* step size increments in composition
*
      del     = 0.01
*
* This routine is complicated because we have to consider the
* H/He interface and He/C interfaces.
* We also have to make sure to avoid non-physical compositions
* such as X<0 or X>1.
*
* H/He interface
*
      if(x(1).gt.0.0 .and. x(2).gt.0.0)then
        if(x(3) .gt. x(1))then
*
* Very thin He buffer layer. Call H/He or He/C transition zone
* depending on case involved.
*
          go to 45
        endif
*
* first, locate T point on grid for Helium
*
        do 992 i=1,lmax
          temp(i) = tk(2,i)
992     continue
        call hunt(temp,lmax,tl,ktlo)
        kthi = ktlo + 1
*
* set size of temporary arrays
*
        kmax1 = ml(2,ktlo)
        kmax2 = ml(2,kthi)
*
* fill up temporary 1-D arrays for rho and T
* He arrays
*
        do 10 k=1,kmax1
          rhot1(k) = rhok(2,ktlo,k)
          ptmp1(k) =   pp(2,ktlo,k)
10      continue
        do 20 k=1,kmax2
          rhot2(k) = rhok(2,kthi,k)
          ptmp2(k) =   pp(2,kthi,k)
20      continue
*
* first, locate T point on grid for Hydrogen
*
        do 993 i=1,lmax
          temp(i) = tk(1,i)
993     continue
        call hunt(temp,lmax,tl,ktlo)
        kthi = ktlo + 1
*
* set size of temporary arrays
*
        kmax3 = ml(1,ktlo)
        kmax4 = ml(1,kthi)
*
* H arrays
*
        do 30 k=1,kmax3
          rhot3(k) = rhok(1,ktlo,k)
          ptmp3(k) =   pp(1,ktlo,k)
30      continue
        do 40 k=1,kmax4
          rhot4(k) = rhok(1,kthi,k)
          ptmp4(k) =   pp(1,kthi,k)
40      continue
*
* set BC's so second der is zero there
*
        yp1 = 1.0e+40
        ypn = 1.0e+40
*
* call for spline coefficients
*
        call spline(rhot1,ptmp1,kmax1,yp1,ypn,y21)
        call spline(rhot2,ptmp2,kmax2,yp1,ypn,y22)
        call spline(rhot3,ptmp3,kmax3,yp1,ypn,y23)
        call spline(rhot4,ptmp4,kmax4,yp1,ypn,y24)
*
* now compute P for each isotherm at a given rho,comp point
*
        call splint(rhot1,ptmp1,y21,kmax1,rho,p11)
        call splint(rhot2,ptmp2,y22,kmax2,rho,p12)
        call splint(rhot3,ptmp3,y23,kmax3,rho,p13)
        call splint(rhot4,ptmp4,y24,kmax4,rho,p14)
*
* linear interpolation for P at the actual temp.
*
        frac=(tl-tk(je,ktlo))/(tk(je,kthi)-tk(je,ktlo))
        pone=p11+frac*(p12-p11)
        ptwo=p13+frac*(p14-p13)
        pone=10.**(pone)
        ptwo=10.**(ptwo)
*
* now to consider specific composition mixtures!
* CASE 1: x(1) or x(2) is less than del (use a 3 point formula)
* CASE 1.1: Nearly pure helium: Treat first point as pure He.
*           treat second point as H / He mixture 
*
           if(x(1).lt.del)then
             diff = 1.0-x(2)
             xone = x(2) - diff
             xtwo = 1.0
*
* compute P for mixture by partial pressures.
* pmix1 is for x(2)-diff; pmix2 is for x(2)+diff
*
             pmix1 = xone*pone + (diff)*ptwo
             pmix1 = dlog(pmix1)
             pmix2 = pone + (diff)*ptwo
             pmix2 = dlog(pmix2)
*
* compute d ln P/d ln Y and return
*
             bled=(x(2)/(2.*diff))*(pmix2-pmix1)
             return
*
* CASE 1.2: Nearly pure hydrogen: Treat first point as pure H.
*           treat second point as H / He mixture 
*
           else if(x(2).lt.del)then
             diff=x(2)
             xone = 0.0
             xtwo = x(2) + diff
*
* compute P for mixture by summing partial pressures.
* pmix1 is for x(2)-diff; pmix2 is for x(2)+diff
*
             pmix1 = x(1)*ptwo
             pmix1 = dlog(pmix1)
             pmix2 = xtwo*pone + x(1)*ptwo
             pmix2 = dlog(pmix2)
             bled=(x(2)/(2.*diff))*(pmix2-pmix1)
             return
           endif
*
* CASE 2: We are not near 0 or 1 in composition
* we can use 3-point differencing formula with impunity.
*
           xone = x(2) - del
           xtwo = x(2) + del
           pmix1 = xone*pone + x(1)*ptwo
           pmix2 = xtwo*pone + x(1)*ptwo
           pmix1 = dlog(pmix1)
           pmix2 = dlog(pmix2)
*
* compute d ln P/d ln Y and return
*
           bled=(x(2)/(2.*del))*(pmix2-pmix1)
           return
         endif
*
* finished with H/He mixtures
*
* He/C interface
*
      if(x(2).gt.0.0 .and. x(3).gt.0.0)then
*
* first, locate T point on grid for Helium
*
45      do 994 i=1,lmax
          temp(i) = tk(2,i)
994     continue
        call hunt(temp,lmax,tl,ktlo)
        kthi = ktlo + 1
*
* set size of temporary arrays
*
        kmax1 = ml(2,ktlo)
        kmax2 = ml(2,kthi)
*
* fill up temporary 1-D arrays for rho and T
* He arrays
*
        do 50 k=1,kmax1
          rhot1(k) = rhok(2,ktlo,k)
          ptmp1(k) =   pp(2,ktlo,k)
50      continue
        do 60 k=1,kmax2
          rhot2(k) = rhok(2,kthi,k)
          ptmp2(k) =   pp(2,kthi,k)
60      continue
*
* first, locate T point on grid for Carbon
*
        do 995 i=1,lmax
          temp(i) = tk(3,i)
995     continue
        call hunt(temp,lmax,tl,ktlo)
        kthi = ktlo + 1
*
* set size of temporary arrays
*
        kmax3 = ml(3,ktlo)
        kmax4 = ml(3,kthi)
*
* C arrays
*
        do 70 k=1,kmax3
          rhot3(k) = rhok(3,ktlo,k)
          ptmp3(k) =   pp(3,ktlo,k)
70      continue
        do 80 k=1,kmax4
          rhot4(k) = rhok(3,kthi,k)
          ptmp4(k) =   pp(3,kthi,k)
80      continue
*
* set BC's so second der is zero there
*
        yp1 = 1.0e+40
        ypn = 1.0e+40
*
* call for spline coefficients
*
        call spline(rhot1,ptmp1,kmax1,yp1,ypn,y21)
        call spline(rhot2,ptmp2,kmax2,yp1,ypn,y22)
        call spline(rhot3,ptmp3,kmax3,yp1,ypn,y23)
        call spline(rhot4,ptmp4,kmax4,yp1,ypn,y24)
*
* now compute P for each isotherm at a given rho,comp point
*
        call splint(rhot1,ptmp1,y21,kmax1,rho,p11)
        call splint(rhot2,ptmp2,y22,kmax2,rho,p12)
        call splint(rhot3,ptmp3,y23,kmax3,rho,p13)
        call splint(rhot4,ptmp4,y24,kmax4,rho,p14)
*
* linear interpolation for P at the actual temp.
*
        frac=(tl-tk(je,ktlo))/(tk(je,kthi)-tk(je,ktlo))
        pone=p11+frac*(p12-p11)
        ptwo=p13+frac*(p14-p13)
        pone=10.**(pone)
        ptwo=10.**(ptwo)
*
* now to consider specific composition mixtures!
* CASE 3: x(3) or x(2) is less than del (use a 3 point formula)
* CASE 3.1: Nearly pure helium: Treat first point as pure He.
*           treat second point as He / C mixture 
*
         if(x(3).lt.del)then
           diff = 1.0-x(2)
           xone = x(2) - diff
           xtwo = 1.0
*
* compute P for mixture by summing partial pressures.
* pmix1 is for x(2)-diff; pmix2 is for x(2)+diff
*
           pmix1 = xone*pone + (1.0-xone)*ptwo
           pmix2 = xtwo*pone + (1.0-xtwo)*ptwo
           pmix1 = dlog(pmix1)
           pmix2 = dlog(pmix2)
           bled=(x(2)/(2.*diff))*(pmix1-pmix2)
           return
*
* CASE 3.2: Nearly pure Carbon: Treat first point as pure C.
*           treat second point as He / C mixture 
*
         else if(x(2).lt.del)then
           diff=x(2)
           xone = 0.0
           xtwo = x(2) + diff
*
* compute P for mixture by additive volume method.
* pmix1 is for x(2)-diff; pmix2 is for x(2)+diff
*
           pmix1 = xone*pone + (1.0-xone)*ptwo
           pmix2 = xtwo*pone + (1.0-xtwo)*ptwo
           pmix1 = dlog(pmix1)
           pmix2 = dlog(pmix2)
*
* compute d ln P/d ln Y and return.
*
           bled=(x(2)/(2.*diff))*(pmix1-pmix2)
           return
         endif
*
* CASE 2: We are not near 0 or 1 in composition
* we can use 3-point differencing formula with impunity.
*
        xone = x(2) - del
        xtwo = x(2) + del
        pmix1 = xone*pone + (1.0-xone)*ptwo
        pmix2 = xtwo*pone + (1.0-xtwo)*ptwo
        pmix1 = dlog(pmix1)
        pmix2 = dlog(pmix2)
*
* compute d ln P/d ln Y and return.
*
        bled=(x(2)/(2.*del))*(pmix1-pmix2)
        return
      endif
*
* finished with He/C mixtures.
*
* If no transition zone, then set bled to 0.0
*
      bled = 0.0

      return
      end

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

