import numpy as np

from checkfield import checkfield
from fielddisplay import fielddisplay
from project3d import project3d
from WriteData import WriteData


class SMBgemb(object):
    """
    SMBgemb Class definition

       Usage:
          SMB = SMBgemb()
    """

    def __init__(self, *args):  # {{{
        #each one of these properties is a transient forcing to the GEMB model, loaded from meteorological data derived
        #from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number
        #of time steps. )

        #solution choices
        #check these:
        #self.isgraingrowth = 0
        #self.issmbgradients = 0
        #isalbedo
        #isshortwave
        #isthermal
        #isaccumulation
        #ismelt
        #isdensification
        #isturbulentflux

        #inputs:
        self.Ta                     = np.nan    # 2 m air temperature, in Kelvin
        self.V                      = np.nan    # wind speed (m/s-1)
        self.dswrf                  = np.nan    # downward shortwave radiation flux [W/m^2]
        self.dlwrf                  = np.nan    # downward longwave radiation flux [W/m^2]
        self.P                      = np.nan    # precipitation [mm w.e. / m^2]
        self.eAir                   = np.nan    # screen level vapor pressure [Pa]
        self.pAir                   = np.nan    # surface pressure [Pa]
        self.Tmean                  = np.nan    # mean annual temperature [K]
        self.Vmean                  = np.nan    # mean annual wind velocity [m s-1]
        self.C                      = np.nan    # mean annual snow accumulation [kg m-2 yr-1]
        self.Tz                     = np.nan    # height above ground at which temperature (T) was sampled [m]
        self.Vz                     = np.nan    # height above ground at which wind (V) was sampled [m]

        #optional inputs:
        self.aValue                 = np.nan    # Albedo forcing at every element.  Used only if aIdx == 0.
        self.teValue                = np.nan    # Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)

        # Initialization of snow properties
        self.Dzini                  = np.nan    # cell depth (m)
        self.Dini                   = np.nan    # snow density (kg m-3)
        self.Reini                  = np.nan    # effective grain size (mm)
        self.Gdnini                 = np.nan    # grain dricity (0-1)
        self.Gspini                 = np.nan    # grain sphericity (0-1)
        self.ECini                  = np.nan    # evaporation/condensation (kg m-2)
        self.Wini                   = np.nan    # Water content (kg m-2)
        self.Aini                   = np.nan    # albedo (0-1)
        self.Tini                   = np.nan    # snow temperature (K)
        self.Sizeini                = np.nan    # Number of layers

        #settings:
        self.aIdx                   = np.nan    # method for calculating albedo and subsurface absorption (default is 1)
            # 0: direct input from aValue parameter
            # 1: effective grain radius [Gardner & Sharp, 2009]
            # 2: effective grain radius [Brun et al., 1992; LeFebre et al., 2003], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992)
            # 3: density and cloud amount [Greuell & Konzelmann, 1994]
            # 4: exponential time decay & wetness [Bougamont & Bamber, 2005]

        self.swIdx                  = np.nan    # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1, with snow density (taken from Bassford, 2004))
        self.denIdx                 = np.nan    # densification model to use (default is 2):
            # 1 = emperical model of Herron and Langway (1980)
            # 2 = semi-emperical model of Anthern et al. (2010)
            # 3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010)
            # 4 = DO NOT USE: emperical model of Li and Zwally (2004)
            # 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)
            # 6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)
            # 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)

        self.dsnowIdx               = np.nan    # model for fresh snow accumulation density (default is 1):
            # 0 = Original GEMB value, 150 kg/m^3
            # 1 = Antarctica value of fresh snow density, 350 kg/m^3
            # 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)
            # 3 = Antarctica model of Kaspers et al. (2004)
            # 4 = Greenland model of Kuipers Munneke et al. (2015)

        self.zTop                   = np.nan    # depth over which grid length is constant at the top of the snopack (default 10) [m]
        self.dzTop                  = np.nan    # initial top vertical grid spacing (default .05) [m]
        self.dzMin                  = np.nan    # initial min vertical allowable grid spacing (default dzMin/2) [m]
        self.zY                     = np.nan    # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]
        self.zMax                   = np.nan    # initial max model depth (default is min(thickness, 250)) [m]
        self.zMin                   = np.nan    # initial min model depth (default is min(thickness, 130)) [m]
        self.outputFreq             = np.nan    # output frequency in days (default is monthly, 30)

        #specific albedo parameters:
        #Method 1 and 2:
        self.aSnow                  = np.nan    # new snow albedo (0.64 - 0.89)
        self.aIce                   = np.nan    # range 0.27-0.58 for old snow
        #Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo
        self.cldFrac                = np.nan    # average cloud amount
        #Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005)
        self.t0wet                  = np.nan    # time scale for wet snow (15-21.9)
        self.t0dry                  = np.nan    # warm snow timescale (30)
        self.K                      = np.nan    # time scale temperature coef. (7)
        self.adThresh               = np.nan    # Apply aIdx method to all areas with densities below this value,
        # or else apply direct input value from aValue, allowing albedo to be altered.
        # Default value is rho water (1023 kg m-3).

        #densities:
        self.InitDensityScaling     = np.nan    # initial scaling factor multiplying the density of ice, which describes the density of the snowpack.

        #thermo:
        self.ThermoDeltaTScaling    = np.nan    # scaling factor to multiply the thermal diffusion timestep (delta t)

        self.steps_per_step         = 1
        self.averaging              = 0
        self.requested_outputs      = []

        #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM.
        #dateN: that's the last row of the above fields.
        #dt:    included in dateN. Not an input.
        #elev:  this is taken from the ISSM surface itself.

        nargin = len(args)
        if nargin:
            if nargin == 2:
                mesh = args[0]
                geometry = args[1]
                self.setdefaultparameters(mesh, geometry)
            else:
                raise Exception('constructor not supported: need mesh and geometry to set defaults')
        #}}}

    def __repr__(self):  # {{{
        #string = "   surface forcings parameters:"
        #string = "#s\n#s"%(string, fielddisplay(self, 'mass_balance', 'surface mass balance [m/yr ice eq]'))
        #string = "#s\n#s"%(string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
        string = '   surface forcings for SMB GEMB model :'
        string = "%s\n%s" % (string, fielddisplay(self, 'issmbgradients', 'is smb gradients method activated (0 or 1, default is 0)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'isgraingrowth', 'run grain growth module (default true)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'isalbedo', 'run albedo module (default true)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'isshortwave', 'run short wave module (default true)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'isthermal', 'run thermal module (default true)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'isaccumulation', 'run accumulation module (default true)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'ismelt', 'run melting  module (default true)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'isdensification', 'run densification module (default true)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'isturbulentflux', 'run turbulant heat fluxes module (default true)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'isconstrainsurfaceT', 'constrain surface temperatures to air temperature, turn off EC and surface flux contribution to surface temperature change'))
        string = "%s\n%s" % (string, fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Ta', '2 m air temperature, in Kelvin'))
        string = "%s\n%s" % (string, fielddisplay(self, 'V', 'wind speed (m s-1)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'dlwrf', 'downward shortwave radiation flux [W/m^2]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'dswrf', 'downward longwave radiation flux [W/m^2]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'P', 'precipitation [mm w.e. / m^2]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'eAir', 'screen level vapor pressure [Pa]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'pAir', 'surface pressure [Pa]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Tmean', 'mean annual temperature [K]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'C', 'mean annual snow accumulation [kg m-2 yr-1]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Vmean', 'mean annual temperature [m s-1] (default 10 m/s)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Tz', 'height above ground at which temperature (T) was sampled [m]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Vz', 'height above ground at which wind (V) eas sampled [m]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'zTop', 'depth over which grid length is constant at the top of the snopack (default 10) [m]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'dzTop', 'initial top vertical grid spacing (default .05) [m] '))
        string = "%s\n%s" % (string, fielddisplay(self, 'dzMin', 'initial min vertical allowable grid spacing (default dzMin/2) [m] '))
        string = "%s\n%s" % (string, fielddisplay(self, 'zMax', 'initial max model depth (default is min(thickness, 500)) [m]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'zMin', 'initial min model depth (default is min(thickness, 30)) [m]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'zY', 'strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'InitDensityScaling', ['initial scaling factor multiplying the density of ice', 'which describes the density of the snowpack.']))
        string = "%s\n%s" % (string, fielddisplay(self, 'ThermoDeltaTScaling', 'scaling factor to multiply the thermal diffusion timestep (delta t)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'outputFreq', 'output frequency in days (default is monthly, 30)'))
        string = "%s\n%s" % (string, fielddisplay(self, 'adThresh', 'Apply aIdx method to all areas with densities below this value, ', 'or else apply direct input value from aValue, allowing albedo to be altered.'))
        string = "%s\n%s" % (string, fielddisplay(self, 'aIdx', ['method for calculating albedo and subsurface absorption (default is 1)',
                                                                 '0: direct input from aValue parameter',
                                                                 '1: effective grain radius [Gardner & Sharp, 2009]',
                                                                 '2: effective grain radius [Brun et al., 1992; LeFebre et al., 2003], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992)',
                                                                 '3: density and cloud amount [Greuell & Konzelmann, 1994]',
                                                                 '4: exponential time decay & wetness [Bougamont & Bamber, 2005]']))
        string = "%s\n%s" % (string, fielddisplay(self, 'teValue', 'Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)'))
        #snow properties init
        string = "%s\n%s" % (string, fielddisplay(self, 'Dzini', 'Initial cell depth when restart [m]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Dini', 'Initial snow density when restart [kg m-3]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Reini', 'Initial grain size when restart [mm]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Gdnini', 'Initial grain dricity when restart [-]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Gspini', 'Initial grain sphericity when restart [-]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'ECini', 'Initial evaporation/condensation when restart [kg m-2]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Wini', 'Initial snow water content when restart [kg m-2]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Aini', 'Initial albedo when restart [-]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Tini', 'Initial snow temperature when restart [K]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [-]'))

        #additional albedo parameters:
        if isinstance(self.aIdx, (list, type(np.array([1, 2])))) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)):
            string = "%s\n%s" % (string, fielddisplay(self, 'aSnow', 'new snow albedo (0.64 - 0.89)'))
            string = "%s\n%s" % (string, fielddisplay(self, 'aIce', 'albedo of ice (0.27-0.58)'))
        elif self.aIdx == 0:
            string = "%s\n%s" % (string, fielddisplay(self, 'aValue', 'Albedo forcing at every element.  Used only if aIdx == {0, 5}'))
        elif self.aIdx == 5:
            string = "%s\n%s" % (string, fielddisplay(self, 'aValue', 'Albedo forcing at every element.  Used only if aIdx == {0, 5}'))
        elif self.aIdx == 3:
            string = "%s\n%s" % (string, fielddisplay(self, 'cldFrac', 'average cloud amount'))
        elif self.aIdx == 4:
            string = "%s\n%s" % (string, fielddisplay(self, 't0wet', 'time scale for wet snow (15-21.9) [d]'))
            string = "%s\n%s" % (string, fielddisplay(self, 't0dry', 'warm snow timescale (30) [d]'))
            string = "%s\n%s" % (string, fielddisplay(self, 'K', 'time scale temperature coef. (7) [d]'))

        string = "%s\n%s" % (string, fielddisplay(self, 'swIdx', 'apply all SW to top grid cell (0) or allow SW to penetrate surface (1) [default 1, with snow density (taken from Bassford, 2004)]'))
        string = "%s\n%s" % (string, fielddisplay(self, 'denIdx', ['densification model to use (default is 2):',
                                                                   '1 = emperical model of Herron and Langway (1980)',
                                                                   '2 = semi-emperical model of Anthern et al. (2010)',
                                                                   '3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010)',
                                                                   '4 = DO NOT USE: emperical model of Li and Zwally (2004)',
                                                                   '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)',
                                                                   '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)',
                                                                   '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)']))
        string = "%s\n%s" % (string, fielddisplay(self, 'dsnowIdx', ['model for fresh snow accumulation density (default is 1):',
                                                                     '0 = Original GEMB value, 150 kg/m^3',
                                                                     '1 = Antarctica value of fresh snow density, 350 kg/m^3',
                                                                     '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)',
                                                                     '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately',
                                                                     '4 = Greenland model of Kuipers Munneke et al. (2015)']))
        string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
        string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
        string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)')
        string = "%s\n\t\t%s" % (string, '1: Geometric')
        string = "%s\n\t\t%s" % (string, '2: Harmonic')
        string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
        return string
    #}}}

    def extrude(self, md):  # {{{
        self.Ta = project3d(md, 'vector', self.Ta, 'type', 'node')
        self.V = project3d(md, 'vector', self.V, 'type', 'node')
        self.dswrf = project3d(md, 'vector', self.dswrf, 'type', 'node')
        self.dswrf = project3d(md, 'vector', self.dswrf, 'type', 'node')
        self.P = project3d(md, 'vector', self.P, 'type', 'node')
        self.eAir = project3d(md, 'vector', self.eAir, 'type', 'node')
        self.pAir = project3d(md, 'vector', self.pAir, 'type', 'node')

        if (self.aIdx == 0) and np.isnan(self.aValue):
            self.aValue = project3d(md, 'vector', self.aValue, 'type', 'node')
        if np.isnan(self.teValue):
            self.teValue = project3d(md, 'vector', self.teValue, 'type', 'node')

        return self
    #}}}

    def defaultoutputs(self, md):  # {{{
        return ['SmbMassBalance']
    #}}}

    def setdefaultparameters(self, mesh, geometry):  # {{{
        self.isgraingrowth = 1
        self.issmbgradients = 0
        self.isalbedo = 1
        self.isshortwave = 1
        self.isthermal = 1
        self.isaccumulation = 1
        self.ismelt = 1
        self.isdensification = 1
        self.isturbulentflux = 1
        self.isclimatology = 0
        self.isconstrainsurfaceT = 0

        self.aIdx = 1
        self.swIdx = 1
        self.denIdx = 2
        self.dsnowIdx = 1
        self.zTop = 10 * np.ones((mesh.numberofelements,))
        self.dzTop = .05 * np.ones((mesh.numberofelements,))
        self.dzMin = self.dzTop / 2
        self.InitDensityScaling = 1.0
        self.ThermoDeltaTScaling = 1 / 11.0

        self.Vmean = 10 * np.ones((mesh.numberofelements,))

        self.zMax = 250 * np.ones((mesh.numberofelements,))
        self.zMin = 130 * np.ones((mesh.numberofelements,))
        self.zY = 1.025 * np.ones((mesh.numberofelements,))
        self.outputFreq = 30

        #additional albedo parameters
        self.aSnow = 0.85
        self.aIce = 0.48
        self.cldFrac = 0.1
        self.t0wet = 15
        self.t0dry = 30
        self.K = 7
        self.adThresh = 1023

        self.teValue = np.ones((mesh.numberofelements,))
        self.aValue = self.aSnow * np.ones(mesh.numberofelements,)

        self.Dzini = 0.05 * np.ones((mesh.numberofelements, 2))
        self.Dini = 910.0 * np.ones((mesh.numberofelements, 2))
        self.Reini = 2.5 * np.ones((mesh.numberofelements, 2))
        self.Gdnini = 0.0 * np.ones((mesh.numberofelements, 2))
        self.Gspini = 0.0 * np.ones((mesh.numberofelements, 2))
        self.ECini = 0.0 * np.ones((mesh.numberofelements,))
        self.Wini = 0.0 * np.ones((mesh.numberofelements, 2))
        self.Aini = self.aSnow * np.ones((mesh.numberofelements, 2))
        self.Tini = 273.15 * np.ones((mesh.numberofelements, 2))
#       /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh).
#       If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp
        self.Sizeini = 2 * np.ones((mesh.numberofelements,))
    #}}}

    def checkconsistency(self, md, solution, analyses):    # {{{

        md = checkfield(md, 'fieldname', 'smb.isgraingrowth', 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.issmbgradients', 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.isalbedo', 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.isshortwave', 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.isthermal', 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.isaccumulation', 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.ismelt', 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.isdensification', 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.isturbulentflux', 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.isclimatology', 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.isconstrainsurfaceT', 'values', [0, 1])

        md = checkfield(md, 'fieldname', 'smb.Ta', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100)  #-100/100 celsius min/max value
        md = checkfield(md, 'fieldname', 'smb.V', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '<', 45, 'size', np.shape(self.Ta))  #max 500 km/h
        md = checkfield(md, 'fieldname', 'smb.dswrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1400, 'size', np.shape(self.Ta))
        md = checkfield(md, 'fieldname', 'smb.dlwrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, 'size', np.shape(self.Ta))
        md = checkfield(md, 'fieldname', 'smb.P', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 200, 'size', np.shape(self.Ta))
        md = checkfield(md, 'fieldname', 'smb.eAir', 'timeseries', 1, 'NaN', 1, 'Inf', 1, 'size', np.shape(self.Ta))

        if (self.isclimatology > 0):
            md = checkfield(md, 'fieldname', 'smb.Ta', 'size', [md.mesh.numberofelements + 1], 'message', 'Ta must have md.mesh.numberofelements+1 rows in order to force a climatology')
            md = checkfield(md, 'fieldname', 'smb.V', 'size', [md.mesh.numberofelements + 1], 'message', 'V must have md.mesh.numberofelements+1 rows in order to force a climatology')
            md = checkfield(md, 'fieldname', 'smb.dswrf', 'size', [md.mesh.numberofelements + 1], 'message', 'dswrf must have md.mesh.numberofelements+1 rows in order to force a climatology')
            md = checkfield(md, 'fieldname', 'smb.dlwrf', 'size', [md.mesh.numberofelements + 1], 'message', 'dlwrf must have md.mesh.numberofelements+1 rows in order to force a climatology')
            md = checkfield(md, 'fieldname', 'smb.P', 'size', [md.mesh.numberofelements + 1], 'message', 'P must have md.mesh.numberofelements+1 rows in order to force a climatology')
            md = checkfield(md, 'fieldname', 'smb.eAir', 'size', [md.mesh.numberofelements + 1], 'message', 'eAir must have md.mesh.numberofelements+1 rows in order to force a climatology')

        md = checkfield(md, 'fieldname', 'smb.Tmean', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100)  #-100/100 celsius min/max value
        md = checkfield(md, 'fieldname', 'smb.C', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0)
        md = checkfield(md, 'fieldname', 'smb.Vmean', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0)
        md = checkfield(md, 'fieldname', 'smb.Tz', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 5000)
        md = checkfield(md, 'fieldname', 'smb.Vz', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 5000)

        md = checkfield(md, 'fieldname', 'smb.teValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)

        md = checkfield(md, 'fieldname', 'smb.aIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1, 2, 3, 4])
        md = checkfield(md, 'fieldname', 'smb.swIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
        md = checkfield(md, 'fieldname', 'smb.denIdx', 'NaN', 1, 'Inf', 1, 'values', [1, 2, 3, 4, 5, 6, 7])
        md = checkfield(md, 'fieldname', 'smb.dsnowIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1, 2, 3, 4])

        md = checkfield(md, 'fieldname', 'smb.zTop', 'NaN', 1, 'Inf', 1, '> = ', 0)
        md = checkfield(md, 'fieldname', 'smb.dzTop', 'NaN', 1, 'Inf', 1, '>', 0)
        md = checkfield(md, 'fieldname', 'smb.dzMin', 'NaN', 1, 'Inf', 1, '>', 0)
        md = checkfield(md, 'fieldname', 'smb.zY', 'NaN', 1, 'Inf', 1, '> = ', 1)
        md = checkfield(md, 'fieldname', 'smb.outputFreq', 'NaN', 1, 'Inf', 1, '>', 0, '<', 10 * 365)  #10 years max
        md = checkfield(md, 'fieldname', 'smb.InitDensityScaling', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1)
        md = checkfield(md, 'fieldname', 'smb.ThermoDeltaTScaling', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1)
        md = checkfield(md, 'fieldname', 'smb.adThresh', 'NaN', 1, 'Inf', 1, '>=', 0)

        if isinstance(self.aIdx, (list, type(np.array([1, 2])))) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)):
            md = checkfield(md, 'fieldname', 'smb.aSnow', 'NaN', 1, 'Inf', 1, '> = ', .64, '< = ', .89)
            md = checkfield(md, 'fieldname', 'smb.aIce', 'NaN', 1, 'Inf', 1, '> = ', .27, '< = ', .58)
        elif self.aIdx == 0:
            md = checkfield(md, 'fieldname', 'smb.aValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
        elif self.aIdx == 3:
            md = checkfield(md, 'fieldname', 'smb.cldFrac', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1)
        elif self.aIdx == 4:
            md = checkfield(md, 'fieldname', 'smb.t0wet', 'NaN', 1, 'Inf', 1, '> = ', 15, '< = ', 21.9)
            md = checkfield(md, 'fieldname', 'smb.t0dry', 'NaN', 1, 'Inf', 1, '> = ', 30, '< = ', 30)
            md = checkfield(md, 'fieldname', 'smb.K', 'NaN', 1, 'Inf', 1, '> = ', 7, '< = ', 7)

        #check zTop is < local thickness:
        he = np.sum(md.geometry.thickness[md.mesh.elements - 1], axis=1) / np.size(md.mesh.elements, 1)
        if np.any(he < self.zTop):
            raise IOError('SMBgemb consistency check error: zTop should be smaller than local ice thickness')
        md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
        md = checkfield(md, 'fieldname', 'smb.requested_outputs', 'stringrow', 1)
        return md
    # }}}

    def marshall(self, prefix, md, fid):    # {{{
        yts = md.constants.yts

        WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 8, 'format', 'Integer')

        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isgraingrowth', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'issmbgradients', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isalbedo', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isshortwave', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isthermal', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isaccumulation', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ismelt', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isdensification', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isturbulentflux', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isclimatology', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isconstrainsurfaceT', 'format', 'Boolean')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Ta', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'V', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dswrf', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dlwrf', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'P', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'eAir', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'pAir', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)

        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Tmean', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'C', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Vmean', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Tz', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Vz', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'zTop', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dzTop', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dzMin', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'zY', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'zMax', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'zMin', 'format', 'DoubleMat', 'mattype', 2)

        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'aIdx', 'format', 'Integer')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'swIdx', 'format', 'Integer')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'denIdx', 'format', 'Integer')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dsnowIdx', 'format', 'Integer')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'InitDensityScaling', 'format', 'Double')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ThermoDeltaTScaling', 'format', 'Double')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'outputFreq', 'format', 'Double')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'aSnow', 'format', 'Double')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'aIce', 'format', 'Double')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'cldFrac', 'format', 'Double')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 't0wet', 'format', 'Double')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 't0dry', 'format', 'Double')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'K', 'format', 'Double')
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'adThresh', 'format', 'Double')

        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'aValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'teValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)

        #snow properties init
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Dzini', 'format', 'DoubleMat', 'mattype', 3)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Dini', 'format', 'DoubleMat', 'mattype', 3)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Reini', 'format', 'DoubleMat', 'mattype', 3)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Gdnini', 'format', 'DoubleMat', 'mattype', 3)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Gspini', 'format', 'DoubleMat', 'mattype', 3)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ECini', 'format', 'DoubleMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Wini', 'format', 'DoubleMat', 'mattype', 3)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Aini', 'format', 'DoubleMat', 'mattype', 3)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Tini', 'format', 'DoubleMat', 'mattype', 3)
        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Sizeini', 'format', 'IntMat', 'mattype', 2)
        WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
        #figure out dt from forcings:
        if (np.any(self.P[-1] - self.Ta[-1] != 0) | np.any(self.V[-1] - self.Ta[-1] != 0) | np.any(self.dswrf[-1] - self.Ta[-1] != 0) | np.any(self.dlwrf[-1] - self.Ta[-1] != 0) | np.any(self.eAir[-1] - self.Ta[-1] != 0) | np.any(self.pAir[-1] - self.Ta[-1] != 0)):
            raise IOError('All GEMB forcings (Ta, P, V, dswrf, dlwrf, eAir, pAir) must have the same time steps in the final row!')

        time = self.Ta[-1]  #assume all forcings are on the same time step
        dtime = np.diff(time, n=1, axis=0)
        dt = min(dtime)

        WriteData(fid, prefix, 'data', dt, 'name', 'md.smb.dt', 'format', 'Double', 'scale', yts)

        # Check if smb_dt goes evenly into transient core time step
        if (md.timestepping.time_step % dt >= 1e-10):
            raise IOError('smb_dt/dt = #f. The number of SMB time steps in one transient core time step has to be an an integer', md.timestepping.time_step / dt)

        #process requested outputs
        outputs = self.requested_outputs
        indices = [i for i, x in enumerate(outputs) if x == 'default']
        if len(indices) > 0:
            outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
            outputs = outputscopy

        WriteData(fid, prefix, 'data', outputs, 'name', 'md.smb.requested_outputs', 'format', 'StringArray')
    # }}}
