source: issm/oecreview/Archive/18296-19100/ISSM-18967-18968.diff@ 19102

Last change on this file since 19102 was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 92.1 KB
RevLine 
[19102]1Index: ../trunk-jpl/test/NightlyRun/test236.py
2===================================================================
3--- ../trunk-jpl/test/NightlyRun/test236.py (revision 18967)
4+++ ../trunk-jpl/test/NightlyRun/test236.py (revision 18968)
5@@ -15,6 +15,7 @@
6 # Use of ispdd and isdelta18o methods
7 md.surfaceforcings = SMBpdd();
8 md.surfaceforcings.isdelta18o=1
9+md.surfaceforcings.ismungsm=0
10
11 # Add temperature, precipitation and delta18o needed to measure the surface mass balance
12 # creating delta18o
13@@ -49,14 +50,20 @@
14 for imonth in xrange(0,12):
15 md.surfaceforcings.precipitations_presentday[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5
16 md.surfaceforcings.precipitations_presentday[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.)
17+ md.surfaceforcings.precipitations_lgm[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5
18+ md.surfaceforcings.precipitations_lgm[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.)
19
20+# Interpolation factors
21+md.surfaceforcings.Tdiff[1,1:md.timestepping.final_time]=0.5;
22+md.surfaceforcings.sealev[1,1:md.timestepping.final_time]=0.5;
23+# Year of each data point
24+md.surfaceforcings.Tdiff[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time;
25+md.surfaceforcings.sealev[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time;
26+
27 # time steps and resolution
28 md.timestepping.time_step=20.
29 md.timestepping.final_time=60.
30
31-md.surfaceforcings.Pfac=[1,1]
32-md.surfaceforcings.Tdiff=[1,1]
33-md.surfaceforcings.sealev=[1,1]
34
35 #
36 md.transient.requested_outputs=['default','SurfaceforcingsMonthlytemperatures']
37Index: ../trunk-jpl/test/NightlyRun/test237.py
38===================================================================
39--- ../trunk-jpl/test/NightlyRun/test237.py (revision 18967)
40+++ ../trunk-jpl/test/NightlyRun/test237.py (revision 18968)
41@@ -11,20 +11,21 @@
42
43 md=triangle(model(),'../Exp/Square.exp',600000) #180000
44 md=setmask(md,'all','')
45+md=parameterize(md,'../Par/SquareShelf.py')
46
47 # Use of ispdd and isdelta18o methods
48 md.surfaceforcings = SMBpdd();
49-md.surfaceforcings.isdelta18o=1
50+md.surfaceforcings.isdelta18o=0
51+md.surfaceforcings.ismungsm=1
52
53-md=parameterize(md,'../Par/SquareShelf.py')
54
55-# Add temperature, precipitation and delta18o needed to measure the surface mass balance
56-# creating delta18o
57-delta18o=numpy.loadtxt('../Data/delta18o.data')
58-md.surfaceforcings.delta18o=delta18o
59-# creating delta18oSurface
60-md.surfaceforcings.delta18o_surface = numpy.zeros((2,numpy.size(delta18o,axis=1)))
61-md.surfaceforcings.delta18o_surface[1,:] = delta18o[1,:]
62+# # Add temperature, precipitation and delta18o needed to measure the surface mass balance
63+# # creating delta18o
64+# delta18o=numpy.loadtxt('../Data/delta18o.data')
65+# md.surfaceforcings.delta18o=delta18o
66+# # creating delta18oSurface
67+# md.surfaceforcings.delta18o_surface = numpy.zeros((2,numpy.size(delta18o,axis=1)))
68+# md.surfaceforcings.delta18o_surface[1,:] = delta18o[1,:]
69
70 # creating Present day and lgm temperatures
71 # Same temperature over the all region:
72@@ -51,16 +52,23 @@
73 for imonth in xrange(0,12):
74 md.surfaceforcings.precipitations_presentday[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5
75 md.surfaceforcings.precipitations_presentday[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.)
76+ md.surfaceforcings.precipitations_lgm[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5
77+ md.surfaceforcings.precipitations_lgm[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.)
78
79+# Interpolation factors
80+md.surfaceforcings.Pfac[1,1:md.timestepping.final_time]=0.5;
81+md.surfaceforcings.Tdiff[1,1:md.timestepping.final_time]=0.5;
82+md.surfaceforcings.sealev[1,1:md.timestepping.final_time]=0.5;
83+# Year of each data point
84+md.surfaceforcings.Pfac[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time;
85+md.surfaceforcings.Tdiff[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time;
86+md.surfaceforcings.sealev[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time;
87+
88 # time steps and resolution
89 md.timestepping.time_step=20.
90 md.settings.output_frequency=1
91 md.timestepping.final_time=60.
92
93-md.surfaceforcings.Pfac=[1,1]
94-md.surfaceforcings.Tdiff=[1,1]
95-md.surfaceforcings.sealev=[1,1]
96-
97 #
98 md.transient.requested_outputs=['default','SurfaceforcingsMonthlytemperatures']
99 md.extrude(3,1.)
100Index: ../trunk-jpl/test/NightlyRun/test236.m
101===================================================================
102--- ../trunk-jpl/test/NightlyRun/test236.m (revision 18967)
103+++ ../trunk-jpl/test/NightlyRun/test236.m (revision 18968)
104@@ -2,12 +2,18 @@
105 md=setmask(md,'all','');
106 md=parameterize(md,'../Par/SquareShelf.par');
107
108+%md.verbose=verbose('all');
109+
110 % Use of ispdd and isdelta18o methods
111 md.surfaceforcings = SMBpdd();
112 md.surfaceforcings.isdelta18o=1;
113+md.surfaceforcings.ismungsm=0;
114
115+%md.surfaceforcings.precipitation(1:md.mesh.numberofvertices,1:12)=0;
116+%md.surfaceforcings.monthlytemperatures(1:md.mesh.numberofvertices,1:12)=273;
117+
118 % Add temperature, precipitation and delta18o needed to measure the surface mass balance
119-% creating delta18o
120+% creating delta18o
121 load '../Data/delta18o.data'
122 md.surfaceforcings.delta18o=delta18o;
123 % creating delta18oSurface
124@@ -25,7 +31,8 @@
125 md.surfaceforcings.temperatures_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
126 end
127
128-% creating initialization and spc temperatures initialization and spc
129+% creating initialization and spc temperatures initialization and
130+% spc
131 md.thermal.spctemperature=mean(md.surfaceforcings.temperatures_lgm(1:md.mesh.numberofvertices,1:12),2); %-10*ones(md.mesh.numberofvertices,1);
132 md.thermal.spctemperature=repmat(md.thermal.spctemperature,1,md.timestepping.final_time/md.timestepping.time_step);
133 itemp=0:md.timestepping.time_step:md.timestepping.final_time-md.timestepping.time_step;
134@@ -36,21 +43,28 @@
135 % creating precipitation
136 for imonth=0:11
137 md.surfaceforcings.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5;
138+ md.surfaceforcings.precipitations_lgm(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5;
139+ % Time for the last line:
140 md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
141+ md.surfaceforcings.precipitations_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
142 end
143
144+% Interpolation factors
145+md.surfaceforcings.Tdiff(1,1:md.timestepping.final_time)=0.5;
146+md.surfaceforcings.sealev(1,1:md.timestepping.final_time)=0.5;
147+% Year of each data point
148+md.surfaceforcings.Tdiff(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time;
149+md.surfaceforcings.sealev(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time;
150+
151 % time steps and resolution
152 md.timestepping.time_step=20;
153+md.settings.output_frequency=1;
154 md.timestepping.final_time=60;
155
156-md.surfaceforcings.Pfac=[1;1];
157-md.surfaceforcings.Tdiff=[1;1];
158-md.surfaceforcings.sealev=[1;1];
159-
160 %
161 md.transient.requested_outputs={'default','SurfaceforcingsMonthlytemperatures'};
162 md=setflowequation(md,'SSA','all');
163-md.cluster=generic('name',oshostname(),'np',3);
164+md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster
165 md=solve(md,TransientSolutionEnum());
166
167 %Fields and tolerances to track changes
168Index: ../trunk-jpl/test/NightlyRun/test237.m
169===================================================================
170--- ../trunk-jpl/test/NightlyRun/test237.m (revision 18967)
171+++ ../trunk-jpl/test/NightlyRun/test237.m (revision 18968)
172@@ -1,19 +1,27 @@
173 md=triangle(model(),'../Exp/Square.exp',600000.);%180000
174 md=setmask(md,'all','');
175+md=parameterize(md,'../Par/SquareShelf.par');
176
177-% Use of ispdd and isdelta18o methods
178+%md.verbose=verbose('all');
179+
180+% Use of ispdd methods
181 md.surfaceforcings = SMBpdd();
182-md.surfaceforcings.isdelta18o=1;
183+md.surfaceforcings.isdelta18o=0;
184+md.surfaceforcings.ismungsm=1;
185
186-md=parameterize(md,'../Par/SquareShelf.par');
187+if md.surfaceforcings.isdelta18o==0 & md.surfaceforcings.ismungsm==0
188+ md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1);
189+ md.surfaceforcings.monthlytemperatures=273*ones(md.mesh.numberofvertices,1);
190+end
191+
192
193 % Add temperature, precipitation and delta18o needed to measure the surface mass balance
194-% creating delta18o
195-load '../Data/delta18o.data'
196-md.surfaceforcings.delta18o=delta18o;
197-% creating delta18oSurface
198-md.surfaceforcings.delta18o_surface(1,1:(length(delta18o))) = 0;
199-md.surfaceforcings.delta18o_surface(2,:) = delta18o(2,:);
200+% % creating delta18o
201+% load '../Data/delta18o.data'
202+% md.surfaceforcings.delta18o=delta18o;
203+% % creating delta18oSurface
204+% md.surfaceforcings.delta18o_surface(1,1:(length(delta18o))) = 0;
205+% md.surfaceforcings.delta18o_surface(2,:) = delta18o(2,:);
206
207 % creating Present day and lgm temperatures
208 % Same temperature over the all region:
209@@ -37,12 +45,19 @@
210 % creating precipitation
211 for imonth=0:11
212 md.surfaceforcings.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5;
213+ md.surfaceforcings.precipitations_lgm(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5;
214+ % Time for the last line:
215 md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
216+ md.surfaceforcings.precipitations_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
217 end
218
219-md.surfaceforcings.Pfac=[1;1];
220-md.surfaceforcings.Tdiff=[1;1];
221-md.surfaceforcings.sealev=[1;1];
222+md.surfaceforcings.Pfac(1,1:md.timestepping.final_time)=0.5;
223+md.surfaceforcings.Tdiff(1,1:md.timestepping.final_time)=0.5;
224+md.surfaceforcings.sealev(1,1:md.timestepping.final_time)=0.5;
225+% Year of each data point
226+md.surfaceforcings.Pfac(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time;
227+md.surfaceforcings.Tdiff(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time;
228+md.surfaceforcings.sealev(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time;
229
230 % time steps and resolution
231 md.timestepping.time_step=20;
232@@ -51,10 +66,11 @@
233
234 %
235 md.transient.requested_outputs={'default','SurfaceforcingsMonthlytemperatures'};
236-md=extrude(md,3,1.);
237+md=extrude(md,3,1);
238+
239 md=setflowequation(md,'SSA','all');
240-md.cluster=generic('name',oshostname(),'np',1);
241-md=solve(md,TransientSolutionEnum());
242+md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster
243+md=solve(md,TransientSolutionEnum);
244
245 %Fields and tolerances to track changes
246 field_names ={'Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Temperature1','BasalforcingsGroundediceMeltingRate1', ...
247Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
248===================================================================
249--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 18967)
250+++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 18968)
251@@ -119,7 +119,7 @@
252
253 int stabilization,finiteelement,smb_model;
254 bool dakota_analysis;
255- bool isdelta18o;
256+ bool isdelta18o,ismungsm;
257 bool isgroundingline;
258 bool islevelset;
259
260@@ -178,17 +178,19 @@
261 break;
262 case SMBpddEnum:
263 iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
264+ iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum);
265 iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
266- if(isdelta18o){
267+ if(isdelta18o || ismungsm){
268 iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesLgmEnum);
269 iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesPresentdayEnum);
270 iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsPresentdayEnum);
271 iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsLgmEnum);
272 }
273 else{
274- iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
275+ iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
276 iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum);
277 }
278+
279 break;
280 case SMBgradientsEnum:
281 iomodel->FetchDataToInput(elements,SurfaceforcingsHrefEnum);
282Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
283===================================================================
284--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18967)
285+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18968)
286@@ -343,6 +343,7 @@
287 SurfaceforcingsDelta18oEnum,
288 SurfaceforcingsDelta18oSurfaceEnum,
289 SurfaceforcingsIsdelta18oEnum,
290+ SurfaceforcingsIsmungsmEnum,
291 SurfaceforcingsPrecipitationsPresentdayEnum,
292 SurfaceforcingsPrecipitationsLgmEnum,
293 SurfaceforcingsTemperaturesPresentdayEnum,
294Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
295===================================================================
296--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18967)
297+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18968)
298@@ -349,6 +349,7 @@
299 case SurfaceforcingsDelta18oEnum : return "SurfaceforcingsDelta18o";
300 case SurfaceforcingsDelta18oSurfaceEnum : return "SurfaceforcingsDelta18oSurface";
301 case SurfaceforcingsIsdelta18oEnum : return "SurfaceforcingsIsdelta18o";
302+ case SurfaceforcingsIsmungsmEnum : return "SurfaceforcingsIsmungsm";
303 case SurfaceforcingsPrecipitationsPresentdayEnum : return "SurfaceforcingsPrecipitationsPresentday";
304 case SurfaceforcingsPrecipitationsLgmEnum : return "SurfaceforcingsPrecipitationsLgm";
305 case SurfaceforcingsTemperaturesPresentdayEnum : return "SurfaceforcingsTemperaturesPresentday";
306Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
307===================================================================
308--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18967)
309+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18968)
310@@ -355,6 +355,7 @@
311 else if (strcmp(name,"SurfaceforcingsDelta18o")==0) return SurfaceforcingsDelta18oEnum;
312 else if (strcmp(name,"SurfaceforcingsDelta18oSurface")==0) return SurfaceforcingsDelta18oSurfaceEnum;
313 else if (strcmp(name,"SurfaceforcingsIsdelta18o")==0) return SurfaceforcingsIsdelta18oEnum;
314+ else if (strcmp(name,"SurfaceforcingsIsmungsm")==0) return SurfaceforcingsIsmungsmEnum;
315 else if (strcmp(name,"SurfaceforcingsPrecipitationsPresentday")==0) return SurfaceforcingsPrecipitationsPresentdayEnum;
316 else if (strcmp(name,"SurfaceforcingsPrecipitationsLgm")==0) return SurfaceforcingsPrecipitationsLgmEnum;
317 else if (strcmp(name,"SurfaceforcingsTemperaturesPresentday")==0) return SurfaceforcingsTemperaturesPresentdayEnum;
318@@ -381,11 +382,11 @@
319 else if (strcmp(name,"SurfaceforcingsRunoff")==0) return SurfaceforcingsRunoffEnum;
320 else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
321 else if (strcmp(name,"SurfaceforcingsMelt")==0) return SurfaceforcingsMeltEnum;
322- else if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum;
323 else stage=4;
324 }
325 if(stage==4){
326- if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
327+ if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum;
328+ else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
329 else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum;
330 else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
331 else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
332@@ -504,11 +505,11 @@
333 else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
334 else if (strcmp(name,"Masscon")==0) return MassconEnum;
335 else if (strcmp(name,"MassconName")==0) return MassconNameEnum;
336- else if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum;
337 else stage=5;
338 }
339 if(stage==5){
340- if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
341+ if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum;
342+ else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum;
343 else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
344 else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum;
345 else if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum;
346@@ -627,11 +628,11 @@
347 else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
348 else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
349 else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
350- else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
351 else stage=6;
352 }
353 if(stage==6){
354- if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
355+ if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
356+ else if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
357 else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum;
358 else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum;
359 else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum;
360@@ -750,11 +751,11 @@
361 else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
362 else if (strcmp(name,"Gset")==0) return GsetEnum;
363 else if (strcmp(name,"Index")==0) return IndexEnum;
364- else if (strcmp(name,"Indexed")==0) return IndexedEnum;
365 else stage=7;
366 }
367 if(stage==7){
368- if (strcmp(name,"Intersect")==0) return IntersectEnum;
369+ if (strcmp(name,"Indexed")==0) return IndexedEnum;
370+ else if (strcmp(name,"Intersect")==0) return IntersectEnum;
371 else if (strcmp(name,"Nodal")==0) return NodalEnum;
372 else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
373 else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
374Index: ../trunk-jpl/src/c/shared/Elements/elements.h
375===================================================================
376--- ../trunk-jpl/src/c/shared/Elements/elements.h (revision 18967)
377+++ ../trunk-jpl/src/c/shared/Elements/elements.h (revision 18968)
378@@ -12,18 +12,21 @@
379 IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n);
380 IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat);
381 // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n);
382-IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec,
383- IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
384- IssmDouble* PrecipitationsLgm, IssmDouble* PrecipitationsPresentday,
385- IssmDouble* pdds, IssmDouble* pds,IssmDouble signorm, IssmDouble yts,
386- IssmDouble h, IssmDouble s, IssmDouble desfac,
387- IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm,
388- IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime);
389+IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec,
390+ IssmDouble* pdds, IssmDouble* pds,IssmDouble signorm, IssmDouble yts,
391+ IssmDouble h, IssmDouble s, IssmDouble desfac,IssmDouble s0t,
392+ IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm,
393+ IssmDouble TdiffTime,IssmDouble sealevTime,
394+ IssmDouble rho_water, IssmDouble rho_ice);
395 void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime,
396- IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime,
397- IssmDouble* PrecipitationsPresentday,
398- IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
399- IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
400+ IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime,
401+ IssmDouble* PrecipitationsPresentday,
402+ IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
403+ IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
404+void ComputeMungsmTemperaturePrecipitation(IssmDouble TdiffTime, IssmDouble PfacTime,
405+ IssmDouble* PrecipitationsLgm,IssmDouble* PrecipitationsPresentday,
406+ IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
407+ IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
408 IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.);
409 IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness);
410
411Index: ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp
412===================================================================
413--- ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 18967)
414+++ ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 18968)
415@@ -1,22 +1,18 @@
416 /* file: PddSurfaceMassBlance.cpp
417 Calculating the surface mass balance using the positive degree day method.
418+ Updating the precipitation and temperature to the new elevation
419 */
420
421 #include "./elements.h"
422 #include "../Numerics/numerics.h"
423
424-IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec,
425- IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
426- IssmDouble* PrecipitationsLgm, IssmDouble* PrecipitationsPresentday,
427- IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm,
428- IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac,
429- IssmDouble s0t,
430- IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm,
431- IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime){
432+IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec,
433+ IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm,
434+ IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac,
435+ IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps,IssmDouble rlapslgm,
436+ IssmDouble TdiffTime,IssmDouble sealevTime,
437+ IssmDouble rho_water,IssmDouble rho_ice){
438
439- // CURENTLY monthlytemperatures and monthlyprec ARE NOT USED HERE.
440- // THEY ARE REPLACE BY tdiffh No usage of deltO18 anymore
441-
442 // output:
443 IssmDouble B; // surface mass balance, melt+accumulation
444 int iqj,imonth;
445@@ -27,15 +23,15 @@
446 IssmDouble prect; // total precipitation during 1 year taking into account des. ef.
447 IssmDouble water; //water=rain + snowmelt
448 IssmDouble runoff; //meltwater only, does not include rain
449+ IssmDouble sconv; //rhow_rain/rhoi / 12 months
450
451 //IssmDouble sealev=0.; // degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
452 //IssmDouble Pfac=0.5,Tdiff=0.5;
453- IssmDouble rtlaps;
454+ IssmDouble rtlaps;
455 // IssmDouble lapser=6.5 // lapse rate
456 // IssmDouble desfac = 0.3; // desert elevation factor
457 // IssmDouble s0p=0.; // should be set to elevation from precip source
458 // IssmDouble s0t=0.; // should be set to elevation from temperature source
459- IssmDouble tdiffh;
460 IssmDouble st; // elevation between altitude of the temp record and current altitude
461 IssmDouble sp; // elevation between altitude of the prec record and current altitude
462 IssmDouble deselcut=1.0;
463@@ -48,7 +44,6 @@
464 IssmDouble DT = 0.02;
465 IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
466
467- IssmDouble premt; // monthly precipitation (m of ice equivalent)
468 IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
469 IssmDouble qm = 0.; // snow part of the precipitation
470 IssmDouble qmt = 0.; // precipitation without desertification effect adjustment
471@@ -72,6 +67,7 @@
472 IssmDouble fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice
473 IssmDouble pddtj, hmx2;
474
475+ sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
476
477 /*PDD constant*/
478 siglim = 2.5*signorm;
479@@ -89,8 +85,8 @@
480 rtlaps=TdiffTime*rlapslgm + (1.-TdiffTime)*rlaps; // lapse rate
481
482 /*********compute Surface temperature *******/
483- tdiffh = TdiffTime*( TemperaturesLgm[imonth] - TemperaturesPresentday[imonth] );
484- tstar = tdiffh + TemperaturesPresentday[imonth] - rtlaps *max(st,sealevTime*0.001);
485+ monthlytemperatures[imonth]=monthlytemperatures[imonth] - rtlaps *max(st,sealevTime*0.001);
486+ tstar = monthlytemperatures[imonth];
487 Tsurf = tstar*deltm+Tsurf;
488
489 /*********compute PD ****************/
490@@ -105,10 +101,8 @@
491 if (sp>0.0){q = exp(-desfac*sp);}
492 else {q = 1.0;}
493
494- premt =min(1.5, PrecipitationsPresentday[imonth] * pow(PrecipitationsLgm[imonth],PfacTime)); // [m/month]
495-
496- qmt= qmt + premt;
497- qmpt= q*premt;
498+ qmt= qmt + monthlyprec[imonth]*sconv; //*sconv to convert in m of ice equivalent per month
499+ qmpt= q*monthlyprec[imonth]*sconv;
500 qmp= qmp + qmpt;
501 qm= qm + qmpt*pd;
502
503@@ -123,6 +117,10 @@
504 pdd = pdd + pddsig*deltm;
505 frzndd = frzndd - (tstar-pddsig)*deltm;}
506 else{frzndd = frzndd - tstar*deltm; }
507+
508+ /*Assign output pointer*/
509+ *(monthlytemperatures+imonth) = monthlytemperatures[imonth];
510+ *(monthlyprec+imonth) = monthlyprec[imonth];
511 } // end of seasonal loop
512 //******************************************************************
513
514Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
515===================================================================
516--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 18967)
517+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 18968)
518@@ -21,7 +21,7 @@
519 int numoutputs,materialtype,smb_model,basalforcing_model;
520 char** requestedoutputs = NULL;
521 IssmDouble time;
522- bool isdelta18o;
523+ bool isdelta18o,ismungsm;
524
525 /*parameters for mass flux:*/
526 int mass_flux_num_profiles = 0;
527@@ -107,27 +107,31 @@
528 break;
529 case SMBpddEnum:
530 parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum));
531+ parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsmungsmEnum));
532 parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsDesfacEnum));
533 parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsS0pEnum));
534 parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsS0tEnum));
535 parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsRlapsEnum));
536 parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsRlapslgmEnum));
537 iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
538+ iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum);
539
540- iomodel->FetchData(&temp,&N,&M,SurfaceforcingsPfacEnum); _assert_(N==2);
541- for(i=0;i<M;i++) temp[M+i]=temp[M+i];
542- parameters->AddObject(new TransientParam(SurfaceforcingsPfacEnum,&temp[0],&temp[M],M));
543- iomodel->DeleteData(temp,SurfaceforcingsPfacEnum);
544+ if(ismungsm){
545+ iomodel->FetchData(&temp,&N,&M,SurfaceforcingsPfacEnum); _assert_(N==2);
546+ for(i=0;i<M;i++) temp[M+i]=temp[M+i];
547+ parameters->AddObject(new TransientParam(SurfaceforcingsPfacEnum,&temp[0],&temp[M],M));
548+ iomodel->DeleteData(temp,SurfaceforcingsPfacEnum);
549
550- iomodel->FetchData(&temp,&N,&M,SurfaceforcingsTdiffEnum); _assert_(N==2);
551- for(i=0;i<M;i++) temp[M+i]=temp[M+i];
552- parameters->AddObject(new TransientParam(SurfaceforcingsTdiffEnum,&temp[0],&temp[M],M));
553- iomodel->DeleteData(temp,SurfaceforcingsTdiffEnum);
554+ iomodel->FetchData(&temp,&N,&M,SurfaceforcingsTdiffEnum); _assert_(N==2);
555+ for(i=0;i<M;i++) temp[M+i]=temp[M+i];
556+ parameters->AddObject(new TransientParam(SurfaceforcingsTdiffEnum,&temp[0],&temp[M],M));
557+ iomodel->DeleteData(temp,SurfaceforcingsTdiffEnum);
558
559- iomodel->FetchData(&temp,&N,&M,SurfaceforcingsSealevEnum); _assert_(N==2);
560- for(i=0;i<M;i++) temp[M+i]=temp[M+i];
561- parameters->AddObject(new TransientParam(SurfaceforcingsSealevEnum,&temp[0],&temp[M],M));
562- iomodel->DeleteData(temp,SurfaceforcingsSealevEnum);
563+ iomodel->FetchData(&temp,&N,&M,SurfaceforcingsSealevEnum); _assert_(N==2);
564+ for(i=0;i<M;i++) temp[M+i]=temp[M+i];
565+ parameters->AddObject(new TransientParam(SurfaceforcingsSealevEnum,&temp[0],&temp[M],M));
566+ iomodel->DeleteData(temp,SurfaceforcingsSealevEnum);
567+ }
568
569 if(isdelta18o){
570 iomodel->Constant(&yts,ConstantsYtsEnum);
571Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h
572===================================================================
573--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 18967)
574+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 18968)
575@@ -11,6 +11,7 @@
576 void SurfaceMassBalancex(FemModel* femmodel);
577 void SmbGradientsx(FemModel* femmodel);
578 void Delta18oParameterizationx(FemModel* femmodel);
579+void MungsmtpParameterizationx(FemModel* femmodel);
580 void PositiveDegreeDayx(FemModel* femmodel);
581 void SmbHenningx(FemModel* femmodel);
582 void SmbComponentsx(FemModel* femmodel);
583Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
584===================================================================
585--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 18967)
586+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 18968)
587@@ -10,7 +10,7 @@
588
589 /*Intermediaties*/
590 int smb_model;
591- bool isdelta18o;
592+ bool isdelta18o,ismungsm;
593
594 /*First, get SMB model from parameters*/
595 femmodel->parameters->FindParam(&smb_model,SurfaceforcingsEnum);
596@@ -22,10 +22,15 @@
597 break;
598 case SMBpddEnum:
599 femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
600+ femmodel->parameters->FindParam(&ismungsm,SurfaceforcingsIsmungsmEnum);
601 if(isdelta18o){
602- if(VerboseSolution()) _printf0_(" call Delta18oParametrization module\n");
603+ if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n");
604 Delta18oParameterizationx(femmodel);
605 }
606+ if(ismungsm){
607+ if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n");
608+ MungsmtpParameterizationx(femmodel);
609+ }
610 if(VerboseSolution()) _printf0_(" call positive degree day module\n");
611 PositiveDegreeDayx(femmodel);
612 break;
613@@ -117,6 +122,14 @@
614 }
615
616 }/*}}}*/
617+void MungsmtpParameterizationx(FemModel* femmodel){/*{{{*/
618+
619+ for(int i=0;i<femmodel->elements->Size();i++){
620+ Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
621+ element->MungsmtpParameterization();
622+ }
623+
624+}/*}}}*/
625 void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/
626
627 // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
628@@ -143,6 +156,8 @@
629 IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
630 IssmDouble tstar; // monthly mean surface temp
631
632+ bool ismungsm;
633+
634 IssmDouble *pdds = NULL;
635 IssmDouble *pds = NULL;
636 Element *element = NULL;
637@@ -150,6 +165,9 @@
638 pdds=xNew<IssmDouble>(NPDMAX+1);
639 pds=xNew<IssmDouble>(NPDCMAX+1);
640
641+ // Get ismungsm parameter
642+ femmodel->parameters->FindParam(&ismungsm,SurfaceforcingsIsmungsmEnum);
643+
644 /* initialize PDD (creation of a lookup table)*/
645 tstep = 0.1;
646 tsint = tstep*0.5;
647@@ -204,9 +222,8 @@
648
649 for(i=0;i<femmodel->elements->Size();i++){
650 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
651- element->PositiveDegreeDay(pdds,pds,signorm);
652+ element->PositiveDegreeDay(pdds,pds,signorm,ismungsm);
653 }
654-
655 /*free ressouces: */
656 xDelete<IssmDouble>(pdds);
657 xDelete<IssmDouble>(pds);
658Index: ../trunk-jpl/src/c/Makefile.am
659===================================================================
660--- ../trunk-jpl/src/c/Makefile.am (revision 18967)
661+++ ../trunk-jpl/src/c/Makefile.am (revision 18968)
662@@ -219,6 +219,7 @@
663 ./shared/Elements/PrintArrays.cpp\
664 ./shared/Elements/PddSurfaceMassBalance.cpp\
665 ./shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp\
666+ ./shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp\
667 ./shared/Elements/DrainageFunctionWaterfraction.cpp\
668 ./shared/String/sharedstring.h\
669 ./shared/String/DescriptorIndex.cpp\
670Index: ../trunk-jpl/src/c/classes/Elements/Element.h
671===================================================================
672--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 18967)
673+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 18968)
674@@ -179,6 +179,7 @@
675 virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;
676 virtual void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum)=0;
677 virtual void Delta18oParameterization(void)=0;
678+ virtual void MungsmtpParameterization(void)=0;
679 virtual void ElementResponse(IssmDouble* presponse,int response_enum)=0;
680 virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
681 virtual int FiniteElement(void)=0;
682@@ -250,7 +251,7 @@
683 virtual void NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0;
684 virtual int NumberofNodesPressure(void)=0;
685 virtual int NumberofNodesVelocity(void)=0;
686- virtual void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm)=0;
687+ virtual void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm)=0;
688 virtual void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
689 virtual int PressureInterpolation()=0;
690 virtual void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0;
691Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
692===================================================================
693--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18967)
694+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18968)
695@@ -608,7 +608,6 @@
696 input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
697 input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
698 input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
699- PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]/yts; // converion in m/sec
700 }
701 }
702
703@@ -650,6 +649,69 @@
704 /*clean-up*/
705 delete gauss;
706 }
707+void Tria::MungsmtpParameterization(void){/*{{{*/
708+ /*Are we on the base? If not, return*/
709+ if(!IsOnBase()) return;
710+
711+ int i;
712+ IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
713+ IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
714+ IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
715+ IssmDouble tmp[NUMVERTICES];
716+ IssmDouble TdiffTime,PfacTime;
717+ IssmDouble time,yts;
718+ this->parameters->FindParam(&time,TimeEnum);
719+ this->parameters->FindParam(&yts,ConstantsYtsEnum);
720+
721+ /*Recover present day temperature and precipitation*/
722+ Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input);
723+ Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input2);
724+ Input* input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3);
725+ Input* input4=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input4);
726+ GaussTria* gauss=new GaussTria();
727+ for(int month=0;month<12;month++) {
728+ for(int iv=0;iv<NUMVERTICES;iv++) {
729+ gauss->GaussVertex(iv);
730+ input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
731+ input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
732+ input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
733+ input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
734+ }
735+ }
736+
737+ /*Recover interpolation parameters at time t*/
738+ this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
739+ this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
740+
741+ /*Compute the temperature and precipitation*/
742+ for(int iv=0;iv<NUMVERTICES;iv++){
743+ ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
744+ &PrecipitationsLgm[iv][0],&PrecipitationsPresentday[iv][0],
745+ &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
746+ &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
747+ }
748+
749+ /*Update inputs*/
750+ TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
751+ TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
752+ for (int imonth=0;imonth<12;imonth++) {
753+ for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
754+ TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
755+ NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
756+
757+ for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
758+ TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
759+ NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
760+ }
761+ NewTemperatureInput->Configure(this->parameters);
762+ NewPrecipitationInput->Configure(this->parameters);
763+
764+ this->inputs->AddInput(NewTemperatureInput);
765+ this->inputs->AddInput(NewPrecipitationInput);
766+
767+ /*clean-up*/
768+ delete gauss;
769+}
770 /*}}}*/
771 int Tria::EdgeOnBaseIndex(void){/*{{{*/
772
773@@ -2441,83 +2503,94 @@
774 return TriaRef::NumberofNodes(this->VelocityInterpolation());
775 }
776 /*}}}*/
777-void Tria::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){/*{{{*/
778-
779+void Tria::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/
780+
781+ int i;
782 IssmDouble agd[NUMVERTICES]; // surface mass balance
783 IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
784- IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
785- IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
786+ IssmDouble tmp[NUMVERTICES];
787 IssmDouble h[NUMVERTICES],s[NUMVERTICES];
788 IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
789 IssmDouble PfacTime,TdiffTime,sealevTime;
790- IssmDouble sconv; //rhow_rain/rhoi / 12 months
791
792 /*Get material parameters :*/
793 rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
794 rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
795- //rho_ice=matpar->GetRhoIce();
796- //rho_water=matpar->GetRhoFreshwater();
797
798- sconv=(rho_water/rho_ice)/12.; //to convert monbthly prec in m of ice equivalent per month
799+ /*Get some pdd parameters*/
800+ desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
801+ s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
802+ s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
803+ rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
804+ rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
805
806 /*Recover monthly temperatures and precipitation and Present day and LGm ones*/
807 Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
808 Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
809- Input* input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input3);
810- Input* input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input4);
811- Input* input5=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input5);
812- Input* input6=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input6);
813 GaussTria* gauss=new GaussTria();
814 IssmDouble time,yts;
815 this->parameters->FindParam(&time,TimeEnum);
816 this->parameters->FindParam(&yts,ConstantsYtsEnum);
817+
818+
819 for(int month=0;month<12;month++) {
820 for(int iv=0;iv<NUMVERTICES;iv++) {
821 gauss->GaussVertex(iv);
822 input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
823 monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius
824 input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
825- monthlyprec[iv][month]=monthlyprec[iv][month]*sconv*yts; // convertion in m/y
826- input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
827- TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius
828- input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
829- TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius
830- input5->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
831- PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month];
832- input6->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
833- PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month];
834 }
835 }
836
837- /*Recover Pfac, Tdiff and sealev at time t:*/
838- this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
839- this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
840- this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
841+ /*Recover Pfac, Tdiff and sealev at time t:
842+ This parameters are used to interpolate the temperature
843+ and precipitaton between PD and LGM when ismungsm==1 */
844+ if (ismungsm==1){
845+ this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
846+ this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
847+ }
848+ else {
849+ TdiffTime=0;
850+ sealevTime=0;
851+ }
852
853 /*Recover info at the vertices: */
854 GetInputListOnVertices(&h[0],ThicknessEnum);
855 GetInputListOnVertices(&s[0],SurfaceEnum);
856-
857- /*Get other pdd parameters*/
858- desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
859- s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
860- s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
861- rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
862- rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
863
864 /*measure the surface mass balance*/
865- for (int iv = 0; iv<NUMVERTICES; iv++){
866- agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
867- &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
868- &PrecipitationsLgm[iv][0], &PrecipitationsPresentday[iv][0],
869+ for (int iv = 0; iv<NUMVERTICES; iv++){
870+ agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
871 pdds, pds, signorm, yts, h[iv], s[iv],
872- desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime);
873+ desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime,
874+ rho_water,rho_ice);
875 }
876
877 /*Update inputs*/
878+ // TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
879+ // TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
880+ // for (int imonth=0;imonth<12;imonth++) {
881+ // for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
882+ // TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
883+ // NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
884+ //
885+ // for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
886+ // TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
887+ // NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
888+ // }
889+ // NewTemperatureInput->Configure(this->parameters);
890+ // NewPrecipitationInput->Configure(this->parameters);
891+
892+
893 this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
894+ // this->inputs->AddInput(NewTemperatureInput);
895+ // this->inputs->AddInput(NewPrecipitationInput);
896 // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
897
898+ //this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1);
899+ // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
900+ // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
901+
902 /*clean-up*/
903 delete gauss;
904 }
905Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
906===================================================================
907--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18967)
908+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18968)
909@@ -62,6 +62,7 @@
910 void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
911 void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
912 void Delta18oParameterization(void);
913+ void MungsmtpParameterization(void);
914 int EdgeOnBaseIndex();
915 void EdgeOnBaseIndices(int* pindex1,int* pindex);
916 int EdgeOnSurfaceIndex();
917@@ -108,7 +109,7 @@
918 int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
919 int NumberofNodesPressure(void);
920 int NumberofNodesVelocity(void);
921- void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
922+ void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm);
923 void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
924 int PressureInterpolation();
925 void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
926Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
927===================================================================
928--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18967)
929+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18968)
930@@ -607,7 +607,6 @@
931 input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
932 input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
933 input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
934- PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]/yts; // converion in m/sec
935 }
936 }
937
938@@ -652,6 +651,72 @@
939 /*clean-up*/
940 delete gauss;
941 }
942+void Penta::MungsmtpParameterization(void){/*{{{*/
943+ /*Are we on the base? If not, return*/
944+ if(!IsOnBase()) return;
945+
946+ int i;
947+ IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
948+ IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
949+ IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
950+ IssmDouble tmp[NUMVERTICES];
951+ IssmDouble TdiffTime,PfacTime;
952+ IssmDouble time,yts;
953+ this->parameters->FindParam(&time,TimeEnum);
954+ this->parameters->FindParam(&yts,ConstantsYtsEnum);
955+
956+ /*Recover present day temperature and precipitation*/
957+ Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input);
958+ Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input2);
959+ Input* input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3);
960+ Input* input4=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input4);
961+ GaussPenta* gauss=new GaussPenta();
962+ for(int month=0;month<12;month++) {
963+ for(int iv=0;iv<NUMVERTICES;iv++) {
964+ gauss->GaussVertex(iv);
965+ input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
966+ input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
967+ input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
968+ input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
969+ }
970+ }
971+
972+ /*Recover interpolation parameters at time t*/
973+ this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
974+ this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
975+
976+ /*Compute the temperature and precipitation*/
977+ for(int iv=0;iv<NUMVERTICES;iv++){
978+ ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
979+ &PrecipitationsLgm[iv][0],&PrecipitationsPresentday[iv][0],
980+ &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
981+ &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
982+ }
983+
984+ /*Update inputs*/
985+ TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
986+ TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
987+ for (int imonth=0;imonth<12;imonth++) {
988+ for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
989+ PentaInput* newmonthinput1 = new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
990+ NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
991+
992+ for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
993+ PentaInput* newmonthinput2 = new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
994+ NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
995+ }
996+ NewTemperatureInput->Configure(this->parameters);
997+ NewPrecipitationInput->Configure(this->parameters);
998+
999+ this->inputs->AddInput(NewTemperatureInput);
1000+ this->inputs->AddInput(NewPrecipitationInput);
1001+
1002+ this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
1003+ this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
1004+
1005+ /*clean-up*/
1006+ delete gauss;
1007+}
1008 /*}}}*/
1009 void Penta::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
1010
1011@@ -2081,32 +2146,30 @@
1012
1013 }
1014 /*}}}*/
1015-void Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){/*{{{*/
1016+void Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/
1017
1018+ int i;
1019 IssmDouble agd[NUMVERTICES]; // surface mass balance
1020 IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
1021- IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
1022- IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
1023+ IssmDouble tmp[NUMVERTICES];
1024 IssmDouble h[NUMVERTICES],s[NUMVERTICES];
1025 IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
1026 IssmDouble PfacTime,TdiffTime,sealevTime;
1027- IssmDouble sconv; //rhow_rain/rhoi / 12 months
1028
1029 /*Get material parameters :*/
1030 rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
1031 rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
1032- //rho_ice=matpar->GetRhoIce();
1033- //rho_water=matpar->GetRhoFreshwater();
1034
1035- sconv=(rho_water/rho_ice)/12.; //to convert monbthly prec in m of ice equivalent per month
1036+ /*Get some pdd parameters*/
1037+ desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
1038+ s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
1039+ s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
1040+ rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
1041+ rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
1042
1043 /*Recover monthly temperatures and precipitation*/
1044 Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
1045 Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
1046- Input* input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input3);
1047- Input* input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input4);
1048- Input* input5=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input5);
1049- Input* input6=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input6);
1050 GaussPenta* gauss=new GaussPenta();
1051 IssmDouble time,yts;
1052 this->parameters->FindParam(&time,TimeEnum);
1053@@ -2118,50 +2181,61 @@
1054 input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
1055 monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius
1056 input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
1057- monthlyprec[iv][month]=monthlyprec[iv][month]*sconv; // convertion to m/yr
1058- input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
1059- TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius
1060- input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
1061- TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius
1062- input5->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
1063- PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month];
1064- input6->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
1065- PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month];
1066 }
1067 }
1068
1069- /*Recover Pfac, Tdiff and sealev at time t:*/
1070- this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
1071- this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
1072- this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
1073+ /*Recover Pfac, Tdiff and sealev at time t:
1074+ This parameters are used to interpolate the temperature
1075+ and precipitaton between PD and LGM when ismungsm==1 */
1076+ if (ismungsm==1){
1077+ this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
1078+ this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
1079+ }
1080+ else {
1081+ TdiffTime=0;
1082+ sealevTime=0;
1083+ }
1084
1085 /*Recover info at the vertices: */
1086 GetInputListOnVertices(&h[0],ThicknessEnum);
1087 GetInputListOnVertices(&s[0],SurfaceEnum);
1088
1089- /*Get other pdd parameters*/
1090- desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
1091- s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
1092- s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
1093- rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
1094- rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
1095
1096 /*measure the surface mass balance*/
1097 for (int iv = 0; iv < NUMVERTICES; iv++){
1098- agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
1099- &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
1100- &PrecipitationsLgm[iv][0], &PrecipitationsPresentday[iv][0],
1101+ agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
1102 pdds,pds, signorm, yts, h[iv], s[iv],
1103- desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime);
1104+ desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime,
1105+ rho_water,rho_ice);
1106 }
1107
1108 /*Update inputs*/
1109+ // TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
1110+ // TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
1111+ // for (int imonth=0;imonth<12;imonth++) {
1112+ // for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
1113+ // PentaInput* newmonthinput1 = new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
1114+ // NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
1115+ //
1116+ // for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
1117+ // PentaInput* newmonthinput2 = new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
1118+ // NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
1119+ // }
1120+ // NewTemperatureInput->Configure(this->parameters);
1121+ // NewPrecipitationInput->Configure(this->parameters);
1122+
1123+
1124+
1125 this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
1126- //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
1127- this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1);
1128+ // this->inputs->AddInput(NewTemperatureInput);
1129+ // this->inputs->AddInput(NewPrecipitationInput);
1130+ // //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
1131+ this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1);
1132+ // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
1133+ // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
1134
1135- /*clean-up*/
1136- delete gauss;
1137+ /*clean-up*/
1138+ delete gauss;
1139 }
1140 /*}}}*/
1141 void Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
1142Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
1143===================================================================
1144--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 18967)
1145+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 18968)
1146@@ -58,6 +58,7 @@
1147 void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
1148 ElementMatrix* CreateBasalMassMatrix(void);
1149 void Delta18oParameterization(void);
1150+ void MungsmtpParameterization(void);
1151 void ElementResponse(IssmDouble* presponse,int response_enum);
1152 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
1153 int FiniteElement(void);
1154@@ -136,7 +137,7 @@
1155 int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
1156 int NumberofNodesPressure(void);
1157 int NumberofNodesVelocity(void);
1158- void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
1159+ void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm);
1160 void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
1161 int PressureInterpolation();
1162 void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
1163Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
1164===================================================================
1165--- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 18967)
1166+++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 18968)
1167@@ -53,6 +53,7 @@
1168 void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
1169 void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
1170 void Delta18oParameterization(void){_error_("not implemented yet");};
1171+ void MungsmtpParameterization(void){_error_("not implemented yet");};
1172 void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
1173 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
1174 void FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){_error_("not implemented yet");};
1175@@ -128,7 +129,7 @@
1176 void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
1177 int NumberofNodesPressure(void){_error_("not implemented yet");};
1178 int NumberofNodesVelocity(void){_error_("not implemented yet");};
1179- void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");};
1180+ void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){_error_("not implemented yet");};
1181 void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
1182 int PressureInterpolation(void){_error_("not implemented yet");};
1183 void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){_error_("not implemented yet");};
1184Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
1185===================================================================
1186--- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18967)
1187+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18968)
1188@@ -53,6 +53,7 @@
1189 void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
1190 void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
1191 void Delta18oParameterization(void){_error_("not implemented yet");};
1192+ void MungsmtpParameterization(void){_error_("not implemented yet");};
1193 IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");};
1194 void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
1195 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
1196@@ -134,7 +135,7 @@
1197 void NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
1198 int NumberofNodesPressure(void);
1199 int NumberofNodesVelocity(void);
1200- void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");};
1201+ void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){_error_("not implemented yet");};
1202 void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
1203 int PressureInterpolation(void);
1204 void ResetFSBasalBoundaryCondition(void);
1205Index: ../trunk-jpl/src/m/classes/SMBpdd.py
1206===================================================================
1207--- ../trunk-jpl/src/m/classes/SMBpdd.py (revision 18967)
1208+++ ../trunk-jpl/src/m/classes/SMBpdd.py (revision 18968)
1209@@ -25,11 +25,13 @@
1210 self.Tdiff = float('NaN')
1211 self.sealev = float('NaN')
1212 self.isdelta18o = 0
1213+ self.ismungsm = 0
1214 self.delta18o = float('NaN')
1215 self.delta18o_surface = float('NaN')
1216 self.temperatures_presentday = float('NaN')
1217 self.temperatures_lgm = float('NaN')
1218 self.precipitations_presentday = float('NaN')
1219+ self.precipitations_lgm = float('NaN')
1220
1221 #set defaults
1222 self.setdefaultparameters()
1223@@ -37,46 +39,63 @@
1224 def __repr__(self): # {{{
1225 string=" surface forcings parameters:"
1226
1227- string="%s\n%s"%(string,fielddisplay(self,'precipitation','surface precipitation [m/yr water eq]'))
1228+ string="%s\n%s"%(string,fielddisplay(self,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)'))
1229+ string="%s\n%s"%(string,fielddisplay(self,'ismungsm','is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)'))
1230 string="%s\n%s"%(string,fielddisplay(self,'desfac','desertification elevation factor (between 0 and 1, default is 0.5) [m]'))
1231- string="%s\n%s"%(string,fielddisplay(self,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)'))
1232 string="%s\n%s"%(string,fielddisplay(self,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'))
1233 string="%s\n%s"%(string,fielddisplay(self,'s0t','should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]'))
1234 string="%s\n%s"%(string,fielddisplay(self,'rlaps','present day lapse rate [degree/km]'))
1235 string="%s\n%s"%(string,fielddisplay(self,'rlapslgm','LGM lapse rate [degree/km]'))
1236- string="%s\n%s"%(string,fielddisplay(self,'Pfac','time interpolation parameter for precipitation, 1D (year)'))
1237- string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D (year)'))
1238- string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year)'))
1239- string="%s\n%s"%(string,fielddisplay(self,'monthlytemperatures','monthly surface temperatures [K], required if pdd is activated and delta18o not activated'))
1240- string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if pdd is activated and delta18o activated'))
1241- string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if pdd is activated and delta18o activated'))
1242- string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated'))
1243- string="%s\n%s"%(string,fielddisplay(self,'delta18o','delta18o, required if pdd is activated and delta18o activated'))
1244- string="%s\n%s"%(string,fielddisplay(self,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated'))
1245-
1246+ if not (self.isdelta18o and self.ismungsm):
1247+ string="%s\n%s"%(string,fielddisplay(self,'monthlytemperatures',['monthly surface temperatures [K], required if pdd is activated and delta18o not activated']))
1248+ string="%s\n%s"%(string,fielddisplay(self,'precipitation',['monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated']))
1249+ if self.isdelta18o:
1250+ string="%s\n%s"%(string,fielddisplay(self,'delta18o','delta18o, required if pdd is activated and delta18o activated'))
1251+ string="%s\n%s"%(string,fielddisplay(self,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated'))
1252+ string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated'))
1253+ string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated'))
1254+ string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'))
1255+ string="%s\n%s"%(string,fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'))
1256+ string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated'))
1257+ string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated'))
1258+ if self.ismungsm:
1259+ string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated'))
1260+ string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated'))
1261+ string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'))
1262+ string="%s\n%s"%(string,fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'))
1263+ string="%s\n%s"%(string,fielddisplay(self,'Pfac','time interpolation parameter for precipitation, 1D(year), required if mungsm is activated'))
1264+ string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated'))
1265+ string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated'))
1266 return string
1267 #}}}
1268 def extrude(self,md): # {{{
1269
1270- self.precipitation=project3d(md,'vector',self.precipitation,'type','node');
1271- self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node');
1272+ if not (self.isdelta18o and self.ismungsm):
1273+ self.precipitation=project3d(md,'vector',self.precipitation,'type','node')
1274+ self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node')
1275 if self.isdelta18o: self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node')
1276 if self.isdelta18o: self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node')
1277 if self.isdelta18o: self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node')
1278+ if self.isdelta18o: self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node')
1279+ if self.ismungsm: self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node')
1280+ if self.ismungsm: self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node')
1281+ if self.ismungsm: self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node')
1282+ if self.ismungsm: self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node')
1283 return self
1284 #}}}
1285 def initialize(self,md): # {{{
1286
1287- if numpy.all(numpy.isnan(self.precipitation)):
1288- self.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
1289- print " no SMBpdd.precipitation specified: values set as zero"
1290-
1291- return self
1292+ # if numpy.all(numpy.isnan(self.precipitation)):
1293+ # self.precipitation=numpy.zeros((md.mesh.numberofvertices,1))
1294+ # print " no SMBpdd.precipitation specified: values set as zero"
1295+ #
1296+ # return self
1297 #}}}
1298 def setdefaultparameters(self): # {{{
1299
1300 #pdd method not used in default mode
1301 self.isdelta18o = 0
1302+ self.ismungsm = 0
1303 self.desfac = 0.5
1304 self.s0p = 0.
1305 self.s0t = 0.
1306@@ -88,23 +107,32 @@
1307 def checkconsistency(self,md,solution,analyses): # {{{
1308
1309 if MasstransportAnalysisEnum() in analyses:
1310- md = checkfield(md,'fieldname','surfaceforcings.desfac','<=',1,'numel',[1]);
1311- md = checkfield(md,'fieldname','surfaceforcings.s0p','>=',0,'numel',[1]);
1312- md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',[1]);
1313- md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',[1]);
1314- md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',[1]);
1315- md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,numpy.nan]);
1316- md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan]);
1317- md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan]);
1318- if not self.isdelta18o:
1319- md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1)
1320- md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1)
1321- else:
1322+ md = checkfield(md,'fieldname','surfaceforcings.desfac','<=',1,'numel',[1])
1323+ md = checkfield(md,'fieldname','surfaceforcings.s0p','>=',0,'numel',[1])
1324+ md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',[1])
1325+ md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',[1])
1326+ md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',[1])
1327+
1328+ if not (self.isdelta18o and self.ismungsm):
1329+ md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','NaN',1)
1330+ md = checkfield(md,'fieldname','surfaceforcings.precipitation','NaN',1)
1331+ elif self.isdelta18o:
1332 md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1)
1333 md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1)
1334 md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1)
1335 md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1)
1336 md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1)
1337+ md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1)
1338+ md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan])
1339+ md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan])
1340+ elif self.ismungsm:
1341+ md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1)
1342+ md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1)
1343+ md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1)
1344+ md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1)
1345+ md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,numpy.nan])
1346+ md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan])
1347+ md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan])
1348
1349 return md
1350 # }}}
1351@@ -112,24 +140,34 @@
1352
1353 yts=365.0*24.0*3600.0
1354
1355- WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer');
1356+ WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer')
1357
1358- WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1)
1359 WriteData(fid,'object',self,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean')
1360- WriteData(fid,'object',self,'class','surfaceforcings','fieldname','desfac','format','Double');
1361- WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0p','format','Double');
1362- WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0t','format','Double');
1363- WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlaps','format','Double');
1364- WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlapslgm','format','Double');
1365- WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1);
1366- WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1);
1367- WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1);
1368- if self.isdelta18o:
1369+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','ismungsm','format','Boolean')
1370+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','desfac','format','Double')
1371+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0p','format','Double')
1372+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0t','format','Double')
1373+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlaps','format','Double')
1374+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlapslgm','format','Double')
1375+
1376+ if not (self.isdelta18o and self.ismungsm):
1377+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1)
1378+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1)
1379+ elif self.isdelta18o:
1380 WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1)
1381 WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1)
1382 WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1)
1383+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1)
1384 WriteData(fid,'object',self,'class','surfaceforcings','fieldname','delta18o_surface','format','DoubleMat','mattype',1)
1385 WriteData(fid,'object',self,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1)
1386- else:
1387- WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)
1388+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1)
1389+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1)
1390+ elif self.ismungsm:
1391+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1)
1392+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1)
1393+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1)
1394+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1)
1395+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1)
1396+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1)
1397+ WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1)
1398 # }}}
1399Index: ../trunk-jpl/src/m/classes/SMBpdd.m
1400===================================================================
1401--- ../trunk-jpl/src/m/classes/SMBpdd.m (revision 18967)
1402+++ ../trunk-jpl/src/m/classes/SMBpdd.m (revision 18968)
1403@@ -16,6 +16,7 @@
1404 Tdiff = NaN;
1405 sealev = NaN;
1406 isdelta18o = 0;
1407+ ismungsm = 0;
1408 delta18o = NaN;
1409 delta18o_surface = NaN;
1410 temperatures_presentday = NaN;
1411@@ -33,27 +34,30 @@
1412 end
1413 end % }}}
1414 function self = extrude(self,md) % {{{
1415-
1416- self.precipitation=project3d(md,'vector',self.precipitation,'type','node');
1417- self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node');
1418+ if(self.isdelta18o==0 & self.ismungsm==0),self.precipitation=project3d(md,'vector',self.precipitation,'type','node');end
1419+ if(self.isdelta18o==0 & self.ismungsm==0),self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node');end
1420 if(self.isdelta18o),self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node');end
1421 if(self.isdelta18o),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end
1422 if(self.isdelta18o),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end
1423 if(self.isdelta18o),self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node');end
1424+ if(self.ismungsm),self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node');end
1425+ if(self.ismungsm),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end
1426+ if(self.ismungsm),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end
1427+ if(self.ismungsm),self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node');end
1428
1429-
1430 end % }}}
1431 function self = initialize(self,md) % {{{
1432+
1433+ % if isnan(self.precipitation),
1434+ % self.precipitation=zeros(md.mesh.numberofvertices,1);
1435+ % disp(' no SMBpdd.precipitation specified: values set as zero');
1436+ % end
1437
1438- if isnan(self.precipitation),
1439- self.precipitation=zeros(md.mesh.numberofvertices,1);
1440- disp(' no SMBpdd.precipitation specified: values set as zero');
1441- end
1442-
1443 end % }}}
1444 function obj = setdefaultparameters(obj) % {{{
1445
1446 obj.isdelta18o = 0;
1447+ obj.ismungsm = 0;
1448 obj.desfac = 0.5;
1449 obj.s0p = 0;
1450 obj.s0t = 0;
1451@@ -69,44 +73,61 @@
1452 md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',1);
1453 md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',1);
1454 md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',1);
1455- md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,NaN]);
1456- md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,NaN]);
1457- md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,NaN]);
1458- if(obj.isdelta18o==0)
1459- md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1);
1460- md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1);
1461- else
1462- md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1);
1463- md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1);
1464- md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1465- md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1466- md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1467- md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1468+ if(obj.isdelta18o==0 & obj.ismungsm==0)
1469+ md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1);
1470+ md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1);
1471+ elseif(obj.isdelta18o==1)
1472+ md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1);
1473+ md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1);
1474+ md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1475+ md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1476+ md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1477+ md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1478+ md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1);
1479+ md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1);
1480+ elseif(obj.ismungsm==1)
1481+ md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1482+ md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1483+ md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1484+ md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1);
1485+ md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,NaN]);
1486+ md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1);
1487+ md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1);
1488 end
1489- end
1490+ end
1491 end % }}}
1492 function disp(obj) % {{{
1493 disp(sprintf(' surface forcings parameters:'));
1494
1495 disp(sprintf('\n PDD and deltaO18 parameters:'));
1496 fielddisplay(obj,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)');
1497+ fielddisplay(obj,'ismungsm','is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)');
1498 fielddisplay(obj,'desfac','desertification elevation factor (between 0 and 1, default is 0.5) [m]');
1499 fielddisplay(obj,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]');
1500 fielddisplay(obj,'s0t','should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]');
1501 fielddisplay(obj,'rlaps','present day lapse rate [degree/km]');
1502 fielddisplay(obj,'rlapslgm','LGM lapse rate [degree/km]');
1503- fielddisplay(obj,'Pfac','time interpolation parameter for precipitation, 1D(year)');
1504- fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year)');
1505- fielddisplay(obj,'sealev','sea level [m], 1D(year)');
1506- fielddisplay(obj,'monthlytemperatures',['CURRENTLY NOT USED monthly surface temperatures [K], required if pdd is activated and delta18o not activated']);
1507- fielddisplay(obj,'precipitation',['CURRENTLY NOT USED monthly surface precipitation [m/yr water eq]']);
1508- fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if pdd is activated and delta18o activated');
1509- fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if pdd is activated and delta18o activated');
1510- fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated');
1511- fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated');
1512- fielddisplay(obj,'delta18o','delta18o, required if pdd is activated and delta18o activated');
1513- fielddisplay(obj,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated');
1514-
1515+ if(obj.isdelta18o==0 & obj.ismungsm==0)
1516+ fielddisplay(obj,'monthlytemperatures',['monthly surface temperatures [K], required if pdd is activated and delta18o not activated']);
1517+ fielddisplay(obj,'precipitation',['monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated']);
1518+ elseif(obj.isdelta18o==1)
1519+ fielddisplay(obj,'delta18o','delta18o, required if pdd is activated and delta18o activated');
1520+ fielddisplay(obj,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated');
1521+ fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated');
1522+ fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated');
1523+ fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
1524+ fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
1525+ fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated');
1526+ fielddisplay(obj,'sealev','sea level [m], 1D(year), required if mungsm is activated');
1527+ elseif(obj.ismungsm==1)
1528+ fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated');
1529+ fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated');
1530+ fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
1531+ fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated');
1532+ fielddisplay(obj,'Pfac','time interpolation parameter for precipitation, 1D(year), required if mungsm is activated');
1533+ fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated');
1534+ fielddisplay(obj,'sealev','sea level [m], 1D(year), required if mungsm is activated');
1535+ end
1536 end % }}}
1537 function marshall(obj,md,fid) % {{{
1538
1539@@ -114,25 +135,35 @@
1540
1541 WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer');
1542
1543- WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
1544+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean');
1545+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','ismungsm','format','Boolean');
1546 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','desfac','format','Double');
1547 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','s0p','format','Double');
1548 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','s0t','format','Double');
1549 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','rlaps','format','Double');
1550 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','rlapslgm','format','Double');
1551- WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1);
1552- WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1);
1553- WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1);
1554- WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean');
1555- if obj.isdelta18o
1556+
1557+ if(obj.isdelta18o==0 & obj.ismungsm==0)
1558+ %WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1);
1559+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
1560+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
1561+ elseif obj.isdelta18o
1562 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1);
1563 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1);
1564 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1);
1565 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1);
1566 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','delta18o_surface','format','DoubleMat','mattype',1);
1567 WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1);
1568- else
1569- WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
1570+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1);
1571+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1);
1572+ elseif obj.ismungsm
1573+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1);
1574+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1);
1575+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1);
1576+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1);
1577+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1);
1578+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1);
1579+ WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1);
1580 end
1581 end % }}}
1582 end
1583Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py
1584===================================================================
1585--- ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18967)
1586+++ ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18968)
1587@@ -341,6 +341,7 @@
1588 def SurfaceforcingsDelta18oEnum(): return StringToEnum("SurfaceforcingsDelta18o")[0]
1589 def SurfaceforcingsDelta18oSurfaceEnum(): return StringToEnum("SurfaceforcingsDelta18oSurface")[0]
1590 def SurfaceforcingsIsdelta18oEnum(): return StringToEnum("SurfaceforcingsIsdelta18o")[0]
1591+def SurfaceforcingsIsmungsmEnum(): return StringToEnum("SurfaceforcingsIsmungsm")[0]
1592 def SurfaceforcingsPrecipitationsPresentdayEnum(): return StringToEnum("SurfaceforcingsPrecipitationsPresentday")[0]
1593 def SurfaceforcingsPrecipitationsLgmEnum(): return StringToEnum("SurfaceforcingsPrecipitationsLgm")[0]
1594 def SurfaceforcingsTemperaturesPresentdayEnum(): return StringToEnum("SurfaceforcingsTemperaturesPresentday")[0]
Note: See TracBrowser for help on using the repository browser.