# Copyright (C) 2011, 2014, 2015, 2016, 2017, 2018 David Maxwell and Constantine Khroulev
#
# 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

"""Module containing classes and functions for managing computational
grids, model physics, and creating standard vectors."""

import PISM

class ModelVecs(object):

    """A dictionary of :cpp:class:`IceModelVec`'s containing model data.
    This plays a similar role as the C++ :cpp:class:`Vars` with the key
    difference that it keeps track of fields that may need to be written
    to a file.::

        # 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

    See also functions :func:`createYieldStressVec` and friends for one
    step creation of standard PISM vecs.

    """

    def __init__(self, variables):
        self.needs_writing = set()
        self._vecs = variables

    def __getattr__(self, name):
        """Allows access to contents via dot notation, e.g. `vecs.thickness`."""
        try:
            return self.get(name)
        except RuntimeError as e:
            print("#### ModelVecs.__getattr__ got a PISM-level exception:", e)
            print("# The contents of this ModelVecs instance:")
            print(self)
            print("#### End of the contents")
            raise AttributeError(name)

    def __repr__(self):
        """:returns: a string representation of the :class:`ModelVecs` listing its contents."""
        s = 'ModelVecs:\n'
        items = [(key, self._vecs.get(key)) for key in list(self._vecs.keys())]
        for (key, value) in items:
            short_name = value.metadata().get_string("short_name")
            long_name = value.metadata().get_string("long_name")
            standard_name = value.metadata().get_string("standard_name")

            if key != short_name or len(standard_name) == 0:
                # alternative name was given, so we can't look up by standard name
                s += "    %s: %s\n" % (key, long_name)
            else:
                s += "    %s or %s: %s\n" % (short_name, standard_name, long_name)
        return s

    def __str__(self):
        return self.__repr__()

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

    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:
            self._vecs.add(var)
        else:
            self._vecs.add(var, name)

        if writing:
            self.needs_writing.add(var)

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

    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(key=lambda v: v.get_name())
        for v in vlist:
            v.write(output_filename)

    def writeall(self, output_filename):
        """Writes all member vectors to the given :file:`.nc` filename."""
        vlist = [self._vecs.get(var) for var in list(self._vecs.keys())]
        vlist.sort(key=lambda v: v.get_name())

        try:
            com = vlist[0].get_grid().ctx().com()
        except:
            com = PISM.PETSc.COMM_WORLD

        f = PISM.PIO(com, "netcdf3", output_filename, PISM.PISM_READWRITE)
        for v in vlist:
            v.define(f, PISM.PISM_DOUBLE)
        for v in vlist:
            v.write(f)
        f.close()

    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)


class ModelData(object):

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

         * ``grid``      the computation grid
         * ``config``    the associated :cpp:class:`Config`
         * ``basal``     a :cpp:class:`YieldStress`.
         * ``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.ctx().config()
        self.config = config

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

        self.vecs = ModelVecs(grid.variables())         #: 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 initGrid(ctx, Lx, Ly, Lz, Mx, My, Mz, r):
    """Initialize a :cpp:class:`IceGrid` intended for 3-d computations.

    :param ctx:  The execution context.
    :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 r:    grid registration (one of ``PISM.CELL_CENTER``, ``PISM.CELL_CORNER``)
    """
    P = PISM.GridParameters(ctx.config)

    P.Lx = Lx
    P.Ly = Ly
    P.Mx = Mx
    P.My = My
    P.registration = r
    z = PISM.IceGrid.compute_vertical_levels(Lz, Mz, PISM.EQUAL)
    P.z = PISM.DoubleVector(z)
    P.horizontal_extent_from_options()
    P.ownership_ranges_from_options(ctx.size)

    return PISM.IceGrid(ctx.ctx, P)

def _stencil_width(grid, ghost_type, width):
    "Compute stencil width for a field."
    if ghost_type == PISM.WITH_GHOSTS:
        if width is None:
            return int(grid.ctx().config().get_double("grid.max_stencil_width"))
        else:
            return width
    else:
        return 0


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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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 createSeaLevelVec(grid, name='sea_level', ghost_type=PISM.WITH_GHOSTS, stencil_width=1):
    """Returns an :cpp:class:`IceModelVec2S` with attributes setup for
    sea level 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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    sea_level = PISM.IceModelVec2S(grid, name, ghost_type, stencil_width)
    sea_level.set_attrs("diagnostic",
                        "sea level elevation above reference ellipsoid",
                        "m",
                        "sea_surface_height_above_reference_ellipsoid")
    return sea_level


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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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.

    """

    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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)
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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)
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

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

    power = 1.0 / grid.ctx().config().get_double("stress_balance.ssa.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()
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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 createStrainHeatingVec(grid, name='strain_heating',
                           ghost_type=PISM.WITH_GHOSTS, stencil_width=None):
    """Returns an :cpp:class:`IceModelVec3` with attributes setup for
    strain heating.

    """

    result = PISM.IceModelVec3()
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    result.create(grid, name, ghost_type, stencil_width)
    result.set_attrs("internal",
                     "rate of strain heating in ice (dissipation heating)",
                     "W m-3", "")
    return result


def createAgeVec(grid, name='age', ghost_type=PISM.WITH_GHOSTS, stencil_width=None):
    """Returns an :cpp:class:`IceModelVec3` with attributes setup for age 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.

    """

    age = PISM.IceModelVec3()
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    age.create(grid, name, ghost_type, stencil_width)
    age.set_attrs("model_state", "age of the ice", "s", "")
    return age


def create3DVelocityVecs(grid, ghost_type=PISM.WITH_GHOSTS, stencil_width=None):
    """Returns a size 3 tuple with instances of :cpp:class:`IceModelVec3`
    with attributes setup for the x, y, and z-components of the 3D
    velocity.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    u = PISM.IceModelVec3()
    u.create(grid, "uvel", ghost_type, stencil_width)
    u.set_attrs("diagnostic", "x-component of the ice velocity", "m s-1", "")

    v = PISM.IceModelVec3()
    v.create(grid, "vvel", ghost_type, stencil_width)
    v.set_attrs("diagnostic", "y-component of the ice velocity", "m s-1", "")

    w = PISM.IceModelVec3()
    w.create(grid, "wvel_rel", ghost_type, stencil_width)
    w.set_attrs("diagnostic",
                "z-component of the ice velocity relative to the base directly below",
                "m s-1", "")

    return u, v, w


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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    bmr = PISM.IceModelVec2S()
    bmr.create(grid, name, 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.write_in_glaciological_units = True
    bmr.metadata().set_string("glaciological_units", "m year-1")
    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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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 createCellAreaVec(grid, name='cell_area', desc="cell area",
                      ghost_type=PISM.WITHOUT_GHOSTS, stencil_width=None):
    """Returns an :cpp:class:`IceModelVec2S` with cell area (dx*dy).

    :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.

    """

    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    cell_area = PISM.IceModelVec2S()
    cell_area.create(grid, name, ghost_type, stencil_width)
    cell_area.set_attrs("model_state", desc, "m2", "")
    cell_area.metadata().set_double("valid_min", 0.0)

    cell_area.set(grid.dx() * grid.dy())

    return cell_area


def createBasalWaterVec(grid, name='tillwat', 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.

    """

    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

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


def createGroundingLineMask(grid, name='gl_mask', desc="fractional floatation mask",
                            ghost_type=PISM.WITH_GHOSTS, stencil_width=None):
    """Returns an :cpp:class:`IceModelVec2S` with attributes setup for the fractional floatation mask.

    :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.

    """

    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    gl_mask = PISM.IceModelVec2S()
    gl_mask.create(grid, name, ghost_type, stencil_width)
    gl_mask.set_attrs("model_state", desc, "1", "")
    gl_mask.metadata().set_double("valid_min", 0.0)
    gl_mask.metadata().set_double("valid_max", 1.0)
    return gl_mask


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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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.metadata(0).set_string("glaciological_units", "m year-1")
    vel.metadata(1).set_string("glaciological_units", "m year-1")
    vel.write_in_glaciological_units = True
    sys = grid.ctx().unit_system()
    huge_vel = PISM.convert(sys, 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.

    """

    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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.

    """

    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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="velbar_mag", 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.

    """

    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    velbar_mag = PISM.IceModelVec2S()

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


def createIceMaskVec(grid, name='mask', ghost_type=PISM.WITH_GHOSTS, stencil_width=None):
    """Returns an :cpp:class:`IceModelVec2CellType` 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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    ice_mask = PISM.IceModelVec2CellType()
    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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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.

    """

    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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.

    """

    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)
    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.

    """
    stencil_width = _stencil_width(grid, ghost_type, stencil_width)

    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
