Ignore:
Timestamp:
10/17/19 06:03:43 (5 years ago)
Author:
bdef
Message:

Adding a substeping framework in hydro and part of smb

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/classes/SMBgemb.py

    r24213 r24240  
    3131
    3232        #inputs:
    33         self.Ta = float('NaN')  #2 m air temperature, in Kelvin
    34         self.V = float('NaN')  #wind speed (m / s - 1)
    35         self.dswrf = float('NaN')  #downward shortwave radiation flux [W / m^2]
    36         self.dlwrf = float('NaN')  #downward longwave radiation flux [W / m^2]
    37         self.P = float('NaN')  #precipitation [mm w.e. / m^2]
    38         self.eAir = float('NaN')  #screen level vapor pressure [Pa]
    39         self.pAir = float('NaN')  #surface pressure [Pa]
    40         self.Tmean = float('NaN')  #mean annual temperature [K]
    41         self.Vmean = float('NaN')  #mean annual wind velocity [m s - 1]
    42         self.C = float('NaN')  #mean annual snow accumulation [kg m - 2 yr - 1]
    43         self.Tz = float('NaN')  #height above ground at which temperature (T) was sampled [m]
    44         self.Vz = float('NaN')  #height above ground at which wind (V) was sampled [m]
     33        self.Ta = float('NaN')       #2 m air temperature, in Kelvin
     34        self.V = float('NaN')        #wind speed (m/s-1)
     35        self.dswrf = float('NaN')    #downward shortwave radiation flux [W/m^2]
     36        self.dlwrf = float('NaN')    #downward longwave radiation flux [W/m^2]
     37        self.P = float('NaN')        #precipitation [mm w.e. / m^2]
     38        self.eAir = float('NaN')     #screen level vapor pressure [Pa]
     39        self.pAir = float('NaN')     #surface pressure [Pa]
     40        self.Tmean = float('NaN')    #mean annual temperature [K]
     41        self.Vmean = float('NaN')    #mean annual wind velocity [m s-1]
     42        self.C = float('NaN')        #mean annual snow accumulation [kg m-2 yr-1]
     43        self.Tz = float('NaN')       #height above ground at which temperature (T) was sampled [m]
     44        self.Vz = float('NaN')       #height above ground at which wind (V) eas sampled [m]
    4545
    4646        #optional inputs:
    … …  
    4949
    5050        # Initialization of snow properties
    51         self.Dzini = float('NaN')  #cell depth (m)
    52         self.Dini = float('NaN')  #snow density (kg m - 3)
    53         self.Reini = float('NaN')  #effective grain size (mm)
    54         self.Gdnini = float('NaN')  #grain dricity (0 - 1)
    55         self.Gspini = float('NaN')  #grain sphericity (0 - 1)
    56         self.ECini = float('NaN')  #evaporation/condensation (kg m - 2)
    57         self.Wini = float('NaN')  #Water content (kg m - 2)
    58         self.Aini = float('NaN')  #albedo (0 - 1)
    59         self.Tini = float('NaN')  #snow temperature (K)
     51        self.Dzini = float('NaN')    #cell depth (m)
     52        self.Dini = float('NaN')     #snow density (kg m-3)
     53        self.Reini = float('NaN')    #effective grain size (mm)
     54        self.Gdnini = float('NaN')   #grain dricity (0-1)
     55        self.Gspini = float('NaN')   #grain sphericity (0-1)
     56        self.ECini = float('NaN')    #evaporation/condensation (kg m-2)
     57        self.Wini = float('NaN')     #Water content (kg m-2)
     58        self.Aini = float('NaN')     #albedo (0-1)
     59        self.Tini = float('NaN')     #snow temperature (K)
    6060        self.Sizeini = float('NaN')  #Number of layers
    6161
    6262        #settings:
    63         self.aIdx = float('NaN')  #method for calculating albedo and subsurface absorption (default is 1)
    64         # 0: direct input from aValue parameter
    65         # 1: effective grain radius [Gardner & Sharp, 2009]
    66         # 2: effective grain radius [Brun et al., 2009]
    67         # 3: density and cloud amount [Greuell & Konzelmann, 1994]
    68         # 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
    69         # 5: ingest MODIS mode, direct input from aValue parameter applied to surface ice only
    70         self.swIdx = float('NaN')  # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1)
    71         self.denIdx = float('NaN')  #densification model to use (default is 2):
    72         # 1 = emperical model of Herron and Langway (1980)
    73         # 2 = semi - emperical model of Anthern et al. (2010)
    74         # 3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010)
    75         # 4 = DO NOT USE: emperical model of Li and Zwally (2004)
    76         # 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)
    77         # 6 = Antarctica semi - emperical model of Ligtenberg et al. (2011)
    78         # 7 = Greenland semi - emperical model of Kuipers Munneke et al. (2015)
     63        self.aIdx = float('NaN')     #method for calculating albedo and subsurface absorption (default is 1)
     64        self.swIdx = float('NaN')    # apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1)
     65        self.denIdx = float('NaN')   #densification model to use (default is 2):
    7966        self.dsnowIdx = float('NaN')  #model for fresh snow accumulation density (default is 1):
    80         # 0 = Original GEMB value, 150 kg / m^3
    81         # 1 = Antarctica value of fresh snow density, 350 kg / m^3
    82         # 2 = Greenland value of fresh snow density, 315 kg / m^3, Fausto et al. (2008)
    83         # 3 = Antarctica model of Kaspers et al. (2004)
    84         # 4 = Greenland model of Kuipers Munneke et al. (2015)
    85 
    86         self.zTop = float('NaN')  # depth over which grid length is constant at the top of the snopack (default 10) [m]
    87         self.dzTop = float('NaN')  # initial top vertical grid spacing (default .05) [m]
    88         self.dzMin = float('NaN')  # initial min vertical allowable grid spacing (default dzTop / 2) [m]
    89         self.zY = float('NaN')  # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]
    90         self.zMax = float('NaN')  #initial max model depth (default is min(thickness, 250)) [m]
    91         self.zMin = float('NaN')  #initial min model depth (default is min(thickness, 130)) [m]
    92         self.outputFreq = float('NaN')  #output frequency in days (default is monthly, 30)
     67        self.zTop = float('NaN')     # depth over which grid length is constant at the top of the snopack (default 10) [m]
     68        self.dzTop = float('NaN')    # initial top vertical grid spacing (default .05) [m]
     69        self.dzMin = float('NaN')    # initial min vertical allowable grid spacing (default dzMin/2) [m]
     70        self.zY = float('NaN')       # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]
     71        self.zMax = float('NaN')     #initial max model depth (default is min(thickness, 500)) [m]
     72        self.zMin = float('NaN')     #initial min model depth (default is min(thickness, 30)) [m]
     73        self.outputFreq = float('NaN')       #output frequency in days (default is monthly, 30)
    9374
    9475        #specific albedo parameters:
    9576        #Method 1 and 2:
    96         self.aSnow = float('NaN')  # new snow albedo (0.64 - 0.89)
    97         self.aIce = float('NaN')  # range 0.27 - 0.58 for old snow
    98         #Method 3: Radiation Correction Factors - > only used for met station data and Greuell & Konzelmann, 1994 albedo
     77        self.aSnow = float('NaN')    # new snow albedo (0.64 - 0.89)
     78        self.aIce = float('NaN')     # range 0.27-0.58 for old snow
     79        #Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo
    9980        self.cldFrac = float('NaN')  # average cloud amount
    10081        #Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005)
    101         self.t0wet = float('NaN')  # time scale for wet snow (15 - 21.9)
    102         self.t0dry = float('NaN')  # warm snow timescale (30)
    103         self.K = float('NaN')  # time scale temperature coef. (7)
     82        self.t0wet = float('NaN')    # time scale for wet snow (15-21.9)
     83        self.t0dry = float('NaN')    # warm snow timescale (30)
     84        self.K = float('NaN')        # time scale temperature coef. (7)
    10485        self.adThresh = float('NaN')  # Apply aIdx method to all areas with densities below this value,
    10586        # or else apply direct input value from aValue, allowing albedo to be altered.
    106         # Default value is rho water (1023 kg m - 3).
     87        # Default value is rho water (1023 kg m-3).
    10788
    10889        #densities:
    109         self.InitDensityScaling = float('NaN')  #initial scaling factor multiplying the density of ice, which describes the density of the snowpack.
     90        self.InitDensityScaling = float('NaN')       #initial scaling factor multiplying the density of ice, which describes the density of the snowpack.
    11091
    11192        #thermo:
    11293        self.ThermoDeltaTScaling = float('NaN')  #scaling factor to multiply the thermal diffusion timestep (delta t)
    11394
     95        self.steps_per_step = 1
    11496        self.requested_outputs = []
    11597
    116     #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM.
    117     #dateN: that's the last row of the above fields.
    118     #dt:    included in dateN. Not an input.
    119     #elev:  this is taken from the ISSM surface itself.
    120 
    121     #}}}
     98        #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM.
     99        #dateN: that's the last row of the above fields.
     100        #dt:    included in dateN. Not an input.
     101        #elev:  this is taken from the ISSM surface itself.
     102
     103        #}}}
     104
    122105
    123106    def __repr__(self):  # {{{
    124107        #string = "   surface forcings parameters:"
    125         #string = "  #s\n  #s" % (string, fielddisplay(self, 'mass_balance', 'surface mass balance [m / yr ice eq]'))
    126         #string = "  #s\n  #s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     108        #string = "#s\n#s"%(string, fielddisplay(self, 'mass_balance', 'surface mass balance [m/yr ice eq]'))
     109        #string = "#s\n#s"%(string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    127110        string = '   surface forcings for SMB GEMB model :'
    128111        string = "%s\n%s" % (string, fielddisplay(self, 'issmbgradients', 'is smb gradients method activated (0 or 1, default is 0)'))
    … …  
    137120        string = "%s\n%s" % (string, fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)'))
    138121        string = "%s\n%s" % (string, fielddisplay(self, 'Ta', '2 m air temperature, in Kelvin'))
    139         string = "%s\n%s" % (string, fielddisplay(self, 'V', 'wind speed (m s - 1)'))
    140         string = "%s\n%s" % (string, fielddisplay(self, 'dlwrf', 'downward shortwave radiation flux [W / m^2]'))
    141         string = "%s\n%s" % (string, fielddisplay(self, 'dswrf', 'downward longwave radiation flux [W / m^2]'))
     122        string = "%s\n%s" % (string, fielddisplay(self, 'V', 'wind speed (m s-1)'))
     123        string = "%s\n%s" % (string, fielddisplay(self, 'dlwrf', 'downward shortwave radiation flux [W/m^2]'))
     124        string = "%s\n%s" % (string, fielddisplay(self, 'dswrf', 'downward longwave radiation flux [W/m^2]'))
    142125        string = "%s\n%s" % (string, fielddisplay(self, 'P', 'precipitation [mm w.e. / m^2]'))
    143126        string = "%s\n%s" % (string, fielddisplay(self, 'eAir', 'screen level vapor pressure [Pa]'))
    144127        string = "%s\n%s" % (string, fielddisplay(self, 'pAir', 'surface pressure [Pa]'))
    145128        string = "%s\n%s" % (string, fielddisplay(self, 'Tmean', 'mean annual temperature [K]'))
    146         string = "%s\n%s" % (string, fielddisplay(self, 'C', 'mean annual snow accumulation [kg m - 2 yr - 1]'))
    147         string = "%s\n%s" % (string, fielddisplay(self, 'Vmean', 'mean annual wind velocity [m s - 1] (default 10 m / s)'))
     129        string = "%s\n%s" % (string, fielddisplay(self, 'C', 'mean annual snow accumulation [kg m-2 yr-1]'))
     130        string = "%s\n%s" % (string, fielddisplay(self, 'Vmean', 'mean annual temperature [m s-1] (default 10 m/s)'))
    148131        string = "%s\n%s" % (string, fielddisplay(self, 'Tz', 'height above ground at which temperature (T) was sampled [m]'))
    149         string = "%s\n%s" % (string, fielddisplay(self, 'Vz', 'height above ground at which wind (V) was sampled [m]'))
     132        string = "%s\n%s" % (string, fielddisplay(self, 'Vz', 'height above ground at which wind (V) eas sampled [m]'))
    150133        string = "%s\n%s" % (string, fielddisplay(self, 'zTop', 'depth over which grid length is constant at the top of the snopack (default 10) [m]'))
    151134        string = "%s\n%s" % (string, fielddisplay(self, 'dzTop', 'initial top vertical grid spacing (default .05) [m] '))
    152         string = "%s\n%s" % (string, fielddisplay(self, 'dzMin', 'initial min vertical allowable grid spacing (default dzMin / 2) [m] '))
     135        string = "%s\n%s" % (string, fielddisplay(self, 'dzMin', 'initial min vertical allowable grid spacing (default dzMin/2) [m] '))
    153136        string = "%s\n%s" % (string, fielddisplay(self, 'zMax', 'initial max model depth (default is min(thickness, 500)) [m]'))
    154137        string = "%s\n%s" % (string, fielddisplay(self, 'zMin', 'initial min model depth (default is min(thickness, 30)) [m]'))
    … …  
    167150        #snow properties init
    168151        string = "%s\n%s" % (string, fielddisplay(self, 'Dzini', 'Initial cell depth when restart [m]'))
    169         string = "%s\n%s" % (string, fielddisplay(self, 'Dini', 'Initial snow density when restart [kg m - 3]'))
     152        string = "%s\n%s" % (string, fielddisplay(self, 'Dini', 'Initial snow density when restart [kg m-3]'))
    170153        string = "%s\n%s" % (string, fielddisplay(self, 'Reini', 'Initial grain size when restart [mm]'))
    171         string = "%s\n%s" % (string, fielddisplay(self, 'Gdnini', 'Initial grain dricity when restart [ - ]'))
    172         string = "%s\n%s" % (string, fielddisplay(self, 'Gspini', 'Initial grain sphericity when restart [ - ]'))
    173         string = "%s\n%s" % (string, fielddisplay(self, 'ECini', 'Initial evaporation / condensation when restart [kg m - 2]'))
    174         string = "%s\n%s" % (string, fielddisplay(self, 'Wini', 'Initial snow water content when restart [kg m - 2]'))
    175         string = "%s\n%s" % (string, fielddisplay(self, 'Aini', 'Initial albedo when restart [ - ]'))
     154        string = "%s\n%s" % (string, fielddisplay(self, 'Gdnini', 'Initial grain dricity when restart [-]'))
     155        string = "%s\n%s" % (string, fielddisplay(self, 'Gspini', 'Initial grain sphericity when restart [-]'))
     156        string = "%s\n%s" % (string, fielddisplay(self, 'ECini', 'Initial evaporation/condensation when restart [kg m-2]'))
     157        string = "%s\n%s" % (string, fielddisplay(self, 'Wini', 'Initial snow water content when restart [kg m-2]'))
     158        string = "%s\n%s" % (string, fielddisplay(self, 'Aini', 'Initial albedo when restart [-]'))
    176159        string = "%s\n%s" % (string, fielddisplay(self, 'Tini', 'Initial snow temperature when restart [K]'))
    177160        string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [K]'))
    178161
    179162        #additional albedo parameters:
    180         if isinstance(self.aIdx, list) or isinstance(self.aIdx, np.ndarray) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)):
     163        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)):
    181164            string = "%s\n%s" % (string, fielddisplay(self, 'aSnow', 'new snow albedo (0.64 - 0.89)'))
    182             string = "%s\n%s" % (string, fielddisplay(self, 'aIce', 'albedo of ice (0.27 - 0.58)'))
    183         elif elf.aIdx == 0:
     165            string = "%s\n%s" % (string, fielddisplay(self, 'aIce', 'albedo of ice (0.27-0.58)'))
     166        elif self.aIdx == 0:
    184167            string = "%s\n%s" % (string, fielddisplay(self, 'aValue', 'Albedo forcing at every element.  Used only if aIdx == {0, 5}'))
    185         elif elf.aIdx == 5:
     168        elif self.aIdx == 5:
    186169            string = "%s\n%s" % (string, fielddisplay(self, 'aValue', 'Albedo forcing at every element.  Used only if aIdx == {0, 5}'))
    187170        elif self.aIdx == 3:
    188171            string = "%s\n%s" % (string, fielddisplay(self, 'cldFrac', 'average cloud amount'))
    189172        elif self.aIdx == 4:
    190             string = "%s\n%s" % (string, fielddisplay(self, 't0wet', 'time scale for wet snow (15 - 21.9) [d]'))
     173            string = "%s\n%s" % (string, fielddisplay(self, 't0wet', 'time scale for wet snow (15-21.9) [d]'))
    191174            string = "%s\n%s" % (string, fielddisplay(self, 't0dry', 'warm snow timescale (30) [d]'))
    192175            string = "%s\n%s" % (string, fielddisplay(self, 'K', 'time scale temperature coef. (7) [d]'))
    … …  
    195178        string = "%s\n%s" % (string, fielddisplay(self, 'denIdx', ['densification model to use (default is 2):',
    196179                                                                   '1 = emperical model of Herron and Langway (1980)',
    197                                                                    '2 = semi - emperical model of Anthern et al. (2010)',
     180                                                                   '2 = semi-emperical model of Anthern et al. (2010)',
    198181                                                                   '3 = DO NOT USE: physical model from Appix B of Anthern et al. (2010)',
    199182                                                                   '4 = DO NOT USE: emperical model of Li and Zwally (2004)',
    200183                                                                   '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)',
    201                                                                    '6 = Antarctica semi - emperical model of Ligtenberg et al. (2011)',
    202                                                                    '7 = Greenland semi - emperical model of Kuipers Munneke et al. (2015)']))
     184                                                                   '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)',
     185                                                                   '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)']))
    203186        string = "%s\n%s" % (string, fielddisplay(self, 'dsnowIdx', ['model for fresh snow accumulation density (default is 1):',
    204                                                                      '0 = Original GEMB value, 150 kg / m^3',
    205                                                                      '1 = Antarctica value of fresh snow density, 350 kg / m^3',
    206                                                                      '2 = Greenland value of fresh snow density, 315 kg / m^3, Fausto et al. (2008)',
     187                                                                     '0 = Original GEMB value, 150 kg/m^3',
     188                                                                     '1 = Antarctica value of fresh snow density, 350 kg/m^3',
     189                                                                     '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)',
    207190                                                                     '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately',
    208191                                                                     '4 = Greenland model of Kuipers Munneke et al. (2015)']))
     192        string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
    209193        string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    210194        return string
    … …  
    247231        self.denIdx = 2
    248232        self.dsnowIdx = 1
    249         self.zTop = 10 * np.ones((mesh.numberofelements, ))
    250         self.dzTop = .05 * np.ones((mesh.numberofelements, ))
     233        self.zTop = 10 * np.ones((mesh.numberofelements,))
     234        self.dzTop = .05 * np.ones((mesh.numberofelements,))
    251235        self.dzMin = self.dzTop / 2
    252236        self.InitDensityScaling = 1.0
    253237        self.ThermoDeltaTScaling = 1 / 11.0
    254238
    255         self.Vmean = 10 * np.ones((mesh.numberofelements, ))
    256 
    257         self.zMax = 250 * np.ones((mesh.numberofelements, ))
    258         self.zMin = 130 * np.ones((mesh.numberofelements, ))
    259         self.zY = 1.10 * np.ones((mesh.numberofelements, ))
     239        self.Vmean = 10 * np.ones((mesh.numberofelements,))
     240
     241        self.zMax = 250 * np.ones((mesh.numberofelements,))
     242        self.zMin = 130 * np.ones((mesh.numberofelements,))
     243        self.zY = 1.10 * np.ones((mesh.numberofelements,))
    260244        self.outputFreq = 30
    261245
    262     #additional albedo parameters
     246        #additional albedo parameters
    263247        self.aSnow = 0.85
    264248        self.aIce = 0.48
    … …  
    269253        self.adThresh = 1023
    270254
    271         self.teValue = np.ones((mesh.numberofelements, ))
    272         self.aValue = self.aSnow * np.ones(mesh.numberofelements, )
     255        self.teValue = np.ones((mesh.numberofelements,))
     256        self.aValue = self.aSnow * np.ones(mesh.numberofelements,)
    273257
    274258        self.Dzini = 0.05 * np.ones((mesh.numberofelements, 2))
    … …  
    277261        self.Gdnini = 0.0 * np.ones((mesh.numberofelements, 2))
    278262        self.Gspini = 0.0 * np.ones((mesh.numberofelements, 2))
    279         self.ECini = 0.0 * np.ones((mesh.numberofelements, ))
     263        self.ECini = 0.0 * np.ones((mesh.numberofelements,))
    280264        self.Wini = 0.0 * np.ones((mesh.numberofelements, 2))
    281265        self.Aini = self.aSnow * np.ones((mesh.numberofelements, 2))
    282266        self.Tini = 273.15 * np.ones((mesh.numberofelements, 2))
    283         #          / !\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh).
    284         #         If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp
    285         self.Sizeini = 2 * np.ones((mesh.numberofelements, ))
     267#       /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh).
     268#       If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp
     269        self.Sizeini = 2 * np.ones((mesh.numberofelements,))
    286270    #}}}
    287271
    288     def checkconsistency(self, md, solution, analyses):  # {{{
     272    def checkconsistency(self, md, solution, analyses):    # {{{
     273
    289274        md = checkfield(md, 'fieldname', 'smb.isgraingrowth', 'values', [0, 1])
    290275        md = checkfield(md, 'fieldname', 'smb.isalbedo', 'values', [0, 1])
    … …  
    297282        md = checkfield(md, 'fieldname', 'smb.isclimatology', 'values', [0, 1])
    298283
    299         md = checkfield(md, 'fieldname', 'smb.Ta', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100)  # - 100 / 100 celsius min / max value
    300         md = checkfield(md, 'fieldname', 'smb.V', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<', 45, 'size', np.shape(self.Ta))  #max 500 km / h
    301         md = checkfield(md, 'fieldname', 'smb.dswrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1400, 'size', np.shape(self.Ta))
    302         md = checkfield(md, 'fieldname', 'smb.dlwrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, 'size', np.shape(self.Ta))
    303         md = checkfield(md, 'fieldname', 'smb.P', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 100, 'size', np.shape(self.Ta))
     284        md = checkfield(md, 'fieldname', 'smb.Ta', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100)  #-100/100 celsius min/max value
     285        md = checkfield(md, 'fieldname', 'smb.V', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '<', 45, 'size', np.shape(self.Ta))  #max 500 km/h
     286        md = checkfield(md, 'fieldname', 'smb.dswrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1400, 'size', np.shape(self.Ta))
     287        md = checkfield(md, 'fieldname', 'smb.dlwrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, 'size', np.shape(self.Ta))
     288        md = checkfield(md, 'fieldname', 'smb.P', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 100, 'size', np.shape(self.Ta))
    304289        md = checkfield(md, 'fieldname', 'smb.eAir', 'timeseries', 1, 'NaN', 1, 'Inf', 1, 'size', np.shape(self.Ta))
    305290
    306291        if (self.isclimatology > 0):
    307             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')
    308             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')
    309             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')
    310             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')
    311             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')
    312             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')
    313 
    314         md = checkfield(md, 'fieldname', 'smb.Tmean', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100)  # - 100 / 100 celsius min / max value
    315         md = checkfield(md, 'fieldname', 'smb.C', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>=', 0)
    316         md = checkfield(md, 'fieldname', 'smb.Vmean', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>=', 0)
    317         md = checkfield(md, 'fieldname', 'smb.Tz', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 5000)
    318         md = checkfield(md, 'fieldname', 'smb.Vz', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 5000)
     292            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')
     293            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')
     294            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')
     295            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')
     296            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')
     297            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')
     298
     299        md = checkfield(md, 'fieldname', 'smb.Tmean', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100)  #-100/100 celsius min/max value
     300        md = checkfield(md, 'fieldname', 'smb.C', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0)
     301        md = checkfield(md, 'fieldname', 'smb.Vmean', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0)
     302        md = checkfield(md, 'fieldname', 'smb.Tz', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 5000)
     303        md = checkfield(md, 'fieldname', 'smb.Vz', 'size', [md.mesh.numberofelements], 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 5000)
    319304
    320305        md = checkfield(md, 'fieldname', 'smb.teValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
    … …  
    325310        md = checkfield(md, 'fieldname', 'smb.dsnowIdx', 'NaN', 1, 'Inf', 1, 'values', [0, 1, 2, 3, 4])
    326311
    327         md = checkfield(md, 'fieldname', 'smb.zTop', 'NaN', 1, 'Inf', 1, '>=', 0)
     312        md = checkfield(md, 'fieldname', 'smb.zTop', 'NaN', 1, 'Inf', 1, '> = ', 0)
    328313        md = checkfield(md, 'fieldname', 'smb.dzTop', 'NaN', 1, 'Inf', 1, '>', 0)
    329314        md = checkfield(md, 'fieldname', 'smb.dzMin', 'NaN', 1, 'Inf', 1, '>', 0)
    330         md = checkfield(md, 'fieldname', 'smb.zY', 'NaN', 1, 'Inf', 1, '>=', 1)
     315        md = checkfield(md, 'fieldname', 'smb.zY', 'NaN', 1, 'Inf', 1, '> = ', 1)
    331316        md = checkfield(md, 'fieldname', 'smb.outputFreq', 'NaN', 1, 'Inf', 1, '>', 0, '<', 10 * 365)  #10 years max
    332         md = checkfield(md, 'fieldname', 'smb.InitDensityScaling', 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
    333         md = checkfield(md, 'fieldname', 'smb.ThermoDeltaTScaling', 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
     317        md = checkfield(md, 'fieldname', 'smb.InitDensityScaling', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1)
     318        md = checkfield(md, 'fieldname', 'smb.ThermoDeltaTScaling', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1)
    334319        md = checkfield(md, 'fieldname', 'smb.adThresh', 'NaN', 1, 'Inf', 1, '>=', 0)
    335320
    336         if isinstance(self.aIdx, list) or isinstance(self.aIdx, np.ndarray) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)):
    337             md = checkfield(md, 'fieldname', 'smb.aSnow', 'NaN', 1, 'Inf', 1, '>=', .64, '<=', .89)
    338             md = checkfield(md, 'fieldname', 'smb.aIce', 'NaN', 1, 'Inf', 1, '>=', .27, '<=', .58)
     321        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)):
     322            md = checkfield(md, 'fieldname', 'smb.aSnow', 'NaN', 1, 'Inf', 1, '> = ', .64, '< = ', .89)
     323            md = checkfield(md, 'fieldname', 'smb.aIce', 'NaN', 1, 'Inf', 1, '> = ', .27, '< = ', .58)
    339324        elif self.aIdx == 0:
    340325            md = checkfield(md, 'fieldname', 'smb.aValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
    341326        elif self.aIdx == 3:
    342             md = checkfield(md, 'fieldname', 'smb.cldFrac', 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
     327            md = checkfield(md, 'fieldname', 'smb.cldFrac', 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1)
    343328        elif self.aIdx == 4:
    344             md = checkfield(md, 'fieldname', 'smb.t0wet', 'NaN', 1, 'Inf', 1, '>=', 15, '<=', 21.9)
    345             md = checkfield(md, 'fieldname', 'smb.t0dry', 'NaN', 1, 'Inf', 1, '>=', 30, '<=', 30)
    346             md = checkfield(md, 'fieldname', 'smb.K', 'NaN', 1, 'Inf', 1, '>=', 7, '<=', 7)
     329            md = checkfield(md, 'fieldname', 'smb.t0wet', 'NaN', 1, 'Inf', 1, '> = ', 15, '< = ', 21.9)
     330            md = checkfield(md, 'fieldname', 'smb.t0dry', 'NaN', 1, 'Inf', 1, '> = ', 30, '< = ', 30)
     331            md = checkfield(md, 'fieldname', 'smb.K', 'NaN', 1, 'Inf', 1, '> = ', 7, '< = ', 7)
    347332
    348333        #check zTop is < local thickness:
    349334        he = np.sum(md.geometry.thickness[md.mesh.elements - 1], axis=1) / np.size(md.mesh.elements, 1)
    350335        if np.any(he < self.zTop):
    351             error('SMBgemb consistency check error: zTop should be smaller than local ice thickness')
    352 
     336            raise IOError('SMBgemb consistency check error: zTop should be smaller than local ice thickness')
     337        md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
    353338        md = checkfield(md, 'fieldname', 'smb.requested_outputs', 'stringrow', 1)
    354339        return md
    355340    # }}}
    356341
    357     def marshall(self, prefix, md, fid):  # {{{
     342    def marshall(self, prefix, md, fid):    # {{{
     343
    358344        yts = md.constants.yts
    359345
    … …  
    408394        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'teValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
    409395
    410     #snow properties init
     396        #snow properties init
    411397        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Dzini', 'format', 'DoubleMat', 'mattype', 3)
    412398        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Dini', 'format', 'DoubleMat', 'mattype', 3)
    … …  
    419405        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Tini', 'format', 'DoubleMat', 'mattype', 3)
    420406        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Sizeini', 'format', 'IntMat', 'mattype', 2)
    421 
    422     #figure out dt from forcings:
     407        WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
     408        #figure out dt from forcings:
    423409        time = self.Ta[-1]  #assume all forcings are on the same time step
    424410        dtime = np.diff(time, n=1, axis=0)
    … …  
    427413        WriteData(fid, prefix, 'data', dt, 'name', 'md.smb.dt', 'format', 'Double', 'scale', yts)
    428414
    429     # Check if smb_dt goes evenly into transient core time step
     415        # Check if smb_dt goes evenly into transient core time step
    430416        if (md.timestepping.time_step % dt >= 1e-10):
    431             error('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)
    432 
    433     #process requested outputs
     417            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)
     418
     419        #process requested outputs
    434420        outputs = self.requested_outputs
    435421        indices = [i for i, x in enumerate(outputs) if x == 'default']
Note: See TracChangeset for help on using the changeset viewer.