      subroutine xbndry
*
* find the xtal boundary and print it's mass fraction
*
      implicit double precision(a-h,o-z)

      common/ac/ac(6)
      common/ao/ao(6)
      common/shells/ sa(400),ra(400),ba(400),pa(400),ta(400),
     1 ea(400),xca(400),fca(400),s(400),r(400),b(400),
     2 p(400),t(400),e(400),xc(400),sk(400),
     3 rk(400),bk(400),pk(400),tk(400)
      common/temp/s1,r1,s2,r2,b2,p2,t2,e2,xc2,xo2,fca2,f2,q2,w2,c
      common/contrl/ds,g,sm,wc,it,nite,ja,jb,j,k,l
      common/xbnd/amxc,amxo,gold
      common/ixswch/ixswch
      common/kill/ killit, nxmods, ikill

      logical killit
      logical oncec, onceo
      dimension y(400,7)
      equivalence (s(1),y(1,1))
*
*   locate xtal boundary between two shells   ....
*
      jx = jb+2
      jxc = jx
      jxo = jx
      oncec = .true.
      onceo = .true.
      do 2 i=1,jb
         temper  = y(i,5)
         xc2 = xc(i)
         xo2 = 1.0-xc2
         if ( ixswch .eq. 1 ) then
            pfreeze = phasec(temper)
         elseif ( ixswch .eq. 2 ) then
            pfreeze = phaseo(temper)
         elseif ( ixswch .eq. 3 ) then
            pfrzc = phasec(temper)
            pfrzo = phaseo(temper)
            pfreeze = xc2 * pfrzc + xo2 * pfrzo
         elseif ( ixswch .eq. 4 ) then
            pfrzc = phasec(temper)
            pfrzo = phaseo(temper)
            pfreeze = dmax1( pfrzc, pfrzo )
         else
            stop
         endif
         press = y(i,4)
 876     format('xbndry: i, xc, xo, t(i), p, px: ',i3,5f9.3)
         if ( ixswch .ne. 4 ) then
            if (press .le. pfreeze) then
               jx=i-1
               goto 3
            endif
         else
            if (press .le. pfrzc .and. oncec) then
               jxc = i-1
               oncec = .false.
            endif
            if (press .le. pfrzo  .and. onceo) then
               jxo = i-1
               onceo = .false.
            endif
            if (press .le. pfrzc .and. press .le. pfrzo ) goto 3
         endif
2     continue

3     continue
      if ( ixswch .ne. 4 ) then
         jxold=jx
         kx = jx + 1
         temp1 = y(jx,5)
         temp2 = y(kx,5)
         pres1 = y(jx,4)
         pres2 = y(kx,4)
         if ( ixswch .eq. 1 ) then
            px1 = phasec(temp1)
            px2 = phasec(temp2)
         elseif ( ixswch .eq. 2 ) then
            px1 = phaseo(temp1)
            px2 = phaseo(temp2)
         elseif ( ixswch .eq. 3 ) then
            pxc = phasec(temp1)
            pxo = phaseo(temp1)
            px1 = xc(jx) * pxc + (1.-xc(jx)) * pxo
            pxc = phasec(temp2)
            pxo = phaseo(temp2)
            px2 = xc(kx) * pxc + (1.-xc(kx)) * pxo
         endif
         sfrac = (px1 - pres1) / (pres2 - px2 + px1 - pres1 )
         am1 = 10.**s(jx)
         am2 = 10.**s(kx)
         amxc = dlog10( am1 + sfrac*(am2 - am1) ) 
      else
         if ( jxc .ne. jb + 2 .and. jxc .ne. 0 ) then
            jxcold = jxc
            pxc=phasec(y(jxc,5))
            kx=jxc+1
            dmdp=(y(kx,1)-y(jxc,1))/(y(kx,4)-y(jxc,4))
            amxc= y(jxc,1)+dmdp*(pxc-y(jxc,4))
         endif
         if ( jxo .ne. jb + 2 .and. jxo .ne. 0 ) then
            jxoold = jxo
            pxo=phaseo(y(jxo,5))
            kx=jxo+1
            dmdp=(y(kx,1)-y(jxo,1))/(y(kx,4)-y(jxo,4))
            amxo= y(jxo,1)+dmdp*(pxo-y(jxo,4))
         endif
      endif
*
*   find mx/mstar for this model              
*
91    format(1h ,'px,jx,g  ',e18.8,i5,e18.8)
      if ( jx .eq. jb-1 ) then
	 ikill = ikill+1
         if ( ikill .gt. nxmods ) killit = .true.
      endif

      return
      end 

************************************************************************
* FUNCTIONS
************************************************************************

      function cndlc(dens,t,gamma)
*
* this function returns the (carbon) thermal conductivity computed for a
* specific dens, temp, and gamma using the results if itoh et al (1983)
* for the liquid metal regime.
*
      implicit double precision (a-z)
      common /cdatl/ a(4), b(3), c(3), d(4), e(3), f(3)
 
      data at /12./
*
*  the following data from itoh et al. are for carbon (z=6) 
*
      t8 = 10.**(t -8)
      d6 = 10.**(dens - 6 )
      d6mu = d6 / 2.
      d6mu23 = d6mu**(2./3.)
 
      rs    = 1.388e-02 / d6mu**(1./3.) 
 
      x = 0.45641 * dlog( gamma ) - 1.31636
      x2 = x*x
      x3 = x2 * x
 
      asum = a(1) + a(2) * x + a(3) * x2 + a(4) * x3
      bsum = b(1) + b(2) * x + b(3) * x2
      csum = c(1) + c(2) * x + c(3) * x2
      dsum = d(1) + d(2) * x + d(3) * x2 + d(4) * x3
      esum = e(1) + e(2) * x + e(3) * x2
      fsum = f(1) + f(2) * x + f(3) * x2
 
      sminus1 = asum * ( 1. + bsum * rs + csum * rs*rs )
      splus1  = dsum * ( 1. + esum * rs + fsum * rs*rs )
 
      cstd6  = 1.018 * d6mu23 
      tmpsum = 1. + cstd6
      s = sminus1 - splus1 * cstd6 / tmpsum
 
      opacity = 2.363e+17 * d6 * t8 / ( at * tmpsum * s )
      cndlc  = dlog10( opacity )
 
      return
      end 

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

      function cndlhe(dens,t,gamma)
*
* this function returns the (helium) thermal conductivity computed for a
* specific dens, temp, and gamma using the results if itoh et al (1983)
* for the liquid metal regime.
*
      implicit double precision (a-z)
 
      common /hedatl/ a(4), b(3), c(3), d(4), e(3), f(3)
 
      data at /4./
*
*  the following data from itoh et al. are for helium (z=2) 
*
      t8 = 10.**(t - 8.)
      d6 = 10.**(dens - 6. )
      d6mu = d6 / 2.
      d6mu23 = d6mu**(2./3.)
 
      rs    = 1.388e-02 / d6mu**(1./3.) 
 
      x = 0.45641 * dlog( gamma ) - 1.31636
      x2 = x*x
      x3 = x2 * x
 
      asum = a(1) + a(2) * x + a(3) * x2 + a(4) * x3
      bsum = b(1) + b(2) * x + b(3) * x2
      csum = c(1) + c(2) * x + c(3) * x2
      dsum = d(1) + d(2) * x + d(3) * x2 + d(4) * x3
      esum = e(1) + e(2) * x + e(3) * x2
      fsum = f(1) + f(2) * x + f(3) * x2
 
      sminus1 = asum * ( 1. + bsum * rs + csum * rs*rs )
      splus1  = dsum * ( 1. + esum * rs + fsum * rs*rs )
 
      cstd6  = 1.018 * d6mu23 
      tmpsum = 1. + cstd6
      s = sminus1 - splus1 * cstd6 / tmpsum
 
      opacity = 2.363e+17 * d6 * t8 / ( at * tmpsum * s )
      cndlhe = dlog10( opacity )
 
      return
      end 

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

      function cndlhy(dens,t,gamma)
*
* this function returns the (hydrogen) thermal conductivity computed for a
* specific dens, temp, and gamma using the results if itoh et al (1983)
* for the liquid metal regime.
*
      implicit double precision (a-z)

      common /hydatl/ a(4), b(3), c(3), d(4), e(3), f(3)
      data at /1./
*
*  the following data from itoh et al. are for hydrogen (z=1) 
*
      t8 = 10.**(t -8.)
      d6 = 10.**(dens - 6. )
      d6mu = d6
      d6mu23 = d6mu**(2./3.)
 
      rs    = 1.388e-02 / d6mu**(1./3.) 
 
      x = 0.45641 * dlog( gamma ) - 1.31636
      x2 = x*x
      x3 = x2 * x
 
      asum = a(1) + a(2) * x + a(3) * x2 + a(4) * x3
      bsum = b(1) + b(2) * x + b(3) * x2
      csum = c(1) + c(2) * x + c(3) * x2
      dsum = d(1) + d(2) * x + d(3) * x2 + d(4) * x3
      esum = e(1) + e(2) * x + e(3) * x2
      fsum = f(1) + f(2) * x + f(3) * x2
 
      sminus1 = asum * ( 1. + bsum * rs + csum * rs*rs )
      splus1  = dsum * ( 1. + esum * rs + fsum * rs*rs )
 
      cstd6  = 1.018 * d6mu23 
      tmpsum = 1. + cstd6
      s = sminus1 - splus1 * cstd6 / tmpsum
 
      opacity = 2.363e+17 * d6 * t8 / ( at * tmpsum * s )
      cndlhy = dlog10( opacity )
 
      return
      end 

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

      function cndlo(dens,t,gamma)
*
* this function returns the (oxygen) thermal conductivity computed for a
* specific dens, temp, and gamma using the results if itoh et al (1983)
* for the liquid metal regime.
*
      implicit double precision (a-z)
 
      common /odatl/ a(4), b(3), c(3), d(4), e(3), f(3)
 
      data at /16./
*
*  the following data from itoh et al. are for oxygen (z=8) 
*
      t8 = 10.**(t -8)
      d6 = 10.**(dens - 6 )
      d6mu = d6 / 2.
      d6mu23 = d6mu**(2./3.)
 
      rs    = 1.388e-02 / d6mu**(1./3.) 
 
      x = 0.45641 * dlog( gamma ) - 1.31636
      x2 = x*x
      x3 = x2 * x
 
      asum = a(1) + a(2) * x + a(3) * x2 + a(4) * x3
      bsum = b(1) + b(2) * x + b(3) * x2
      csum = c(1) + c(2) * x + c(3) * x2
      dsum = d(1) + d(2) * x + d(3) * x2 + d(4) * x3
      esum = e(1) + e(2) * x + e(3) * x2
      fsum = f(1) + f(2) * x + f(3) * x2
 
      sminus1 = asum * ( 1. + bsum * rs + csum * rs*rs )
      splus1  = dsum * ( 1. + esum * rs + fsum * rs*rs )
 
      cstd6  = 1.018 * d6mu23 
      tmpsum = 1. + cstd6
      s = sminus1 - splus1 * cstd6 / tmpsum
 
      opacity = 2.363e+17 * d6 * t8 / ( at * tmpsum * s )
      cndlo  = dlog10( opacity )
 
      return
      end 

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

      function cndsc(dens,t,gamma)
*
* this function returns the thermal conductivity computed for a specific
* dens, temp, and gamma using the results if itoh et al (1984) for carbon
* in the xtal regime.
*
      implicit double precision (a-z)
*
*  declare variables (find the values in the block data
*  piece of code).
*
      common/cdats/ a(5), b(3), e(5), f(3), i(5), j(3), p(5), q(3),
     1          alpha(4), beta(4), c, d, g, h, k, l, r, s
      data at, z, twopi, pisqr /12.,6., 6.283185308, 9.869604404/
 
      x13 = -1./3.
      x23 = -2./3.
      t8 = 10.**t / 1.e+08
      d6 = 10.**dens / 1.e+06 
      d6mu = d6 * z / at
      d6mu23 = d6mu**(2./3.)
*
* compute the interpolation coefficients v and w
*
      u = twopi * ( dens - 3. ) / 9.
      v = alpha(1) + alpha(2) * gamma**x13 + alpha(3) * gamma**x23
     1    + alpha(4) / gamma
      w = beta(1) + beta(2) * gamma**x13 + beta(3) * gamma**x23
     1    + beta(4) / gamma
*
* set up the variables (sines and cosines) for the loops
*
      su  = dsin(u)
      su2 = su * su 
      cu2 = 1. - su2
      cu  = dsqrt( cu2 )
      s2u = 2. * su * cu
      c2u = cu2 - su2
      s3u = su * ( 3. - 4. * su2 )
      c3u = cu * ( 4. * cu2 - 3. )
      c4u = -8. * cu2 * su2 + 1.
*
*  here's the loops, unwrapped
*
      asum = a(1) + a(2) * cu    + a(3) * c2u
     1            + a(4) * c3u + a(5) * c4u
      bsum = b(1) * su + b(2) * s2u + b(3) * s3u
     1       + c*u + d
      esum = e(1) + e(2) * cu    + e(3) * c2u
     1            + e(4) * c3u + e(5) * c4u
      fsum = f(1) * su + f(2) * s2u + f(3) * s3u
     1       + g*u + h
      isum = i(1) + i(2) * cu    + i(3) * c2u
     1            + i(4) * c3u + i(5) * c4u
      jsum = j(1) * su + j(2) * s2u + j(3) * s3u
     1       + k*u + l
      psum = p(1) + p(2) * cu    + p(3) * c2u
     1            + p(4) * c3u + p(5) * c4u
      qsum = q(1) * su + q(2) * s2u + q(3) * s3u
     1       + r*u + s
 
      isig171 = asum + bsum
      isig5k  = esum + fsum
      ikap171 = isum + jsum
      ikap5k  = psum + qsum
 
      isig    = (1.-v) * isig171 + v * isig5k
      ikap    = (1.-w) * ikap171 + w * ikap5k
 
      smgamma = .07832 * (z/at) * dsqrt(d6) / t8
      sg2     = smgamma * smgamma
 
      g0      = 13.00 / dsqrt( 1. + 0.0174 * sg2 ) 
      g2   = (sg2 / pisqr) * (1. + 0.0118 * sg2)**(-3./2.)
 
      fkap    = isig * g0 + ikap * g2
      cstd6   = 1.018 * d6mu23
      nukap   = 9.554e+16 * t8 * dsqrt( 1. + 1./cstd6 ) * fkap
 
      opacity = 4.146d+15 * d6mu / dsqrt( 1. + cstd6 )
     1          * t8 * 1.e+18 / nukap
 
      cndsc  = dlog10( opacity )
 
      return
      end 

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

      function cndshe(dens,t,gamma)
*
* this function returns the thermal conductivity computed for a specific
* dens, temp, and gamma using the results if itoh et al (1984) for helium
* in the xtal regime.
*
      implicit double precision (a-z)
*
*  declare variables (find the values in the block data
*  piece of code).
*
      common/hedats/ a(5), b(3), e(5), f(3), i(5), j(3), p(5), q(3),
     1          alpha(4), beta(4), c, d, g, h, k, l, r, s
      data at, z, twopi, pisqr /4.,2., 6.283185308, 9.869604404/
 
      x13 = -1./3.
      x23 = -2./3.
      t8 = 10.**t / 1.e+08
      d6 = 10.**dens / 1.e+06 
      d6mu = d6 * z / at
      d6mu23 = d6mu**(2./3.)
*
* compute the interpolation coefficients v and w
*
      u = twopi * ( dens - 3. ) / 9.
      v = alpha(1) + alpha(2) * gamma**x13 + alpha(3) * gamma**x23
     1    + alpha(4) / gamma
      w = beta(1) + beta(2) * gamma**x13 + beta(3) * gamma**x23
     1    + beta(4) / gamma
*
* set up the variables (sines and cosines) for the loops
*
      su  = dsin(u)
      su2 = su * su 
      cu2 = 1. - su2
      cu  = dsqrt( cu2 )
      s2u = 2. * su * cu
      c2u = cu2 - su2
      s3u = su * ( 3. - 4. * su2 )
      c3u = cu * ( 4. * cu2 - 3. )
      c4u = -8. * cu2 * su2 + 1.
*
*  here's the loops, unwrapped
*
      asum = a(1) + a(2) * cu    + a(3) * c2u
     1            + a(4) * c3u + a(5) * c4u
      bsum = b(1) * su + b(2) * s2u + b(3) * s3u
     1       + c*u + d
      esum = e(1) + e(2) * cu    + e(3) * c2u
     1            + e(4) * c3u + e(5) * c4u
      fsum = f(1) * su + f(2) * s2u + f(3) * s3u
     1       + g*u + h
      isum = i(1) + i(2) * cu    + i(3) * c2u
     1            + i(4) * c3u + i(5) * c4u
      jsum = j(1) * su + j(2) * s2u + j(3) * s3u
     1       + k*u + l
      psum = p(1) + p(2) * cu    + p(3) * c2u
     1            + p(4) * c3u + p(5) * c4u
      qsum = q(1) * su + q(2) * s2u + q(3) * s3u
     1       + r*u + s
 
      isig171 = asum + bsum
      isig5k  = esum + fsum
      ikap171 = isum + jsum
      ikap5k  = psum + qsum
 
      isig    = (1.-v) * isig171 + v * isig5k
      ikap    = (1.-w) * ikap171 + w * ikap5k
 
      smgamma = .07832 * (z/at) * dsqrt(d6) / t8
      sg2     = smgamma * smgamma
 
      g0      = 13.00 / dsqrt( 1. + 0.0174 * sg2 ) 
      g2   = (sg2 / pisqr) * (1. + 0.0118 * sg2)**(-3./2.)
 
      fkap    = isig * g0 + ikap * g2
      cstd6   = 1.018 * d6mu23
      nukap   = 9.554e+16 * t8 * dsqrt( 1. + 1./cstd6 ) * fkap
 
      opacity = 4.146d+15 * d6mu / dsqrt( 1. + cstd6 )
     1          * t8 * 1.e+18 / nukap
 
      cndshe = dlog10( opacity )
 
      return
      end 

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

      function cndso(dens,t,gamma)
*
* this function returns the thermal conductivity computed for a specific
* dens, temp, and gamma using the results if itoh et al (1984) for oxygen
* in the xtal regime.
*
      implicit double precision (a-z)
*
*  declare variables (find the values in the block data
*  piece of code).
*
      common/odats/ a(5), b(3), e(5), f(3), i(5), j(3), p(5), q(3),
     1          alpha(4), beta(4), c, d, g, h, k, l, r, s
      data at, z, twopi, pisqr /16.,8., 6.283185308, 9.869604404/
 
      x13 = -1./3.
      x23 = -2./3.
      t8 = 10.**t / 1.e+08
      d6 = 10.**dens / 1.e+06 
      d6mu = d6 / 2.
      d6mu23 = d6mu**(2./3.)
*
* compute the interpolation coefficients v and w
*
      u = twopi * ( dens - 3. ) / 9.
      v = alpha(1) + alpha(2) * gamma**x13 + alpha(3) * gamma**x23
     1    + alpha(4) / gamma
      w = beta(1) + beta(2) * gamma**x13 + beta(3) * gamma**x23
     1    + beta(4) / gamma
*
* set up the variables (sines and cosines) for the loops
*
      su  = dsin(u)
      su2 = su * su 
      cu2 = 1. - su2
      cu  = dsqrt( cu2 )
      s2u = 2. * su * cu
      c2u = cu2 - su2
      s3u = su * ( 3. - 4. * su2 )
      c3u = cu * ( 4. * cu2 - 3. )
      c4u = -8. * cu2 * su2 + 1.
*
*  here's the loops, unwrapped
*
      asum = a(1) + a(2) * cu    + a(3) * c2u
     1            + a(4) * c3u + a(5) * c4u
      bsum = b(1) * su + b(2) * s2u + b(3) * s3u
     1       + c*u + d
      esum = e(1) + e(2) * cu    + e(3) * c2u
     1            + e(4) * c3u + e(5) * c4u
      fsum = f(1) * su + f(2) * s2u + f(3) * s3u
     1       + g*u + h
      isum = i(1) + i(2) * cu    + i(3) * c2u
     1            + i(4) * c3u + i(5) * c4u
      jsum = j(1) * su + j(2) * s2u + j(3) * s3u
     1       + k*u + l
      psum = p(1) + p(2) * cu    + p(3) * c2u
     1            + p(4) * c3u + p(5) * c4u
      qsum = q(1) * su + q(2) * s2u + q(3) * s3u
     1       + r*u + s
 
      isig171 = asum + bsum
      isig5k  = esum + fsum
      ikap171 = isum + jsum
      ikap5k  = psum + qsum
 
      isig    = (1.-v) * isig171 + v * isig5k
      ikap    = (1.-w) * ikap171 + w * ikap5k
 
      smgamma = .07832 * (z/at) * dsqrt(d6) / t8
      sg2     = smgamma * smgamma
 
      g0      = 13.00 / dsqrt( 1. + 0.0174 * sg2 ) 
      g2   = (sg2 / pisqr) * (1. + 0.0118 * sg2)**(-3./2.)
 
      fkap    = isig * g0 + ikap * g2
      cstd6   = 1.018 * d6mu23
      nukap   = 9.554e+16 * t8 * dsqrt( 1. + 1./cstd6 ) * fkap
 
      opacity = 4.146d+15 * d6mu / dsqrt( 1. + cstd6 )
     1          * t8 * 1.e+18 / nukap
 
      cndso  = dlog10( opacity )
 
      return
      end 

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

      function f12(x)       
*
* Compute Fermi-Dirac integral F1/2 to accuracy of 10^-4 using 
* rational function approximation of Anita (1993, ApJS, 84, 101).
*
      implicit double precision(a-h,o-z)
      dimension a1(8),b1(8),a2(11),b2(12)
*
* high precision formulae
*
      data an,m1,k1,m2,k2/0.5d0,7,7,10,11/
      data a1/5.75834152995465d+6,1.30964880355883d+7,
     1        1.07608632249013d+7,3.93536421893014d+6,
     2        6.42493233715640d+5,4.16031909245777d+4,
     3        7.77238678539648d+2,1.00000000000000d+0/
      data b1/6.49759261942269d+6,1.70750501625775d+7,
     1        1.69288134856160d+7,7.95192647756086d+6,
     2        1.83167424554505d+6,1.95155948326832d+5,
     3        8.17922106644547d+3,9.02129136642157d+1/
      data a2/4.85378381173415d-14,1.64429113030738d-11,
     1        3.76794942277806d-09,4.69233883900644d-07,
     2        3.40679845803144d-05,1.32212995937796d-03,
     3        2.60768398973913d-02,2.48653216266227d-01,
     4        1.08037861921488d+00,1.91247528779676d+00,
     5        1.00000000000000d+00/
      data b2/7.28067571760518d-14,2.45745452167585d-11,
     1        5.62152894375277d-09,6.96888634549649d-07,
     2        5.02360015186394d-05,1.92040136756592d-03,
     3        3.66887808002874d-02,3.24095226486468d-01,
     4        1.16434871200131d+00,1.34981244060549d+00,
     5        2.01311836975930d-01,-2.14562434782759d-02/
*
* small eta rational function       
*
      if(x .lt. 2.d0)then
        xx = dexp(x)
        rn = xx + a1(m1)
        do 100 i=m1-1,1,-1
          rn = rn*xx + a1(i)
100     continue
        den = b1(k1+1)
        do 200 i=k1,1,-1
          den = den*xx + b1(i)
200     continue
        f12 = xx * rn / den
       else
*
* large eta rational function       
*
        xx = 1.d0 / x**2
        rn = xx + a2(m2)
        do 300 i=m2-1,1,-1
          rn = rn*xx + a2(i)
300     continue
        den = b2(k2+1)
        do 400 i=k2,1,-1
          den = den*xx + b2(i)
400     continue
        f12 = ((x**(an+1.d0)) * rn) / den
      endif

      return      
      end 

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

      function finv(f)       
*
* Compute Fermi-Dirac integral F1/2 to accuracy of 10^-4 using 
* rational function approximation of Anita (1993, ApJS, 84, 101).
*
      implicit double precision(a-h,o-z)
      dimension a1(3),b1(3),a2(3),b2(3)
*
* low precision formulae
*
      data an,m1,k1,m2,k2/0.5d0,2,2,2,2/
      data a1/4.4593646d+01,+1.1288764d+01,1.0000000d+00/
      data b1/3.9519346d+01,-5.7517464d+00,2.6594291d-01/
      data a2/3.4873722d+01,-2.6922515d+01,1.0000000d+00/
      data b2/2.6612832d+01,-2.0452930d+01,1.1808945d+01/
*
* small eta rational function       
*
      if(f .lt. 4.d0)then
        rn = f + a1(m1)
        do 100 i=m1-1,1,-1
          rn = rn*f + a1(i)
100     continue
        den = b1(k1+1)
        do 200 i=k1,1,-1
          den = den*f + b1(i)
200     continue
        finv = dlog(f * rn / den)
       else
*
* large eta rational function       
*
        ff = 1.d0 / f**(1.d0/(1.d0+an))
        rn = ff + a2(m2)
        do 300 i=m2-1,1,-1
          rn = rn*ff + a2(i)
300     continue
        den = b2(k2+1)
        do 400 i=k2,1,-1
          den = den*ff + b2(i)
400     continue
        finv =  rn / (den * ff)
      endif

      return      
      end 

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

        function func1(x)
*
*  func1 (for diffusion exponents) Top of H/He transition zone.
*
        implicit double precision(a-h,o-z)
        common/alpha/alph(4),qmix1,qmix2,zone(2,2)
 
        x5 = x**alph(1)
        result = ( 2. - x5 ) / ( 2. + 3. * x5 )
        func1 = result
 
        return
        end

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

        function func2(x)
*
*  func2 (for diffusion exponents) Bottom of H/He transition zone.
*
        implicit double precision(a-h,o-z)
        common/alpha/alph(4),qmix1,qmix2,zone(2,2)
 
        xx = x**alph(2)
        result =  xx / ( 8. - 3. * xx ) 
        func2 = result
 
        return
        end

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

        function func3(x)
*
*  func3 (for diffusion exponents) Top of He/C transition zone.
*
        implicit double precision(a-h,o-z)
        common/alpha/alph(4),qmix1,qmix2,zone(2,2)
 
        x2 = x**alph(3)
        result = ( 2. - x2 ) / ( 2. + 2. * x2 )
        func3 = result
 
        return
        end

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

        function func4(x)
*
*  func4 (for diffusion exponents) Bottom of He/C transition zone.
*
        implicit double precision(a-h,o-z)
        common/alpha/alph(4),qmix1,qmix2,zone(2,2)
 
        xx = x**alph(4)
        result = xx / ( 6. - 2. * xx )
        func4 = result
 
        return
        end

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

      function hlcnc(d,t)
*
* Don Lamb's fit to the Hubbard-Lampe carbon conductive opacities
*
      implicit double precision(a-h,o-z)
      dimension a(25),w(25)

      data a/-5.18321e0,2.88323e-2,-1.68700e0,-3.80388e-1,-1.42808e0, 
     1       1.42606e0,-6.77026e-2,2.99633e-1,3.02198e-1,-1.94700e-1, 
     2       -2.72927e-3,4.24382e-2,-9.77295e-2,-1.50302e-2,8.16693e-3,
     3       4.80447e-4,1.04512e-3,-8.95778e-3,1.27411e-2,-3.37754e-4,
     4       4.17955e-5,-7.52181e-5,-1.38498e-4,5.92069e-4,-5.46719e-4/

      w(1)=1.
      m=0 
      n=1 
      do 2 jmax=1,4 
        n=n+1
        do 1 j=1,jmax 
          m=m+1
          w(n)=d*w(m)
          n=n+1
1       continue
        w(n)=t*w(m)
2     continue
      do 3 i=16,25
        w(i)=d*w(i-5) 
3     continue
      hlcnc=0
      do 4 i=1,25
        hlcnc=hlcnc+a(i)*w(i) 
4     continue
      return
      end 

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

      function hlcnh(d,t)
*
* Gilles Fontaine's fit to the Hubbard-Lampe Hydrogen conductive opacities
*
      implicit double precision(a-h,o-z)

      dimension a(5,5),temp(5), rho(5)
      data a/12.40928,     -2.545420,   -3.019045,    6.746942e-2,
     1    4.686838e-2,-15.29732,      .8154491,   2.387849, 
     2   -6.939705e-2, -4.129124e-2, 5.162517,   -3.839103e-1,
     3   -6.899530e-1,  2.568522e-2, 1.281765e-2, -.6347533,
     4    7.203151e-2,  8.405360e-2,-3.965926e-3,-1.652045e-3,
     5    2.699066e-2, -4.243979e-3,-3.648879e-3, 
     6    2.127303e-4,  7.509063e-5/

      sum = 0.
      temp(1) = 1.
      rho(1)  = 1.
      do 10 i = 2,5 
        temp(i) = temp(i-1) * t
        rho(i)  = rho(i-1) * d
10    continue
      do 2 i=1,5
        do 1 j=1,5
           sum = sum + a(i,j)*temp(j)*rho(i)
1       continue
2     continue
      hlcnh = sum
      return
      end 

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

      function hlcnhe(d,t)
*
* Gilles Fontaine's fit to the Hubbard-Lampe Helium conductive opacities
*
      implicit double precision(a-h,o-z)

      dimension a(5,5),temp(5), rho(5)
      data a /9.246777e+0, -2.097036,   -1.267003,  5.490678e-2,
     1    2.463979e-3, -1.265387e+1, 1.716563e-1,   1.088986e+0,
     2   -4.576479e-2, -6.966255e-3, 4.539412e+0, -1.427574e-1,
     3     -3.499339e-1,  1.589662e-2, 3.421386e-3, -5.842722e-1,
     4    3.930197e-2,  4.702301e-2,-2.529960e-3, -5.814372e-4,
     5    2.617848e-2, -2.791279e-3,-2.238884e-3, 
     6    1.442103e-4,  3.234929e-5/

      sum = 0.
      temp(1) = 1.
      rho(1)  = 1.
      do 10 i = 2,5 
        temp(i) = temp(i-1) * t
        rho(i)  = rho(i-1) * d
10    continue
      do 2 i=1,5
        do 1 j=1,5
           sum = sum + a(i,j)*temp(j)*rho(i)
1       continue
2     continue
      hlcnhe = sum 
      return
      end 

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

      function oconhl(d,t)
*
* get Hubbard-Lampe opacities for the appropriate composition
*
      implicit double precision(a-h,o-z)

      common/je/je
      common/comp/amhyhe,amheca,x(4)

      if ( je .eq. 3 ) then
         polytmp = hlcnc(d,t)
      elseif ( je .eq. 2 ) then
         polytmp = hlcnhe(d,t)
      elseif ( je .eq. 1 ) then
         polytmp = hlcnh(d,t)
      endif

      oconhl = polytmp

      return
      end 

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

      function oconit(d,t)
*
*  compute the conductive opacities using the method of itoh et al.
*  see top of program for references.
*
      implicit double precision(a-h,o-z)

      common/je/je
      common/comp/amhyhe,amheca,x(4)
      common/temp/s1,r1,s2,r2,b2,p2,t2,e2,xc2,xo2,fca2,f2,q2,w2,c
      data alac43 /-3.519188/ 
*
*  gamma constants
*
      data gamc3c, gamc3o, gamc2, gamc1 
     1     /3.5772e+06, 5.7782e+06, 573264., 2.275e+05/
*
* compute gamma, and call the appropriate routine (liquid metal or
* crystalline).  if gamma > 5000, then set gamma = 5000.
* This is the suggested proceedure from itoh et al (1984), since
* the results barely change when gamma > 5000.
* NOTE: there is no fit to xtal hydrogen.
*
      if ( je .eq. 3 ) then
         gammac = 10.**(d/3. - t) * gamc3c
         gammao = gammac / gamc3c * gamc3o
         gamma = xc2 * gammac + xo2 * gammao
*
* Carbon contribution
*
         if ( xc2 .gt. 0. ) then
            if ( gamma .lt. 171. ) then 
               condc = cndlc(d,t,gammac)
            elseif ( gamma .ge. 171 ) then
               if ( gammac .gt. 5000.) gammac = 5000.
               condc = cndsc(d,t,gammac)
            endif
         else
            condc = 0.
         endif
*
* Oxygen contribution
*
         if ( xo2 .gt. 0. ) then
            if ( gamma .lt. 171. ) then 
               condo = cndlo(d,t,gammao)
            elseif ( gamma .ge. 171 ) then
               if ( gammao .gt. 5000.) gammao = 5000.
               condo = cndso(d,t,gammao)
            endif
         else
            condo = 0.
         endif
         oconit = alac43 + (t*3. - d - (xc2*condc + xo2*condo) )
*
* Helium contribution
*
      elseif ( je .eq. 2 ) then
         gamma = gamc2 * 10.**(d/3. - t)
         if ( gamma .lt. 171. ) then
            conduct = cndlhe(d,t,gamma)
         elseif ( gamma .ge. 171 ) then 
            if ( gamma .gt. 5000.) gamma = 5000.
            conduct = cndshe(d,t,gamma)
         endif
         oconit = alac43 + (t*3. - d - conduct)
*
* Hydrogen contribution
*
      elseif ( je .eq. 1 ) then
         gamma = gamc1 * 10.**(d/3. - t)
         conduct = cndlhy(d,t,gamma)
         oconit = alac43 + (t*3. - d - conduct)
      endif

      return
      end 

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

      function ocon(d,t)
*
*  compute the conductive opacity of a specific rho-t point. 
*  note:  in this version, the hubbard lampe opacities have been replaced
*  by the conductive opacities of itoh et al (see top of code for
*  references.)
*
      implicit double precision(a-h,o-z)

      common/je/je
      common/comp/amhyhe,amheca,x(4)
*
*  the itoh opacities are only good down to log(rho) = 1.6
*  therefore, maw put in this kludge where if the density is
*  less than 1.8 and greater than 1.5, the interpolate between
*  opacity tables.  if d < 1.5, then use hubbard-lampe conductive
*  opacities.
*
      if (d .ge. 1.8) then
         ocon = oconit(d,t)
      elseif ( d .ge. 1.5 ) then
         dfrac = (d - 1.5) / 0.3
         ocon = dfrac * oconit(d,t) + (1.-dfrac) * oconhl(d,t)
      else
         ocon = oconhl(d,t)
      endif
      if(dabs(ocon).gt.40.)then
        ocon=dsign(40.d0,ocon)
      endif

      return
      end 

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

      function opac(d,t)
*
*  compute the total (rad. + cond.) opacity of a specific rho-t point. 
*  This function assumes a pure composition only. (See opac3 for mixtures.)
*  NOTE:  in this version, the Hubbard-Lampe opacities have been replaced
*  by the conductive opacities of itoh et al (see top of code for
*  references.)
*
      implicit double precision(a-h,o-z)

      common/opacs/cso(3,29,8),csr(29),cst(29)
      common/je/je
      common/comp/amhyhe,amheca,x(4)
      dimension two(2)
      data jrow,krow/8,29/
*
* trap point on staggered rho,T grid
*
      j=1 
      k=krow
6     i=(j+k)/2
      if ( t - cst(i) ) 7,7,8 
7     k=i 
      go to 9
8     j=i 
9     if(iabs(j-k)-1) 10,10,6 
10    do 14 i=1,2
      l=k+i-2
      delta=d-csr(l)+1
      j=delta+1
      if(j-1) 11,11,12
11    j=2 
      go to 14
12    if(j-jrow) 14,14,13
13    j=jrow
14    two(i)=cso(je,l,j-1)+(cso(je,l,j)-cso(je,l,j-1))*(delta-j+1)

      orad=two(1)+(two(2)-two(1))*(t-cst(k-1))/(cst(k)-cst(k-1))
      if(dabs(orad).gt.40.) orad=dsign(40.d0,orad)
      ocond = ocon(d,t)
      opac=dlog10(1./(10.**(-orad)+10.**(-ocond)))

      return
      end 

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

      function opac3(rhol,tl)
*
*  this acts like opac, but interpolates between opacities.
*  the vector x(1-3) must be set for this to work!
*  13 Nov 89 MAW
*  Added in calls to opacder under certain conditions. 4-91 PAB
*
      implicit double precision(a-h,o-z)
      common/opacfin/ifinal
      common/je/je
      common/comp/amhyhe,amheca,x(4)
      real*8 z(3),az(3)
*
* initialize constants
*
      az(1)=1.0
      az(2)=4.0
      az(3)=12.0
      z(1)=1.0
      z(2)=2.0
      z(3)=6.0
      opl = 0.
      jesav = je
*
* pure composition (comp. > 0.999) 
*
      do 50 i = 1,3
        xx = x(i)
        if(xx .gt. 0.999)then
          je = i
          if(ifinal.eq.0)then
 	    opac3 = opac(rhol,tl)
           else
            call opacder(rhol,tl,opcs,dlklr,dlklt)
            opac3 = opcs
          endif
          je = jesav
          return
        endif
50	continue
*
* mixtures
*
	do 100 i = 1,3
	   xje = x(i)
	   if ( xje.gt. 0.001) then
	     je = i
             opl = opl + xje * 10.**opac(rhol,tl)
	   endif
100	continue
	je = jesav
        opac3   = dlog10(opl)

999     return
	end

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

      function phasec(t)
*
*  is it liquid metal, or is it xtal lattice?
*  ac(6) contains the melting coefficients for carbon.
*
      implicit double precision(a-h,o-z)

      common/ac/ac(6)

      phasec = ac(6)
      do 1 i=1,5
         phasec = phasec * t + ac(6-i)
1     continue

      return
      end 

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

      function phaseo(t)
*
*  is it liquid metal, or is it xtal lattice?
*  ao(6) contains the melting coefficients for oxygen.
*
      implicit double precision(a-h,o-z)

      common/ao/ao(6)

      phaseo=ao(6)
      do 1 i=1,5
         phaseo = phaseo * t + ao(6-i)
1     continue

      return
      end 

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

      function rjay(zetfac,gnoo)       
*
* Function to determine value of J for each element
*
      implicit double precision(a-h,o-z)

      if(gnoo.ge.-20.0 .and. gnoo.lt.0.d0)then
	zfac = zetfac
        rjnum = (0.0151/zfac)+(0.242/zfac**2.25)+(1.21/zfac**4.55)
        rjden = 1.d0 + 0.0371*dexp(0.906*gnoo)*(1.d0+0.928*zfac)
	rjay = rjnum*dexp(gnoo)/rjden
      elseif(gnoo .le. 10.0)then
        den = 1.d0+(-0.0120*gnoo)+(0.029*gnoo*gnoo)+(-0.00104*gnoo**3)
	zfac = zetfac/den
        rjnum = (0.0123/zfac)+(0.266/zfac**2.25)+(1.30/zfac**4.55)
        rjden = 1.d0 + 0.117*dexp(0.897*gnoo)*(1.d0+0.177*zfac)
	rjay = rjnum*dexp(gnoo)/rjden
      endif

      return      
      end 

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