subroutine istat1(p2,t2,iw,iswch,converg) * * compute the thermodynamic variables for a specific P, T * point. This version searches either the C or O eos table, depending * on the composition and returns the thermodynamic properties. * 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/neut/en(5),fn common/iii/ iii common/ixswch/ixswch common/je/je logical converg dimension porder(4),torder(4),iorder(4),jorder(4), 1 xorder(4), yorder(4) dimension thermo(24) equivalence (d2,thermo(1)) data delta/2.e-3/ data fnat/4.342945e-1/ * * initialize some variables * xc = xc2 xo = xo2 je=3 iii=0 p=p2 t=t2 if ( iswch .eq. 12 ) then ci=1./(xc/12.+xo/16.) else if ( iswch .eq. 16 ) then ci=1./(xo/16.+(1-xo)/21.45) else stop endif * * find eos box and interpolate. * We search carbon if ISWCH is 12 and oxygen if ISWCH is 16. * if( iswch .eq. 12 ) then call wdselect(t,tminc,deltc,torder,iorder,1,itc,converg) kmax=6 if( iw .ge. 1 ) kmax=1 do 5 k=1,kmax do 4 i=1,4 indexc=iorder(i) pfrzc = phasec(t) 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 call wdselect(p,pstartc,delpc,porder,jorder,jminc,jmaxc, 1 converg) do 3 j=1,4 xorder(j)=tablec(k,indexc,jorder(j)) 3 continue call dali(p,porder,xorder,yorder(i)) 4 continue l=k if(k.gt.4) l=k+1 call dali(t,torder,yorder,thermo(l)) 5 continue * * or if oxygen, then come here * else if ( iswch .eq. 16 ) then call wdselect(t,tmino,delto,torder,iorder,1,ito,converg) kmax=6 if ( iw .ge. 1 ) kmax=1 do 35 k=1,kmax do 34 i=1,4 indexo=iorder(i) pfrzo = phaseo(t) pfrzc = pfrzo 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,pstarto,delpo,porder,jorder,jmino,jmaxo, 1 converg) do 33 j=1,4 xorder(j)=tableo(k,indexo,jorder(j)) 33 continue call dali(p,porder,xorder,yorder(i)) 34 continue l=k if(k.gt.4) l=k+1 call dali(t,torder,yorder,thermo(l)) 35 continue endif 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 xtmp = .999999 fcc=4.368*xtmp*xtmp*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) e = 10.**e et = e * et ep = (dt/fnat)*10.**(p-d-t) 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 ************************************************************************