Index: /issm/trunk-jpl/src/m/classes/SMBsemic.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBsemic.m (revision 27512) +++ /issm/trunk-jpl/src/m/classes/SMBsemic.m (revision 27513) @@ -32,6 +32,6 @@ % albedo - albedo = 0; % required for first energy balance calculation of SEMIC. - albedo_snow = 0; % required for ISBA method. + albedo = 0; % required for first energy balance calculation of SEMIC + albedo_snow = 0; % required for ISBA method albedo_scheme = 0; alb_smax = NaN; @@ -138,5 +138,5 @@ self.ismethod = 0; self.requested_outputs={'default'}; - end % }}}zo + end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -159,5 +159,5 @@ md = checkfield(md,'fieldname','smb.ismethod','numel',1,'values',[0,1]); - if self.ismethod == 1 % transient mode. + if self.ismethod == 1 % transient mode md = checkfield(md,'fieldname','smb.desfacElevation','>=',0,'numel',1); @@ -204,5 +204,5 @@ fielddisplay(self,'ismethod','method for calculating SMB with SEMIC. Default version of SEMIC is really slow. 0: steady, 1: transient (default: 0)'); - if self.ismethod == 1 % tranisent mode + if self.ismethod == 1 % transient mode fielddisplay(self,'desfacElevation','desertification elevation (default is 2000 m; Vizcaino et al. 2010)'); fielddisplay(self,'Tamp','amplitude of diurnal cycle [K]'); @@ -219,8 +219,8 @@ fielddisplay(self,'alb_smax','maximum snow albedo (default: 0.79)'); fielddisplay(self,'alb_smin','minimum snow albedo (default: 0.6)'); - fielddisplay(self,'albi','background albeod for bare ice (default: 0.41)'); - fielddisplay(self,'albl','background albeod for bare land (default: 0.07)'); - end - % albeod_scheme - 0: none, 1: slater, 2: isba, 3: denby, 4: alex. + fielddisplay(self,'albi','background albedo for bare ice (default: 0.41)'); + fielddisplay(self,'albl','background albedo for bare land (default: 0.07)'); + end + % albedo_scheme - 0: none, 1: slater, 2: isba, 3: denby, 4: alex. if self.albedo_scheme == 1 disp(sprintf('\n\tSEMIC snow albedo parameters for Slater et al, (1998).')); @@ -234,11 +234,11 @@ elseif self.albedo_scheme == 2, disp(sprintf('\n\tSEMIC snow albedo parameters for ISBA.? where is citation?')); - fielddisplay(self,'mcrit','critical melt rate (defaut: 6e-8) [unit: m/sec]'); - fielddisplay(self,'wcrit','critical liquid water content (defaut: 15) [unit: kg/m2]'); + fielddisplay(self,'mcrit','critical melt rate (default: 6e-8) [unit: m/sec]'); + fielddisplay(self,'wcrit','critical liquid water content (default: 15) [unit: kg/m2]'); fielddisplay(self,'tau_a','dry albedo decline [unit: 1/day]'); fielddisplay(self,'tau_f','wet albedo decline [unit: 1/day]'); elseif self.albedo_scheme == 3, disp(sprintf('\n\tSEMIC snow albedo parameters for Denby et al. (2002 Tellus)')); - fielddisplay(self,'mcrit','critical melt rate (defaut: 6e-8) [unit: m/sec]'); + fielddisplay(self,'mcrit','critical melt rate (default: 6e-8) [unit: m/sec]'); elseif self.albedo_scheme == 4, disp(sprintf('\n\tSEMIC snow albedo parameters for Alex.?')); @@ -278,5 +278,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailytemperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); % TODO: transient mode should be merged with SEMIC developed by Ruckamp et al. (2018). - if self.ismethod == 1% transient mode. + if self.ismethod == 1% transient mode WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tamp','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'class','smb','fieldname','mask','format','DoubleMat','mattype',1); @@ -287,5 +287,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','rcrit','format','Double'); - %albedo... + %albedo WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo_snow','format','DoubleMat','mattype',1); @@ -296,5 +296,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','alb_smax','format','Double'); - %abledo parameters for ? + %albedo parameters for ? %for slater WriteData(fid,prefix,'object',self,'class','smb','fieldname','tmin','format','Double'); Index: /issm/trunk-jpl/src/m/classes/SMBsemic.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBsemic.py (revision 27512) +++ /issm/trunk-jpl/src/m/classes/SMBsemic.py (revision 27513) @@ -25,5 +25,10 @@ self.dailyairhumidity = np.nan self.dailytemperature = np.nan + self.Tamp = np.nan + self.mask = np.nan + self.hice = np.nan + self.hsnow = np.nan self.desfac = 0 + self.desfacElevation = 0 self.rlaps = 0 self.rdl = 0 @@ -32,4 +37,36 @@ self.averaging = 0 self.requested_outputs = [] + + self.hcrit = 0 + self.rcrit = 0 + + # albedo + self.albedo = 0 # required for first energy balance calculation of SEMIC + self.albedo_snow = 0 # required for ISBA method + self.albedo_scheme = 0 + self.alb_smax = np.nan + self.alb_smin = np.nan + self.albi = np.nan + self.albl = np.nan + + # albedo parameters depending on albedo_scheme + # for slater + self.tmin = np.nan + self.tmax = np.nan + + # for isba & denby method + self.mcrit = np.nan + + # for isba + self.tau_a = np.nan + self.tau_f = np.nan + self.wcrit = np.nan + + # for alex + self.tmid = np.nan + self.afac = np.nan + + # method + self.ismethod = 0 if len(args) == 0: @@ -57,4 +94,48 @@ s += '{}\n'.format(fielddisplay(self, 'rdl', 'longwave downward radiation decrease (default is 0.29 [W/m^2/km]; Marty et al. 2002)')) s += '{}\n'.format(fielddisplay(self, 's0gcm', 'GCM reference elevation; (default is 0) [m]')) + s += '{}\n'.format(fielddisplay(self,'ismethod','method for calculating SMB with SEMIC. Default version of SEMIC is really slow. 0: steady, 1: transient (default: 0)')) + if self.ismethod: # transient mode + s += '{}\n'.format(fielddisplay(self,'desfacElevation','desertification elevation (default is 2000 m; Vizcaino et al. 2010)')) + s += '{}\n'.format(fielddisplay(self,'Tamp','amplitude of diurnal cycle [K]')) + s += '{}\n'.format(fielddisplay(self,'albedo','initial albedo [no unit]')) + s += '{}\n'.format(fielddisplay(self,'albedo_snow','initial albedo for snow [no unit]')) + s += '{}\n'.format(fielddisplay(self,'hice','initial thickness of ice [unit: m]')) + s += '{}\n'.format(fielddisplay(self,'hsnow','initial thickness of snow [unit: m]')) + s += '{}\n'.format(fielddisplay(self,'mask','masking for albedo. 0: ocean, 1: land, 2: ice (default: 2)')) + s += '{}\n'.format(fielddisplay(self,'hcrit','critical snow height for albedo [unit: m]')) + s += '{}\n'.format(fielddisplay(self,'rcrit','critical refreezing height for albedo [no unit]')) + + s += '\nSEMIC albedo parameters.\n' + s += '{}\n'.format(fielddisplay(self,'albedo_scheme','albedo scheme for SEMIC. 0: none, 1: slater, 2: isba, 3: denby, 4: alex (default is 0)')) + s += '{}\n'.format(fielddisplay(self,'alb_smax','maximum snow albedo (default: 0.79)')) + s += '{}\n'.format(fielddisplay(self,'alb_smin','minimum snow albedo (default: 0.6)')) + s += '{}\n'.format(fielddisplay(self,'albi','background albedo for bare ice (default: 0.41)')) + s += '{}\n'.format(fielddisplay(self,'albl','background albedo for bare land (default: 0.07)')) + # albedo_scheme - 0: none, 1: slater, 2: isba, 3: denby, 4: alex. + if self.albedo_scheme == 1: + s += '\n\tSEMIC snow albedo parameters for Slater et al, (1998).\n' + s += '\t alb = alb_smax - (alb_smax - alb_smin)*tm^(3.0)\n' + s += '\t tm = 1 (tsurf > 273.15 K)\n' + s += '\t tm = f*(tsurf-tmin) (tmin <= tsurf < 273.15)\n' + s += '\t 0 (tsurf < tmin)\n' + s += '\t f = 1/(273.15-tmin)\n' + s += '{}\n'.format(fielddisplay(self, 'tmin', 'minimum temperature for which albedo decline become effective. (default: 263.15 K)[unit: K])')) + s += '{}\n'.format(fielddisplay(self, 'tmax', 'maxmium temperature for which albedo decline become effective. This value should be fixed. (default: 273.15 K)[unit: K])')) + elif self.albedo_scheme == 2: + s += '\n\tSEMIC snow albedo parameters for ISBA.? where is citation?\n' + s += '{}\n'.format(fielddisplay(self, 'mcrit', 'critical melt rate (default: 6e-8) [unit: m/sec]')) + s += '{}\n'.format(fielddisplay(self, 'wcrit', 'critical liquid water content (default: 15) [unit: kg/m2]')) + s += '{}\n'.format(fielddisplay(self, 'tau_a', 'dry albedo decline [unit: 1/day]')) + s += '{}\n'.format(fielddisplay(self, 'tau_f', 'wet albedo decline [unit: 1/day]')) + elif self.albedo_scheme == 3: + s += '\n\tSEMIC snow albedo parameters for Denby et al. (2002 Tellus)\n' + s += '{}\n'.format(fielddisplay(self,'mcrit','critical melt rate (defaut: 6e-8) [unit: m/sec]')) + elif self.albedo_scheme == 4: + s += '\n\tSEMIC snow albedo parameters for Alex.?\n' + s += '{}\n'.format(fielddisplay(self,'afac','[unit: ?]')) + s += '{}\n'.format(fielddisplay(self,'tmid','[unit: ?]')) + else: + raise Exception('ERROR: {} is not supported albedo scheme.'.format(self.albedo_scheme)) + s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step')) s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) @@ -84,15 +165,60 @@ # }}} + def outputlists(self, md): # {{{ + if self.ismethod: + list = ['SmbMassBalance', 'SmbMassBalanceSnow', 'SmbMelt', 'SmbAccumulation', 'SmbHIce', 'SmbHSnow', 'SmbAlbedo', 'SmbAlbedoSnow', 'TemperatureSEMIC'] + else: + list = ['SmbMassBalance'] + return list + # }}} + def initialize(self, md): # {{{ if np.all(np.isnan(self.s0gcm)): self.s0gcm = np.zeros((md.mesh.numberofvertices)) - print(" no SMBsemic.s0gcm specified: values set as zero") + print(' no SMBsemic.s0gcm specified: values set as zero') + + self.Tamp = 3 * np.ones((md.mesh.numberofvertices,)) + #self.albedo = 0.8 * np.ones((md.mesh.numberofvertices,)) + #self.albedo_snow = 0.5 * np.ones((md.mesh.numberofvertices,)) + self.hice = np.zeros((md.mesh.numberofvertices,)) + self.hsnow = 5 * np.ones((md.mesh.numberofvertices,)) + return self # }}} def setdefaultparameters(self): #{{{ + # albedo parameters + self.albedo_scheme = 0 + self.alb_smax = 0.79 + self.alb_smin = 0.6 + self.albi = 0.41 + self.albl = 0.07 + + # albedo parameters for? + # for slater + self.tmin = 263.15 + self.tmax = 273.15 + + # for isba & denby + self.mcrit = 6e-8 + + # for isba + self.tau_a = 0.008 + self.tau_f = 0.24 + self.wcrit = 15.0 + + # for alex + self.tmid = 273.35 + self.afac = -0.18 + + self.hcrit = 0.028 # from Krapp et al. (2017) + self.rcrit = 0.85 # from Krapp et al. (2017) + self.desfac = -log(2.0) / 1000 + self.desfacElevation = 2000 self.rlaps = 7.4 self.rdl = 0.29 + + self.ismethod = 0 self.requested_outputs = ['default'] return self @@ -114,4 +240,22 @@ md = checkfield(md, 'fieldname', 'smb.dailyairhumidity', 'timeseries', 1, 'NaN', 1, 'Inf', 1) md = checkfield(md, 'fieldname', 'smb.dailytemperature', 'timeseries', 1, 'NaN', 1, 'Inf', 1) + + # TODO: transient model should be merged with SEMIC developed by Ruckamp et al. (2018) + md = checkfield(md, 'fieldname', 'smb.ismethod', 'numel', 1, 'values', [0, 1]) + if self.ismethod: # transient mode + md = checkfield(md, 'fieldname', 'smb.desfacElevation', '>=', 0, 'numel', 1) + md = checkfield(md, 'fieldname', 'smb.albedo_scheme', 'NaN', 1, 'Inf', 1, 'numel', 1, 'values', [0, 1, 2, 3, 4]) + md = checkfield(md, 'fieldname', 'smb.alb_smax', '>=', 0, 'NaN', 1, 'Inf', 1, 'numel', 1) + md = checkfield(md, 'fieldname', 'smb.mask', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 1], 'values', [0, 1, 2]) + + # initial values + md = checkfield(md, 'fieldname', 'smb.albedo', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 1]) + md = checkfield(md, 'fieldname', 'smb.albedo_snow', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 1]) + md = checkfield(md, 'fieldname', 'smb.alb_smax', '>=', 0, '<=', 1, 'NaN', 1, 'Inf', 1, 'numel', 1) + md = checkfield(md, 'fieldname', 'smb.alb_smin', '>=', 0, '<=', 1, 'NaN', 1, 'Inf', 1, 'numel', 1) + md = checkfield(md, 'fieldname', 'smb.albi', '>=', 0, '<=', 1, 'NaN', 1, 'Inf', 1, 'numel', 1) + md = checkfield(md, 'fieldname', 'smb.albl', '>=', 0, '<=', 1, 'NaN', 1, 'Inf', 1, 'numel', 1) + md = checkfield(md, 'fieldname', 'smb.hice', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 1]) + md = checkfield(md, 'fieldname', 'smb.hsnow', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 1]) md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2]) @@ -123,5 +267,7 @@ yts = md.constants.yts WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 12, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ismethod', 'format', 'Integer', 'values', [0, 1]) WriteData(fid, prefix, 'object', self, 'class' ,'smb', 'fieldname', 'desfac', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'desfacElevation', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class' ,'smb', 'fieldname', 's0gcm', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'class' ,'smb', 'fieldname', 'rlaps', 'format', 'Double') @@ -136,4 +282,36 @@ WriteData(fid, prefix, 'object', self, 'class' ,'smb', 'fieldname', 'dailyairhumidity', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class' ,'smb', 'fieldname', 'dailytemperature', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) + # TODO: transient mode should be merged with SEMIC developed by Ruckamp et al. (2018). + if self.ismethod: # transient mode + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Tamp', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'mask', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'hice', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'hsnow', 'format', 'DoubleMat', 'mattype', 1) + + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'hcrit', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'rcrit', 'format', 'Double') + + # albedo + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'albedo', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'albedo_snow', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'albedo_scheme', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'albi', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'albl', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'alb_smin', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'alb_smax', 'format', 'Double') + + # albedo parameters for ? + # for slater + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'tmin', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'tmax', 'format', 'Double') + # for isba & denby + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'mcrit', 'format', 'Double') + # for isba + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'wcrit', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'tau_a', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'tau_f', 'format', 'Double') + # for alex + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'tmid', 'format', 'Double') + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'afac', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')