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

Last change on this file since 27498 was 27498, checked in by inwoo, 2 years ago

CHG: update SEMIC - enable transient mode.
File size: 13.3 KB
RevLine 
[23540]1%SMBsemic Class definition
2%
3% Usage:
4% SMBsemic=SMBsemic();
5
6classdef SMBsemic
[24240]7 properties (SetAccess=public)
[23540]8 dailysnowfall = NaN;
9 dailyrainfall = NaN;
10 dailydsradiation = NaN;
11 dailydlradiation = NaN;
12 dailywindspeed = NaN;
13 dailypressure = NaN;
[27453]14 dailyairdensity = NaN;
[23540]15 dailyairhumidity = NaN;
16 dailytemperature = NaN;
[27498]17 Tamp = NaN;
18 mask = NaN;
19 hice = NaN;
20 hsnow = NaN;
[23540]21 desfac = 0;
[27498]22 desfacElevation = 0;
[27453]23 rlaps = 0;
[23540]24 rdl = 0;
[27453]25 s0gcm = NaN;
26 steps_per_step = 1;
27 averaging = 0;
28 requested_outputs = {};
[27498]29
30 hcrit = 0;
31 rcrit = 0;
32
33 % albedo
34 albedo = 0;
35 albedo_snow = 0;
36 albedo_scheme = 0;
37 alb_smax = NaN;
38 alb_smin = NaN;
39 albi = NaN;
40 albl = NaN;
41
42 % method
43 ismethod = 0;
[23540]44 end
45 methods
46 function self = SMBsemic(varargin) % {{{
47 switch nargin
48 case 0
49 self=setdefaultparameters(self);
50 otherwise
51 error('constructor not supported');
52 end
53 end % }}}
54 function self = extrude(self,md) % {{{
55 self.dailysnowfall=project3d(md,'vector',self.dailysnowfall,'type','node');
56 self.dailyrainfall=project3d(md,'vector',self.dailyrainfall,'type','node');
57 self.dailydsradiation=project3d(md,'vector',self.dailydsradiation,'type','node');
58 self.dailydlradiation=project3d(md,'vector',self.dailydlradiation,'type','node');
59 self.dailywindspeed=project3d(md,'vector',self.dailywindspeed,'type','node');
60 self.dailypressure=project3d(md,'vector',self.dailypressure,'type','node');
61 self.dailyairdensity=project3d(md,'vector',self.dailyairdensity,'type','node');
62 self.dailyairhumidity=project3d(md,'vector',self.dailyairhumidity,'type','node');
63 self.dailytemperature=project3d(md,'vector',self.dailytemperature,'type','node');
64 self.s0gcm=project3d(md,'vector',self.s0gcm,'type','node');
65
66 end % }}}
67 function list = defaultoutputs(self,md) % {{{
[27403]68 list = {'SmbMassBalance'};
[23540]69 end % }}}
70 function self = initialize(self,md) % {{{
71
72 if isnan(self.s0gcm),
73 self.s0gcm=zeros(md.mesh.numberofvertices,1);
74 disp(' no SMBsemic.s0gcm specified: values set as zero');
75 end
[27498]76 self.Tamp = 3*ones(md.mesh.numberofvertices,1);
77 self.albedo = 0.8*ones(md.mesh.numberofvertices,1);
78 self.albedo_snow= 0.5*ones(md.mesh.numberofvertices,1);
79 self.hice = zeros(md.mesh.numberofvertices,1);
80 self.hsnow = 5*ones(md.mesh.numberofvertices,1);
[23540]81 end % }}}
82 function self = setdefaultparameters(self) % {{{
83
[27498]84 self.albedo_scheme = 0;
85 self.alb_smax = 0.79;
86 self.alb_smin = 0.6;
87 self.albi = 0.41;
88 self.albl = 0.07;
89
90 self.hcrit = 0.028;% from Krapp et al. (2017)
91 self.rcrit = 0.85; % from Krapp et al. (2017)
92
93 self.desfac = -log(2.0)/1000;
94 self.desfacElevation = 2000;
95 self.rlaps = 7.4;
96 self.rdl = 0.29;
97
98 self.ismethod = 0;
[27409]99 self.requested_outputs={'default'};
[23540]100 end % }}}zo
101 function md = checkconsistency(self,md,solution,analyses) % {{{
102
103 if ismember('MasstransportAnalysis',analyses),
104 md = checkfield(md,'fieldname','smb.desfac','<=',1,'numel',1);
105 md = checkfield(md,'fieldname','smb.s0gcm','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
106 md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',1);
107 md = checkfield(md,'fieldname','smb.rdl','>=',0,'numel',1);
108 md = checkfield(md,'fieldname','smb.dailysnowfall','timeseries',1,'NaN',1,'Inf',1);
109 md = checkfield(md,'fieldname','smb.dailyrainfall','timeseries',1,'NaN',1,'Inf',1);
110 md = checkfield(md,'fieldname','smb.dailydsradiation','timeseries',1,'NaN',1,'Inf',1);
111 md = checkfield(md,'fieldname','smb.dailydlradiation','timeseries',1,'NaN',1,'Inf',1);
112 md = checkfield(md,'fieldname','smb.dailywindspeed','timeseries',1,'NaN',1,'Inf',1);
113 md = checkfield(md,'fieldname','smb.dailypressure','timeseries',1,'NaN',1,'Inf',1);
114 md = checkfield(md,'fieldname','smb.dailyairdensity','timeseries',1,'NaN',1,'Inf',1);
115 md = checkfield(md,'fieldname','smb.dailyairhumidity','timeseries',1,'NaN',1,'Inf',1);
116 md = checkfield(md,'fieldname','smb.dailytemperature','timeseries',1,'NaN',1,'Inf',1);
[27498]117
118 % TODO: transient model should be merged with SEMIC developed by Ruckamp et al. (2018)
119
120 md = checkfield(md,'fieldname','smb.ismethod','numel',1,'values',[0,1]);
121 if self.ismethod == 1 % transient mode.
122 md = checkfield(md,'fieldname','smb.desfacElevation','>=',0,'numel',1);
123
124 md = checkfield(md,'fieldname','smb.albedo_scheme','NaN',1,'Inf',1,'numel',1,'values',[0,1,2,3,4]);
125 md = checkfield(md,'fieldname','smb.alb_smax','>=',0,'NaN',1,'Inf',1,'numel',1);
126 md = checkfield(md,'fieldname','smb.mask','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1],'values',[0,1,2]);
127
128 % initial values
129 md = checkfield(md,'fieldname','smb.albedo','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
130 md = checkfield(md,'fieldname','smb.albedo_snow','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
131 md = checkfield(md,'fieldname','smb.alb_smax','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
132 md = checkfield(md,'fieldname','smb.alb_smin','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
133 md = checkfield(md,'fieldname','smb.albi','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
134 md = checkfield(md,'fieldname','smb.albl','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
135 md = checkfield(md,'fieldname','smb.hice','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
136 md = checkfield(md,'fieldname','smb.hsnow','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
137 end
[23540]138 end
[24240]139 md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
[24800]140 md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]);
[23540]141 md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
142
143 end % }}}
144 function disp(self) % {{{
145 disp(sprintf(' surface forcings parameters:'));
[24240]146
[23540]147 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).'));
[24240]148 disp(sprintf(' The implemented coupling uses daily mean GCM input to calculate yearly mean smb, accumulation, ablation, and surface temperature.'));
[23540]149 disp(sprintf(' smb and temperatures are updated every year'));
150 disp(sprintf('\n SEMIC parameters:'));
[24240]151 fielddisplay(self,'dailysnowfall','daily surface dailysnowfall [m/s]');
[23540]152 fielddisplay(self,'dailyrainfall','daily surface dailyrainfall [m/s]');
153 fielddisplay(self,'dailydsradiation','daily downwelling shortwave radiation [W/m2]');
154 fielddisplay(self,'dailydlradiation','daily downwelling longwave radiation [W/m2]');
[24240]155 fielddisplay(self,'dailywindspeed','daily surface wind speed [m/s]');
[23540]156 fielddisplay(self,'dailypressure','daily surface pressure [Pa]');
[24240]157 fielddisplay(self,'dailyairdensity','daily air density [kg/m3]');
[23540]158 fielddisplay(self,'dailyairhumidity','daily air specific humidity [kg/kg]');
159 fielddisplay(self,'dailytemperature','daily surface air temperature [K]');
[27454]160 fielddisplay(self,'rlaps','present day lapse rate (default is 7.4 [degree/km]; Erokhina et al. 2017)');
[23540]161 fielddisplay(self,'desfac','desertification elevation factor (default is -log(2.0)/1000 [1/km]; Vizcaino et al. 2010)');
162 fielddisplay(self,'rdl','longwave downward radiation decrease (default is 0.29 [W/m^2/km]; Marty et al. 2002)');
163 fielddisplay(self,'s0gcm','GCM reference elevation; (default is 0) [m]');
[27498]164
165 fielddisplay(self,'ismethod','method for calculating SMB with SEMIC. Default version of SEMIC is really slow. 0: steady, 1: transient (default: 0)');
166 if self.ismethod == 1 % tranisent mode
167 fielddisplay(self,'desfacElevation','desertification elevation (default is 2000 m; Vizcaino et al. 2010)');
168 fielddisplay(self,'Tamp','amplitude of diurnal cycle [K]');
169 fielddisplay(self,'albedo','initial albedo [no unit]');
170 fielddisplay(self,'albedo_snow','initial albedo for snow [no unit]');
171 fielddisplay(self,'hice','initial thickness of ice [unit: m]');
172 fielddisplay(self,'hsnow','initial thickness of snow [unit: m]');
173 fielddisplay(self,'mask','masking for albedo. 0: ocean, 1: land, 2: ice (default: 2)');
174 fielddisplay(self,'hcrit','critical snow height for albedo [unit: m]');
175 fielddisplay(self,'rcrit','critical refreezing height for albedo [no unit]');
176
177 disp(sprintf('\nSEMIC albedo parameters.'));
178 fielddisplay(self,'albedo_scheme','albedo scheme for SEMIC. 0: none, 1: slater, 2: isba, 3: denby, 4: alex (default is 0)');
179 fielddisplay(self,'abl_smax','maximum snow albedo (default: 0.79)');
180 fielddisplay(self,'abl_smin','minimum snow albedo (default: 0.6)');
181 fielddisplay(self,'albi','background albeod for bare ice (default: 0.41)');
182 fielddisplay(self,'albl','background albeod for bare land (default: 0.07)');
183 end
184
[24240]185 fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
[24793]186 fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
[24806]187 disp(sprintf('%51s 0: Arithmetic (default)',' '));
188 disp(sprintf('%51s 1: Geometric',' '));
189 disp(sprintf('%51s 2: Harmonic',' '));
[23540]190 fielddisplay(self,'requested_outputs','additional outputs requested');
191 end % }}}
192 function marshall(self,prefix,md,fid) % {{{
[24240]193
[27453]194 yts=md.constants.yts;
195
[23540]196 WriteData(fid,prefix,'name','md.smb.model','data',12,'format','Integer');
197
[27498]198 WriteData(fid,prefix,'object',self,'class','smb','fieldname','ismethod','format','Integer','values',[0, 1]);
[23540]199 WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfac','format','Double');
[27498]200 WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfacElevation','format','Double');
[23540]201 WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0gcm','format','DoubleMat','mattype',1);
202 WriteData(fid,prefix,'object',self,'class','smb','fieldname','rlaps','format','Double');
203 WriteData(fid,prefix,'object',self,'class','smb','fieldname','rdl','format','Double');
[27453]204 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailysnowfall','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
205 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyrainfall','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
206 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailydsradiation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
207 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailydlradiation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
208 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailywindspeed','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
209 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailypressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
210 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyairdensity','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
211 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyairhumidity','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
212 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailytemperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
[27498]213 % TODO: transient mode should be merged with SEMIC developed by Ruckamp et al. (2018).
214 if self.ismethod == 1% transient mode.
215 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tamp','format','DoubleMat','mattype',1);
216 WriteData(fid,prefix,'object',self,'class','smb','fieldname','mask','format','DoubleMat','mattype',1);
217 WriteData(fid,prefix,'object',self,'class','smb','fieldname','hice','format','DoubleMat','mattype',1);
218 WriteData(fid,prefix,'object',self,'class','smb','fieldname','hsnow','format','DoubleMat','mattype',1);
219
220 WriteData(fid,prefix,'object',self,'class','smb','fieldname','hcrit','format','Double');
221 WriteData(fid,prefix,'object',self,'class','smb','fieldname','rcrit','format','Double');
222
223 %albedo...
224 WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo','format','DoubleMat','mattype',1);
225 WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo_snow','format','DoubleMat','mattype',1);
226 WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo_scheme','format','Integer');
227 WriteData(fid,prefix,'object',self,'class','smb','fieldname','albi','format','Double');
228 WriteData(fid,prefix,'object',self,'class','smb','fieldname','albl','format','Double');
229 WriteData(fid,prefix,'object',self,'class','smb','fieldname','alb_smin','format','Double');
230 WriteData(fid,prefix,'object',self,'class','smb','fieldname','alb_smax','format','Double');
231 end
232
[27453]233 WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer');
234 WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer');
[27498]235
[23540]236 %process requested outputs
237 outputs = self.requested_outputs;
238 pos = find(ismember(outputs,'default'));
239 if ~isempty(pos),
240 outputs(pos) = []; %remove 'default' from outputs
241 outputs = [outputs defaultoutputs(self,md)]; %add defaults
242 end
243 WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray');
244
245 end % }}}
246 end
247end
Note: See TracBrowser for help on using the repository browser.