source: issm/oecreview/Archive/25834-26739/ISSM-26058-26059.diff@ 26740

Last change on this file since 26740 was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 119.2 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m
2===================================================================
3--- ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m (revision 26058)
4+++ ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m (revision 26059)
5@@ -30,10 +30,10 @@
6
7 %Initialize surface and basal forcings
8 md.smb = initialize(md.smb,md);
9-md.basalforcings = initialize(md.basalforcings,md);
10+md.basalforcings = initialize(md.basalforcings,md);
11
12 %Initialize ocean forcings and sealevel
13-md.dsl= initialize(md.dsl,md);
14+md.dsl = initialize(md.dsl,md);
15
16 %Deal with other boundary conditions
17 if isnan(md.balancethickness.thickening_rate),
18Index: ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py
19===================================================================
20--- ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py (revision 26058)
21+++ ../trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py (revision 26059)
22@@ -2,13 +2,12 @@
23
24
25 def SetIceSheetBC(md):
26- """
27- SETICESHEETBC - Create the boundary conditions for stressbalance and thermal models for an IceSheet with no Ice Front
28+ """SETICESHEETBC - Create the boundary conditions for stressbalance and thermal models for an IceSheet with no Ice Front
29
30- Usage:
31- md = SetIceSheetBC(md)
32+ Usage:
33+ md = SetIceSheetBC(md)
34
35- See also: SETICESHELFBC, SETMARINEICESHEETBC
36+ See also: SETICESHELFBC, SETMARINEICESHEETBC
37 """
38
39 #node on Dirichlet
40@@ -32,10 +31,13 @@
41
42 #No ice front -> do nothing
43
44- #Create zeros basalforcings and smb
45+ #Initialize surface and basal forcings
46 md.smb.initialize(md)
47 md.basalforcings.initialize(md)
48
49+ #Initialize ocean forcings and sealevel
50+ md.dsl.initialize(md)
51+
52 #Deal with other boundary conditions
53 if np.all(np.isnan(md.balancethickness.thickening_rate)):
54 md.balancethickness.thickening_rate = np.zeros((md.mesh.numberofvertices))
55Index: ../trunk-jpl/src/m/classes/hydrologydc.py
56===================================================================
57--- ../trunk-jpl/src/m/classes/hydrologydc.py (revision 26058)
58+++ ../trunk-jpl/src/m/classes/hydrologydc.py (revision 26059)
59@@ -6,11 +6,10 @@
60
61
62 class hydrologydc(object):
63- """
64- Hydrologydc class definition
65+ """HYDROLOGYDC class definition
66
67 Usage:
68- hydrologydc = hydrologydc()
69+ hydrologydc = hydrologydc()
70 """
71
72 def __init__(self): # {{{
73@@ -48,11 +47,14 @@
74 self.epl_conductivity = 0
75 self.eplflip_lock = 0
76
77- #set defaults
78 self.setdefaultparameters()
79 # }}}
80
81 def __repr__(self): # {{{
82+ # TODO:
83+ # - Convert all formatting to calls to <string>.format (see any
84+ # already converted <class>.__repr__ method for examples)
85+ #
86 string = ' hydrology Dual Porous Continuum Equivalent parameters:'
87 string = ' - general parameters'
88 string = "%s\n%s" % (string, fielddisplay(self, 'water_compressibility', 'compressibility of water [Pa^ - 1]'))
89Index: ../trunk-jpl/src/m/classes/hydrologytws.m
90===================================================================
91--- ../trunk-jpl/src/m/classes/hydrologytws.m (revision 26058)
92+++ ../trunk-jpl/src/m/classes/hydrologytws.m (revision 26059)
93@@ -23,7 +23,7 @@
94 end % }}}
95 function list = defaultoutputs(self,md) % {{{
96 list = {''};
97- end % }}}
98+ end % }}}
99
100 function self = setdefaultparameters(self) % {{{
101 self.requested_outputs = {'default'};
102@@ -45,16 +45,16 @@
103 function marshall(self,prefix,md,fid) % {{{
104 WriteData(fid,prefix,'name','md.hydrology.model','data',2,'format','Integer');
105 WriteData(fid,prefix,'object',self,'fieldname','spcwatercolumn','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
106- outputs = self.requested_outputs;
107- pos = find(ismember(outputs,'default'));
108- if ~isempty(pos),
109- outputs(pos) = []; %remove 'default' from outputs
110- outputs = [outputs defaultoutputs(self,md)]; %add defaults
111- end
112- WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
113+ outputs = self.requested_outputs;
114+ pos = find(ismember(outputs,'default'));
115+ if ~isempty(pos),
116+ outputs(pos) = []; %remove 'default' from outputs
117+ outputs = [outputs defaultoutputs(self,md)]; %add defaults
118+ end
119+ WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
120 end % }}}
121 function savemodeljs(self,fid,modelname) % {{{
122-
123+
124 writejs1Darray(fid,[modelname '.hydrology.spcwatercolumn'],self.spcwatercolumn);
125
126 end % }}}
127Index: ../trunk-jpl/src/m/classes/hydrologytws.py
128===================================================================
129--- ../trunk-jpl/src/m/classes/hydrologytws.py (nonexistent)
130+++ ../trunk-jpl/src/m/classes/hydrologytws.py (revision 26059)
131@@ -0,0 +1,64 @@
132+import numpy as np
133+
134+from structtoobj import structtoobj
135+
136+class hydrologytws(object):
137+ """HYDROLOGYTWS class definition
138+
139+ Usage:
140+ hydrologytws = hydrologytws()
141+ """
142+
143+ def __init__(self): # {{{
144+ self.spcwatercolumn = np.nan
145+ self.requested_outputs = np.nan
146+
147+ nargs = len(args)
148+ if nargs == 0:
149+ self.setdefaultparameters()
150+ elif nargs == 1:
151+ self = structtoobj(self, args[0])
152+ else:
153+ raise RuntimeError('constructor not supported')
154+ #}}}
155+
156+ def __repr__(self): # {{{
157+ s = ' hydrologytws solution parameters:\n'
158+ s += '{}\n'.format(fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
159+ s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
160+ return s
161+ #}}}
162+
163+ def defaultoutputs(self, md): # {{{
164+ return ['']
165+ #}}}
166+
167+ def setdefaultparameters(self): # {{{
168+ self.requested_outputs = ['defualt']
169+ return self
170+ #}}}
171+
172+ def extrude(self, md): # {{{
173+ return self
174+ #}}}
175+
176+ def checkconsistency(self, md, solution, analyses): # {{{
177+ # Early return
178+ if 'HydrologyTwsAnalysis' not in analyses:
179+ return
180+ md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1)
181+ #}}}
182+
183+ def marshall(self, prefix, md, fid): # {{{
184+ WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer')
185+ WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
186+ outputs = self.requested_outputs
187+ pos = find(ismember(outputs,'default'))
188+ if not len(pos),
189+ outputs[pos] = []; # remove 'default' from outputs
190+ outputs.extend(defaultoutputs(self, md)) # add defaults
191+ end
192+ WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray')
193+ # }}}
194+
195+
196Index: ../trunk-jpl/src/m/classes/matenhancedice.m
197===================================================================
198--- ../trunk-jpl/src/m/classes/matenhancedice.m (revision 26058)
199+++ ../trunk-jpl/src/m/classes/matenhancedice.m (revision 26059)
200@@ -5,26 +5,26 @@
201
202 classdef matenhancedice
203 properties (SetAccess=public)
204- rho_ice = 0.;
205- rho_water = 0.;
206- rho_freshwater = 0.;
207- mu_water = 0.;
208- heatcapacity = 0.;
209- latentheat = 0.;
210- thermalconductivity = 0.;
211- temperateiceconductivity = 0.;
212+ rho_ice = 0.;
213+ rho_water = 0.;
214+ rho_freshwater = 0.;
215+ mu_water = 0.;
216+ heatcapacity = 0.;
217+ latentheat = 0.;
218+ thermalconductivity = 0.;
219+ temperateiceconductivity = 0.;
220 effectiveconductivity_averaging = 0.;
221- meltingpoint = 0.;
222- beta = 0.;
223- mixed_layer_capacity = 0.;
224- thermal_exchange_velocity = 0.;
225- rheology_E = NaN;
226- rheology_B = NaN;
227- rheology_n = NaN;
228- rheology_law = '';
229+ meltingpoint = 0.;
230+ beta = 0.;
231+ mixed_layer_capacity = 0.;
232+ thermal_exchange_velocity = 0.;
233+ rheology_E = NaN;
234+ rheology_B = NaN;
235+ rheology_n = NaN;
236+ rheology_law = '';
237
238- %slr
239- earth_density = 0;
240+ %SLC
241+ earth_density = 0;
242
243 end
244 methods
245Index: ../trunk-jpl/src/m/classes/matestar.m
246===================================================================
247--- ../trunk-jpl/src/m/classes/matestar.m (revision 26058)
248+++ ../trunk-jpl/src/m/classes/matestar.m (revision 26059)
249@@ -5,26 +5,26 @@
250
251 classdef matestar
252 properties (SetAccess=public)
253- rho_ice = 0.;
254- rho_water = 0.;
255- rho_freshwater = 0.;
256- mu_water = 0.;
257- heatcapacity = 0.;
258- latentheat = 0.;
259- thermalconductivity = 0.;
260- temperateiceconductivity = 0.;
261+ rho_ice = 0.;
262+ rho_water = 0.;
263+ rho_freshwater = 0.;
264+ mu_water = 0.;
265+ heatcapacity = 0.;
266+ latentheat = 0.;
267+ thermalconductivity = 0.;
268+ temperateiceconductivity = 0.;
269 effectiveconductivity_averaging = 0;
270- meltingpoint = 0.;
271- beta = 0.;
272- mixed_layer_capacity = 0.;
273- thermal_exchange_velocity = 0.;
274- rheology_B = NaN;
275- rheology_Ec = NaN;
276- rheology_Es = NaN;
277- rheology_law = '';
278+ meltingpoint = 0.;
279+ beta = 0.;
280+ mixed_layer_capacity = 0.;
281+ thermal_exchange_velocity = 0.;
282+ rheology_B = NaN;
283+ rheology_Ec = NaN;
284+ rheology_Es = NaN;
285+ rheology_law = '';
286
287 %slc
288- earth_density = 0;
289+ earth_density = 0;
290
291 end
292 methods
293Index: ../trunk-jpl/src/m/classes/matice.py
294===================================================================
295--- ../trunk-jpl/src/m/classes/matice.py (revision 26058)
296+++ ../trunk-jpl/src/m/classes/matice.py (revision 26059)
297@@ -7,12 +7,11 @@
298
299
300 class matice(object):
301- '''
302- MATICE class definition
303+ """MATICE class definition
304
305 Usage:
306 matice = matice()
307- '''
308+ """
309
310 def __init__(self): #{{{
311 self.rho_ice = 0.
312@@ -32,13 +31,7 @@
313 self.rheology_n = np.nan
314 self.rheology_law = ''
315
316- #giaivins
317- self.lithosphere_shear_modulus = 0.
318- self.lithosphere_density = 0.
319- self.mantle_shear_modulus = 0.
320- self.mantle_density = 0.
321-
322- #slr
323+ #slc
324 self.earth_density = 0
325 self.setdefaultparameters()
326 #}}}
327@@ -61,10 +54,6 @@
328 s = "%s\n%s" % (s, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]"))
329 s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
330 s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', or 'NyeH2O'"))
331- s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
332- s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]"))
333- s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
334- s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]"))
335 s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]"))
336
337 return s
338@@ -107,11 +96,9 @@
339 #available: none, paterson and arrhenius
340 self.rheology_law = 'Paterson'
341
342- # GIA:
343- self.lithosphere_shear_modulus = 6.7e10 # (Pa)
344- self.lithosphere_density = 3.32 # (g/cm^-3)
345- self.mantle_shear_modulus = 1.45e11 # (Pa)
346- self.mantle_density = 3.34 # (g/cm^-3)
347+ # Rheology for ice
348+ self.rheology_B = 2.1 * 1e8
349+ self.rheology_n = 3
350
351 # SLR
352 self.earth_density = 5512 # average density of the Earth, (kg/m^3)
353@@ -118,25 +105,18 @@
354 #}}}
355
356 def checkconsistency(self, md, solution, analyses): #{{{
357- if solution != 'SealevelriseSolution':
358+ if solution == 'TransientSolution' and md.transient.isslc:
359+ md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
360+ else:
361 md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
362 md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0)
363 md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0)
364 md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0)
365 md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'universal', 1, 'NaN', 1, 'Inf', 1)
366- md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements])
367+ md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'universal',1, 'NaN', 1, 'Inf', 1)
368 md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O'])
369- md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
370+ md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
371
372- if 'GiaAnalysis' in analyses:
373- md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', [1])
374- md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', [1])
375- md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', [1])
376- md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
377-
378- if 'SealevelriseAnalysis' in analyses:
379- md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
380-
381 return md
382 #}}}
383
384@@ -158,9 +138,5 @@
385 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
386 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
387 WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
388- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
389- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.)
390- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
391- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.)
392 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
393 #}}}
394Index: ../trunk-jpl/src/m/classes/surfaceload.m
395===================================================================
396--- ../trunk-jpl/src/m/classes/surfaceload.m (revision 26058)
397+++ ../trunk-jpl/src/m/classes/surfaceload.m (revision 26059)
398@@ -19,13 +19,13 @@
399 end
400 end % }}}
401 function self = setdefaultparameters(self) % {{{
402-
403+
404 icethicknesschange=[];
405 waterheightchange=[];
406 otherchange=[];
407-
408+
409 end % }}}
410- function md = checkconsistency(self,md,solution,analyses) % {{{
411+ function md = checkconsistency(self,md,solution,analyses) % {{{
412
413 if ~ismember('SealevelchangeAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslc==0),
414 return;
415Index: ../trunk-jpl/src/m/classes/transient.m
416===================================================================
417--- ../trunk-jpl/src/m/classes/transient.m (revision 26058)
418+++ ../trunk-jpl/src/m/classes/transient.m (revision 26059)
419@@ -7,7 +7,7 @@
420 properties (SetAccess=public)
421 issmb = 0;
422 ismasstransport = 0;
423- isoceantransport = 0;
424+ isoceantransport = 0;
425 isstressbalance = 0;
426 isthermal = 0;
427 isgroundingline = 0;
428@@ -15,7 +15,7 @@
429 isdamageevolution = 0;
430 ismovingfront = 0;
431 ishydrology = 0;
432- issampling = 0;
433+ issampling = 0;
434 isslc = 0;
435 amr_frequency = 0;
436 isoceancoupling = 0;
437@@ -43,7 +43,7 @@
438 self.isdamageevolution = 0;
439 self.ismovingfront =0;
440 self.ishydrology = 0;
441- self.issampling = 0;
442+ self.issampling = 0;
443 self.isslc = 0;
444 self.isoceancoupling = 0;
445 self.amr_frequency = 0;
446@@ -64,7 +64,7 @@
447 self.isdamageevolution = 0;
448 self.ismovingfront = 0;
449 self.ishydrology = 0;
450- self.issampling = 0;
451+ self.issampling = 0;
452 self.isslc = 0;
453 self.isoceancoupling = 0;
454 self.amr_frequency = 0;
455@@ -97,7 +97,7 @@
456 md = checkfield(md,'fieldname','transient.requested_outputs','stringrow',1);
457 md = checkfield(md,'fieldname','transient.isslc','numel',[1],'values',[0 1]);
458 md = checkfield(md,'fieldname','transient.isoceancoupling','numel',[1],'values',[0 1]);
459- md = checkfield(md,'fieldname','transient.issampling','numel',[1],'values',[0 1]);
460+ md = checkfield(md,'fieldname','transient.issampling','numel',[1],'values',[0 1]);
461 md = checkfield(md,'fieldname','transient.amr_frequency','numel',[1],'>=',0,'NaN',1,'Inf',1);
462
463 if (~strcmp(solution,'TransientSolution') & md.transient.iscoupling==1),
464@@ -120,7 +120,7 @@
465 fielddisplay(self,'isdamageevolution','indicates whether damage evolution is used in the transient');
466 fielddisplay(self,'ismovingfront','indicates whether a moving front capability is used in the transient');
467 fielddisplay(self,'ishydrology','indicates whether an hydrology model is used');
468- fielddisplay(self,'issampling','indicates whether sampling is used in the transient')
469+ fielddisplay(self,'issampling','indicates whether sampling is used in the transient')
470 fielddisplay(self,'isslc','indicates whether a sea-level change solution is used in the transient');
471 fielddisplay(self,'isoceancoupling','indicates whether a coupling with an ocean model is used in the transient');
472 fielddisplay(self,'amr_frequency','frequency at which mesh is refined in simulations with multiple time_steps');
473@@ -138,7 +138,7 @@
474 WriteData(fid,prefix,'object',self,'fieldname','isdamageevolution','format','Boolean');
475 WriteData(fid,prefix,'object',self,'fieldname','ishydrology','format','Boolean');
476 WriteData(fid,prefix,'object',self,'fieldname','ismovingfront','format','Boolean');
477- WriteData(fid,prefix,'object',self,'fieldname','issampling','format','Boolean');
478+ WriteData(fid,prefix,'object',self,'fieldname','issampling','format','Boolean');
479 WriteData(fid,prefix,'object',self,'fieldname','isslc','format','Boolean');
480 WriteData(fid,prefix,'object',self,'fieldname','isoceancoupling','format','Boolean');
481 WriteData(fid,prefix,'object',self,'fieldname','amr_frequency','format','Integer');
482@@ -164,7 +164,7 @@
483 writejsdouble(fid,[modelname '.trans.isdamageevolution'],self.isdamageevolution);
484 writejsdouble(fid,[modelname '.trans.ismovingfront'],self.ismovingfront);
485 writejsdouble(fid,[modelname '.trans.ishydrology'],self.ishydrology);
486- writejsdouble(fid,[modelname '.trans.issampling'],self.issampling);
487+ writejsdouble(fid,[modelname '.trans.issampling'],self.issampling);
488 writejsdouble(fid,[modelname '.trans.isslc'],self.isslc);
489 writejsdouble(fid,[modelname '.trans.isoceancoupling'],self.isoceancoupling);
490 writejsdouble(fid,[modelname '.trans.amr_frequency'],self.amr_frequency);
491Index: ../trunk-jpl/src/m/consistency/ismodelselfconsistent.m
492===================================================================
493--- ../trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 26058)
494+++ ../trunk-jpl/src/m/consistency/ismodelselfconsistent.m (revision 26059)
495@@ -2,7 +2,7 @@
496 %ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
497 %
498 % Usage:
499-% ismodelselfconsistent(md),
500+% ismodelselfconsistent(md);
501
502 %initialize consistency as true
503 md.private.isconsistent=true;
504Index: ../trunk-jpl/src/m/plot/radarpower.py
505===================================================================
506--- ../trunk-jpl/src/m/plot/radarpower.py (nonexistent)
507+++ ../trunk-jpl/src/m/plot/radarpower.py (revision 26059)
508@@ -0,0 +1,16 @@
509+def solveslm(slm, solutionstringi, *args):
510+ """RADARPOWER - overlay a power radar image on an existing mesh
511+
512+ This routine will overlay a power radar image on an existing mesh.
513+ The power amplitude will be output to vel for now.
514+ In the future, think about a field to hold this value.
515+
516+ Usage:
517+ md=radarpower(md,options);
518+ md=radarpower(md)
519+
520+ TODO:
521+ - Translate from MATLAB API as we bring Python plotting capabilities online
522+ """
523+
524+ return md
525Index: ../trunk-jpl/src/m/solve/solveslm.m
526===================================================================
527--- ../trunk-jpl/src/m/solve/solveslm.m (revision 26058)
528+++ ../trunk-jpl/src/m/solve/solveslm.m (revision 26059)
529@@ -1,5 +1,5 @@
530 function slm=solveslm(slm,solutionstringi,varargin)
531-%SOLVESLR - apply solution sequence for this sealevel model
532+%SOLVESLM - apply solution sequence for this sealevel model
533 %
534 % Usage:
535 % slm=solve(slm,solutionstring,varargin)
536@@ -42,7 +42,7 @@
537 slm.private.solution=solutionstring;
538 cluster=slm.cluster;
539 batch=0;
540-%now, go through icecaps, glacies and earth, and upload all the data independently:
541+%now, go through icecaps, glaciers and earth, and upload all the data independently:
542 disp('solving ice caps first');
543 for i=1:length(slm.icecaps),
544 slm.icecaps{i}=solve(slm.icecaps{i},solutionstringi,'batch','yes');
545@@ -50,7 +50,7 @@
546 disp('solving earth now');
547 slm.earth=solve(slm.earth,solutionstringi,'batch','yes');
548
549-%Firs, build a runtime name that is unique
550+%First, build a runtime name that is unique
551 c=clock;
552 slm.private.runtimename=sprintf('%s-%02i-%02i-%04i-%02i-%02i-%02i-%i',slm.miscellaneous.name,c(2),c(3),c(1),c(4),c(5),floor(c(6)),feature('GetPid'));
553
554Index: ../trunk-jpl/src/m/solve/solveslm.py
555===================================================================
556--- ../trunk-jpl/src/m/solve/solveslm.py (nonexistent)
557+++ ../trunk-jpl/src/m/solve/solveslm.py (revision 26059)
558@@ -0,0 +1,95 @@
559+from datetime import datetime
560+import os
561+
562+import numpy as np
563+
564+from loadresultsfromcluster import loadresultsfromcluster
565+from pairoptions import pairoptions
566+from waitonlock import waitonlock
567+
568+
569+def solveslm(slm, solutionstringi, *args):
570+ """SOLVESLM - apply solution sequence for this sealevel model
571+
572+ Usage:
573+ slm=solve(slm,solutionstring,varargin)
574+ where varargin is a lit of paired arguments of string OR enums
575+
576+ solution types available comprise:
577+ - 'Transient'
578+
579+ extra options:
580+
581+ Examples:
582+ slm=solve(slm,'Transient');
583+ """
584+
585+ # Recover and process solve options
586+ if solutionstringi.lower() == 'tr' or solutionstringi.lower() == 'transient':
587+ solutionstring = 'TransientSolution'
588+ else:
589+ raise RuntimeError('solutionstring {} not supported!'.format(solutionstringi))
590+
591+ # Default settings for debugging
592+ valgrind = 0
593+ #slm.cluster.interactive = 0
594+ #valgrind = 1
595+
596+ # Check consistency
597+ slm.checkconsistency(solutionstring)
598+
599+ # Process options
600+ options = pairoptions('solutionstring', solutionstring, *args)
601+
602+ # Make sure we request sum of cluster processors
603+ totalnp = 0
604+ for i in range(len(slm.icecaps)):
605+ totalnp = totalnp + slm.icecaps[i].cluster.np
606+ totalnp = totalnp + slm.earth.cluster.np
607+ if totalnp != slm.cluster.np:
608+ raise RuntimeError('sum of all icecaps and earch cluster processors requestes should be equal to slm.cluster.np')
609+
610+ # Recover some fields
611+ slm.private.solution = solutionstring
612+ cluster = slm.cluster
613+ batch = 0
614+ # Now, go through icecaps, glaciers and earth, and upload all the data independently
615+ print('solving ice caps first')
616+ for i in range(len(slm.icecaps)):
617+ slm.icecaps[i] = solve(slm.icecaps[i], solutionastringi,'batch','yes')
618+ print('solving earth now')
619+ slm.earth = solve(slm.earth, solutionstringi, 'batch', 'yes')
620+
621+ # First, build a runtime name that is unique
622+ c = datetime.now()
623+ md.private.runtimename = "%s-%02i-%02i-%04i-%02i-%02i-%02i-%i" % (md.miscellaneous.name, c.month, c.day, c.year, c.hour, c.minute, c.second, os.getpid())
624+
625+ # Write all input files
626+ privateruntimenames = []
627+ miscellaneousnames = []
628+ nps = []
629+ for i in range(len(slm.icecaps)):
630+ privateruntimenames.append(slm.icecaps[i],private.runtimename)
631+ miscellaneousnames.append(slm.earth.miscellaneous.name)
632+ nps.append(slm.earth.cluster.np)
633+
634+ BuildQueueScriptMultipleModels(cluster, slm.private.runtimename, slm.miscellaneous.name, slm.private.solution, valgrind, privateruntimenames, miscellaneousnames, nps)
635+
636+ # Upload all required files, given that each individual solution for icecaps and earth model already did
637+ filelist = [slm.miscellaneous.name + '.queue']
638+ UploadQueueJob(cluster, slm.miscellaneous.name, slm.private.runtimename, filelist)
639+
640+ # Launch queue job
641+ LaunchQueueJob(cluster, slm.miscellaneous.name, slm.private.runtimename, filelist, '', batch)
642+
643+ # Wait on lock
644+ if slm.settings.waitonlock > 0:
645+ islock = waitonlock(slm)
646+ if islock == 0: # no results to be loaded
647+ print('The results must be loaded manually with md = loadresultsfromcluster(md).')
648+ else: # load results
649+ if slm.verbose.solution:
650+ print('loading results from cluster')
651+ slm = loadresultsfromcluster(slm)
652+
653+ return slm
654Index: ../trunk-jpl/src/m/classes/dsl.m
655===================================================================
656--- ../trunk-jpl/src/m/classes/dsl.m (revision 26058)
657+++ ../trunk-jpl/src/m/classes/dsl.m (revision 26059)
658@@ -6,9 +6,9 @@
659 classdef dsl
660 properties (SetAccess=public)
661
662- global_average_thermosteric_sea_level; %corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m)
663- sea_surface_height_above_geoid; %corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable (in m)
664- sea_water_pressure_at_sea_floor; %corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!)
665+ global_average_thermosteric_sea_level; %Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).
666+ sea_surface_height_above_geoid; %Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).
667+ sea_water_pressure_at_sea_floor; %Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).
668
669 end
670 methods
671@@ -50,9 +50,9 @@
672 function disp(self) % {{{
673
674 disp(sprintf(' dsl parameters:'));
675- fielddisplay(self,'global_average_thermosteric_sea_level','corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m)');
676- fielddisplay(self,'sea_surface_height_above_geoid','corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m)');
677- fielddisplay(self,'sea_water_pressure_at_sea_floor','corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent)');
678+ fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).');
679+ fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).');
680+ fielddisplay(self,'sea_water_pressure_at_sea_floor','Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).');
681
682 end % }}}
683 function marshall(self,prefix,md,fid) % {{{
684Index: ../trunk-jpl/src/m/classes/dsl.py
685===================================================================
686--- ../trunk-jpl/src/m/classes/dsl.py (revision 26058)
687+++ ../trunk-jpl/src/m/classes/dsl.py (revision 26059)
688@@ -14,24 +14,22 @@
689 """
690
691 def __init__(self): #{{{
692- self.global_average_thermosteric_sea_level_change = 0 #corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)
693- self.sea_surface_height_change_above_geoid = float('NaN') #corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)
694- self.sea_water_pressure_change_at_sea_floor = float('NaN') #corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)
695- self.compute_fingerprints = 0; #do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid
696+ self.global_average_thermosteric_sea_level = np.nan # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).
697+ self.sea_surface_height_above_geoid = np.nan # Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).
698+ self.sea_water_pressure_at_sea_floor = np.nan # Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).
699 #}}}
700
701 def __repr__(self): #{{{
702 s = ' dsl parameters:\n'
703- s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)'))
704- s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)'))
705- s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)'))
706- s += '{}\n'.format(fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
707+ s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level', 'Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).'))
708+ s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_above_geoid', 'Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).'))
709+ s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_at_sea_floor', 'Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).'))
710 return s
711 #}}}
712
713 def extrude(self, md): #{{{
714- self.sea_surface_height_change_above_geoid = project3d(md, 'vector', self.sea_surface_height_change_above_geoid, 'type', 'node')
715- self.sea_water_pressure_change_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor, 'type', 'node')
716+ self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node')
717+ self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node')
718 return self
719 #}}}
720
721@@ -41,18 +39,16 @@
722
723 def checkconsistency(self, md, solution, analyses): #{{{
724 # Early return
725- if 'SealevelriseAnalysis' not in analyses:
726+ if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
727 return md
728- if solution == 'TransientSolution' and md.transient.isslc == 0:
729- return md
730- md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level_change', 'NaN', 1, 'Inf', 1)
731- md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_change_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
732- md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_change_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
733- md = checkfield(md, 'fieldname', 'dsl.compute_fingerprints', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
734- if self.compute_fingerprints:
735- # Check if geodetic flag of slr is on
736- if not md.slr.geodetic:
737- raise RuntimeError('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on')
738+
739+ md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level', 'NaN', 1, 'Inf', 1)
740+ md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
741+ md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
742+
743+ if md.solidearth.settings.compute_bp_grd:
744+ md = checkfield(md, 'fieldname', dsl.sea_water_pressure_at_sea_floor, 'empty', 1)
745+
746 return md
747 # }}}
748
749@@ -59,8 +55,7 @@
750 def marshall(self, prefix, md, fid): #{{{
751 yts = md.constants.yts
752 WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')
753- WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'compute_fingerprints', 'format', 'Integer')
754- WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts)
755- WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts)
756- WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts)
757+ WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 2, because we are sending a GMSL value identical everywhere on each element.
758+ WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 1 because we specify DSL at vertex locations.
759+ WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts) # mattype 1 because we specify bottom pressure at vertex locations.
760 # }}}
761Index: ../trunk-jpl/src/m/classes/dslmme.m
762===================================================================
763--- ../trunk-jpl/src/m/classes/dslmme.m (revision 26058)
764+++ ../trunk-jpl/src/m/classes/dslmme.m (revision 26059)
765@@ -7,9 +7,9 @@
766 properties (SetAccess=public)
767
768 modelid; %index into the multi-model ensemble, determine which field will be used.
769- global_average_thermosteric_sea_level; %corresponds to zostoga fields in CMIP5 archives. Specified as a temporally quantity (in m) for each ensemble.
770- sea_surface_height_above_geoid; %corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m) for each ensemble
771- sea_water_pressure_at_sea_floor; %corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble
772+ global_average_thermosteric_sea_level; %Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.
773+ sea_surface_height_above_geoid; %Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m) for each ensemble.
774+ sea_water_pressure_at_sea_floor; %Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble.
775
776 end
777 methods
778@@ -52,9 +52,9 @@
779
780 disp(sprintf(' dsl mme parameters:'));
781 fielddisplay(self,'modelid','index into the multi-model ensemble, determine which field will be used.');
782- fielddisplay(self,'global_average_thermosteric_sea_level','corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.');
783- fielddisplay(self,'sea_surface_height_above_geoid','corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally quantity (in m) for each ensemble.');
784- fielddisplay(self,'sea_water_pressure_at_sea_floor','corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally quantity (in m) for each ensemble.');
785+ fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.');
786+ fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m) for each ensemble.');
787+ fielddisplay(self,'sea_water_pressure_at_sea_floor','Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble.');
788 end % }}}
789 function marshall(self,prefix,md,fid) % {{{
790
791Index: ../trunk-jpl/src/m/classes/dslmme.py
792===================================================================
793--- ../trunk-jpl/src/m/classes/dslmme.py (revision 26058)
794+++ ../trunk-jpl/src/m/classes/dslmme.py (revision 26059)
795@@ -14,15 +14,14 @@
796 """
797
798 def __init__(self, *args): #{{{
799- self.modelid = 0 #index into the multi-model ensemble
800- self.global_average_thermosteric_sea_level_change = [] #corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble.
801- self.sea_surface_height_change_above_geoid = [] #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble
802- self.sea_water_pressure_change_at_sea_floor = [] #corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr equivalent, not in Pa/yr!) for each ensemble
803- self.compute_fingerprints = 0 #corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble
804-
805- nargin = len(args)
806+ self.modelid = 0 # Index into the multi-model ensemble
807+ self.global_average_thermosteric_sea_level = [] # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.
808+ self.sea_surface_height_above_geoid = [] # Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m) for each ensemble.
809+ self.sea_water_pressure_at_sea_floor = [] #Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble.
810
811- if nargin == 0:
812+ nargs = len(args)
813+
814+ if nargs == 0:
815 self.setdefaultparameters()
816 else:
817 raise Exception('constructor not supported')
818@@ -31,10 +30,9 @@
819 def __repr__(self): # {{{
820 s = ' dsl mme parameters:\n'
821 s += '{}\n'.format(fielddisplay(self, 'modelid', 'index into the multi-model ensemble, determines which field will be used.'))
822- s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level_change', 'corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble.'))
823- s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_change_above_geoid', 'corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble.'))
824- s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_change_at_sea_floor', 'corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr) for each ensemble.'))
825- s += '{}\n'.format(fielddisplay(self, 'compute_fingerprints', 'do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'))
826+ s += '{}\n'.format(fielddisplay(self, 'global_average_thermosteric_sea_level', 'Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m) for each ensemble.'))
827+ s += '{}\n'.format(fielddisplay(self, 'sea_surface_height_above_geoid', 'Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m) for each ensemble.'))
828+ s += '{}\n'.format(fielddisplay(self, 'sea_water_pressure_at_sea_floor', 'Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!) for each ensemble.'))
829 return s
830 #}}}
831
832@@ -43,34 +41,35 @@
833 #}}}
834
835 def checkconsistency(self, md, solution, analyses): # {{{
836- if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslr):
837+ # Early return
838+ if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
839 return md
840- for i in range(len(self.global_average_thermosteric_sea_level_change)):
841- md = checkfield(md, 'field', self.global_average_thermosteric_sea_level_change[i], 'NaN', 1, 'Inf', 1)
842- md = checkfield(md, 'field', self.sea_surface_height_change_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
843- md = checkfield(md, 'field', self.sea_water_pressure_change_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
844- md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level_change))
845- if self.compute_fingerprints:
846- #check geodetic flag of slr is on
847- if not md.solidearth.settings.computesealevelchange:
848- raise Exception('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on')
849+
850+ for i in range(len(self.global_average_thermosteric_sea_level)):
851+ md = checkfield(md, 'field', self.global_average_thermosteric_sea_level[i], 'NaN', 1, 'Inf', 1)
852+ md = checkfield(md, 'field', self.sea_surface_height_above_geoid[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
853+ md = checkfield(md, 'field', self.sea_water_pressure_at_sea_floor[i], 'NaN', 1, 'Inf', 1, 'timeseries', 1)
854+ md = checkfield(md, 'field', self.modelid, 'NaN', 1, 'Inf', 1, '>=', 1, '<=',len(self.global_average_thermosteric_sea_level))
855+
856+ if self.solidearth.settings.compute_bp_grd:
857+ md = checkfield(md, 'field', self.sea_water_pressure_at_sea_floor, 'empty', 1)
858+
859 return md
860 #}}}
861
862 def marshall(self, prefix, md, fid): #{{{
863 WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 2, 'format', 'Integer')
864- WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_fingerprints', 'format', 'Integer')
865 WriteData(fid, prefix, 'object', self, 'fieldname', 'modelid', 'format', 'Double')
866- WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level_change), 'format', 'Integer')
867- WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level_change', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
868- WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_change_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3)
869- WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_change_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
870+ WriteData(fid, prefix, 'name', 'md.dsl.nummodels', 'data', len(self.global_average_thermosteric_sea_level), 'format', 'Integer')
871+ WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'MatArray', 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
872+ WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3)
873+ WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_above_geoid', 'format', 'MatArray', 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts, 'scale', 1e-3 / md.constants.yts)
874 #}}}
875
876 def extrude(self, md): #{{{
877- for i in range(len(self.global_average_thermosteric_sea_level_change)):
878- self.sea_surface_height_change_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_change_above_geoid[i], 'type', 'node', 'layer', 1)
879- self.sea_water_pressure_change_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_change_at_sea_floor[i], 'type', 'node', 'layer', 1)
880+ for i in range(len(self.global_average_thermosteric_sea_level)):
881+ self.sea_surface_height_above_geoid[i] = project3d(md, 'vector', self.self.sea_surface_height_above_geoid[i], 'type', 'node', 'layer', 1)
882+ self.sea_water_pressure_at_sea_floor[i] = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor[i], 'type', 'node', 'layer', 1)
883
884 return self
885 #}}}
886Index: ../trunk-jpl/src/m/classes/fourierlove.m
887===================================================================
888--- ../trunk-jpl/src/m/classes/fourierlove.m (revision 26058)
889+++ ../trunk-jpl/src/m/classes/fourierlove.m (revision 26059)
890@@ -79,12 +79,8 @@
891 end
892
893 %need 'litho' material:
894- if ~isa(md.materials,'materials')
895+ if ~isa(md.materials,'materials') | ~sum(strcmpi(md.materials.nature,'litho'))
896 error('Need a ''litho'' material to run a Fourier Love number analysis');
897- else
898- if ~sum(strcmpi(md.materials.nature,'litho')),
899- error('Need a ''litho'' material to run a Fourier Love number analysis');
900- end
901 end
902
903 end % }}}
904Index: ../trunk-jpl/src/m/classes/fourierlove.py
905===================================================================
906--- ../trunk-jpl/src/m/classes/fourierlove.py (revision 26058)
907+++ ../trunk-jpl/src/m/classes/fourierlove.py (revision 26059)
908@@ -4,11 +4,10 @@
909
910
911 class fourierlove(object):
912- """
913- Fourier Love Number class definition
914+ """FOURIERLOVE - Fourier Love Number class definition
915
916- Usage:
917- fourierlove = fourierlove()
918+ Usage:
919+ fourierlove = fourierlove()
920 """
921 def __init__(self): # {{{
922 self.nfreq = 0
923@@ -22,11 +21,14 @@
924 self.love_kernels = 0
925 self.forcing_type = 0
926
927- #set defaults
928 self.setdefaultparameters()
929 #}}}
930
931 def __repr__(self): # {{{
932+ # TODO:
933+ # - Convert all formatting to calls to <string>.format (see any
934+ # already converted <class>.__repr__ method for examples)
935+ #
936 string = ' Fourier Love class:'
937 string = "%s\n%s" % (string, fielddisplay(self, 'nfreq', 'number of frequencies sampled (default 1, elastic) [Hz]'))
938 string = "%s\n%s" % (string, fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'))
939@@ -75,6 +77,8 @@
940 #}}}
941
942 def checkconsistency(self, md, solution, analyses): # {{{
943+ if 'LoveAnalysis' not in analyses:
944+ return md
945 md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
946 md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', [md.love.nfreq])
947 md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
948@@ -88,6 +92,9 @@
949 if md.love.sh_nmin <= 1 and md.love.forcing_type == 9:
950 raise RuntimeError("Degree 1 not supported for Volumetric Potential forcing. Use sh_min >= 2 for this kind of calculation.")
951
952+ # need 'litho' material
953+ if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature:
954+ raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis')
955 return md
956 # }}}
957
958Index: ../trunk-jpl/src/m/classes/geometry.m
959===================================================================
960--- ../trunk-jpl/src/m/classes/geometry.m (revision 26058)
961+++ ../trunk-jpl/src/m/classes/geometry.m (revision 26059)
962@@ -84,10 +84,10 @@
963
964 end % }}}
965 function marshall(self,prefix,md,fid) % {{{
966-
967- if (size(self.thickness==md.mesh.numberofvertices) | (self.thickness==md.mesh.numberofvertices+1)),
968+
969+ if (size(self.thickness)==md.mesh.numberofvertices) | (size(self.thickness)==md.mesh.numberofvertices+1)),
970 WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
971- elseif (size(self.thickness==md.mesh.numberofelements) | (self.thickness==md.mesh.numberofelements+1)),
972+ elseif (size(self.thickness)==md.mesh.numberofelements) | (size(self.thickness)==md.mesh.numberofelements+1)),
973 WriteData(fid,prefix,'object',self,'fieldname','thickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
974 else
975 error('geometry thickness time series should be a vertex or element time series');
976Index: ../trunk-jpl/src/m/classes/geometry.py
977===================================================================
978--- ../trunk-jpl/src/m/classes/geometry.py (revision 26058)
979+++ ../trunk-jpl/src/m/classes/geometry.py (revision 26059)
980@@ -50,12 +50,8 @@
981 #}}}
982
983 def checkconsistency(self, md, solution, analyses): #{{{
984- if (solution == 'TransientSolution' and md.transient.isgia) or (solution == 'GiaSolution'):
985- md = checkfield(md, 'fieldname', 'geometry.thickness', 'timeseries', 1, 'NaN', 1, 'Inf', 1)
986- elif solution == 'SealevelriseSolution':
987- md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
988- elif solution == 'LoveSolution':
989- return
990+ if solution == 'LoveSolution':
991+ return md
992 else:
993 md = checkfield(md, 'fieldname', 'geometry.surface', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
994 md = checkfield(md, 'fieldname', 'geometry.base', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
995@@ -70,13 +66,18 @@
996 if np.any(np.abs(self.bed[pos] - self.base[pos]) > 10**-9):
997 md.checkmessage('equality base = bed on grounded ice violated')
998 md = checkfield(md, 'fieldname', 'geometry.bed', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
999-
1000 return md
1001 # }}}
1002
1003 def marshall(self, prefix, md, fid): #{{{
1004+ if (len(self.thickness) == md.mesh.numberofvertices) or (len(self.thickness) == md.mesh.numberofvertices + 1):
1005+ WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
1006+ elif (len(self.thickness) == md.mesh.numberofelements) or (len(self.thickness) == md.mesh.numberofelements + 1):
1007+ WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
1008+ else:
1009+ raise RuntimeError('geometry thickness time series should be a vertex or element time series')
1010+
1011 WriteData(fid, prefix, 'object', self, 'fieldname', 'surface', 'format', 'DoubleMat', 'mattype', 1)
1012- WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
1013 WriteData(fid, prefix, 'object', self, 'fieldname', 'base', 'format', 'DoubleMat', 'mattype', 1)
1014 WriteData(fid, prefix, 'object', self, 'fieldname', 'bed', 'format', 'DoubleMat', 'mattype', 1)
1015 WriteData(fid, prefix, 'object', self, 'fieldname', 'hydrostatic_ratio', 'format', 'DoubleMat', 'mattype', 1)
1016Index: ../trunk-jpl/src/m/classes/hydrologyshreve.m
1017===================================================================
1018--- ../trunk-jpl/src/m/classes/hydrologyshreve.m (revision 26058)
1019+++ ../trunk-jpl/src/m/classes/hydrologyshreve.m (revision 26059)
1020@@ -30,7 +30,7 @@
1021
1022 %Type of stabilization to use 0:nothing 1:artificial_diffusivity
1023 self.stabilization = 1;
1024- self.requested_outputs = {'default'};
1025+ self.requested_outputs = {'default'};
1026 end % }}}
1027 function md = checkconsistency(self,md,solution,analyses) % {{{
1028
1029Index: ../trunk-jpl/src/m/classes/hydrologyshreve.py
1030===================================================================
1031--- ../trunk-jpl/src/m/classes/hydrologyshreve.py (revision 26058)
1032+++ ../trunk-jpl/src/m/classes/hydrologyshreve.py (revision 26059)
1033@@ -4,11 +4,10 @@
1034
1035
1036 class hydrologyshreve(object):
1037- """
1038- HYDROLOGYSHREVE class definition
1039+ """HYDROLOGYSHREVE class definition
1040
1041- Usage:
1042- hydrologyshreve = hydrologyshreve()
1043+ Usage:
1044+ hydrologyshreve = hydrologyshreve()
1045 """
1046
1047 def __init__(self): # {{{
1048@@ -15,13 +14,15 @@
1049 self.spcwatercolumn = float('NaN')
1050 self.stabilization = 0
1051 self.requested_outputs = []
1052- #set defaults
1053+
1054 self.setdefaultparameters()
1055-
1056 #}}}
1057
1058 def __repr__(self): # {{{
1059-
1060+ # TODO:
1061+ # - Convert all formatting to calls to <string>.format (see any
1062+ # already converted <class>.__repr__ method for examples)
1063+ #
1064 string = ' hydrologyshreve solution parameters:'
1065 string = "%s\n%s" % (string, fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
1066 string = "%s\n%s" % (string, fielddisplay(self, 'stabilization', 'artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.'))
1067@@ -47,7 +48,7 @@
1068
1069 def checkconsistency(self, md, solution, analyses): # {{{
1070 #Early return
1071- if 'HydrologyShreveAnalysis' not in analyses:
1072+ if 'HydrologyShreveAnalysis' not in analyses or (solution == 'TransientSolution' and not md.transient.ishydrology):
1073 return md
1074
1075 md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1)
1076@@ -60,7 +61,7 @@
1077 WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer')
1078 WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
1079 WriteData(fid, prefix, 'object', self, 'fieldname', 'stabilization', 'format', 'Double')
1080- #process requested outputs
1081+ #process requested outputs
1082 outputs = self.requested_outputs
1083 indices = [i for i, x in enumerate(outputs) if x == 'default']
1084 if len(indices) > 0:
1085Index: ../trunk-jpl/src/m/classes/initialization.m
1086===================================================================
1087--- ../trunk-jpl/src/m/classes/initialization.m (revision 26058)
1088+++ ../trunk-jpl/src/m/classes/initialization.m (revision 26059)
1089@@ -21,7 +21,7 @@
1090 channelarea = NaN;
1091 sealevel = NaN;
1092 bottompressure = NaN;
1093- sample = NaN;
1094+ sample = NaN;
1095 end
1096 methods
1097 function self = extrude(self,md) % {{{
1098@@ -51,7 +51,6 @@
1099 end
1100 end % }}}
1101 function self = setdefaultparameters(self) % {{{
1102-
1103 end % }}}
1104 function md = checkconsistency(self,md,solution,analyses) % {{{
1105 if ismember('StressbalanceAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isstressbalance == 0),
1106@@ -97,12 +96,9 @@
1107 end
1108 if ismember('HydrologyShreveAnalysis',analyses),
1109 if isa(md.hydrology,'hydrologyshreve'),
1110- if strcmp(solution,'TransientSolution') & md.transient.ishydrology,
1111+ if (strcmp(solution,'TransientSolution') & md.transient.ishydrology) | strcmp(solution,'HydrologySolution'),
1112 md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
1113 end
1114- if strcmp(solution,'HydrologySolution'),
1115- md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
1116- end
1117 end
1118 end
1119 if ismember('HydrologyTwsAnalysis',analyses),
1120@@ -110,7 +106,7 @@
1121 md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
1122 end
1123 end
1124- if ismember('SealevelchangeAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isslc == 0),
1125+ if ismember('SealevelchangeAnalysis',analyses),
1126 if strcmp(solution,'TransientSolution') & md.transient.isslc,
1127 md = checkfield(md,'fieldname','initialization.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
1128 end
1129@@ -133,13 +129,13 @@
1130 md = checkfield(md,'fieldname','initialization.epl_head','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
1131 md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
1132 end
1133- end
1134- end
1135- if ismember('SamplingAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.issampling == 0),
1136- if ~isnan(md.initialization.sample)
1137- md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
1138- end
1139+ end
1140 end
1141+ if ismember('SamplingAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.issampling == 0),
1142+ if ~isnan(md.initialization.sample)
1143+ md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
1144+ end
1145+ end
1146 end % }}}
1147 function disp(self) % {{{
1148 disp(sprintf(' initial field values:'));
1149@@ -158,7 +154,7 @@
1150 fielddisplay(self,'watercolumn','subglacial water sheet thickness (for Shreve and GlaDS) [m]');
1151 fielddisplay(self,'hydraulic_potential','Hydraulic potential (for GlaDS) [Pa]');
1152 fielddisplay(self,'channelarea','subglacial water channel area (for GlaDS) [m2]');
1153- fielddisplay(self,'sample','Realization of a Gaussian random field');
1154+ fielddisplay(self,'sample','Realization of a Gaussian random field');
1155 end % }}}
1156 function marshall(self,prefix,md,fid) % {{{
1157
1158@@ -178,8 +174,8 @@
1159 WriteData(fid,prefix,'object',self,'fieldname','watercolumn','format','DoubleMat','mattype',1);
1160 WriteData(fid,prefix,'object',self,'fieldname','channelarea','format','DoubleMat','mattype',1);
1161 WriteData(fid,prefix,'object',self,'fieldname','hydraulic_potential','format','DoubleMat','mattype',1);
1162- WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1);
1163-
1164+ WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1);
1165+
1166 if md.thermal.isenthalpy,
1167 if numel(self.enthalpy) <= 1,
1168 %reconstruct enthalpy
1169@@ -194,7 +190,7 @@
1170 end
1171 end % }}}
1172 function savemodeljs(self,fid,modelname) % {{{
1173-
1174+
1175 writejs1Darray(fid,[modelname '.initialization.vx'],self.vx);
1176 writejs1Darray(fid,[modelname '.initialization.vy'],self.vy);
1177 writejs1Darray(fid,[modelname '.initialization.vz'],self.vz);
1178@@ -209,8 +205,8 @@
1179 writejs1Darray(fid,[modelname '.initialization.watercolumn'],self.watercolumn);
1180 writejs1Darray(fid,[modelname '.initialization.hydraulic_potential'],self.hydraulic_potential);
1181 writejs1Darray(fid,[modelname '.initialization.channel'],self.channelarea);
1182- writejs1Darray(fid,[modelname '.initialization.sample'],self.sample);
1183-
1184+ writejs1Darray(fid,[modelname '.initialization.sample'],self.sample);
1185+
1186 end % }}}
1187 end
1188 end
1189Index: ../trunk-jpl/src/m/classes/initialization.py
1190===================================================================
1191--- ../trunk-jpl/src/m/classes/initialization.py (revision 26058)
1192+++ ../trunk-jpl/src/m/classes/initialization.py (revision 26059)
1193@@ -28,6 +28,8 @@
1194 self.epl_thickness = np.nan
1195 self.hydraulic_potential = np.nan
1196 self.channelarea = np.nan
1197+ self.sealevel = np.nan
1198+ self.bottompressure = np.nan
1199 self.sample = np.nan
1200
1201 #set defaults
1202@@ -66,6 +68,8 @@
1203 self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1)
1204 self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1)
1205 self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1)
1206+ self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1)
1207+ self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1)
1208
1209 #Lithostatic pressure by default
1210 if np.ndim(md.geometry.surface) == 2:
1211@@ -89,6 +93,9 @@
1212 if 'MasstransportAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.ismasstransport:
1213 md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1214 md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1215+ if 'OceantransportAnalysis' in analyses:
1216+ if solution == 'TransientSolution' and md.transient.isslc and md.transient.isoceantransport:
1217+ md = checkfield(md, 'fieldname', 'initialization.bottompressure', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1218 if 'BalancethicknessAnalysis' in analyses and solution == 'BalancethicknessSolution':
1219 md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1220 md = checkfield(md, 'fieldname', 'initialization.vy', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1221@@ -111,15 +118,14 @@
1222 md = checkfield(md, 'fieldname', 'delta Tpmp', 'field', np.absolute(md.initialization.temperature[pos] - (md.materials.meltingpoint - md.materials.beta * md.initialization.pressure[pos])), '<', 1e-11, 'message', 'set temperature to pressure melting point at locations with waterfraction > 0')
1223 if 'HydrologyShreveAnalysis' in analyses:
1224 if hasattr(md.hydrology, 'hydrologyshreve'):
1225+ if (solution == 'TransientSolution' and md.transient.ishydrology) or solution == 'HydrologySolution':
1226+ md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1227+ if 'HydrologyTwsAnalysis' in analyses:
1228+ if hasattr(md.hydrology, 'hydrologytws'):
1229 md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1230- if 'HydrologyDCInefficientAnalysis' in analyses:
1231- if hasattr(md.hydrology, 'hydrologydc'):
1232- md = checkfield(md, 'fieldname', 'initialization.sediment_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1233- if 'HydrologyDCEfficientAnalysis' in analyses:
1234- if hasattr(md.hydrology, 'hydrologydc'):
1235- if md.hydrology.isefficientlayer == 1:
1236- md = checkfield(md, 'fieldname', 'initialization.epl_head', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1237- md = checkfield(md, 'fieldname', 'initialization.epl_thickness', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1238+ if 'SealevelchangeAnalysis' in analyses:
1239+ if solution == 'TransientSolution' and md.transient.isslc:
1240+ md = checkfield(md, 'fieldname', 'initialization.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1241 if 'HydrologyGlaDSAnalysis' in analyses:
1242 if hasattr(md.hydrology, 'hydrologyglads'):
1243 md = checkfield(md, 'fieldname', 'initialization.watercolumn', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
1244@@ -137,6 +143,8 @@
1245 WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
1246 WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
1247 WriteData(fid, prefix, 'object', self, 'fieldname', 'pressure', 'format', 'DoubleMat', 'mattype', 1)
1248+ WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1)
1249+ WriteData(fid, prefix, 'object', self, 'fieldname', 'bottompressure', 'format', 'DoubleMat', 'mattype', 1)
1250 WriteData(fid, prefix, 'object', self, 'fieldname', 'temperature', 'format', 'DoubleMat', 'mattype', 1)
1251 WriteData(fid, prefix, 'object', self, 'fieldname', 'waterfraction', 'format', 'DoubleMat', 'mattype', 1)
1252 WriteData(fid, prefix, 'object', self, 'fieldname', 'sediment_head', 'format', 'DoubleMat', 'mattype', 1)
1253Index: ../trunk-jpl/src/m/consistency/ismodelselfconsistent.py
1254===================================================================
1255--- ../trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 26058)
1256+++ ../trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 26059)
1257@@ -1,10 +1,9 @@
1258 def ismodelselfconsistent(md): #{{{
1259- '''
1260- ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
1261+ """ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
1262
1263- Usage:
1264- ismodelselfconsistent(md),
1265- '''
1266+ Usage:
1267+ ismodelselfconsistent(md)
1268+ """
1269
1270 #initialize consistency as true
1271 md.private.isconsistent = True
1272@@ -52,6 +51,8 @@
1273 analyses = ['EnthalpyAnalysis', 'ThermalAnalysis', 'MeltingAnalysis']
1274 elif solutiontype == 'MasstransportSolution':
1275 analyses = ['MasstransportAnalysis']
1276+ elif solutiontype == 'OceantransportSolution':
1277+ analyses = ['OceantransportAnalysis']
1278 elif solutiontype == 'BalancethicknessSolution':
1279 analyses = ['BalancethicknessAnalysis']
1280 elif solutiontype == 'Balancethickness2Solution':
1281@@ -71,11 +72,11 @@
1282 elif solutiontype == 'EsaSolution':
1283 analyses = ['EsaAnalysis']
1284 elif solutiontype == 'TransientSolution':
1285- analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis']
1286+ analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'OceantransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyShreveAnalysis', 'HydrologyTwsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis']
1287 elif solutiontype == 'SealevelriseSolution':
1288 analyses = ['SealevelriseAnalysis']
1289 elif solutiontype == 'HydrologySolution':
1290- analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis']
1291+ analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'HydrologyGladsAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyTwsAnalysis']
1292 elif 'DamageEvolutionSolution':
1293 analyses = ['DamageEvolutionAnalysis']
1294 elif 'SamplingSolution':
1295Index: ../trunk-jpl/src/m/solve/solve.m
1296===================================================================
1297--- ../trunk-jpl/src/m/solve/solve.m (revision 26058)
1298+++ ../trunk-jpl/src/m/solve/solve.m (revision 26059)
1299@@ -9,7 +9,7 @@
1300 % Solution types available comprise:
1301 % - 'Stressbalance' or 'sb'
1302 % - 'Masstransport' or 'mt'
1303-% - 'Oceantransport' or 'oceant'
1304+% - 'Oceantransport' or 'oceant'
1305 % - 'Thermal' or 'th'
1306 % - 'Steadystate' or 'ss'
1307 % - 'Transient' or 'tr'
1308Index: ../trunk-jpl/test/NightlyRun/test2001.m
1309===================================================================
1310--- ../trunk-jpl/test/NightlyRun/test2001.m (revision 26058)
1311+++ ../trunk-jpl/test/NightlyRun/test2001.m (revision 26059)
1312@@ -21,8 +21,8 @@
1313 %Loading history
1314 md.timestepping.start_time=-2400000; %4,800 kyr :: EVALUATION TIME
1315 md.timestepping.time_step= 2400000; %2,400 kyr :: EVALUATION TIME
1316-% to get rid of default final_time: make sure final_time>start_time
1317-md.timestepping.final_time=2400000; %2,500 kyr
1318+% to get rid of default final_time, make sure final_time > start_time
1319+md.timestepping.final_time=2400000; %2,400 kyr
1320 md.masstransport.spcthickness=[...
1321 [md.geometry.thickness; 0],...
1322 [md.geometry.thickness; 2400000]...
1323@@ -45,7 +45,7 @@
1324 %md.cluster=generic('name',oshostname(),'np',3);
1325 md.verbose=verbose('11111111111');
1326 md.verbose.solver=0;
1327-md=solve(md,'tr');
1328+md=solve(md,'Transient');
1329
1330 %Fields and tolerances to track changes
1331 field_names ={'UGrd'};
1332Index: ../trunk-jpl/src/m/classes/matdamageice.m
1333===================================================================
1334--- ../trunk-jpl/src/m/classes/matdamageice.m (revision 26058)
1335+++ ../trunk-jpl/src/m/classes/matdamageice.m (revision 26059)
1336@@ -5,25 +5,25 @@
1337
1338 classdef matdamageice
1339 properties (SetAccess=public)
1340- rho_ice = 0.;
1341- rho_water = 0.;
1342- rho_freshwater = 0.;
1343- mu_water = 0.;
1344- heatcapacity = 0.;
1345- latentheat = 0.;
1346- thermalconductivity = 0.;
1347- temperateiceconductivity = 0.;
1348+ rho_ice = 0.;
1349+ rho_water = 0.;
1350+ rho_freshwater = 0.;
1351+ mu_water = 0.;
1352+ heatcapacity = 0.;
1353+ latentheat = 0.;
1354+ thermalconductivity = 0.;
1355+ temperateiceconductivity = 0.;
1356 effectiveconductivity_averaging = 0.;
1357- meltingpoint = 0.;
1358- beta = 0.;
1359- mixed_layer_capacity = 0.;
1360- thermal_exchange_velocity = 0.;
1361- rheology_B = NaN;
1362- rheology_n = NaN;
1363- rheology_law = '';
1364+ meltingpoint = 0.;
1365+ beta = 0.;
1366+ mixed_layer_capacity = 0.;
1367+ thermal_exchange_velocity = 0.;
1368+ rheology_B = NaN;
1369+ rheology_n = NaN;
1370+ rheology_law = '';
1371
1372 %slc
1373- earth_density = 0;
1374+ earth_density = 0;
1375
1376 end
1377 methods
1378Index: ../trunk-jpl/src/m/classes/matdamageice.py
1379===================================================================
1380--- ../trunk-jpl/src/m/classes/matdamageice.py (revision 26058)
1381+++ ../trunk-jpl/src/m/classes/matdamageice.py (revision 26059)
1382@@ -30,18 +30,17 @@
1383 self.rheology_n = float('NaN')
1384 self.rheology_law = ''
1385
1386- #giaivins:
1387- self.lithosphere_shear_modulus = 0.
1388- self.lithosphere_density = 0.
1389- self.mantle_shear_modulus = 0.
1390- self.mantle_density = 0.
1391+ #SLC
1392+ self.earth_density = 5512 # average density of the Earth, (kg / m^3)
1393
1394- #SLR
1395- self.earth_density = 5512 # average density of the Earth, (kg / m^3)
1396 self.setdefaultparameters()
1397 #}}}
1398
1399 def __repr__(self): # {{{
1400+ # TODO:
1401+ # - Convert all formatting to calls to <string>.format (see any
1402+ # already converted <class>.__repr__ method for examples)
1403+ #
1404 string = " Materials:"
1405 string = "%s\n%s" % (string, fielddisplay(self, "rho_ice", "ice density [kg / m^3]"))
1406 string = "%s\n%s" % (string, fielddisplay(self, "rho_water", "water density [kg / m^3]"))
1407@@ -59,10 +58,6 @@
1408 string = "%s\n%s" % (string, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1 / n)]"))
1409 string = "%s\n%s" % (string, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
1410 string = "%s\n%s" % (string, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'"))
1411- string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
1412- string = "%s\n%s" % (string, fielddisplay(self, "lithosphere_density", "Lithosphere density [g / cm^ - 3]"))
1413- string = "%s\n%s" % (string, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
1414- string = "%s\n%s" % (string, fielddisplay(self, "mantle_density", "Mantle density [g / cm^ - 3]"))
1415 string = "%s\n%s" % (string, fielddisplay(self, "earth_density", "Mantle density [kg / m^ - 3]"))
1416 return string
1417 #}}}
1418@@ -104,13 +99,7 @@
1419 #available: none, paterson and arrhenius
1420 self.rheology_law = 'Paterson'
1421
1422- # GIA:
1423- self.lithosphere_shear_modulus = 6.7e10 # (Pa)
1424- self.lithosphere_density = 3.32 # (g / cm^ - 3)
1425- self.mantle_shear_modulus = 1.45e11 # (Pa)
1426- self.mantle_density = 3.34 # (g / cm^ - 3)
1427-
1428- #SLR
1429+ #SLC
1430 self.earth_density = 5512 #average density of the Earth, (kg / m^3)
1431 return self
1432 #}}}
1433@@ -130,12 +119,6 @@
1434 md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', [1])
1435 md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', [1])
1436
1437- if 'GiaAnalysis' in analyses:
1438- md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
1439- md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
1440- md = checkfield(md,'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
1441- md = checkfield(md,'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
1442-
1443 if 'SealevelriseAnalysis' in analyses:
1444 md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1)
1445
1446@@ -160,10 +143,6 @@
1447 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1)
1448 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
1449 WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
1450- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
1451- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10.**3.)
1452- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
1453- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10.**3.)
1454 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
1455
1456 # }}}
1457Index: ../trunk-jpl/src/m/classes/matenhancedice.py
1458===================================================================
1459--- ../trunk-jpl/src/m/classes/matenhancedice.py (revision 26058)
1460+++ ../trunk-jpl/src/m/classes/matenhancedice.py (revision 26059)
1461@@ -1,15 +1,14 @@
1462+from checkfield import checkfield
1463 from fielddisplay import fielddisplay
1464 from project3d import project3d
1465-from checkfield import checkfield
1466 from WriteData import WriteData
1467
1468
1469 class matenhancedice(object):
1470- """
1471- MATICE class definition
1472+ """MATICE class definition
1473
1474- Usage:
1475- matenhancedice = matenhancedice()
1476+ Usage:
1477+ matenhancedice = matenhancedice()
1478 """
1479
1480 def __init__(self): #{{{
1481@@ -31,18 +30,17 @@
1482 self.rheology_n = float('NaN')
1483 self.rheology_law = ''
1484
1485- #GIA
1486- self.lithosphere_shear_modulus = 0.
1487- self.lithosphere_density = 0.
1488- self.mantle_shear_modulus = 0.
1489- self.mantle_density = 0.
1490+ #SLC
1491+ self.earth_density = 0 # average density of the Earth, (kg/m^3)
1492
1493- #SLR
1494- self.earth_density = 0 # average density of the Earth, (kg/m^3)
1495 self.setdefaultparameters()
1496 #}}}
1497
1498 def __repr__(self): #{{{
1499+ # TODO:
1500+ # - Convert all formatting to calls to <string>.format (see any
1501+ # already converted <class>.__repr__ method for examples)
1502+ #
1503 s = " Materials:"
1504 s = "%s\n%s" % (s, fielddisplay(self, "rho_ice", "ice density [kg/m^3]"))
1505 s = "%s\n%s" % (s, fielddisplay(self, "rho_water", "water density [kg/m^3]"))
1506@@ -61,10 +59,6 @@
1507 s = "%s\n%s" % (s, fielddisplay(self, "rheology_B", "flow law parameter [Pa s^(1/n)]"))
1508 s = "%s\n%s" % (s, fielddisplay(self, "rheology_n", "Glen's flow law exponent"))
1509 s = "%s\n%s" % (s, fielddisplay(self, "rheology_law", "law for the temperature dependance of the rheology: 'None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius' or 'LliboutryDuval'"))
1510- s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_shear_modulus", "Lithosphere shear modulus [Pa]"))
1511- s = "%s\n%s" % (s, fielddisplay(self, "lithosphere_density", "Lithosphere density [g/cm^-3]"))
1512- s = "%s\n%s" % (s, fielddisplay(self, "mantle_shear_modulus", "Mantle shear modulus [Pa]"))
1513- s = "%s\n%s" % (s, fielddisplay(self, "mantle_density", "Mantle density [g/cm^-3]"))
1514 s = "%s\n%s" % (s, fielddisplay(self, "earth_density", "Mantle density [kg/m^-3]"))
1515
1516 return s
1517@@ -108,13 +102,13 @@
1518 #available: none, paterson and arrhenius
1519 self.rheology_law = 'Paterson'
1520
1521- # GIA:
1522+ #GIA
1523 self.lithosphere_shear_modulus = 6.7 * 10**10 # (Pa)
1524 self.lithosphere_density = 3.32 # (g / cm^ - 3)
1525 self.mantle_shear_modulus = 1.45 * 10**11 # (Pa)
1526 self.mantle_density = 3.34 # (g / cm^ - 3)
1527
1528- #SLR
1529+ #SLC
1530 self.earth_density = 5512 #average density of the Earth, (kg / m^3)
1531
1532 return self
1533@@ -131,12 +125,6 @@
1534 md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval'])
1535 md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
1536
1537- if 'GiaAnalysis' in analyses:
1538- md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
1539- md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
1540- md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
1541- md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
1542-
1543 if 'SealevelriseAnalysis' in analyses:
1544 md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1)
1545 return md
1546@@ -161,9 +149,5 @@
1547 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
1548 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2)
1549 WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
1550- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
1551- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 10**3)
1552- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
1553- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 10**3)
1554 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
1555 # }}}
1556Index: ../trunk-jpl/src/m/classes/materials.m
1557===================================================================
1558--- ../trunk-jpl/src/m/classes/materials.m (revision 26058)
1559+++ ../trunk-jpl/src/m/classes/materials.m (revision 26059)
1560@@ -122,7 +122,6 @@
1561 self.rheology_B = 1*1e8;
1562 self.rheology_n = 3;
1563
1564-
1565 case 'litho'
1566 %we default to a configuration that enables running GIA solutions using giacaron and/or giaivins.
1567 self.numlayers=2;
1568Index: ../trunk-jpl/src/m/classes/materials.py
1569===================================================================
1570--- ../trunk-jpl/src/m/classes/materials.py (revision 26058)
1571+++ ../trunk-jpl/src/m/classes/materials.py (revision 26059)
1572@@ -7,12 +7,11 @@
1573
1574
1575 class materials(object):
1576- '''
1577- MATERIALS class definition
1578+ """MATERIALS class definition
1579
1580- Usage:
1581- materials = materials()
1582- '''
1583+ Usage:
1584+ materials = materials()
1585+ """
1586
1587 def __init__(self, *args): #{{{
1588 self.nature = []
1589@@ -37,6 +36,7 @@
1590 setattr(self, 'latentheat', 0)
1591 setattr(self, 'thermalconductivity', 0)
1592 setattr(self, 'temperateiceconductivity', 0)
1593+ setattr(self, 'effectiveconductivity_averaging', 0)
1594 setattr(self, 'meltingpoint', 0)
1595 setattr(self, 'beta', 0)
1596 setattr(self, 'mixed_layer_capacity', 0)
1597@@ -59,9 +59,9 @@
1598 setattr(self, 'rho_ice', 0)
1599 setattr(self, 'rho_water', 0)
1600 setattr(self, 'rho_freshwater', 0)
1601- setattr(self, 'earth_density', 0)
1602 else:
1603 raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
1604+ setattr(self, 'earth_density', 0)
1605
1606 #set default parameters:
1607 self.setdefaultparameters()
1608@@ -73,37 +73,58 @@
1609 if nat == 'ice':
1610 #ice density (kg/m^3)
1611 self.rho_ice = 917.
1612+
1613 #ocean water density (kg/m^3)
1614 self.rho_water = 1023.
1615+
1616 #fresh water density (kg/m^3)
1617 self.rho_freshwater = 1000.
1618+
1619 #water viscosity (N.s/m^2)
1620 self.mu_water = 0.001787
1621+
1622 #ice heat capacity cp (J/kg/K)
1623 self.heatcapacity = 2093.
1624+
1625 #ice latent heat of fusion L (J / kg)
1626 self.latentheat = 3.34e5
1627+
1628 #ice thermal conductivity (W/m/K)
1629 self.thermalconductivity = 2.4
1630+
1631 #wet ice thermal conductivity (W/m/K)
1632 self.temperateiceconductivity = 0.24
1633+
1634+ #computation of effective conductivity
1635+ self.effectiveconductivity_averaging = 1
1636+
1637 #the melting point of ice at 1 atmosphere of pressure in K
1638 self.meltingpoint = 273.15
1639+
1640 #rate of change of melting point with pressure (K/Pa)
1641 self.beta = 9.8e-8
1642+
1643 #mixed layer (ice-water interface) heat capacity (J/kg/K)
1644 self.mixed_layer_capacity = 3974.
1645+
1646 #thermal exchange velocity (ice-water interface) (m/s)
1647 self.thermal_exchange_velocity = 1.00e-4
1648+
1649 #Rheology law: what is the temperature dependence of B with T
1650 #available: none, paterson and arrhenius
1651 self.rheology_law = 'Paterson'
1652+
1653+ #Rheology fields default
1654+ self.rheology_B = 1e8
1655+ self.rheology_n = 3
1656 elif nat == 'litho':
1657- #we default to a configuration that enables running GIA solutions using giacaron and / or giaivins.
1658+ #we default to a configuration that enables running GIA solutions using giacaron and/or giaivins.
1659 self.numlayers = 2
1660+
1661 #center of the earth (approximation, must not be 0), then the lab (lithosphere / asthenosphere boundary) then the surface
1662 #(with 1d3 to avoid numerical singularities)
1663 self.radius = [1e3, 6278e3, 6378e3]
1664+
1665 self.viscosity = [1e21, 1e40] #mantle and lithosphere viscosity (respectively) [Pa.s]
1666 self.lame_mu = [1.45e11, 6.7e10] # (Pa) #lithosphere and mantle shear modulus (respectively) [Pa]
1667 self.lame_lambda = self.lame_mu # (Pa) #mantle and lithosphere lamba parameter (respectively) [Pa]
1668@@ -115,14 +136,17 @@
1669 elif nat == 'hydro':
1670 #ice density (kg/m^3)
1671 self.rho_ice = 917.
1672+
1673 #ocean water density (kg/m^3)
1674 self.rho_water = 1023.
1675- #average density of the Earth (kg/m^3)
1676- self.earth_density = 5512
1677+
1678 #fresh water density (kg/m^3)
1679 self.rho_freshwater = 1000.
1680 else:
1681 raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
1682+
1683+ #average density of the Earth (kg/m^3)
1684+ self.earth_density = 5512
1685 #}}}
1686
1687 def __repr__(self): #{{{
1688@@ -222,7 +246,6 @@
1689 #1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
1690 WriteData(fid, prefix, 'name', 'md.materials.nature', 'data', naturetointeger(self.nature), 'format', 'IntMat', 'mattype', 3)
1691 WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER, this can evolve if you have classes
1692-
1693 for i in range(len(self.nature)):
1694 nat = self.nature[i]
1695 if nat == 'ice':
1696@@ -235,6 +258,7 @@
1697 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'latentheat', 'format', 'Double')
1698 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'thermalconductivity', 'format', 'Double')
1699 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'temperateiceconductivity', 'format', 'Double')
1700+ WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'effectiveconductivity_averaging', 'format', 'Integer')
1701 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'meltingpoint', 'format', 'Double')
1702 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'beta', 'format', 'Double')
1703 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mixed_layer_capacity', 'format', 'Double')
1704@@ -253,13 +277,19 @@
1705 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'isburgers', 'format', 'DoubleMat', 'mattype', 3)
1706 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_viscosity', 'format', 'DoubleMat', 'mattype', 3)
1707 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'burgers_mu', 'format', 'DoubleMat', 'mattype', 3)
1708+ # Compute earth density compatible with our layer density distribution
1709+ earth_density = 0
1710+ for i in range(len(self.numlayers)):
1711+ earth_density = earth_density + (self.radius[i + 1] ** 3 - self.radius[i] ** 3) * self.density[i]
1712+ earth_density = earth_density / self.radius[self.numlayers + 1] ** 3
1713+ self.earth_density = earth_density
1714 elif nat == 'hydro':
1715 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')
1716 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_water', 'format', 'Double')
1717- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
1718 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_freshwater', 'format', 'Double')
1719 else:
1720 raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
1721+ WriteData(fid, prefix, 'data', self.earth_density, 'name', 'md.materials.earth_density', 'format', 'Double')
1722 #}}}
1723
1724 def extrude(self, md): #{{{
1725Index: ../trunk-jpl/src/m/classes/matestar.py
1726===================================================================
1727--- ../trunk-jpl/src/m/classes/matestar.py (revision 26058)
1728+++ ../trunk-jpl/src/m/classes/matestar.py (revision 26059)
1729@@ -33,13 +33,7 @@
1730 self.rheology_Es = np.nan
1731 self.rheology_law = ''
1732
1733- #GIA
1734- self.lithosphere_shear_modulus = 0.
1735- self.lithosphere_density = 0.
1736- self.mantle_shear_modulus = 0.
1737- self.mantle_density = 0.
1738-
1739- #SLR
1740+ #slc
1741 self.earth_density = 0
1742
1743 #set default parameters
1744@@ -65,10 +59,6 @@
1745 s = "%s\n%s" % (s, fielddisplay(self, 'rheology_Ec', 'compressive enhancement factor'))
1746 s = "%s\n%s" % (s, fielddisplay(self, 'rheology_Es', 'shear enhancement factor'))
1747 s = "%s\n%s" % (s, fielddisplay(self, 'rheology_law', ['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', ''Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']))
1748- s = "%s\n%s" % (s, fielddisplay(self, 'lithosphere_shear_modulus', 'Lithosphere shear modulus [Pa]'))
1749- s = "%s\n%s" % (s, fielddisplay(self, 'lithosphere_density', 'Lithosphere density [g/cm^-3]'))
1750- s = "%s\n%s" % (s, fielddisplay(self, 'mantle_shear_modulus', 'Mantle shear modulus [Pa]'))
1751- s = "%s\n%s" % (s, fielddisplay(self, 'mantle_density', 'Mantle density [g/cm^-3]'))
1752 s = "%s\n%s" % (s, fielddisplay(self, 'earth_density', 'Mantle density [kg/m^-3]'))
1753
1754 return s
1755@@ -112,12 +102,7 @@
1756 #Rheology law: what is the temperature dependence of B with T
1757 #available: none, paterson and arrhenius
1758 self.rheology_law = 'Paterson'
1759- # GIA
1760- self.lithosphere_shear_modulus = 6.7e10 # (Pa)
1761- self.lithosphere_density = 3.32 # (g / cm^ - 3)
1762- self.mantle_shear_modulus = 1.45e11 # (Pa)
1763- self.mantle_density = 3.34 # (g / cm^ - 3)
1764- #SLR
1765+ #slc
1766 self.earth_density = 5512 # average density of the Earth, (kg / m^3)
1767
1768 return self
1769@@ -134,12 +119,6 @@
1770 md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval'])
1771 md = checkfield(md, 'fieldname', 'materials.effectiveconductivity_averaging', 'numel', [1], 'values', [0, 1, 2])
1772
1773- if 'GiaAnalysis' in analyses:
1774- md = checkfield(md, 'fieldname', 'materials.lithosphere_shear_modulus', '>', 0, 'numel', 1)
1775- md = checkfield(md, 'fieldname', 'materials.lithosphere_density', '>', 0, 'numel', 1)
1776- md = checkfield(md, 'fieldname', 'materials.mantle_shear_modulus', '>', 0, 'numel', 1)
1777- md = checkfield(md, 'fieldname', 'materials.mantle_density', '>', 0, 'numel', 1)
1778-
1779 if 'SealevelriseAnalysis' in analyses:
1780 md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1)
1781
1782@@ -165,9 +144,5 @@
1783 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_Ec', 'format', 'DoubleMat', 'mattype', 1)
1784 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rheology_Es', 'format', 'DoubleMat', 'mattype', 1)
1785 WriteData(fid, prefix, 'data', self.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String')
1786- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_shear_modulus', 'format', 'Double')
1787- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'lithosphere_density', 'format', 'Double', 'scale', 1.0e3)
1788- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_shear_modulus', 'format', 'Double')
1789- WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'mantle_density', 'format', 'Double', 'scale', 1.0e3)
1790 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'earth_density', 'format', 'Double')
1791 # }}}
1792Index: ../trunk-jpl/src/m/classes/matice.m
1793===================================================================
1794--- ../trunk-jpl/src/m/classes/matice.m (revision 26058)
1795+++ ../trunk-jpl/src/m/classes/matice.m (revision 26059)
1796@@ -103,20 +103,18 @@
1797
1798 end % }}}
1799 function md = checkconsistency(self,md,solution,analyses) % {{{
1800-
1801+
1802 if strcmpi(solution,'TransientSolution') & md.transient.isslc,
1803 md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1);
1804- return;
1805- end
1806-
1807- md = checkfield(md,'fieldname','materials.rho_ice','>',0);
1808- md = checkfield(md,'fieldname','materials.rho_water','>',0);
1809- md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
1810- md = checkfield(md,'fieldname','materials.mu_water','>',0);
1811- md = checkfield(md,'fieldname','materials.rheology_B','>',0,'universal',1,'NaN',1,'Inf',1);
1812- md = checkfield(md,'fieldname','materials.rheology_n','>',0,'universal',1,'NaN',1,'Inf',1);
1813- md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'});
1814- md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
1815+ else
1816+ md = checkfield(md,'fieldname','materials.rho_ice','>',0);
1817+ md = checkfield(md,'fieldname','materials.rho_water','>',0);
1818+ md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
1819+ md = checkfield(md,'fieldname','materials.mu_water','>',0);
1820+ md = checkfield(md,'fieldname','materials.rheology_B','>',0,'universal',1,'NaN',1,'Inf',1);
1821+ md = checkfield(md,'fieldname','materials.rheology_n','>',0,'universal',1,'NaN',1,'Inf',1);
1822+ md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval' 'NyeCO2' 'NyeH2O'});
1823+ md = checkfield(md,'fieldname','materials.effectiveconductivity_averaging','numel',[1],'values',[0 1 2]);
1824
1825 end % }}}
1826 function disp(self) % {{{
1827Index: ../trunk-jpl/src/m/classes/model.m
1828===================================================================
1829--- ../trunk-jpl/src/m/classes/model.m (revision 26058)
1830+++ ../trunk-jpl/src/m/classes/model.m (revision 26059)
1831@@ -43,7 +43,7 @@
1832 frontalforcings = 0;
1833 love = 0;
1834 esa = 0;
1835- sampling = 0;
1836+ sampling = 0;
1837
1838 autodiff = 0;
1839 inversion = 0;
1840@@ -155,7 +155,7 @@
1841 %2019 Jan..
1842 if isa(md.frontalforcings,'double');
1843 if(isprop('meltingrate',md.calving) & ~isnan(md.calving.meltingrate))
1844- disp('Warning: md.calving.meltingrate is now in md.frontalforcings');
1845+ gia disp('Warning: md.calving.meltingrate is now in md.frontalforcings');
1846 end
1847 md.frontalforcings=frontalforcings(md.calving);
1848 end
1849@@ -186,8 +186,8 @@
1850 md.smb.isconstrainsurfaceT = 0;
1851 end
1852 end
1853- end
1854- %2021 February 17
1855+ end
1856+ %2021 February 17
1857 if isa(md.sampling,'double'); md.sampling=sampling(); end
1858 end% }}}
1859 end
1860@@ -241,7 +241,7 @@
1861 md.frontalforcings = frontalforcings();
1862 md.love = fourierlove();
1863 md.esa = esa();
1864- md.sampling = sampling();
1865+ md.sampling = sampling();
1866 md.autodiff = autodiff();
1867 md.inversion = inversion();
1868 md.qmu = qmu();
1869@@ -614,8 +614,8 @@
1870 %size = number of elements * n
1871 elseif fieldsize(1)==numberofelements1
1872 md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:);
1873- elseif (fieldsize(1)==numberofelements1+1)
1874- md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)];
1875+ elseif (fieldsize(1)==numberofelements1+1)
1876+ md2.(model_fields{i}).(object_fields{j})=[field(pos_elem,:); field(end,:)];
1877 end
1878 end
1879 else
1880@@ -627,7 +627,7 @@
1881 %size = number of elements * n
1882 elseif fieldsize(1)==numberofelements1
1883 md2.(model_fields{i})=field(pos_elem,:);
1884- elseif (fieldsize(1)==numberofelements1+1)
1885+ elseif (fieldsize(1)==numberofelements1+1)
1886 md2.(model_fields{i})=[field(pos_elem,:); field(end,:)];
1887 end
1888 end
1889@@ -774,18 +774,18 @@
1890 %get subfields
1891 % loop over time steps
1892 for p=1:length(md1.results.(solutionfields{i}))
1893- current = md1.results.(solutionfields{i})(p);
1894- solutionsubfields=fields(current);
1895- for j=1:length(solutionsubfields),
1896+ current = md1.results.(solutionfields{i})(p);
1897+ solutionsubfields=fields(current);
1898+ for j=1:length(solutionsubfields),
1899 field=md1.results.(solutionfields{i})(p).(solutionsubfields{j});
1900 if length(field)==numberofvertices1,
1901- md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node);
1902+ md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_node);
1903 elseif length(field)==numberofelements1,
1904- md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem);
1905+ md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field(pos_elem);
1906 else
1907- md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field;
1908+ md2.results.(solutionfields{i})(p).(solutionsubfields{j})=field;
1909 end
1910- end
1911+ end
1912 end
1913 else
1914 field=md1.results.(solutionfields{i});
1915@@ -1500,7 +1500,7 @@
1916 disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings'));
1917 disp(sprintf('%19s: %-22s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
1918 disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
1919- disp(sprintf('%19s: %-22s -- %s','sampling' ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
1920+ disp(sprintf('%19s: %-22s -- %s','sampling' ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
1921 disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
1922 disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
1923 disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(self.qmu) ']'],'Dakota properties'));
1924Index: ../trunk-jpl/src/m/classes/model.py
1925===================================================================
1926--- ../trunk-jpl/src/m/classes/model.py (revision 26058)
1927+++ ../trunk-jpl/src/m/classes/model.py (revision 26059)
1928@@ -54,8 +54,6 @@
1929 from thermal import thermal
1930 from steadystate import steadystate
1931 from transient import transient
1932-from giaivins import giaivins
1933-from giamme import giamme
1934 from esa import esa
1935 from autodiff import autodiff
1936 from inversion import inversion
1937@@ -113,7 +111,6 @@
1938 self.levelset = None
1939 self.calving = None
1940 self.frontalforcings = None
1941- self.gia = None
1942 self.love = None
1943 self.esa = None
1944 self.sampling = None
1945@@ -171,7 +168,6 @@
1946 'levelset',
1947 'calving',
1948 'frontalforcings',
1949- 'gia',
1950 'love',
1951 'esa',
1952 'sampling',
1953@@ -224,7 +220,6 @@
1954 s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("levelset", "[%s %s]" % ("1x1", obj.levelset.__class__.__name__), "parameters for moving boundaries (level - set method)"))
1955 s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("calving", "[%s %s]" % ("1x1", obj.calving.__class__.__name__), "parameters for calving"))
1956 s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("frontalforcings", "[%s %s]" % ("1x1", obj.frontalforcings.__class__.__name__), "parameters for frontalforcings"))
1957- s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("gia", "[%s %s]" % ("1x1", obj.gia.__class__.__name__), "parameters for gia solution"))
1958 s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("esa", "[%s %s]" % ("1x1", obj.esa.__class__.__name__), "parameters for elastic adjustment solution"))
1959 s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("sampling", "[%s %s]" % ("1x1", obj.sampling.__class__.__name__), "parameters for stochastic sampler"))
1960 s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("love", "[%s %s]" % ("1x1", obj.love.__class__.__name__), "parameters for love solution"))
1961@@ -271,7 +266,6 @@
1962 self.levelset = levelset()
1963 self.calving = calving()
1964 self.frontalforcings = frontalforcings()
1965- self.gia = giamme()
1966 self.love = fourierlove()
1967 self.esa = esa()
1968 self.sampling = sampling()
1969@@ -851,13 +845,6 @@
1970 if not np.isnan(md.initialization.watercolumn).all():
1971 md.initialization.watercolumn = project2d(md, md.initialization.watercolumn, 1)
1972
1973- # giaivins
1974- if md.gia.__class__.__name__ == 'giaivins':
1975- if not np.isnan(md.gia.mantle_viscosity).all():
1976- md.gia.mantle_viscosity = project2d(md, md.gia.mantle_viscosity, 1)
1977- if not np.isnan(md.gia.lithosphere_thickness).all():
1978- md.gia.lithosphere_thickness = project2d(md, md.gia.lithosphere_thickness, 1)
1979-
1980 # elementstype
1981 if not np.isnan(md.flowequation.element_equation).all():
1982 md.flowequation.element_equation = project2d(md, md.flowequation.element_equation, 1)
1983Index: ../trunk-jpl/src/m/classes/sealevelmodel.m
1984===================================================================
1985--- ../trunk-jpl/src/m/classes/sealevelmodel.m (revision 26058)
1986+++ ../trunk-jpl/src/m/classes/sealevelmodel.m (revision 26059)
1987@@ -90,15 +90,15 @@
1988 for i=1:length(slm.icecaps),
1989 md= slm.icecaps{i};
1990 if ~isempty(md.solidearth.external),
1991- error(sprintf('cannot run external forcings on an ice sheet when runing a coupling earth/ice sheet model');
1992+ error('cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model');
1993 end
1994
1995 end
1996- %make sure that we have the rigth grd model for computing ouf sealevel patterns:
1997+ %make sure that we have the right grd model for computing out sealevel patterns:
1998 for i=1:length(slm.icecaps),
1999 md= slm.icecaps{i};
2000 if md.solidearth.settings.grdmodel~=0
2001- error(sprintf('sealevelmodel checkconsistenty error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i));
2002+ error(sprintf('sealevelmodel checkconsistency error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap %i',i));
2003 end
2004 end
2005
2006Index: ../trunk-jpl/src/m/classes/sealevelmodel.py
2007===================================================================
2008--- ../trunk-jpl/src/m/classes/sealevelmodel.py (revision 26058)
2009+++ ../trunk-jpl/src/m/classes/sealevelmodel.py (revision 26059)
2010@@ -110,6 +110,22 @@
2011 md = slm.icecaps[i]
2012 if np.nonzero(md.slr.steric_rate - slm.earth.slr.steric_rate[slm.earth.slr.transitions[i]]) != []:
2013 raise Exception('steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name))
2014+
2015+ # Make sure grd is the same everywhere
2016+ for i in range(len(slm.icecaps)):
2017+ md = slm.icecaps[i]
2018+ if md.solidearthsettings.isgrd != slm.earth.solidearthsettings.isgrd:
2019+ raise RuntimeError('isgrd on ice cap {} is not the same as for the earth\n'.format(md.miscellaneous.name))
2020+ # Make sure that there is no solid earth external forcing on the basins
2021+ for i in range(len(slm.icecaps)):
2022+ md = slm.icecaps[i]
2023+ if md.solidearth.external:
2024+ raise RuntimeError('cannot run external forcings on an ice sheet when running a coupling earth/ice sheet model')
2025+ # Make sure that we have the right grd model for computing our sealevel patterns
2026+ for i in range(len(slm.icecaps)):
2027+ md = slm.icecaps[i]
2028+ if md.solidearth.settings.grdmodel != 0:
2029+ raise RuntimeError('sealevelmodel checkconsistency error message: ice sheets do not run GRD module, specify solidearth.settings.grdmodel=0 on ice cap {}'.format(i))
2030 #}}}
2031
2032 def mergeresults(self): # {{{
2033Index: ../trunk-jpl/src/m/classes/transient.py
2034===================================================================
2035--- ../trunk-jpl/src/m/classes/transient.py (revision 26058)
2036+++ ../trunk-jpl/src/m/classes/transient.py (revision 26059)
2037@@ -13,10 +13,10 @@
2038 def __init__(self): # {{{
2039 self.issmb = 0
2040 self.ismasstransport = 0
2041+ self.isoceantransport = 0
2042 self.isstressbalance = 0
2043 self.isthermal = 0
2044 self.isgroundingline = 0
2045- self.isgia = 0
2046 self.isesa = 0
2047 self.isdamageevolution = 0
2048 self.ismovingfront = 0
2049@@ -23,7 +23,6 @@
2050 self.ishydrology = 0
2051 self.issampling = 0
2052 self.isslc = 0
2053- self.iscoupler = 0
2054 self.amr_frequency = 0
2055 self.isoceancoupling = 0
2056 self.requested_outputs = []
2057@@ -35,10 +34,10 @@
2058 s = ' transient solution parameters:\n'
2059 s += '{}\n'.format(fielddisplay(self, 'issmb', 'indicates if a surface mass balance solution is used in the transient'))
2060 s += '{}\n'.format(fielddisplay(self, 'ismasstransport', 'indicates if a masstransport solution is used in the transient'))
2061+ s += '{}\n'.format(fielddisplay(self, 'isoceantransport', 'indicates whether an ocean masstransport solution is used in the transient'))
2062 s += '{}\n'.format(fielddisplay(self, 'isstressbalance', 'indicates if a stressbalance solution is used in the transient'))
2063 s += '{}\n'.format(fielddisplay(self, 'isthermal', 'indicates if a thermal solution is used in the transient'))
2064 s += '{}\n'.format(fielddisplay(self, 'isgroundingline', 'indicates if a groundingline migration is used in the transient'))
2065- s += '{}\n'.format(fielddisplay(self, 'isgia', 'indicates if a postglacial rebound is used in the transient'))
2066 s += '{}\n'.format(fielddisplay(self, 'isesa', 'indicates whether an elastic adjustment model is used in the transient'))
2067 s += '{}\n'.format(fielddisplay(self, 'isdamageevolution', 'indicates whether damage evolution is used in the transient'))
2068 s += '{}\n'.format(fielddisplay(self, 'ismovingfront', 'indicates whether a moving front capability is used in the transient'))
2069@@ -46,7 +45,6 @@
2070 s += '{}\n'.format(fielddisplay(self, 'issampling', 'indicates whether sampling is used in the transient'))
2071 s += '{}\n'.format(fielddisplay(self, 'isslc', 'indicates if a sea level change solution is used in the transient'))
2072 s += '{}\n'.format(fielddisplay(self, 'isoceancoupling', 'indicates whether coupling with an ocean model is used in the transient'))
2073- s += '{}\n'.format(fielddisplay(self, 'iscoupler', 'indicates whether different models are being run with need for coupling'))
2074 s += '{}\n'.format(fielddisplay(self, 'amr_frequency', 'frequency at which mesh is refined in simulations with multiple time_steps'))
2075 s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'list of additional outputs requested'))
2076 return s
2077@@ -62,10 +60,10 @@
2078 def deactivateall(self): #{{{
2079 self.issmb = 0
2080 self.ismasstransport = 0
2081+ self.isoceantransport = 0
2082 self.isstressbalance = 0
2083 self.isthermal = 0
2084 self.isgroundingline = 0
2085- self.isgia = 0
2086 self.isesa = 0
2087 self.isdamageevolution = 0
2088 self.ismovingfront = 0
2089@@ -73,7 +71,6 @@
2090 self.issampling = 0
2091 self.isslc = 0
2092 self.isoceancoupling = 0
2093- self.iscoupler = 0
2094 self.amr_frequency = 0
2095
2096 # Default output
2097@@ -85,10 +82,10 @@
2098 #full analysis: Stressbalance, Masstransport and Thermal but no groundingline migration for now
2099 self.issmb = 1
2100 self.ismasstransport = 1
2101+ self.isoceantransport = 0
2102 self.isstressbalance = 1
2103 self.isthermal = 1
2104 self.isgroundingline = 0
2105- self.isgia = 0
2106 self.isesa = 0
2107 self.isdamageevolution = 0
2108 self.ismovingfront = 0
2109@@ -96,7 +93,6 @@
2110 self.issampling = 0
2111 self.isslc = 0
2112 self.isoceancoupling = 0
2113- self.iscoupler = 0
2114 self.amr_frequency = 0
2115
2116 # Default output
2117@@ -111,10 +107,10 @@
2118
2119 md = checkfield(md, 'fieldname', 'transient.issmb', 'numel', [1], 'values', [0, 1])
2120 md = checkfield(md, 'fieldname', 'transient.ismasstransport', 'numel', [1], 'values', [0, 1])
2121+ md = checkfield(md, 'fieldname', 'transient.isocealtransport', 'numel', [1], 'values', [0, 1])
2122 md = checkfield(md, 'fieldname', 'transient.isstressbalance', 'numel', [1], 'values', [0, 1])
2123 md = checkfield(md, 'fieldname', 'transient.isthermal', 'numel', [1], 'values', [0, 1])
2124 md = checkfield(md, 'fieldname', 'transient.isgroundingline', 'numel', [1], 'values', [0, 1])
2125- md = checkfield(md, 'fieldname', 'transient.isgia', 'numel', [1], 'values', [0, 1])
2126 md = checkfield(md, 'fieldname', 'transient.isesa', 'numel', [1], 'values', [0, 1])
2127 md = checkfield(md, 'fieldname', 'transient.isdamageevolution', 'numel', [1], 'values', [0, 1])
2128 md = checkfield(md, 'fieldname', 'transient.ishydrology', 'numel', [1], 'values', [0, 1])
2129@@ -122,7 +118,6 @@
2130 md = checkfield(md, 'fieldname', 'transient.ismovingfront', 'numel', [1], 'values', [0, 1])
2131 md = checkfield(md, 'fieldname', 'transient.isslc', 'numel', [1], 'values', [0, 1])
2132 md = checkfield(md, 'fieldname', 'transient.isoceancoupling', 'numel', [1], 'values', [0, 1])
2133- md = checkfield(md, 'fieldname', 'transient.iscoupler', 'numel', [1], 'values', [0, 1])
2134 md = checkfield(md, 'fieldname', 'transient.amr_frequency', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
2135
2136 if solution != 'TransientSolution' and md.transient.iscoupling:
2137@@ -135,10 +130,10 @@
2138 def marshall(self, prefix, md, fid): # {{{
2139 WriteData(fid, prefix, 'object', self, 'fieldname', 'issmb', 'format', 'Boolean')
2140 WriteData(fid, prefix, 'object', self, 'fieldname', 'ismasstransport', 'format', 'Boolean')
2141+ WriteData(fid, prefix, 'object', self, 'fieldname', 'isoceantransport', 'format', 'Boolean')
2142 WriteData(fid, prefix, 'object', self, 'fieldname', 'isstressbalance', 'format', 'Boolean')
2143 WriteData(fid, prefix, 'object', self, 'fieldname', 'isthermal', 'format', 'Boolean')
2144 WriteData(fid, prefix, 'object', self, 'fieldname', 'isgroundingline', 'format', 'Boolean')
2145- WriteData(fid, prefix, 'object', self, 'fieldname', 'isgia', 'format', 'Boolean')
2146 WriteData(fid, prefix, 'object', self, 'fieldname', 'isesa', 'format', 'Boolean')
2147 WriteData(fid, prefix, 'object', self, 'fieldname', 'isdamageevolution', 'format', 'Boolean')
2148 WriteData(fid, prefix, 'object', self, 'fieldname', 'ishydrology', 'format', 'Boolean')
2149@@ -146,7 +141,6 @@
2150 WriteData(fid, prefix, 'object', self, 'fieldname', 'issampling', 'format', 'Boolean')
2151 WriteData(fid, prefix, 'object', self, 'fieldname', 'isslc', 'format', 'Boolean')
2152 WriteData(fid, prefix, 'object', self, 'fieldname', 'isoceancoupling', 'format', 'Boolean')
2153- WriteData(fid, prefix, 'object', self, 'fieldname', 'iscoupler', 'format', 'Boolean')
2154 WriteData(fid, prefix, 'object', self, 'fieldname', 'amr_frequency', 'format', 'Integer')
2155
2156 # Process requested outputs
2157Index: ../trunk-jpl/src/m/solve/solve.py
2158===================================================================
2159--- ../trunk-jpl/src/m/solve/solve.py (revision 26058)
2160+++ ../trunk-jpl/src/m/solve/solve.py (revision 26059)
2161@@ -1,6 +1,5 @@
2162 from datetime import datetime
2163 import os
2164-import shutil
2165
2166 from ismodelselfconsistent import ismodelselfconsistent
2167 from loadresultsfromcluster import loadresultsfromcluster
2168@@ -21,6 +20,7 @@
2169 Solution types available comprise:
2170 - 'Stressbalance' or 'sb'
2171 - 'Masstransport' or 'mt'
2172+ - 'Oceantransport' or 'oceant'
2173 - 'Thermal' or 'th'
2174 - 'Steadystate' or 'ss'
2175 - 'Transient' or 'tr'
2176@@ -31,9 +31,8 @@
2177 - 'Hydrology' or 'hy'
2178 - 'DamageEvolution' or 'da'
2179 - 'Gia' or 'gia'
2180+ - 'Love' or 'lv'
2181 - 'Esa' or 'esa'
2182- - 'Sealevelchange' or 'slc'
2183- - 'Love' or 'lv'
2184 - 'Sampling' or 'smp'
2185
2186 Extra options:
2187@@ -54,6 +53,8 @@
2188 solutionstring = 'StressbalanceSolution'
2189 elif solutionstring.lower() == 'mt' or solutionstring.lower() == 'masstransport':
2190 solutionstring = 'MasstransportSolution'
2191+ elif solutionstring.lower() == 'oceant' or solutionstring.lower() == 'oceantransport':
2192+ solutionstring = 'OceantransportSolution'
2193 elif solutionstring.lower() == 'th' or solutionstring.lower() == 'thermal':
2194 solutionstring = 'ThermalSolution'
2195 elif solutionstring.lower() == 'st' or solutionstring.lower() == 'steadystate':
2196@@ -78,8 +79,6 @@
2197 solutionstring = 'LoveSolution'
2198 elif solutionstring.lower() == 'esa':
2199 solutionstring = 'EsaSolution'
2200- elif solutionstring.lower() == 'slc' or solutionstring.lower() == 'sealevelchange':
2201- solutionstring = 'SealevelchangeSolution'
2202 elif solutionstring.lower() == 'smp' or solutionstring.lower() == 'sampling':
2203 solutionstring = 'SamplingSolution'
2204 else:
2205@@ -150,9 +149,9 @@
2206 # Wait on lock
2207 if md.settings.waitonlock > 0:
2208 islock = waitonlock(md)
2209- if islock == 0: # no results to be loaded
2210+ if islock == 0: # no results to be loaded
2211 print('The results must be loaded manually with md = loadresultsfromcluster(md).')
2212- else: # load results
2213+ else: # load results
2214 if md.verbose.solution:
2215 print('loading results from cluster')
2216 md = loadresultsfromcluster(md)
2217Index: ../trunk-jpl/test/NightlyRun/test2001.py
2218===================================================================
2219--- ../trunk-jpl/test/NightlyRun/test2001.py (revision 26058)
2220+++ ../trunk-jpl/test/NightlyRun/test2001.py (revision 26059)
2221@@ -3,7 +3,7 @@
2222
2223 import numpy as np
2224
2225-from giaivins import giaivins
2226+from materials import *
2227 from model import *
2228 from parameterize import *
2229 from setflowequation import *
2230@@ -16,35 +16,50 @@
2231 md = setmask(md, '', '')
2232 md = parameterize(md, '../Par/SquareSheetConstrained.py')
2233
2234-#GIA
2235-md.gia = giaivins()
2236-md.gia.lithosphere_thickness = 100. * np.ones(md.mesh.numberofvertices) # in km
2237-md.gia.mantle_viscosity = 1.0e21 * np.ones(md.mesh.numberofvertices) # in Pa.s
2238-md.materials.lithosphere_shear_modulus = 6.7e10 # in Pa
2239-md.materials.lithosphere_density = 3.32 # in g/cm^3
2240-md.materials.mantle_shear_modulus = 1.45e11 # in Pa
2241-md.materials.mantle_density = 3.34 # in g/cm^3
2242+# GIA Ivins, 2 layer model
2243+md.solidearth.settings.grdmodel = 2
2244
2245-#Indicate what you want to compute
2246-md.gia.cross_section_shape = 1 # for square-edged x-section
2247+md.materials = materials('litho','ice')
2248+md.materials.numlayers = 2;
2249+md.materials.radius = [10, 6271e3, 6371e3]
2250+md.materials.density = [3.34e3, 3.32e3]
2251+md.materials.lame_mu = [1.45e11, 6.7e10]
2252+md.materials.viscosity = [1e21, 0]
2253+md.initialization.sealevel = np.zeros(md.mesh.numberofvertices)
2254+md.solidearth.settings.cross_section_shape = 1 # for square-edged x-section
2255+md.solidearth.settings.computesealevelchange = 0 # do not compute sea level, only deformation
2256+md.solidearth.requested_outputs = ['Sealevel', 'SealevelUGrd']
2257
2258-#Define loading history (see test2001.m for the description)
2259-md.timestepping.start_time = 2400000 # 2, 400 kyr
2260-md.timestepping.final_time = 2500000 # 2, 500 kyr
2261-md.geometry.thickness = np.array([
2262- np.append(md.geometry.thickness * 0.0, 0.0),
2263- np.append(md.geometry.thickness / 2.0, 0.1),
2264- np.append(md.geometry.thickness, 0.2),
2265- np.append(md.geometry.thickness, 1.0),
2266- np.append(md.geometry.thickness, md.timestepping.start_time)
2267+# Loading history
2268+md.timestepping.start_time = -2400000 # 4,800 kyr :: EVALUATION TIME
2269+md.timestepping.time_step = 2400000 # 2,400 kyr :: EVALUATION TIME
2270+# To get rid of default final_time, make sure final_time > start_time
2271+md.timestepping.final_time = 2400000 # 2,400 kyr
2272+md.masstransport.spcthickness np.array([
2273+ np.append(md.geometry.thickness, 0),
2274+ np.append(md.geometry.thickness, 2400000)
2275 ]).T
2276
2277+# Geometry at 0 initially
2278+md.geometry.thickness = np.zeros(md.mesh.numberofvertices)
2279+md.geometry.surface = np.zeros(md.mesh.numberofvertices)
2280+md.geometry.base = np.zeros(md.mesh.numberofvertices)
2281+
2282+# Physics
2283+md.transient.issmb = 0
2284+md.transient.isstressbalance = 0
2285+md.transient.isthermal = 0
2286+md.transient.ismasstransport = 0
2287+md.transient.isslc = 0
2288+
2289 #Solve for GIA deflection
2290+md.cluster = generic('name', gethostname(), 'np', 1)
2291 md.cluster = generic('name', gethostname(), 'np', 3)
2292-md.verbose = verbose('1111111')
2293-md = solve(md, 'Gia')
2294+md.verbose = verbose('11111111111')
2295+md.verbose.solver = 0
2296+md = solve(md, 'Transient')
2297
2298 #Fields and tolerances to track changes
2299-field_names = ['UGia', 'UGiaRate']
2300-field_tolerances = [1e-13, 1e-13]
2301-field_values = [md.results.GiaSolution.UGia, md.results.GiaSolution.UGiaRate]
2302+field_names = ['UGrd']
2303+field_tolerances = [1e-13]
2304+field_values = [md.results.TransientSolution[0].SealevelUGrd]
Note: See TracBrowser for help on using the repository browser.