Changeset 27513


Ignore:
Timestamp:
01/12/23 14:28:07 (2 years ago)
Author:
jdquinn
Message:

CHG: MATLAB -> Python

Location:
issm/trunk-jpl/src/m/classes
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/classes/SMBsemic.m

    r27511 r27513  
    3232
    3333                % albedo
    34                 albedo            = 0; % required for first energy balance calculation of SEMIC.
    35                 albedo_snow       = 0; % required for ISBA method.
     34                albedo            = 0; % required for first energy balance calculation of SEMIC
     35                albedo_snow       = 0; % required for ISBA method
    3636                albedo_scheme     = 0;
    3737                alb_smax = NaN;
    … …  
    138138                        self.ismethod        = 0;
    139139                        self.requested_outputs={'default'};
    140                 end % }}}zo
     140                end % }}}
    141141                function md = checkconsistency(self,md,solution,analyses) % {{{
    142142
    … …  
    159159
    160160                                md = checkfield(md,'fieldname','smb.ismethod','numel',1,'values',[0,1]);
    161                                 if self.ismethod == 1 % transient mode.
     161                                if self.ismethod == 1 % transient mode
    162162                                        md = checkfield(md,'fieldname','smb.desfacElevation','>=',0,'numel',1);
    163163
    … …  
    204204
    205205                        fielddisplay(self,'ismethod','method for calculating SMB with SEMIC. Default version of SEMIC is really slow. 0: steady, 1: transient (default: 0)');
    206                         if self.ismethod == 1 % tranisent mode
     206                        if self.ismethod == 1 % transient mode
    207207                                fielddisplay(self,'desfacElevation','desertification elevation (default is 2000 m; Vizcaino et al. 2010)');
    208208                                fielddisplay(self,'Tamp','amplitude of diurnal cycle [K]');
    … …  
    219219                                fielddisplay(self,'alb_smax','maximum snow albedo (default: 0.79)');
    220220                                fielddisplay(self,'alb_smin','minimum snow albedo (default: 0.6)');
    221                                 fielddisplay(self,'albi','background albeod for bare ice (default: 0.41)');
    222                                 fielddisplay(self,'albl','background albeod for bare land (default: 0.07)');
    223                         end
    224                         % albeod_scheme - 0: none, 1: slater, 2: isba, 3: denby, 4: alex.
     221                                fielddisplay(self,'albi','background albedo for bare ice (default: 0.41)');
     222                                fielddisplay(self,'albl','background albedo for bare land (default: 0.07)');
     223                        end
     224                        % albedo_scheme - 0: none, 1: slater, 2: isba, 3: denby, 4: alex.
    225225                        if self.albedo_scheme == 1
    226226                                disp(sprintf('\n\tSEMIC snow albedo parameters for Slater et al, (1998).'));
    … …  
    234234                        elseif self.albedo_scheme == 2,
    235235                                disp(sprintf('\n\tSEMIC snow albedo parameters for ISBA.? where is citation?'));
    236                                 fielddisplay(self,'mcrit','critical melt rate (defaut: 6e-8) [unit: m/sec]');
    237                                 fielddisplay(self,'wcrit','critical liquid water content (defaut: 15) [unit: kg/m2]');
     236                                fielddisplay(self,'mcrit','critical melt rate (default: 6e-8) [unit: m/sec]');
     237                                fielddisplay(self,'wcrit','critical liquid water content (default: 15) [unit: kg/m2]');
    238238                                fielddisplay(self,'tau_a','dry albedo decline [unit: 1/day]');
    239239                                fielddisplay(self,'tau_f','wet albedo decline [unit: 1/day]');
    240240                        elseif self.albedo_scheme == 3,
    241241                                disp(sprintf('\n\tSEMIC snow albedo parameters for Denby et al. (2002 Tellus)'));
    242                                 fielddisplay(self,'mcrit','critical melt rate (defaut: 6e-8) [unit: m/sec]');
     242                                fielddisplay(self,'mcrit','critical melt rate (default: 6e-8) [unit: m/sec]');
    243243                        elseif self.albedo_scheme == 4,
    244244                                disp(sprintf('\n\tSEMIC snow albedo parameters for Alex.?'));
    … …  
    278278                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailytemperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
    279279                        % TODO: transient mode should be merged with SEMIC developed by Ruckamp et al. (2018).
    280                         if self.ismethod == 1% transient mode.
     280                        if self.ismethod == 1% transient mode
    281281                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tamp','format','DoubleMat','mattype',1);
    282282                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','mask','format','DoubleMat','mattype',1);
    … …  
    287287                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','rcrit','format','Double');
    288288
    289                                 %albedo...
     289                                %albedo
    290290                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo','format','DoubleMat','mattype',1);
    291291                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo_snow','format','DoubleMat','mattype',1);
    … …  
    296296                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','alb_smax','format','Double');
    297297
    298                                 %abledo parameters for ?
     298                                %albedo parameters for ?
    299299                                %for slater
    300300                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','tmin','format','Double');

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

    r27458 r27513  
    2525        self.dailyairhumidity = np.nan
    2626        self.dailytemperature = np.nan
     27        self.Tamp = np.nan
     28        self.mask = np.nan
     29        self.hice = np.nan
     30        self.hsnow = np.nan
    2731        self.desfac = 0
     32        self.desfacElevation = 0
    2833        self.rlaps = 0
    2934        self.rdl = 0
    … …  
    3237        self.averaging = 0
    3338        self.requested_outputs = []
     39
     40        self.hcrit = 0
     41        self.rcrit = 0
     42
     43        # albedo
     44        self.albedo = 0 # required for first energy balance calculation of SEMIC
     45        self.albedo_snow = 0 # required for ISBA method
     46        self.albedo_scheme = 0
     47        self.alb_smax = np.nan
     48        self.alb_smin = np.nan
     49        self.albi = np.nan
     50        self.albl = np.nan
     51
     52        # albedo parameters depending on albedo_scheme
     53        # for slater
     54        self.tmin = np.nan
     55        self.tmax = np.nan
     56
     57        # for isba & denby method
     58        self.mcrit = np.nan
     59
     60        # for isba
     61        self.tau_a = np.nan
     62        self.tau_f = np.nan
     63        self.wcrit = np.nan
     64
     65        # for alex
     66        self.tmid = np.nan
     67        self.afac = np.nan
     68
     69        # method
     70        self.ismethod = 0
    3471
    3572        if len(args) == 0:
    … …  
    5794        s += '{}\n'.format(fielddisplay(self, 'rdl', 'longwave downward radiation decrease (default is 0.29 [W/m^2/km]; Marty et al. 2002)'))
    5895        s += '{}\n'.format(fielddisplay(self, 's0gcm', 'GCM reference elevation; (default is 0) [m]'))
     96        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)'))
     97        if self.ismethod: # transient mode
     98            s += '{}\n'.format(fielddisplay(self,'desfacElevation','desertification elevation (default is 2000 m; Vizcaino et al. 2010)'))
     99            s += '{}\n'.format(fielddisplay(self,'Tamp','amplitude of diurnal cycle [K]'))
     100            s += '{}\n'.format(fielddisplay(self,'albedo','initial albedo [no unit]'))
     101            s += '{}\n'.format(fielddisplay(self,'albedo_snow','initial albedo for snow [no unit]'))
     102            s += '{}\n'.format(fielddisplay(self,'hice','initial thickness of ice [unit: m]'))
     103            s += '{}\n'.format(fielddisplay(self,'hsnow','initial thickness of snow [unit: m]'))
     104            s += '{}\n'.format(fielddisplay(self,'mask','masking for albedo. 0: ocean, 1: land, 2: ice (default: 2)'))
     105            s += '{}\n'.format(fielddisplay(self,'hcrit','critical snow height for albedo [unit: m]'))
     106            s += '{}\n'.format(fielddisplay(self,'rcrit','critical refreezing height for albedo [no unit]'))
     107
     108            s += '\nSEMIC albedo parameters.\n'
     109            s += '{}\n'.format(fielddisplay(self,'albedo_scheme','albedo scheme for SEMIC. 0: none, 1: slater, 2: isba, 3: denby, 4: alex (default is 0)'))
     110            s += '{}\n'.format(fielddisplay(self,'alb_smax','maximum snow albedo (default: 0.79)'))
     111            s += '{}\n'.format(fielddisplay(self,'alb_smin','minimum snow albedo (default: 0.6)'))
     112            s += '{}\n'.format(fielddisplay(self,'albi','background albedo for bare ice (default: 0.41)'))
     113            s += '{}\n'.format(fielddisplay(self,'albl','background albedo for bare land (default: 0.07)'))
     114        # albedo_scheme - 0: none, 1: slater, 2: isba, 3: denby, 4: alex.
     115        if self.albedo_scheme == 1:
     116            s += '\n\tSEMIC snow albedo parameters for Slater et al, (1998).\n'
     117            s += '\t   alb = alb_smax - (alb_smax - alb_smin)*tm^(3.0)\n'
     118            s += '\t   tm  = 1 (tsurf > 273.15 K)\n'
     119            s += '\t         tm = f*(tsurf-tmin) (tmin <= tsurf < 273.15)\n'
     120            s += '\t         0 (tsurf < tmin)\n'
     121            s += '\t   f = 1/(273.15-tmin)\n'
     122            s += '{}\n'.format(fielddisplay(self, 'tmin', 'minimum temperature for which albedo decline become effective. (default: 263.15 K)[unit: K])'))
     123            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])'))
     124        elif self.albedo_scheme == 2:
     125            s += '\n\tSEMIC snow albedo parameters for ISBA.? where is citation?\n'
     126            s += '{}\n'.format(fielddisplay(self, 'mcrit', 'critical melt rate (default: 6e-8) [unit: m/sec]'))
     127            s += '{}\n'.format(fielddisplay(self, 'wcrit', 'critical liquid water content (default: 15) [unit: kg/m2]'))
     128            s += '{}\n'.format(fielddisplay(self, 'tau_a', 'dry albedo decline [unit: 1/day]'))
     129            s += '{}\n'.format(fielddisplay(self, 'tau_f', 'wet albedo decline [unit: 1/day]'))
     130        elif self.albedo_scheme == 3:
     131            s += '\n\tSEMIC snow albedo parameters for Denby et al. (2002 Tellus)\n'
     132            s += '{}\n'.format(fielddisplay(self,'mcrit','critical melt rate (defaut: 6e-8) [unit: m/sec]'))
     133        elif self.albedo_scheme == 4:
     134            s += '\n\tSEMIC snow albedo parameters for Alex.?\n'
     135            s += '{}\n'.format(fielddisplay(self,'afac','[unit: ?]'))
     136            s += '{}\n'.format(fielddisplay(self,'tmid','[unit: ?]'))
     137        else:
     138            raise Exception('ERROR: {} is not supported albedo scheme.'.format(self.albedo_scheme))
     139
    59140        s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
    60141        s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps'))
    … …  
    84165    # }}}
    85166
     167    def outputlists(self, md):  # {{{
     168        if self.ismethod:
     169            list = ['SmbMassBalance', 'SmbMassBalanceSnow', 'SmbMelt', 'SmbAccumulation', 'SmbHIce', 'SmbHSnow', 'SmbAlbedo', 'SmbAlbedoSnow', 'TemperatureSEMIC']
     170        else:
     171            list = ['SmbMassBalance']
     172        return list
     173    # }}}
     174
    86175    def initialize(self, md):  # {{{
    87176        if np.all(np.isnan(self.s0gcm)):
    88177            self.s0gcm = np.zeros((md.mesh.numberofvertices))
    89             print("      no SMBsemic.s0gcm specified: values set as zero")
     178            print('      no SMBsemic.s0gcm specified: values set as zero')
     179
     180        self.Tamp = 3 * np.ones((md.mesh.numberofvertices,))
     181        #self.albedo = 0.8 * np.ones((md.mesh.numberofvertices,))
     182        #self.albedo_snow = 0.5 * np.ones((md.mesh.numberofvertices,))
     183        self.hice = np.zeros((md.mesh.numberofvertices,))
     184        self.hsnow = 5 * np.ones((md.mesh.numberofvertices,))
     185
    90186        return self
    91187    # }}}
    92188
    93189    def setdefaultparameters(self):  #{{{
     190        # albedo parameters
     191        self.albedo_scheme = 0
     192        self.alb_smax = 0.79
     193        self.alb_smin = 0.6
     194        self.albi = 0.41
     195        self.albl = 0.07
     196
     197        # albedo parameters for?
     198        # for slater
     199        self.tmin = 263.15
     200        self.tmax = 273.15
     201
     202        # for isba & denby
     203        self.mcrit = 6e-8
     204
     205        # for isba
     206        self.tau_a = 0.008
     207        self.tau_f = 0.24
     208        self.wcrit = 15.0
     209
     210        # for alex
     211        self.tmid = 273.35
     212        self.afac = -0.18
     213
     214        self.hcrit = 0.028 # from Krapp et al. (2017)
     215        self.rcrit = 0.85 # from Krapp et al. (2017)
     216
    94217        self.desfac = -log(2.0) / 1000
     218        self.desfacElevation = 2000
    95219        self.rlaps = 7.4
    96220        self.rdl = 0.29
     221
     222        self.ismethod = 0
    97223        self.requested_outputs = ['default']
    98224        return self
    … …  
    114240            md = checkfield(md, 'fieldname', 'smb.dailyairhumidity', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
    115241            md = checkfield(md, 'fieldname', 'smb.dailytemperature', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
     242
     243            # TODO: transient model should be merged with SEMIC developed by Ruckamp et al. (2018)
     244            md = checkfield(md, 'fieldname', 'smb.ismethod', 'numel', 1, 'values', [0, 1])
     245            if self.ismethod: # transient mode
     246                md = checkfield(md, 'fieldname', 'smb.desfacElevation', '>=', 0, 'numel', 1)
     247                md = checkfield(md, 'fieldname', 'smb.albedo_scheme', 'NaN', 1, 'Inf', 1, 'numel', 1, 'values', [0, 1, 2, 3, 4])
     248                md = checkfield(md, 'fieldname', 'smb.alb_smax', '>=', 0, 'NaN', 1, 'Inf', 1, 'numel', 1)
     249                md = checkfield(md, 'fieldname', 'smb.mask', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 1], 'values', [0, 1, 2])
     250
     251                # initial values
     252                md = checkfield(md, 'fieldname', 'smb.albedo', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 1])
     253                md = checkfield(md, 'fieldname', 'smb.albedo_snow', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 1])
     254                md = checkfield(md, 'fieldname', 'smb.alb_smax', '>=', 0, '<=', 1, 'NaN', 1, 'Inf', 1, 'numel', 1)
     255                md = checkfield(md, 'fieldname', 'smb.alb_smin', '>=', 0, '<=', 1, 'NaN', 1, 'Inf', 1, 'numel', 1)
     256                md = checkfield(md, 'fieldname', 'smb.albi', '>=', 0, '<=', 1, 'NaN', 1, 'Inf', 1, 'numel', 1)
     257                md = checkfield(md, 'fieldname', 'smb.albl', '>=', 0, '<=', 1, 'NaN', 1, 'Inf', 1, 'numel', 1)
     258                md = checkfield(md, 'fieldname', 'smb.hice', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 1])
     259                md = checkfield(md, 'fieldname', 'smb.hsnow', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices, 1])
    116260        md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
    117261        md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2])
    … …  
    123267        yts = md.constants.yts
    124268        WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 12, 'format', 'Integer')
     269        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ismethod', 'format', 'Integer', 'values', [0, 1])
    125270        WriteData(fid, prefix, 'object', self, 'class' ,'smb', 'fieldname', 'desfac', 'format', 'Double')
     271        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'desfacElevation', 'format', 'Double')
    126272        WriteData(fid, prefix, 'object', self, 'class' ,'smb', 'fieldname', 's0gcm', 'format', 'DoubleMat', 'mattype', 1)
    127273        WriteData(fid, prefix, 'object', self, 'class' ,'smb', 'fieldname', 'rlaps', 'format', 'Double')
    … …  
    136282        WriteData(fid, prefix, 'object', self, 'class' ,'smb', 'fieldname', 'dailyairhumidity', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
    137283        WriteData(fid, prefix, 'object', self, 'class' ,'smb', 'fieldname', 'dailytemperature', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
     284        # TODO: transient mode should be merged with SEMIC developed by Ruckamp et al. (2018).
     285        if self.ismethod: # transient mode
     286            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Tamp', 'format', 'DoubleMat', 'mattype', 1)
     287            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'mask', 'format', 'DoubleMat', 'mattype', 1)
     288            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'hice', 'format', 'DoubleMat', 'mattype', 1)
     289            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'hsnow', 'format', 'DoubleMat', 'mattype', 1)
     290
     291            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'hcrit', 'format', 'Double')
     292            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'rcrit', 'format', 'Double')
     293
     294            # albedo
     295            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'albedo', 'format', 'DoubleMat', 'mattype', 1)
     296            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'albedo_snow', 'format', 'DoubleMat', 'mattype', 1)
     297            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'albedo_scheme', 'format', 'Integer')
     298            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'albi', 'format', 'Double')
     299            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'albl', 'format', 'Double')
     300            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'alb_smin', 'format', 'Double')
     301            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'alb_smax', 'format', 'Double')
     302
     303            # albedo parameters for ?
     304            # for slater
     305            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'tmin', 'format', 'Double')
     306            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'tmax', 'format', 'Double')
     307            # for isba & denby
     308            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'mcrit', 'format', 'Double')
     309            # for isba
     310            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'wcrit', 'format', 'Double')
     311            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'tau_a', 'format', 'Double')
     312            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'tau_f', 'format', 'Double')
     313            # for alex
     314            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'tmid', 'format', 'Double')
     315            WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'afac', 'format', 'Double')
    138316        WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer')
    139317        WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer')
Note: See TracChangeset for help on using the changeset viewer.