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 ************************************************************************