# Copyright (C) 2011, 2014 David Maxwell
# 
# This file is part of PISM.
# 
# PISM is free software; you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation; either version 3 of the License, or (at your option) any later
# version.
# 
# PISM is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
# 
# You should have received a copy of the GNU General Public License
# along with PISM; if not, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

import PISM

class ModelVecs:
  """A dictionary of :cpp:class:`IceModelVec`\s containing model data.  This plays
  a similar role as the C++ :cpp:class:`PISMVars` with the key difference that it keeps
  references to its contents.  By contrast, if an :cpp:class:`IceModelVec` is created in python
  and added to a :cpp:class:`PISMVars`, and if no other reference to the python object is kept, the
  :cpp:class:`IceModelVec` will be deleted, leaving a null pointer lurking in the :cpp:class:`PISMVars`.  
  See method :meth:`asPISMVars` to generate a :cpp:class:`PISMVars` from the contents of a :class:`ModelVecs`

  The lookup name associated with a vector in a :class:`ModelVecs` has no meaning outside the context of the :class:`ModelVecs`.
  In particular, when converting to a :cpp:class:`PISMVars`, the lookup names are by default simply the intrinsic
  names of the vectors given at the time of their creation::
  
    # Construct a tau_c vector with intrinsic name 'tauc'
    yieldstress = PISM.IceModelVec2S()
    yieldstress.create(grid, 'tauc', PISM.WITHOUT_GHOSTS);
    yieldstress.set_attrs("diagnostic", desc, "Pa", ""); 

    # Construct a thickness vector with intrinsic name 'thk'
    thickness = PISM.IceModelVec2S();
    thickness.create(grid, 'thk', PISM.WITHOUT_GHOSTS);
    thickness.set_attrs("model_state", "land ice thickness", "m", "land_ice_thickness");
    thickness.metadata().set_double("valid_min", 0.0);
    
    vecs = PISM.model.ModelVecs()
    
    # Add 'yieldstress' to be accessed via the given name 'basal_yield_stress'
    vecs.add(yieldstress,'basal_yield_stress')

    # Add 'thickness' to be accessed by its intrinsic name, (i.e. 'thk')
    vecs.add(thickness)

    # Access the vectors
    v1 = vecs.basal_yield_stress
    v2 = vecs.thk
    
    # Convert to a PISMVars
    pismVars = vecs.asPISMVars()
    
    # At this point, pismVars contains two vectors that can be accessed by their intrinsic names.
    w1 = vecs.get('tauc')
    w2 = vecs.get('thk')

  See also :meth:`setPISMVarsName` for associating an alternate name for a member vector to be given when 
  converting to a :cpp:class:`PISMVars`.  And see also functions :func:`createBasalYieldStressVec` and friends
  for one step creation of standard PISM vecs.
  """

  def __init__(self):
    self.needs_writing = set()
    self._vecs = {}
    self.pvarnames = {}

  def __getattr__(self,key):
    """Allows access to contents via dot notation, e.g. `vecs.thickness`."""
    v = self._vecs.get(key)
    if v is None:
      raise AttributeError(key)
    return v

  def __repr__(self):
    """:returns: a string representation of the :class:`ModelVecs` listing its contents."""
    s= 'ModelVecs:\n'
    for (key,value) in self._vecs.items():
      s += "    %s: %s %s\n" %(key,value,value.string_attr("name"))
    return s

  def get(self,name):
    """:returns: the vector with the given *name*, raising an :exc:AttributeError if one is not found."""
    v = self._vecs.get(name)
    if v is None:
      raise AttributeError(name)
    return v

  def add(self,var,name=None,writing=False):
    """Adds a vector to the collection.
    
    :param var: The :cpp:class:`IceModelVec` to add.
    :param name: The name to associate with `var`.  If `None`, then
                 the `var`'s `name` attribute is used.
    :param writing: A boolean. `True` indicates the vector is among those that should be written to a file if
                    :meth:`write` is called."""
              
    if name is None:
      name = var.name()
    if self._vecs.has_key(name):
      raise RuntimeError("Cannot add variable %s to model.ModelVecs; it is already present" % name)
    self._vecs[name] = var
    if writing:
      self.needs_writing.add(var)
  
  def rename(self,old_name,new_name):
    var = self._vecs.pop(old_name)
    self.add(var,new_name)

  def remove(self,name):
    v = self._vecs.pop(name)
    if v in self.needs_writing:
      self.needs_writing.remove(v)
    if self.pvarnames.has_key(name):
      self.pvarnames.remove(name)

  def has(self,var):
    """:returns: `True` if the dictionary contains a vector with the given name."""
    return self._vecs.has_key(var)

  def items( self ):
    """Iterator returning key-value pairs of the contents."""
    for (name,var) in self._vecs.items():
      yield (name,var)

  def asPISMVars(self):
    """
    :returns: a :cpp:class:`PISMVars` containing all the vecs in the :class:`ModelVecs`.
    
    The name associated with each vector in the :cpp:class:`PISMVars` will be its intrinsic name
    unless this has been overridden with :meth:`setPISMVarsName`::

      # Construct a tau_c vector
      yieldstress = PISM.IceModelVec2S()
      yieldstress.create(grid, 'tauc', PISM.WITHOUT_GHOSTS);
      yieldstress.set_attrs("diagnostic", desc, "Pa", ""); 

      # Add it to a ModelVecs
      vecs = PISM.model.ModelVecs()
      vecs.add(yieldstress,'basal_yield_stress')

      # Access the vector
      v = vecs.basal_yield_stress
      
      pismVars = vecs.asPISMVars()
      # Access the vector from the pismVars
      v2 = pismVars.get('tauc')

    Note that the key in ``pismVars`` is ``tauc``, not ``basal_yield_stress``.
    """
    pismVars = PISM.PISMVars()
    for (name,var) in self.items():
      pv_name = self.pvarnames.get(name)
      if pv_name is not None:
        pismVars.add(var,pv_name)
      else:
        pismVars.add(var)
    return pismVars

  def setPISMVarsName(self,varname,pvarname):
    """
    Associates a name that a member variable should have if the :class:`ModelVecs` is converted into
    a :cpp:class:`PISMVars`.  For example::

      # Construct a tau_c vector with intrinsic name 'tauc'
      yieldstress = PISM.IceModelVec2S()
      yieldstress.create(grid, 'tauc', PISM.WITHOUT_GHOSTS);
      yieldstress.set_attrs("diagnostic", desc, "Pa", ""); 

      vecs = PISM.model.ModelVecs()

      # Add 'yieldstress' to be accessed via the given name 'basal_yield_stress'
      vecs.add(yieldstress,'basal_yield_stress')
      vecs.setPISMVarsName('basal_yield_stress','foo')

      # Access
      v = vecs.basal_yield_stress

      # Convert to a PISMVars
      pismVars = vecs.asPISMVars()

      # Access from the PISMVars
      w = pismVars.get("foo") 
    """
    if not self.has(varname):
      raise RuntimeError("Cannot set export name for %s to %s: variable %s does not exist." %(varname,pvarname,varname))
    self.pvarnames[varname]=pvarname

  def write(self,output_filename):
    """Writes any member vectors that have been flagged as need writing to the given :file:`.nc` filename."""
    vlist = [ v for v in self.needs_writing ]
    vlist.sort(cmp=var_cmp)
    for v in vlist:
      v.write(output_filename)

  def writeall(self,output_filename):
    """Writes all member vectors to the given :file:`.nc` filename."""
    vlist = [ var for (name,var) in self.items() ]
    vlist.sort(cmp=var_cmp)
    for v in vlist:
      v.write(output_filename)

  def markForWriting(self,var):
    """Marks a given vector as needing writing.
    
    :param `var`: Either the name of a member vector, or the vector itself."""
  
    if isinstance(var,str):
      var = self.get(var)
    self.needs_writing.add(var)

# There must be a better way to do this.
def var_cmp(v1,v2):
  ":cpp:class:`IceModelVec` name comparison function used for sorting a list of vecs."
  n1 = v1.name()
  n2 = v2.name()
  if n1 < n2:
    return -1
  if n1 > n2:
    return 1
  return 0


class ModelData:
  """
    A collection of data and physics comprising a PISM model.  The collection consists of

       * ``grid``      the computation grid
       * ``config``    the associated :cpp:class:`PISMConfig`
       * ``basal``     a :cpp:class:`PISMYieldStress`.
       * ``enthalpy``  an :cpp:class:`EnthalpyConverter`.
       * ``vecs``      a :class:`ModelVecs` containing vector data.
    
    These member variables are intended to be accessed as public attributes."""
  
  def __init__(self,grid,config=None):
    self.grid = grid
    if config is None: config = grid.config
    self.config = config

    # Components for the physics
    self.basal = None               #: Attribute for a  :cpp:class:`PISMYieldStress`
    self.enthalpyconverter = None   #: Attribute for an :cpp:class:`EnthalpyConverter`

    self.vecs = ModelVecs()         #: The list 

  def setPhysics(self,enthalpyconverter):
    """Sets the physics components of the model.
    
    :param enthalpyconverter: An instance of :cpp:class:`EnthalpyConverter`    
    """

    self.enthalpyconverter = enthalpyconverter;

def initShallowGrid(grid,Lx,Ly,Mx,My,p):
  """
  Initialize a :cpp:class:`IceGrid` intended for 2-d computations.
  
  :param grid: The grid to initialize.
  :param Lx:   Half width of the ice model grid in x-direction (m)
  :param Ly:   Half width of the ice model grid in y-direction (m)
  :param Mx:   number of grid points in the x-direction
  :param My:   number of grid points in the y-direction
  :param p:    grid periodicity ( one of ``PISM.NOT_PERIODIC``, ``PISM.X_PERIODIC``, 
                ``PISM.Y_PERIODIC``, or ``PISM.XY_PERIODIC``)  
  """ 
  grid.Lx = Lx;
  grid.Ly = Ly;
  grid.periodicity = p;
  grid.Mx = Mx; grid.My=My; grid.Mz=3;

  grid.compute_nprocs();
  grid.compute_ownership_ranges();
  grid.compute_vertical_levels()
  grid.compute_horizontal_spacing();
  grid.allocate()

def initGrid(grid,Lx,Ly,Lz,Mx,My,Mz,p):
  """Initialize a :cpp:class:`IceGrid` intended for 3-d computations.
  
  :param grid: The grid to initialize.
  :param Lx:   Half width of the ice model grid in x-direction (m)
  :param Ly:   Half width of the ice model grid in y-direction (m)
  :param Lz:   Extent of the grid in the vertical direction (m)
  :param Mx:   number of grid points in the x-direction
  :param My:   number of grid points in the y-direction
  :param Mz:   number of grid points in the z-direction
  :param p:    grid periodicity ( one of ``PISM.NOT_PERIODIC``, ``PISM.X_PERIODIC``, 
                ``PISM.Y_PERIODIC``, or ``PISM.XY_PERIODIC``)  
  """ 
  grid.Lx = Lx;
  grid.Ly = Ly;
  grid.Lz = Lz;
  grid.periodicity = p;
  grid.Mx = Mx; grid.My=My; grid.Mz = Mz;

  grid.compute_nprocs();
  grid.compute_ownership_ranges();
  grid.compute_vertical_levels()
  grid.compute_horizontal_spacing();
  grid.allocate()

def initGridFromFile(grid,input_file,periodicity,variable='enthalpy'):
  """Initialize a :cpp:class:`IceGrid` to match the variable dimensions in an :file:`.nc` file.
  
  :param grid:         The grid to initialize.
  :param input_file:   Name of ``.nc`` file.
  :param periodicity:  Desired periodicity ( one of ``PISM.NOT_PERIODIC``, ``PISM.X_PERIODIC``, 
                ``PISM.Y_PERIODIC``, or ``PISM.XY_PERIODIC``).  This parameter cannot be inferred from
                `input_file` and must be specified.
  :param variable:     Name of variable in `input_file` used to infer grid dimensions.
  """ 
  pio = PISM.PIO(grid.com, "netcdf3", grid.get_unit_system())

  pio.open(input_file, PISM.PISM_NOWRITE)
  pio.inq_grid(variable, grid, periodicity)
  pio.close()

  grid.compute_nprocs();
  grid.compute_ownership_ranges();
  grid.allocate()



def createIceSurfaceVec(grid,name='usurf',ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for surface elevation data.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0

  surface = PISM.IceModelVec2S();
  surface.create(grid, name, ghost_type, stencil_width)
  surface.set_attrs("diagnostic", "ice upper surface elevation", "m", "surface_altitude");
  return surface;

def createIceThicknessVec(grid,name='thk',ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for ice thickness data.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0

  thickness = PISM.IceModelVec2S();
  thickness.create(grid, name, ghost_type, stencil_width);
  thickness.set_attrs("model_state", "land ice thickness", "m", "land_ice_thickness");
  thickness.metadata().set_double("valid_min", 0.0);
  return thickness

def createIceSurfaceStoreVec(grid,name='usurfstore',ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for saved surface elevation data in regional models.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0

  usurfstore = PISM.IceModelVec2S()
  usurfstore.create(grid, name, ghost_type, stencil_width)
  usurfstore.set_attrs("model_state", 
    "saved surface elevation for use to keep surface gradient constant in no_model strip",
    "m", "");
  return usurfstore

def createIceThicknessStoreVec(grid,name='thkstore',ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for saved ice thickness data in regional models.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0

  thkstore = PISM.IceModelVec2S();
  thkstore.create(grid, name, ghost_type, stencil_width)
  thkstore.set_attrs("model_state", 
    "saved ice thickness for use to keep driving stress constant in no_model strip",
    "m", "")
  return thkstore

def createBedrockElevationVec(grid,name='topg',desc="bedrock surface elevation",
                                ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for bedrock elevation data.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""

  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0

  bed = PISM.IceModelVec2S()
  bed.create(grid, name, ghost_type, stencil_width);
  bed.set_attrs("model_state", desc, "m", "bedrock_altitude");
  return bed

def createYieldStressVec(grid,name='tauc',desc="yield stress for basal till (plastic or pseudo-plastic model)",
                          ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for basal yield stress data.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""
  # yield stress for basal till (plastic or pseudo-plastic model)
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0
  
  tauc = PISM.IceModelVec2S()
  tauc.create(grid, name, ghost_type, stencil_width);
  tauc.set_attrs("diagnostic", desc, "Pa", ""); 
  return tauc;


def createAveragedHardnessVec(grid,name='hardav',desc="vertical average of ice hardness",
                          ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for vertically-averaged ice-hardness data.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""
  # yield stress for basal till (plastic or pseudo-plastic model)
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0

  hardav = PISM.IceModelVec2S()
  hardav.create(grid, name, ghost_type, stencil_width);

  power = 1.0 / grid.config.get("Glen_exponent");
  units = "Pa s%f" % power

  hardav.set_attrs("diagnostic", desc, units, ""); 

  hardav.metadata().set_double("valid_min", 0.0);

  return hardav;

def createEnthalpyVec(grid,name='enthalpy',ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec3` with attributes setup for enthalpy data.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""


  enthalpy = PISM.IceModelVec3()
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  
  enthalpy.create(grid, name,  ghost_type, stencil_width);
  enthalpy.set_attrs("model_state", "ice enthalpy (includes sensible heat, latent heat, pressure)", "J kg-1", "");
  return enthalpy;

def createBasalMeltRateVec(grid,name='bmelt',ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for basal melt rate data.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  
  bmr = PISM.IceModelVec2S()
  bmr.create(grid, "bmelt", ghost_type, stencil_width);
  bmr.set_attrs("model_state",
                 "ice basal melt rate in ice thickness per time",
                 "m s-1", "land_ice_basal_melt_rate")
  bmr.set_glaciological_units("m year-1")
  bmr.write_in_glaciological_units = True;
  bmr.metadata().set_string("comment", "positive basal melt rate corresponds to ice loss");
  return bmr

def createTillPhiVec(grid,name='tillphi',desc="friction angle for till under grounded ice sheet",
                       ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for till friction angle data.

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.
    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
    :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                          If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                          width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  

  tillphi = PISM.IceModelVec2S()
  tillphi.create(grid, name, ghost_type, stencil_width)
  # // ghosted to allow the "redundant" computation of tauc
  # // PROPOSED standard_name = land_ice_basal_material_friction_angle
  tillphi.set_attrs("climate_steady", desc,
                    "degrees", "")
  tillphi.time_independent = True
  return tillphi

def createBasalWaterVec(grid,name='bwat',desc="effective thickness of subglacial melt water",
                          ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for basal melt water thickness data.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""
                          
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  

  bwat = PISM.IceModelVec2S()
  bwat.create(grid, name, ghost_type, stencil_width)
  bwat.set_attrs("model_state", desc,
                 "m", "")
  #// NB! Effective thickness of subglacial melt water *does* vary from 0 to hmelt_max meters only.
  bwat.metadata().set_double("valid_min", 0.0)
  valid_max = PISM.Context().config.get("bwat_max")
  bwat.metadata().set_double("valid_max", valid_max )
  return bwat

def create2dVelocityVec(grid,name="",desc="",intent="",ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2V` with attributes setup for horizontal velocity data.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  
  if name is None: name="" # FIXME

  vel = PISM.IceModelVec2V();
  vel.create(grid, name, ghost_type, stencil_width)    
  vel.set_attrs(intent, "%s%s" %("X-component of the ",desc), "m s-1", "", 0);
  vel.set_attrs(intent, "%s%s" %("Y-component of the ",desc), "m s-1", "", 1);
  vel.set_glaciological_units("m year-1");
  vel.write_in_glaciological_units = True
  huge_vel = grid.convert(1e10, "m/year", "m/second")
  attrs = [ ("valid_min", -huge_vel), ("valid_max", huge_vel), ("_FillValue", 2*huge_vel) ]
  for a in attrs: 
    for component in [0,1]:
      vel.metadata(component).set_double(a[0],a[1])
  vel.set(2*huge_vel)
  return vel

def createDrivingStressXVec(grid,name="ssa_driving_stress_x",desc="driving stress",intent="model_state",
                           ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for driving stresses int the X-direction.

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.
    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
    :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                          If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                          width is used."""

  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  

  stress = PISM.IceModelVec2S();
  stress.create(grid, name, ghost_type, stencil_width)
  stress.set_attrs(intent, "%s%s" %("X-component of the ",desc), "Pa", "");
  return stress

def createDrivingStressYVec(grid,name="ssa_driving_stress_y",desc="driving stress",intent="model_state",
                           ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for driving stresses int the Y-direction.

     :param grid: The grid to associate with the vector.
     :param name: The intrinsic name to give the vector.
     :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
     :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                           If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                           width is used."""
                           
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  

  stress = PISM.IceModelVec2S();
  stress.create(grid, name, ghost_type, stencil_width)
  stress.set_attrs(intent, "%s%s" %("X-component of the ",desc), "Pa", "");
  return stress

def createVelocityMisfitWeightVec(grid,name="vel_misfit_weight",ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for weights for inversion velocity misfit functionals.

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.
    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
    :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                          If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                          width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0

  vel_misfit_weight = PISM.IceModelVec2S();
  vel_misfit_weight.create(grid, name, ghost_type, stencil_width)
  vel_misfit_weight.set_attrs("diagnostic", "weight for surface velocity misfit functional", "", "");
  return vel_misfit_weight

def createCBarVec(grid,name="cbar",ghost_type=PISM.WITHOUT_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2S` with attributes setup for horizontal velocity magnitudes.

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.
    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
    :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                          If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                          width is used."""

  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  
  cbar = PISM.IceModelVec2S()

  cbar.create(grid, name, ghost_type, stencil_width);
  cbar.set_attrs("diagnostic",
                 "magnitude of vertically-integrated horizontal velocity of ice",
                 "m s-1", "")
  cbar.set_glaciological_units("m year-1")
  cbar.metadata().set_double("valid_min", 0.0);
  cbar.write_in_glaciological_units = True
  return cbar


def createIceMaskVec(grid,name='mask',ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2Int` with attributes setup for PISM's mask determining ice
  mode at grid points (grounded vs. floating and icy vs. ice-free).

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.
    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
    :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                          If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                          width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  

  ice_mask = PISM.IceModelVec2Int()
  ice_mask.create(grid, name, ghost_type, stencil_width);
  ice_mask.set_attrs("model_state", "grounded_dragging_floating integer mask", "", "");
  mask_values=[PISM.MASK_ICE_FREE_BEDROCK, PISM.MASK_GROUNDED, PISM.MASK_FLOATING,
               PISM.MASK_ICE_FREE_OCEAN]
  ice_mask.metadata().set_doubles("flag_values", mask_values);
  ice_mask.metadata().set_string("flag_meanings","ice_free_bedrock dragging_sheet floating ice_free_ocean");
  return ice_mask

def createBCMaskVec(grid,name='bc_mask',ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2Int` with attributes setup for masks indicating where Dirichlet data
  should be applied.

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.
    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
    :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                          If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                          width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  

  bc_mask = PISM.IceModelVec2Int()
  bc_mask.create(grid, name, ghost_type, stencil_width);
  bc_mask.set_attrs("model_state", "grounded_dragging_floating integer mask", "", "");
  mask_values=[0,1]
  bc_mask.metadata().set_doubles("flag_values", mask_values);
  bc_mask.metadata().set_string("flag_meanings","no_data ssa_dirichlet_bc_location");

  return bc_mask


def createNoModelMaskVec(grid,name='no_model_mask',ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2Int` with attributes setup for regional model flags indicating 
  grid points outside the primary domain, (i.e. where ice does does not evolve in time).

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.
    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
    :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                          If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                          width is used."""


  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  

  no_model_mask = PISM.IceModelVec2Int()
  no_model_mask.create(grid, name, ghost_type, stencil_width);
  no_model_mask.set_attrs("model_state",
    "mask: zeros (modeling domain) and ones (no-model buffer near grid edges)",
    "", ""); 

  mask_values=[0,1];
  no_model_mask.metadata().set_doubles("flag_values", mask_values)
  no_model_mask.metadata().set_string("flag_meanings", "normal special_treatment")
  no_model_mask.time_independent = True
  no_model_mask.set(0)
  return no_model_mask

def createZetaFixedMaskVec(grid,name='zeta_fixed_mask',ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2s` with attributes for the mask determining where parameterized
design variables have fixed, known values.

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.
    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
    :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                          If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                          width is used."""

  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  

  zeta_fixed_mask = PISM.IceModelVec2Int()
  zeta_fixed_mask.create(grid, name, ghost_type, stencil_width);
  zeta_fixed_mask.set_attrs("model_state", "tauc_unchanging integer mask", "", "");
  mask_values=[0,1]
  zeta_fixed_mask.metadata().set_doubles("flag_values", mask_values);
  zeta_fixed_mask.metadata().set_string("flag_meanings","tauc_changable tauc_unchangeable");

  return zeta_fixed_mask

def createLongitudeVec(grid,name="lon",ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2s` with attributes setup for longitude data.

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.
    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
    :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                          If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                          width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  
  longitude = PISM.IceModelVec2S()
  longitude.create(grid, name, ghost_type, stencil_width)
  longitude.set_attrs("mapping", "longitude", "degree_east", "longitude")
  longitude.time_independent = True
  longitude.metadata().set_string("coordinates", "")
  longitude.metadata().set_string("grid_mapping", "")
  return longitude

def createLatitudeVec(grid,name="lat",ghost_type=PISM.WITH_GHOSTS,stencil_width=None):
  """Returns an :cpp:class:`IceModelVec2s` with attributes setup for latitude data.

    :param grid: The grid to associate with the vector.
    :param name: The intrinsic name to give the vector.
    :param ghost_type: One of ``PISM.WITH_GHOSTS`` or ``PISM.WITHOUT_GHOSTS`` indicating if the vector is ghosted.
    :param stencil_width: The size of the ghost stencil. Ignored if ``ghost_type`` is ``PISM.WITHOUT_GHOSTS``.
                          If ``stencil_width`` is ``None`` and the vector is ghosted, the grid's maximum stencil
                          width is used."""
  if ghost_type==PISM.WITH_GHOSTS:    
    if stencil_width is None:
      stencil_width = grid.max_stencil_width
  else:
    stencil_width = 0  
  latitude = PISM.IceModelVec2S()
  latitude.create(grid, name, ghost_type, stencil_width)
  latitude.set_attrs("mapping", "latitude", "degree_east", "latitude")
  latitude.time_independent = True
  latitude.metadata().set_string("coordinates", "")
  latitude.metadata().set_string("grid_mapping", "")
  return latitude