subroutine estate * * 3-point lagragian interpolation in the pp-tk grid to find thermo * quantities el1,rl1,ol1,atg,rtg,drdtp,cp [PROFILE: 23% runtime] * implicit double precision(a-h,o-z) common/tenv/rhok(3,60,35),eta(3,60,35),pp(3,60,35),datg(3,60,35), 1 od(3,18,35),tk(3,60),uu(3,60,35),xtt(3,60,35),xrt(3,60,35) common/je/je common/comp/amhyhe,amheca,x(4) common/c/az(3),z(3),taux,ts,ps,sm,rs,err,te,gs,rho,ol1x,rl1x,el1x 1 ,atgx,rtgx,tl,pl,xmass,acu,am,bk,an0,pi common/d/ ml(3,60),nmax,kind,lmax,nsuff,nlim,nlim1,j,iter,lmam, 1 nlum,ko,nhp common/grad/ drdtp(400),cp(400),ttg(400),drdtpx,cpx,pecx,ttgx, 1 deptx,dpdr,dtdr,dmdr,dr,xr,xt,cv,drdptx dimension xtl(3),xrl(3),el(3),rl(3),at(3) integer k1(3),k2(3),k3(3) * * trap point in t, after checking if = low boundry or less, high * boundry or more * if ( tl .lt. tk(je,2) ) then l1=1 l2=2 l3=3 else if ( tl .gt. tk(je,lmax-1) ) then l3=lmax l2=lmax-1 l1=lmax-2 else * * The last 6 isotherms are at a different * spacing in log T of 0.25 in the current incarnation of the * 'envelope' EOS extended to high temperatures. * if (tl.lt.7.54d0) then l3=int((tl-tk(je,1))/(tk(je,2)-tk(je,1)))+2 l2=l3-1 l1=l3-2 else l3=int((tl-tk(je,54))/(tk(je,55)-tk(je,54)))+56 l2=l3-1 l1=l3-2 endif endif * * get fractional distance between mesh points in temp * lagrange coefficients. * dtdt=(tl-tk(je,l2))/(tk(je,l3)-tk(je,l2)) dt1=tk(je,l1)-tl dt2=tk(je,l2)-tl dt3=tk(je,l3)-tl dtdt1=dt2*dt3/((dt1-dt2)*(dt1-dt3)) dtdt2=dt1*dt3/((dt2-dt1)*(dt2-dt3)) dtdt3=dt1*dt2/((dt3-dt1)*(dt3-dt2)) * * Now trap point in p * Note, we come loop back to 8 for l = l1, l2, l3 * This is a *slow* (linear) search thru the grid in pressure * and I wouldn't be surprised if the program spends a large * fraction of its time right here * l=l1 8 kmax=ml(je,l) lp=l+1-l1 if ( pl .lt. pp(je,l,2) ) then k1(lp)=1 k2(lp)=2 k3(lp)=3 else do 16 k=3,kmax if ( pl .gt. pp(je,l,k) ) then goto 16 else k3(lp)=k k2(lp)=k-1 k1(lp)=k-2 go to 14 endif 16 continue k3(lp)=kmax k2(lp)=kmax-1 k1(lp)=kmax-2 endif 14 kl1=k1(lp) kl2=k2(lp) kl3=k3(lp) dp1=pp(je,l,kl1)-pl dp2=pp(je,l,kl2)-pl dp3=pp(je,l,kl3)-pl dpdp1=(pl-pp(je,l,kl2))*(pl-pp(je,l,kl3))/(pp(je,l,kl1)- 1 pp(je,l,kl2))/(pp(je,l,kl1)-pp(je,l,kl3)) dpdp2=(pl-pp(je,l,kl1))*(pl-pp(je,l,kl3))/(pp(je,l,kl2)- 1 pp(je,l,kl1))/(pp(je,l,kl2)-pp(je,l,kl3)) dpdp3=(pl-pp(je,l,kl1))*(pl-pp(je,l,kl2))/(pp(je,l,kl3)- 1 pp(je,l,kl2))/(pp(je,l,kl3)-pp(je,l,kl1)) * * Interpolate in P for each T point * el(lp)=eta(je,l,kl1)*dpdp1+eta(je,l,kl2)*dpdp2+eta(je,l,kl3)* 1 dpdp3 at(lp)=datg(je,l,kl1)*dpdp1+datg(je,l,kl2)*dpdp2+datg(je,l,kl3)* 1 dpdp3 rl(lp)=rhok(je,l,kl1)*dpdp1+rhok(je,l,kl2)*dpdp2+rhok(je,l,kl3)* 1 dpdp3 xtl(lp)=xtt(je,l,kl1)*dpdp1+xtt(je,l,kl2)*dpdp2+xtt(je,l,kl3)* 1 dpdp3 xrl(lp)=xrt(je,l,kl1)*dpdp1+xrt(je,l,kl2)*dpdp2+xrt(je,l,kl3)* 1 dpdp3 if ( lp .ge. 3 ) goto 19 l=l+1 go to 8 * * this is the bottom of the loop * Now interpolate in T for actual pressure. * 19 el1x=el(2)+0.5d0*dtdt*(el(3)-el(1))+dtdt*dtdt 1 *(0.5d0*(el(1)+el(3))-el(2)) atgx=at(2)+0.5d0*dtdt*(at(3)-at(1))+dtdt*dtdt 1 *(0.5d0*(at(1)+at(3))-at(2)) rl1x=rl(2)+0.5d0*dtdt*(rl(3)-rl(1))+dtdt*dtdt 1 *(0.5d0*(rl(1)+rl(3))-rl(2)) xt1x=xtl(2)+0.5d0*dtdt*(xtl(3)-xtl(1))+dtdt*dtdt 1 *(0.5d0*(xtl(1)+xtl(3))-xtl(2)) xr1x=xrl(2)+0.5d0*dtdt*(xrl(3)-xrl(1))+dtdt*dtdt 1 *(0.5d0*(xrl(1)+xrl(3))-xrl(2)) * * compute final thermo properties * xt = xt1x xr = xr1x drdtpx = xt/xr drdptx = -1.0d0/xr rho = 10.d0**rl1x cpx = drdtpx/rho/atgx*10.**(pl-tl) cv = cpx*(1.0d0-(xt*atgx)) return end ************************************************************************