![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 25836 for issm/trunk/src/m/classes/SMBgemb.py
- Timestamp:
- 12/08/20 08:45:53 (4 years ago)
- Location:
- issm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk ¶
- Property svn:mergeinfo changed
-
issm/trunk/src ¶
- Property svn:mergeinfo changed
-
TabularUnified issm/trunk/src/m/classes/SMBgemb.py ¶
r24313 r25836 1 1 import numpy as np 2 3 from checkfield import checkfield 2 4 from fielddisplay import fielddisplay 3 from checkfield import checkfield5 from project3d import project3d 4 6 from WriteData import WriteData 5 from project3d import project3d6 7 7 8 … … 14 15 """ 15 16 16 def __init__(self ): # {{{17 def __init__(self, *args): # {{{ 17 18 #each one of these properties is a transient forcing to the GEMB model, loaded from meteorological data derived 18 19 #from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number … … 21 22 #solution choices 22 23 #check these: 23 #isgraingrowth 24 #self.isgraingrowth = 0 25 #self.issmbgradients = 0 24 26 #isalbedo 25 27 #isshortwave … … 31 33 32 34 #inputs: 33 self.Ta = float('NaN') #2 m air temperature, in Kelvin34 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]35 self.Ta = np.nan # 2 m air temperature, in Kelvin 36 self.V = np.nan # wind speed (m/s-1) 37 self.dswrf = np.nan # downward shortwave radiation flux [W/m^2] 38 self.dlwrf = np.nan # downward longwave radiation flux [W/m^2] 39 self.P = np.nan # precipitation [mm w.e. / m^2] 40 self.eAir = np.nan # screen level vapor pressure [Pa] 41 self.pAir = np.nan # surface pressure [Pa] 42 self.Tmean = np.nan # mean annual temperature [K] 43 self.Vmean = np.nan # mean annual wind velocity [m s-1] 44 self.C = np.nan # mean annual snow accumulation [kg m-2 yr-1] 45 self.Tz = np.nan # height above ground at which temperature (T) was sampled [m] 46 self.Vz = np.nan # height above ground at which wind (V) was sampled [m] 45 47 46 48 #optional inputs: 47 self.aValue = float('NaN') #Albedo forcing at every element. Used only if aIdx == 0.48 self.teValue = float('NaN') #Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)49 self.aValue = np.nan # Albedo forcing at every element. Used only if aIdx == 0. 50 self.teValue = np.nan # Outward longwave radiation thermal emissivity forcing at every element (default in code is 1) 49 51 50 52 # 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)60 self.Sizeini = float('NaN') #Number of layers53 self.Dzini = np.nan # cell depth (m) 54 self.Dini = np.nan # snow density (kg m-3) 55 self.Reini = np.nan # effective grain size (mm) 56 self.Gdnini = np.nan # grain dricity (0-1) 57 self.Gspini = np.nan # grain sphericity (0-1) 58 self.ECini = np.nan # evaporation/condensation (kg m-2) 59 self.Wini = np.nan # Water content (kg m-2) 60 self.Aini = np.nan # albedo (0-1) 61 self.Tini = np.nan # snow temperature (K) 62 self.Sizeini = np.nan # Number of layers 61 63 62 64 #settings: 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): 66 self.dsnowIdx = float('NaN') #model for fresh snow accumulation density (default is 1): 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) 65 self.aIdx = np.nan # method for calculating albedo and subsurface absorption (default is 1) 66 # 0: direct input from aValue parameter 67 # 1: effective grain radius [Gardner & Sharp, 2009] 68 # 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) 69 # 3: density and cloud amount [Greuell & Konzelmann, 1994] 70 # 4: exponential time decay & wetness [Bougamont & Bamber, 2005] 71 72 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)) 73 self.denIdx = np.nan # densification model to use (default is 2): 74 # 1 = emperical model of Herron and Langway (1980) 75 # 2 = semi-emperical model of Anthern et al. (2010) 76 # 3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010) 77 # 4 = DO NOT USE: emperical model of Li and Zwally (2004) 78 # 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008) 79 # 6 = Antarctica semi-emperical model of Ligtenberg et al. (2011) 80 # 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015) 81 82 self.dsnowIdx = np.nan # model for fresh snow accumulation density (default is 1): 83 # 0 = Original GEMB value, 150 kg/m^3 84 # 1 = Antarctica value of fresh snow density, 350 kg/m^3 85 # 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008) 86 # 3 = Antarctica model of Kaspers et al. (2004) 87 # 4 = Greenland model of Kuipers Munneke et al. (2015) 88 89 self.zTop = np.nan # depth over which grid length is constant at the top of the snopack (default 10) [m] 90 self.dzTop = np.nan # initial top vertical grid spacing (default .05) [m] 91 self.dzMin = np.nan # initial min vertical allowable grid spacing (default dzMin/2) [m] 92 self.zY = np.nan # strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)] 93 self.zMax = np.nan # initial max model depth (default is min(thickness, 250)) [m] 94 self.zMin = np.nan # initial min model depth (default is min(thickness, 130)) [m] 95 self.outputFreq = np.nan # output frequency in days (default is monthly, 30) 74 96 75 97 #specific albedo parameters: 76 98 #Method 1 and 2: 77 self.aSnow = float('NaN')# new snow albedo (0.64 - 0.89)78 self.aIce = float('NaN')# range 0.27-0.58 for old snow99 self.aSnow = np.nan # new snow albedo (0.64 - 0.89) 100 self.aIce = np.nan # range 0.27-0.58 for old snow 79 101 #Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo 80 self.cldFrac = float('NaN')# average cloud amount102 self.cldFrac = np.nan # average cloud amount 81 103 #Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005) 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)85 self.adThresh = float('NaN')# Apply aIdx method to all areas with densities below this value,104 self.t0wet = np.nan # time scale for wet snow (15-21.9) 105 self.t0dry = np.nan # warm snow timescale (30) 106 self.K = np.nan # time scale temperature coef. (7) 107 self.adThresh = np.nan # Apply aIdx method to all areas with densities below this value, 86 108 # or else apply direct input value from aValue, allowing albedo to be altered. 87 109 # Default value is rho water (1023 kg m-3). 88 110 89 111 #densities: 90 self.InitDensityScaling = float('NaN') #initial scaling factor multiplying the density of ice, which describes the density of the snowpack.112 self.InitDensityScaling = np.nan # initial scaling factor multiplying the density of ice, which describes the density of the snowpack. 91 113 92 114 #thermo: 93 self.ThermoDeltaTScaling = float('NaN') #scaling factor to multiply the thermal diffusion timestep (delta t) 94 95 self.steps_per_step = 1 96 self.requested_outputs = [] 115 self.ThermoDeltaTScaling = np.nan # scaling factor to multiply the thermal diffusion timestep (delta t) 116 117 self.steps_per_step = 1 118 self.averaging = 0 119 self.requested_outputs = [] 97 120 98 121 #Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM. … … 101 124 #elev: this is taken from the ISSM surface itself. 102 125 126 nargin = len(args) 127 if nargin: 128 if nargin == 2: 129 mesh = args[0] 130 geometry = args[1] 131 self.setdefaultparameters(mesh, geometry) 132 else: 133 raise Exception('constructor not supported: need mesh and geometry to set defaults') 103 134 #}}} 104 105 135 106 136 def __repr__(self): # {{{ … … 118 148 string = "%s\n%s" % (string, fielddisplay(self, 'isdensification', 'run densification module (default true)')) 119 149 string = "%s\n%s" % (string, fielddisplay(self, 'isturbulentflux', 'run turbulant heat fluxes module (default true)')) 150 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')) 120 151 string = "%s\n%s" % (string, fielddisplay(self, 'isclimatology', 'repeat all forcings when past last forcing time (default false)')) 121 152 string = "%s\n%s" % (string, fielddisplay(self, 'Ta', '2 m air temperature, in Kelvin')) … … 144 175 '0: direct input from aValue parameter', 145 176 '1: effective grain radius [Gardner & Sharp, 2009]', 146 '2: effective grain radius [Brun et al., 2009]',177 '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)', 147 178 '3: density and cloud amount [Greuell & Konzelmann, 1994]', 148 179 '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'])) … … 158 189 string = "%s\n%s" % (string, fielddisplay(self, 'Aini', 'Initial albedo when restart [-]')) 159 190 string = "%s\n%s" % (string, fielddisplay(self, 'Tini', 'Initial snow temperature when restart [K]')) 160 string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [ K]'))191 string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [-]')) 161 192 162 193 #additional albedo parameters: … … 175 206 string = "%s\n%s" % (string, fielddisplay(self, 'K', 'time scale temperature coef. (7) [d]')) 176 207 177 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 ]'))208 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)]')) 178 209 string = "%s\n%s" % (string, fielddisplay(self, 'denIdx', ['densification model to use (default is 2):', 179 210 '1 = emperical model of Herron and Langway (1980)', … … 191 222 '4 = Greenland model of Kuipers Munneke et al. (2015)'])) 192 223 string = "%s\n%s" % (string, fielddisplay(self, 'steps_per_step', 'number of smb steps per time step')) 224 string = "%s\n%s" % (string, fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) 225 string = "%s\n\t\t%s" % (string, '0: Arithmetic (default)') 226 string = "%s\n\t\t%s" % (string, '1: Geometric') 227 string = "%s\n\t\t%s" % (string, '2: Harmonic') 193 228 string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested')) 194 229 return string … … 218 253 def setdefaultparameters(self, mesh, geometry): # {{{ 219 254 self.isgraingrowth = 1 255 self.issmbgradients = 0 220 256 self.isalbedo = 1 221 257 self.isshortwave = 1 … … 226 262 self.isturbulentflux = 1 227 263 self.isclimatology = 0 264 self.isconstrainsurfaceT = 0 228 265 229 266 self.aIdx = 1 … … 241 278 self.zMax = 250 * np.ones((mesh.numberofelements,)) 242 279 self.zMin = 130 * np.ones((mesh.numberofelements,)) 243 self.zY = 1. 10* np.ones((mesh.numberofelements,))280 self.zY = 1.025 * np.ones((mesh.numberofelements,)) 244 281 self.outputFreq = 30 245 282 … … 273 310 274 311 md = checkfield(md, 'fieldname', 'smb.isgraingrowth', 'values', [0, 1]) 312 md = checkfield(md, 'fieldname', 'smb.issmbgradients', 'values', [0, 1]) 275 313 md = checkfield(md, 'fieldname', 'smb.isalbedo', 'values', [0, 1]) 276 314 md = checkfield(md, 'fieldname', 'smb.isshortwave', 'values', [0, 1]) … … 281 319 md = checkfield(md, 'fieldname', 'smb.isturbulentflux', 'values', [0, 1]) 282 320 md = checkfield(md, 'fieldname', 'smb.isclimatology', 'values', [0, 1]) 321 md = checkfield(md, 'fieldname', 'smb.isconstrainsurfaceT', 'values', [0, 1]) 283 322 284 323 md = checkfield(md, 'fieldname', 'smb.Ta', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>', 273 - 100, '<', 273 + 100) #-100/100 celsius min/max value … … 286 325 md = checkfield(md, 'fieldname', 'smb.dswrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1400, 'size', np.shape(self.Ta)) 287 326 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))327 md = checkfield(md, 'fieldname', 'smb.P', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 200, 'size', np.shape(self.Ta)) 289 328 md = checkfield(md, 'fieldname', 'smb.eAir', 'timeseries', 1, 'NaN', 1, 'Inf', 1, 'size', np.shape(self.Ta)) 290 329 … … 336 375 raise IOError('SMBgemb consistency check error: zTop should be smaller than local ice thickness') 337 376 md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) 377 md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2]) 338 378 md = checkfield(md, 'fieldname', 'smb.requested_outputs', 'stringrow', 1) 339 379 return md … … 341 381 342 382 def marshall(self, prefix, md, fid): # {{{ 343 344 383 yts = md.constants.yts 345 384 … … 347 386 348 387 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isgraingrowth', 'format', 'Boolean') 388 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'issmbgradients', 'format', 'Boolean') 349 389 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isalbedo', 'format', 'Boolean') 350 390 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isshortwave', 'format', 'Boolean') … … 355 395 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isturbulentflux', 'format', 'Boolean') 356 396 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isclimatology', 'format', 'Boolean') 357 397 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'isconstrainsurfaceT', 'format', 'Boolean') 358 398 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Ta', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) 359 399 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'V', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) … … 406 446 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Sizeini', 'format', 'IntMat', 'mattype', 2) 407 447 WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') 448 WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer') 408 449 #figure out dt from forcings: 450 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)): 451 raise IOError('All GEMB forcings (Ta, P, V, dswrf, dlwrf, eAir, pAir) must have the same time steps in the final row!') 452 409 453 time = self.Ta[-1] #assume all forcings are on the same time step 410 454 dtime = np.diff(time, n=1, axis=0) 411 dt = min(dtime) / yts455 dt = min(dtime) 412 456 413 457 WriteData(fid, prefix, 'data', dt, 'name', 'md.smb.dt', 'format', 'Double', 'scale', yts)
Note:
See TracChangeset
for help on using the changeset viewer.
