subroutine isochr (aj,zlog,mlmin,mlmax)
*
*  computes an isochrone for t = aj and log(Z/0.02) = zlog, 
*  between log(mass) = mlmin and mlmax
*
*  Author: O. Pols     First release: 12 Dec 97
*                      Revised: 19 Mar 98
*
      implicit real*8 (a-h,l,m,o-z)
      parameter (nthk=21, ntzahb=71)
      parameter (delm=0.01, dell=0.1, delt=0.02)
*.....adopted solar log Teff, log g, and Mbol
      parameter (tefsun=3.762, gefsun=4.441, mblsun=4.75)
*
      parameter (ndi=500)
      integer jmod(ndi)
      real*8 lum(ndi), tef(ndi), rad(ndi), mass(ndi), mza(ndi), ml(ndi)
      real*8 mv(ndi), umb(ndi), bmv(ndi), vmr(ndi), vmi(ndi)
      common /isodat/ mza, mass, tef, lum, rad, mv, umb, bmv, vmr, vmi, 
     &     jmod, nk
*
      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
*
*.....interpolate evolution grid for Z
      call intplz (zlog)
*
*.....determine log of starting mass, ml(1), and its index in mlg, jm
      kk = 1
      ml(kk) = max(mlmin, mlg(1) - 0.02)
      jm = indx(ml(kk),mlg,nmg)
      jmz = indx(ml(kk),mlz,nmz)
*.....find starting model, jts
      jts = 2
 1    call fmasst (aj,jts,jts-1,mlt,lumt,teft,nt)
      if (mlt.lt.ml(kk) .and. jts.lt.nt) then
         jts = jts + 1
         goto 1
      end if
      lum(kk) = lumt
      tef(kk) = teft
*
*  start calculation of actual isochrone points:
*  for each model number, jt, first find the zero-age mass corresponding
*  to the age, then determine the other quantities by interpolation
* 
      do jt = jts, ndt
*........find log(z.a.mass) = mlt for which t = aj at model jt
         call fmasst (aj,jt,jt,mlt,lumt,teft,nt)
*........if jt exceeds the maximum model nt then finish calculation
         if(jt.gt.nt) goto 3
*........only compute new point if present mass exceeds the previous one
         if(mlt.ne.ml(kk)) then
            mlt = min(mlt,mlmax)
*...........number of masses between this and previous one, determined by
*...........the condition that the spacing in log M, log L and log Te are
*...........less that delm, dell and delt, resp.
            nm = int((mlt - ml(kk))/delm) + 1
            nlt = int(abs(lumt - lum(kk))/dell) + 1
            ntt = int(abs(teft - tef(kk))/delt) + 1
            nm = max(nm, nlt, ntt)
*...........however, do not interpolate between TFGB and ZAHB model
            if (jt.eq.ntzahb) nm = 1
            i1 = 1
            if (kk.eq.1) i1 = 0
            do i = i1, nm
               k = kk + i
               ml(k) = ml(kk) + (mlt - ml(kk))*float(i)/nm
*..............update mass index
 2             if (ml(k).gt.mlg(jm)) then
                  jm = jm + 1
                  goto 2
               end if
               dmh = (ml(k) - mlg(jm-1))/(mlg(jm) - mlg(jm-1))
               dml = 1.0 - dmh
               if(i.lt.nm) then
*.................interpolate in age for this mass
                  ajlo = ajg(jt-1,jm)**dmh * ajg(jt-1,jm-1)**dml
                  ajhi = ajg(jt,jm)**dmh   * ajg(jt,jm-1)**dml
                  dal = (ajhi - aj)/(ajhi - ajlo + 1.0e-6)
                  dah = 1.0 - dal
                  lum(k) = dah*(dmh*lumg(jt,jm)   + dml*lumg(jt,jm-1))
     &                   + dal*(dmh*lumg(jt-1,jm) + dml*lumg(jt-1,jm-1))
                  tef(k) = dah*(dmh*tefg(jt,jm)   + dml*tefg(jt,jm-1))
     &                   + dal*(dmh*tefg(jt-1,jm) + dml*tefg(jt-1,jm-1))
               else
*.................use values already determined in fmasst
                  lum(k) = lumt
                  tef(k) = teft
               end if
               rad(k) = 0.5*lum(k) - 2.0*(tef(k) - tefsun)
               mza(k) = 10.0**(ml(k))
*..............mass loss not included
               mass(k) = mza(k)
*..............correct MS points for detailed ZAMS shape
               if (jt.le.nthk) then
 4                if (ml(k).gt.mlz(jmz)) then
                     jmz = jmz + 1
                     goto 4
                  end if
                  dmzh = (ml(k) - 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)
                  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
*..............compute BC and UBVRI colours from Teff, log g and log Z
               glog = log10(mass(k)) - 2.0*rad(k) + gefsun
               call tgzubv (tef(k),glog,zlog,bc,umb(k),bmv(k),vmr(k),
     &              vmi(k))
               mbol = -2.5*lum(k) + mblsun
               mv(k) = mbol - bc
               jmod(k) = jt
            end do
            kk = kk + nm
*...........if mass exceeds maximum then finish calculation
            if (mlt.ge.min(mlg(nmg),mlmax) .or. jt.eq.nt) goto 3
         end if
      end do
*
 3    nk = kk
      return
      end
*
      subroutine fmasst (aj,jt,jl,ml,lum,tef,nt)
*
*  finds log(mass)=ml, that has t=aj at grid model jt
*  also returns lum=log(L) and tef=log(Te) for this mass, at model jl
*
      implicit real*8 (a-h,l,m,o-z)
      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
*
      jml = 1
      jmh = nmg
 1    if (jmh-jml.gt.1) then
         j = (jmh + jml)/2
         if (aj.lt.ajg(jt,j)) then
            jml = j
         else
            jmh = j
         end if
         goto 1
      end if
      nt = min(ntg(jmh), ntg(jmh-1))
      if(jt.gt.nt) return
*.....if model no. does not exceed maximum for this mass, interpolate:
      dah = log10(aj/ajg(jt,jmh-1))/log10(ajg(jt,jmh)/ajg(jt,jmh-1))
      dah = max(-0.2, min(1.0, dah))
      dal = 1.0 - dah
      ml = dah*mlg(jmh) + dal*mlg(jmh-1)
      lum = dah*lumg(jl,jmh) + dal*lumg(jl,jmh-1)
      tef = dah*tefg(jl,jmh) + dal*tefg(jl,jmh-1)
      end