docs: def headgrid(self, xg, yg, layers=None, printrow=False):
          """Grid of heads
          
          Parameters
          ----------
          xg : array
              x values of grid
          yg : array
              y values of grid
          layers : integer, list or array, optional
              layers for which grid is returned
          printrow : boolean, optional
              prints dot to screen for each row of grid if set to `True`
          
          Returns
          -------
          h : array size `nlayers, ny, nx`
          
          See also
          --------
          
          :func:`~timml.model.Model.headgrid2`
  
          """
          
          nx, ny = len(xg), len(yg)
          if layers is None:
              Nlayers = self.aq.find_aquifer_data(xg[0], yg[0]).naq
          else:
              Nlayers = len(np.atleast_1d(layers))
          h = np.empty((Nlayers, ny, nx))
          for j in range(ny):
              if printrow:
                  print('.', end='', flush=True)
              for i in range(nx):
                  h[:, j, i] = self.head(xg[i], yg[j], layers)
          if printrow:
              print('', flush=True)
          return h