source: issm/trunk-jpl/src/m/classes/SMBsemic.m@ 27513

Last change on this file since 27513 was 27513, checked in by jdquinn, 2 years ago

CHG: MATLAB -> Python

File size: 16.6 KB
Line 
1%SMBsemic Class definition
2%
3% Usage:
4% SMBsemic=SMBsemic();
5
6classdef SMBsemic
7 properties (SetAccess=public)
8 dailysnowfall = NaN;
9 dailyrainfall = NaN;
10 dailydsradiation = NaN;
11 dailydlradiation = NaN;
12 dailywindspeed = NaN;
13 dailypressure = NaN;
14 dailyairdensity = NaN;
15 dailyairhumidity = NaN;
16 dailytemperature = NaN;
17 Tamp = NaN;
18 mask = NaN;
19 hice = NaN;
20 hsnow = NaN;
21 desfac = 0;
22 desfacElevation = 0;
23 rlaps = 0;
24 rdl = 0;
25 s0gcm = NaN;
26 steps_per_step = 1;
27 averaging = 0;
28 requested_outputs = {};
29
30 hcrit = 0;
31 rcrit = 0;
32
33 % albedo
34 albedo = 0; % required for first energy balance calculation of SEMIC
35 albedo_snow = 0; % required for ISBA method
36 albedo_scheme = 0;
37 alb_smax = NaN;
38 alb_smin = NaN;
39 albi = NaN;
40 albl = NaN;
41
42 % albedo parameters depending on albedo_scheme
43 % for slater
44 tmin = NaN;
45 tmax = NaN;
46
47 % for isba & denby method
48 mcrit = NaN;
49
50 % for isba
51 tau_a = NaN;
52 tau_f = NaN;
53 wcrit = NaN;
54
55 % for alex
56 tmid = NaN;
57 afac = NaN;
58
59 % method
60 ismethod = 0;
61 end
62 methods
63 function self = SMBsemic(varargin) % {{{
64 switch nargin
65 case 0
66 self=setdefaultparameters(self);
67 otherwise
68 error('constructor not supported');
69 end
70 end % }}}
71 function self = extrude(self,md) % {{{
72 self.dailysnowfall=project3d(md,'vector',self.dailysnowfall,'type','node');
73 self.dailyrainfall=project3d(md,'vector',self.dailyrainfall,'type','node');
74 self.dailydsradiation=project3d(md,'vector',self.dailydsradiation,'type','node');
75 self.dailydlradiation=project3d(md,'vector',self.dailydlradiation,'type','node');
76 self.dailywindspeed=project3d(md,'vector',self.dailywindspeed,'type','node');
77 self.dailypressure=project3d(md,'vector',self.dailypressure,'type','node');
78 self.dailyairdensity=project3d(md,'vector',self.dailyairdensity,'type','node');
79 self.dailyairhumidity=project3d(md,'vector',self.dailyairhumidity,'type','node');
80 self.dailytemperature=project3d(md,'vector',self.dailytemperature,'type','node');
81 self.s0gcm=project3d(md,'vector',self.s0gcm,'type','node');
82
83 end % }}}
84 function list = defaultoutputs(self,md) % {{{
85 list = {'SmbMassBalance'};
86 end % }}}
87 function list = outputlists(self,md) % {{{
88 if self.ismethod == 1
89 list = {'SmbMassBalance','SmbMassBalanceSnow','SmbMelt','SmbAccumulation',...
90 'SmbHIce','SmbHSnow','SmbAlbedo','SmbAlbedoSnow','TemperatureSEMIC'};
91 else
92 list = {'SmbMassBalance'};
93 end
94 end % }}}
95 function self = initialize(self,md) % {{{
96
97 if isnan(self.s0gcm),
98 self.s0gcm=zeros(md.mesh.numberofvertices,1);
99 disp(' no SMBsemic.s0gcm specified: values set as zero');
100 end
101 self.Tamp = 3*ones(md.mesh.numberofvertices,1);
102 %self.albedo = 0.8*ones(md.mesh.numberofvertices,1);
103 %self.albedo_snow= 0.5*ones(md.mesh.numberofvertices,1);
104 self.hice = zeros(md.mesh.numberofvertices,1);
105 self.hsnow = 5*ones(md.mesh.numberofvertices,1);
106 end % }}}
107 function self = setdefaultparameters(self) % {{{
108
109 % albedo parameters
110 self.albedo_scheme = 0;
111 self.alb_smax = 0.79;
112 self.alb_smin = 0.6;
113 self.albi = 0.41;
114 self.albl = 0.07;
115
116 % albedo parameters for?
117 % for slater
118 self.tmin = 263.15;
119 self.tmax = 273.15;
120 % for isba & denby
121 self.mcrit = 6e-8;
122 % for isba
123 self.tau_a = 0.008;
124 self.tau_f = 0.24;
125 self.wcrit = 15.0;
126 % for alex
127 self.tmid = 273.35;
128 self.afac = -0.18;
129
130 self.hcrit = 0.028;% from Krapp et al. (2017)
131 self.rcrit = 0.85; % from Krapp et al. (2017)
132
133 self.desfac = -log(2.0)/1000;
134 self.desfacElevation = 2000;
135 self.rlaps = 7.4;
136 self.rdl = 0.29;
137
138 self.ismethod = 0;
139 self.requested_outputs={'default'};
140 end % }}}
141 function md = checkconsistency(self,md,solution,analyses) % {{{
142
143 if ismember('MasstransportAnalysis',analyses),
144 md = checkfield(md,'fieldname','smb.desfac','<=',1,'numel',1);
145 md = checkfield(md,'fieldname','smb.s0gcm','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
146 md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',1);
147 md = checkfield(md,'fieldname','smb.rdl','>=',0,'numel',1);
148 md = checkfield(md,'fieldname','smb.dailysnowfall','timeseries',1,'NaN',1,'Inf',1);
149 md = checkfield(md,'fieldname','smb.dailyrainfall','timeseries',1,'NaN',1,'Inf',1);
150 md = checkfield(md,'fieldname','smb.dailydsradiation','timeseries',1,'NaN',1,'Inf',1);
151 md = checkfield(md,'fieldname','smb.dailydlradiation','timeseries',1,'NaN',1,'Inf',1);
152 md = checkfield(md,'fieldname','smb.dailywindspeed','timeseries',1,'NaN',1,'Inf',1);
153 md = checkfield(md,'fieldname','smb.dailypressure','timeseries',1,'NaN',1,'Inf',1);
154 md = checkfield(md,'fieldname','smb.dailyairdensity','timeseries',1,'NaN',1,'Inf',1);
155 md = checkfield(md,'fieldname','smb.dailyairhumidity','timeseries',1,'NaN',1,'Inf',1);
156 md = checkfield(md,'fieldname','smb.dailytemperature','timeseries',1,'NaN',1,'Inf',1);
157
158 % TODO: transient model should be merged with SEMIC developed by Ruckamp et al. (2018)
159
160 md = checkfield(md,'fieldname','smb.ismethod','numel',1,'values',[0,1]);
161 if self.ismethod == 1 % transient mode
162 md = checkfield(md,'fieldname','smb.desfacElevation','>=',0,'numel',1);
163
164 md = checkfield(md,'fieldname','smb.albedo_scheme','NaN',1,'Inf',1,'numel',1,'values',[0,1,2,3,4]);
165 md = checkfield(md,'fieldname','smb.alb_smax','>=',0,'NaN',1,'Inf',1,'numel',1);
166 md = checkfield(md,'fieldname','smb.mask','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1],'values',[0,1,2]);
167
168 % initial values
169 md = checkfield(md,'fieldname','smb.albedo','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
170 md = checkfield(md,'fieldname','smb.albedo_snow','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
171 md = checkfield(md,'fieldname','smb.alb_smax','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
172 md = checkfield(md,'fieldname','smb.alb_smin','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
173 md = checkfield(md,'fieldname','smb.albi','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
174 md = checkfield(md,'fieldname','smb.albl','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
175 md = checkfield(md,'fieldname','smb.hice','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
176 md = checkfield(md,'fieldname','smb.hsnow','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
177 end
178 end
179 md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
180 md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]);
181 md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
182
183 end % }}}
184 function disp(self) % {{{
185 disp(sprintf(' surface forcings parameters:'));
186
187 disp(sprintf(' Interface for coupling GCM data to the energy balance model SEMIC (Krapp et al (2017) https://doi.org/10.5194/tc-11-1519-2017).'));
188 disp(sprintf(' The implemented coupling uses daily mean GCM input to calculate yearly mean smb, accumulation, ablation, and surface temperature.'));
189 disp(sprintf(' smb and temperatures are updated every year'));
190 disp(sprintf('\n SEMIC parameters:'));
191 fielddisplay(self,'dailysnowfall','daily surface dailysnowfall [m/s]');
192 fielddisplay(self,'dailyrainfall','daily surface dailyrainfall [m/s]');
193 fielddisplay(self,'dailydsradiation','daily downwelling shortwave radiation [W/m2]');
194 fielddisplay(self,'dailydlradiation','daily downwelling longwave radiation [W/m2]');
195 fielddisplay(self,'dailywindspeed','daily surface wind speed [m/s]');
196 fielddisplay(self,'dailypressure','daily surface pressure [Pa]');
197 fielddisplay(self,'dailyairdensity','daily air density [kg/m3]');
198 fielddisplay(self,'dailyairhumidity','daily air specific humidity [kg/kg]');
199 fielddisplay(self,'dailytemperature','daily surface air temperature [K]');
200 fielddisplay(self,'rlaps','present day lapse rate (default is 7.4 [degree/km]; Erokhina et al. 2017)');
201 fielddisplay(self,'desfac','desertification elevation factor (default is -log(2.0)/1000 [1/km]; Vizcaino et al. 2010)');
202 fielddisplay(self,'rdl','longwave downward radiation decrease (default is 0.29 [W/m^2/km]; Marty et al. 2002)');
203 fielddisplay(self,'s0gcm','GCM reference elevation; (default is 0) [m]');
204
205 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 % transient mode
207 fielddisplay(self,'desfacElevation','desertification elevation (default is 2000 m; Vizcaino et al. 2010)');
208 fielddisplay(self,'Tamp','amplitude of diurnal cycle [K]');
209 fielddisplay(self,'albedo','initial albedo [no unit]');
210 fielddisplay(self,'albedo_snow','initial albedo for snow [no unit]');
211 fielddisplay(self,'hice','initial thickness of ice [unit: m]');
212 fielddisplay(self,'hsnow','initial thickness of snow [unit: m]');
213 fielddisplay(self,'mask','masking for albedo. 0: ocean, 1: land, 2: ice (default: 2)');
214 fielddisplay(self,'hcrit','critical snow height for albedo [unit: m]');
215 fielddisplay(self,'rcrit','critical refreezing height for albedo [no unit]');
216
217 disp(sprintf('\nSEMIC albedo parameters.'));
218 fielddisplay(self,'albedo_scheme','albedo scheme for SEMIC. 0: none, 1: slater, 2: isba, 3: denby, 4: alex (default is 0)');
219 fielddisplay(self,'alb_smax','maximum snow albedo (default: 0.79)');
220 fielddisplay(self,'alb_smin','minimum snow albedo (default: 0.6)');
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.
225 if self.albedo_scheme == 1
226 disp(sprintf('\n\tSEMIC snow albedo parameters for Slater et al, (1998).'));
227 disp(sprintf('\t alb = alb_smax - (alb_smax - alb_smin)*tm^(3.0)'))
228 disp(sprintf('\t tm = 1 (tsurf > 273.15 K)'));
229 disp(sprintf('\t tm = f*(tsurf-tmin) (tmin <= tsurf < 273.15)'));
230 disp(sprintf('\t 0 (tsurf < tmin)'));
231 disp(sprintf('\t f = 1/(273.15-tmin)'));
232 fielddisplay(self,'tmin','minimum temperature for which albedo decline become effective. (default: 263.15 K)[unit: K])');
233 fielddisplay(self,'tmax','maxmium temperature for which albedo decline become effective. This value should be fixed. (default: 273.15 K)[unit: K])');
234 elseif self.albedo_scheme == 2,
235 disp(sprintf('\n\tSEMIC snow albedo parameters for ISBA.? where is citation?'));
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]');
238 fielddisplay(self,'tau_a','dry albedo decline [unit: 1/day]');
239 fielddisplay(self,'tau_f','wet albedo decline [unit: 1/day]');
240 elseif self.albedo_scheme == 3,
241 disp(sprintf('\n\tSEMIC snow albedo parameters for Denby et al. (2002 Tellus)'));
242 fielddisplay(self,'mcrit','critical melt rate (default: 6e-8) [unit: m/sec]');
243 elseif self.albedo_scheme == 4,
244 disp(sprintf('\n\tSEMIC snow albedo parameters for Alex.?'));
245 fielddisplay(self,'afac','[unit: ?]');
246 fielddisplay(self,'tmid','[unit: ?]');
247 else
248 error(sprintf('ERROR: %d is not supported albedom scheme.',self.albedo_scheme))
249 end
250
251 fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
252 fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
253 disp(sprintf('%51s 0: Arithmetic (default)',' '));
254 disp(sprintf('%51s 1: Geometric',' '));
255 disp(sprintf('%51s 2: Harmonic',' '));
256 fielddisplay(self,'requested_outputs','additional outputs requested');
257 end % }}}
258 function marshall(self,prefix,md,fid) % {{{
259
260 yts=md.constants.yts;
261
262 WriteData(fid,prefix,'name','md.smb.model','data',12,'format','Integer');
263
264 WriteData(fid,prefix,'object',self,'class','smb','fieldname','ismethod','format','Integer','values',[0, 1]);
265 WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfac','format','Double');
266 WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfacElevation','format','Double');
267 WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0gcm','format','DoubleMat','mattype',1);
268 WriteData(fid,prefix,'object',self,'class','smb','fieldname','rlaps','format','Double');
269 WriteData(fid,prefix,'object',self,'class','smb','fieldname','rdl','format','Double');
270 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailysnowfall','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
271 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyrainfall','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
272 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailydsradiation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
273 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailydlradiation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
274 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailywindspeed','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
275 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailypressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
276 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyairdensity','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
277 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyairhumidity','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
278 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailytemperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
279 % TODO: transient mode should be merged with SEMIC developed by Ruckamp et al. (2018).
280 if self.ismethod == 1% transient mode
281 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tamp','format','DoubleMat','mattype',1);
282 WriteData(fid,prefix,'object',self,'class','smb','fieldname','mask','format','DoubleMat','mattype',1);
283 WriteData(fid,prefix,'object',self,'class','smb','fieldname','hice','format','DoubleMat','mattype',1);
284 WriteData(fid,prefix,'object',self,'class','smb','fieldname','hsnow','format','DoubleMat','mattype',1);
285
286 WriteData(fid,prefix,'object',self,'class','smb','fieldname','hcrit','format','Double');
287 WriteData(fid,prefix,'object',self,'class','smb','fieldname','rcrit','format','Double');
288
289 %albedo
290 WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo','format','DoubleMat','mattype',1);
291 WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo_snow','format','DoubleMat','mattype',1);
292 WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo_scheme','format','Integer');
293 WriteData(fid,prefix,'object',self,'class','smb','fieldname','albi','format','Double');
294 WriteData(fid,prefix,'object',self,'class','smb','fieldname','albl','format','Double');
295 WriteData(fid,prefix,'object',self,'class','smb','fieldname','alb_smin','format','Double');
296 WriteData(fid,prefix,'object',self,'class','smb','fieldname','alb_smax','format','Double');
297
298 %albedo parameters for ?
299 %for slater
300 WriteData(fid,prefix,'object',self,'class','smb','fieldname','tmin','format','Double');
301 WriteData(fid,prefix,'object',self,'class','smb','fieldname','tmax','format','Double');
302 %for isba & denby
303 WriteData(fid,prefix,'object',self,'class','smb','fieldname','mcrit','format','Double');
304 %for isba
305 WriteData(fid,prefix,'object',self,'class','smb','fieldname','wcrit','format','Double');
306 WriteData(fid,prefix,'object',self,'class','smb','fieldname','tau_a','format','Double');
307 WriteData(fid,prefix,'object',self,'class','smb','fieldname','tau_f','format','Double');
308 %for alex
309 WriteData(fid,prefix,'object',self,'class','smb','fieldname','tmid','format','Double');
310 WriteData(fid,prefix,'object',self,'class','smb','fieldname','afac','format','Double');
311 end
312
313 WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer');
314 WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer');
315
316 %process requested outputs
317 outputs = self.requested_outputs;
318 pos = find(ismember(outputs,'default'));
319 if ~isempty(pos),
320 outputs(pos) = []; %remove 'default' from outputs
321 outputs = [outputs defaultoutputs(self,md)]; %add defaults
322 end
323 WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray');
324
325 end % }}}
326 end
327end
Note: See TracBrowser for help on using the repository browser.