subroutine track (mlog,zlog)
*
*  interpolate an evolution track for log(M/Msol) = mlog and 
*  log(Z/Zsol) = zlog
*
      implicit real*8 (a-h,l,m,o-z)
      parameter (nthk=21, ntzahb=71)
*.....adopted solar log Teff, log g, and Mbol
      parameter (tefsun=3.762, gefsun=4.441, mblsun=4.75)
*
      parameter (ndmz=500,ndmg=30,ndt=200)
      integer ntg(ndmg)
      real*8 mlz(ndmz), lumz(ndmz), tefz(ndmz)
      real*8 mlg(ndmg), lumg(ndt,ndmg), tefg(ndt,ndmg), ajg(ndt,ndmg)
      common /evgrid/ mlz, lumz, tefz, mlg, lumg, tefg, ajg, ntg, nmz,
     &     nmg
*
      real*8 ag(ndt), lum(ndt), tef(ndt), rad(ndt), mass(ndt)
      real*8 mv(ndt), umb(ndt), bmv(ndt), vmr(ndt), vmi(ndt), glog(ndt)
      common /trkdat/ ag, tef, lum, rad, mass, mv, umb, bmv, vmr, vmi,
     &     nt
*
*.....interpolate evolution grid for Z
      call intplz (zlog)
*
*.....interpolate for the mass
      jm = indx(mlog,mlg,nmg)
      nt = min(ntg(jm), ntg(jm-1))
      dmh = (mlog - mlg(jm-1))/(mlg(jm) - mlg(jm-1))
      dml = 1.0 - dmh
      jmz = indx(mlog,mlz,nmz)
      dmzh = (mlog - mlz(jmz-1))/(mlz(jmz) - mlz(jmz-1))
      dmzl = 1.0 - dmzh
      lum0 = dmzh*lumz(jmz) + dmzl*lumz(jmz-1)
      tef0 = dmzh*tefz(jmz) + dmzl*tefz(jmz-1)
      rad0 = 0.5*lum0 - 2.0*(tef0 - tefsun)
      lum1 = dmh*lumg(1,jm) + dml*lumg(1,jm-1)
      tef1 = dmh*tefg(1,jm) + dml*tefg(1,jm-1)
      rad1 = 0.5*lum1 - 2.0*(tef1 - tefsun)
      lum2 = dmh*lumg(nthk,jm) + dml*lumg(nthk,jm-1)
      tef2 = dmh*tefg(nthk,jm) + dml*tefg(nthk,jm-1)
      rad2 = 0.5*lum2 - 2.0*(tef2 - tefsun)
      do k = 1, nt
         ag(k) = ajg(k,jm)**dmh * ajg(k,jm-1)**dml
         xm = 10.0**mlog
         lum(k) = dmh*lumg(k,jm)   + dml*lumg(k,jm-1)
         tef(k) = dmh*tefg(k,jm)   + dml*tefg(k,jm-1)
         rad(k) = 0.5*lum(k) - 2.0*(tef(k) - tefsun)
c........correct for ZAMS shape
         if (k.le.nthk) then
            hl = (lum(k) - lum2)/(lum1 - lum2)
            hr = (rad(k) - rad2)/(rad1 - rad2)
            lum(k) = hl*lum0 + (1.0 - hl)*lum2
            rad(k) = hr*rad0 + (1.0 - hr)*rad2
            tef(k) = 0.25*lum(k) - 0.5*rad(k) + tefsun
         end if
         mass(k) = xm
         glog(k) = log10(mass(k)) - 2.0*rad(k) + gefsun
         call tgzubv (tef(k),glog(k),zlog,bc,umb(k),bmv(k),vmr(k),
     &        vmi(k))
         mbol = -2.5*lum(k) + mblsun
         mv(k) = mbol - bc
      end do
      end