subroutine istatco(p2,t2,iw,converg)
*
* compute the thermodynamic variables for a specific P, T
* point.  This version searches both the c and o eos tables,
* and then interpolates to the desired composition using the
* additive volume technique
*
      implicit double precision(a-h,o-z)

      common/temp/s1,r1,s2,r2,b2,px2,tx2,e2,xc2,xo2,fca2,fx2,q2,w2,c
      common/tablec/tablec(6,39,87)
      common/tableo/tableo(6,36,87)
      common/sizec/tminc,tmaxc,deltc,pminc,pmaxc,delpc,pmc(50),ipc(50),
     1   iphasec(50),ioverc,itc
      common/sizeo/tmino,tmaxo,delto,pmino,pmaxo,delpo,pmo(50),ipo(50),
     1   iphaseo(50),iovero,ito
      common/thermo/d2,dp,dt,e,ep,et,psi,pg,o2,op,ot,fp,ft, 
     1   en2(5),fn2,fcc2,ce,ci,cif,w,nu 
      common/thermoc/dc,dpc,dtc,ec,epc,etc,psic
      common/thermoo/do,dpo,dto,eo,epo,eto,psio
      common/neut/en(5),fn
      common/iii/ iii
      common/ixswch/ixswch
      common/je/je

      logical converg
      dimension porderc(4),torderc(4),iorderc(4),jorderc(4),
     1          xorderc(4), yorderc(4)
      dimension pordero(4),tordero(4),iordero(4),jordero(4),
     1          xordero(4), yordero(4)
      dimension thermoc(7), thermoo(7)
      equivalence (dc,thermoc(1))
      equivalence (do,thermoo(1))
      data delta/2e-3/
      data fnat/4.342945e-1/
*
* initialize some variables
*
      xc = xc2
      xo = xo2
      je=3
      iii=0
      ci=1/(xc/12.+xo/16.)
      p=p2
      t=t2
*
* find eos box and interpolate.
* note that we're doing both tables in parallel.
*
      call wdselect(t,tminc,deltc,torderc,iorderc,1,itc,converg)
      call wdselect(t,tmino,delto,tordero,iordero,1,ito,converg)
 
      kmax=6
      if ( iw .ge. 1 )  kmax=1
 
      do 5 k=1,kmax 
 
         do 4 i=1,4 
 
            indexc=iorderc(i) 
            indexo=iordero(i) 
 
            if ( ixswch .eq. 1 ) then
               pfrzc = phasec(t)
               pfrzo = pfrzc
            elseif ( ixswch .eq. 2 ) then
               pfrzo = phaseo(t)
               pfrzc = pfrzo
            elseif ( ixswch .eq. 3 ) then
               pfrzc = xc * phasec(t) + xo * phaseo(t)
               pfrzo = pfrzc
            elseif ( ixswch .eq. 4 ) then
               pfrzc = phasec(t)
               pfrzo = phaseo(t)
            else
               stop 
            endif
 
            if ( p .lt. pfrzc ) then
               jminc=iphasec(indexc)+1
               jmaxc=ipc(indexc)
               pstartc=pmaxc-ioverc*delpc
            else
               jminc=1
               jmaxc=iphasec(indexc)
               pstartc=pmaxc
            endif
 
            if ( p .lt. pfrzo ) then
               jmino=iphaseo(indexo)+1
               jmaxo=ipo(indexo)
               pstarto=pmaxo-iovero*delpo
            else
               jmino=1
               jmaxo=iphaseo(indexo)
               pstarto=pmaxo
            endif
 
            call wdselect(p,pstartc,delpc,porderc,jorderc,jminc,jmaxc,
     1          converg)
            call wdselect(p,pstarto,delpo,pordero,jordero,jmino,jmaxo,
     1          converg)
 
            do 3 j=1,4
               xorderc(j)=tablec(k,indexc,jorderc(j))
               xordero(j)=tableo(k,indexo,jordero(j))
3           continue
 
            call dali(p,porderc,xorderc,yorderc(i))
            call dali(p,pordero,xordero,yordero(i))
 
4        continue
 
         l=k
         if(k.gt.4) l=k+1
 
         call dali(t,torderc,yorderc,thermoc(l))
         call dali(t,tordero,yordero,thermoo(l))
 
5     continue
*
* c/o interpolations finished.  at this time we've interpolated
* in the c and o tables to get d, dp, dt, e, et, and psi for each
* composition.
* use additive volume technique to
* interpolate the therm quantities for the mixture.
* see Fontaine, Graboske, & Van Horn 1977, ApJS, 35, 293
* see my log 4 for derivation of these interpolations.
*
      d2 = -dlog10( xc/10.**dc + xo/10.**do )
      d  = d2
*
* Now that we've interpolated the main quantities, we need to
* calculate the opacity, carbon burning and neutrino luminosities
* and their derivitives wrt rho and t
* Note that the carbon burning luminosity is reduced by a factor
* of 5 arbitrarily (otherwise, runaway).
*
      n  = -1
20    x1=49.4+d-2*t/3-36550*(1+7e-1*10.**(t-10))**(1./3)*10.**(-t/3)
 
      if ( x1 .le. -30. ) then
         fcc=0
      else
         fcc=4.368*xc2*xc2*10.**x1
         fcc=fcc/5. 
      endif
 
      if (iw .ge. 1) then
         fcc2=fcc
         return
      endif
 
      if (nu .le. 0) then
         fn=0
      elseif(nu .eq. 1)then
*
* call venerable BPS neutrino rates
*
         call neu(d,t)
      elseif(nu .eq. 2)then
*
* call Munkata etal Plasmon, Photo, and Pair rates
*
         call neutrino(d,t)
      elseif(nu .eq. 3)then
*
* call Itoh etal neutrino rates
*
         call gnutrino(d,t)
      endif
 
      f=fn-fcc*w
      o=10.**opac(d,t)
      if(n) 14,15,16

14    o2=o
      f2=f
      fn2=fn
      fcc2=fcc
      do 17 i=1,5
17    en2(i)=en(i)
      dt  = 1./( d * ( xc/(dc*dtc) + xo/(do*dto) ) )
      dp  = 1./( d * ( xc/(dc*dpc) + xo/(do*dpo) ) )
      e  = xc * 10.**ec + xo * 10.**eo
      et = xc * etc * 10.**ec + xo * eto * 10.**eo
      ep = 10.**(p-t)/fnat * (xc*dtc*10.**(-dc) + xo*dto*10.**(-do))
      psi = xc * psic + xo * psio
      pg=dlog10(10.**p-10.**(4*t-14.59836))
      d=d2+delta
      n=0 
      go to 20

15    od=(o-o2)/delta
      fd=(f-f2)/delta
      d=d2
      t=t2+delta
      n=1 
      go to 20

16    ot=(o-o2)/delta
      ft=(f-f2)/delta
      op=od*dp
      ot=ot+od*dt
      fp=fd*dp
      ft=ft+fd*dt

      return
      end 

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