subroutine ledouxc
*
* compute modified Ledoux B term to include composition gradient
* effects in the Brunt-Vaisala frequency
* This is the version that computes the B term in the core.
* I vary both the C and O abundance at the C/O interface.
* PAB 1 October 1991
*
      implicit double precision(a-h,o-z)
*
* has log T=t2, log Rho=u2, C abun.=xc2, O abun.=xo2, and log P=p2
*
      common/temp/s1,r1,s2,r2,b2,p2,t2,ea2,xc2,xo2,fca2,f2,q2,w2,c
      common/thermo/u2,up2,ut2,e2,ep2,et2,psi2,pg,o2,op2,ot2,fp,ft,
     1           en2(5),fn,fcc,ce,ci,cif,w,nu
*
* ledoux term value
*
      common/bldxc/bledc
      real*8 x(4)
*
* DEL is step size increment in composition
*
      bledc = 0.0
      del  = 0.01
      x(1) = 0.0
      x(2) = 0.0
      x(3) = xc2
      x(4) = xo2
*
* C/O interface
* I am computing partial ln P/partial ln O here.
* now to consider specific composition mixtures!
* CASE 3: x(4) or x(3) is less than del (use a 3 point formula)
* CASE 3.1: Nearly pure oxygen: Treat first point as pure O.
*           treat second point as C/O mixture 
*
         if(x(3).lt.del)then
           diff = 1.0-x(4)
           xone = x(4) - diff
           xc2 = 1.0 - xone
           xo2 = xone
           call istatco(p2,t2,1,.false.)
           dens1 = 10.**(u2)
           xtwo = 1.0
           xc2 = 0.0
           xo2 = xtwo
           call istatco(p2,t2,1,.false.)
           dens2 = 10.**(u2)
           xc2 = x(3)
           xo2 = x(4)
           bledc=(xo2/(2.*diff))*(dlog(dens2)-dlog(dens1))
           return
*
* CASE 3.2: Nearly pure Carbon: Treat first point as pure C.
*           treat second point as C/O mixture 
*
         else if(x(4).lt.del)then
           diff = x(4)
           xone = 0.0
           xc2 = 1.0 - xone
           xo2 = xone
           call istatco(p2,t2,1,.false.)
           dens1 = 10.**(u2)
           xtwo = x(4) + diff
           xc2 = 1.0 - xtwo
           xo2 = xtwo
           call istatco(p2,t2,1,.false.)
           dens2 = 10.**(u2)
           xc2 = x(3)
           xo2 = x(4)
           bledc=(xo2/(2.*diff))*(dlog(dens2)-dlog(dens1))
           return
         endif
*
* CASE 2: We are not near 0 or 1 in composition
* we can use 3-point differencing formula with impunity.
*
           xone = x(4) - del
           xo2 = xone
           xc2 = 1.0 - xone
           call istatco(p2,t2,1,.false.)
           dens1 = 10.**(u2)
           xtwo = x(4) + del
           xo2 = xtwo
           xc2 = 1.0 - xtwo
           call istatco(p2,t2,1,.false.)
           dens2 = 10.**(u2)
           xc2 = x(3)
           xo2 = x(4)
           bledc=(xo2/(2.*del))*(dlog(dens2)-dlog(dens1))
           return
      end

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