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 

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