      subroutine rdgrid (gtype)
*
*  read stellar evolution grids and store in common /evzgrd/
*
      implicit real*8 (a-h,l,m,o-z)
      parameter (ndmz=500,ndmg=30,ndt=200,ndz=7)
      character*3 gtype
      real*8 zgg(ndz)
*
      integer ntgg(ndmg,ndz)
      real*8 zlgg(ndz), mlzg(ndmz), lumzg(ndmz,ndz), tefzg(ndmz,ndz)
      real*8 mlgg(ndmg,ndz), lumgg(ndt,ndmg,ndz), tefgg(ndt,ndmg,ndz)
      real*8 ajgg(ndt,ndmg,ndz)
      common /evzgrd/ zlgg, mlzg, lumzg, tefzg, mlgg, lumgg, tefgg,
     &     ajgg, ntgg, nmzg, nmgg, nz
*
      character*2 direc
      data direc /'./'/
*
      open(10,file=direc//'EVOLGRID',status='old')
      read(10,'(1x,a3,8x,i5)') gtype, nz
      do k = 1, nz
         read(10,1011) zgg(k), nmgg
         zlgg(k) = log10(zgg(k)/0.02)
         do j = 1, nmgg
            read(10,1012) mass, ntgg(j,k)
            mlgg(j,k) = log10(mass)
            imp = 0
 1          read(10,1013) im, aj, tefgg(im,j,k), lumgg(im,j,k)
            ajgg(im,j,k) = aj + 1.0d-6
            do i = imp+1, im-1
               tefgg(i,j,k) = tefgg(im,j,k)
               lumgg(i,j,k) = lumgg(im,j,k)
               ajgg(i,j,k) = ajgg(im,j,k)
            end do
            imp = im
            if(im.lt.ntgg(j,k)) goto 1
            do i = ntgg(j,k)+1, ndt
               ajgg(i,j,k) = ajgg(i-1,j,k)
            end do
         end do
*........fill in some gaps between grid masses by interpolation 
*........(mainly, post-HeF models)
         do j = 2, nmgg-1
            ntgx = min(ntgg(j-1,k), ntgg(j+1,k))
            if(ntgg(j,k).lt.ntgx) then
               dmh = (mlgg(j,k) - mlgg(j-1,k))/(mlgg(j+1,k) -
     &              mlgg(j-1,k))
               dml = 1.0 - dmh
               do i = ntgg(j,k)+1, ntgx
                  tefgg(i,j,k) = dml*tefgg(i,j-1,k) + dmh*tefgg(i,j+1,k)
                  lumgg(i,j,k) = dml*lumgg(i,j-1,k) + dmh*lumgg(i,j+1,k)
                  ajgg(i,j,k) = ajgg(ntgg(j,k),j,k) +
     &                 ((ajgg(i,j-1,k) - ajgg(ntgg(j,k),j-1,k))**dml *
     &                 (ajgg(i,j+1,k) - ajgg(ntgg(j,k),j+1,k))**dmh)
               end do
               ntgg(j,k) = ntgx
               do i = ntgg(j,k)+1, ndt
                  ajgg(i,j,k) = ajgg(i-1,j,k)
               end do
            end if
         end do
      end do
      close(10)
 1011 format (4(/),4x,f8.4,i5)
 1012 format (/,4x,f8.4,i5,/)
 1013 format (i3,2x,1p,e16.9,0p,5f8.4,3f7.4,10f7.5)
*
      open(11,file=direc//'ZAMSGRID',status='old')
      read(11,'(12x,i5)') nzz
      if (nzz.ne.nz) pause '# Z''s unmatched'
      do k = 1, nzz
         read(11,1001) z, nmzg
         if (z.ne.zgg(k)) pause 'Z''s unequal'
         do i = 1, nmzg
            read(11,1002) idum, mass, tefzg(i,k), lumzg(i,k)
            if (k.gt.1 .and. log10(mass).ne.mlzg(i)) 
     &           pause 'masses unequal'
            mlzg(i) = log10(mass)
         end do
      end do
      close(11)
 1001 format (4(/),4x,f8.4,i5,/)
 1002 format (i3,f9.4,3f8.4,3f7.4,5f7.5)
*
      return
      end
