source: issm/oecreview/Archive/24307-24683/ISSM-24682-24683.diff@ 24684

Last change on this file since 24684 was 24684, checked in by Mathieu Morlighem, 5 years ago

CHG: added new review
File size: 35.6 KB
RevLine 
[24684]1Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
2===================================================================
3--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24682)
4+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24683)
5@@ -714,6 +714,7 @@
6 SmbDziniEnum,
7 SmbEAirEnum,
8 SmbECEnum,
9+ SmbECDtEnum,
10 SmbECiniEnum,
11 SmbElaEnum,
12 SmbEvaporationEnum,
13@@ -733,6 +734,7 @@
14 SmbMeanSHFEnum,
15 SmbMeanULWEnum,
16 SmbMeltEnum,
17+ SmbMInitnum,
18 SmbMonthlytemperaturesEnum,
19 SmbNetLWEnum,
20 SmbNetSWEnum,
21@@ -771,6 +773,7 @@
22 SmbVmeanEnum,
23 SmbVzEnum,
24 SmbWEnum,
25+ SmbWAddEnum,
26 SmbWiniEnum,
27 SmbZMaxEnum,
28 SmbZMinEnum,
29Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
30===================================================================
31--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24682)
32+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24683)
33@@ -719,6 +719,7 @@
34 case SmbDziniEnum : return "SmbDzini";
35 case SmbEAirEnum : return "SmbEAir";
36 case SmbECEnum : return "SmbEC";
37+ case SmbECDtEnum : return "SmbECDt";
38 case SmbECiniEnum : return "SmbECini";
39 case SmbElaEnum : return "SmbEla";
40 case SmbEvaporationEnum : return "SmbEvaporation";
41@@ -776,6 +777,7 @@
42 case SmbVmeanEnum : return "SmbVmean";
43 case SmbVzEnum : return "SmbVz";
44 case SmbWEnum : return "SmbW";
45+ case SmbWAddEnum : return "SmbWAdd";
46 case SmbWiniEnum : return "SmbWini";
47 case SmbZMaxEnum : return "SmbZMax";
48 case SmbZMinEnum : return "SmbZMin";
49Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
50===================================================================
51--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24682)
52+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24683)
53@@ -734,6 +734,7 @@
54 else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
55 else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
56 else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
57+ else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
58 else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
59 else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
60 else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
61@@ -750,11 +751,11 @@
62 else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum;
63 else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum;
64 else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum;
65- else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
66 else stage=7;
67 }
68 if(stage==7){
69- if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum;
70+ if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
71+ else if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum;
72 else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
73 else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
74 else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum;
75@@ -794,6 +795,7 @@
76 else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
77 else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
78 else if (strcmp(name,"SmbW")==0) return SmbWEnum;
79+ else if (strcmp(name,"SmbWAdd")==0) return SmbWAddEnum;
80 else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
81 else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum;
82 else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum;
83@@ -872,12 +874,12 @@
84 else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
85 else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
86 else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
87- else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
88- else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
89 else stage=8;
90 }
91 if(stage==8){
92- if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
93+ if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
94+ else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
95+ else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
96 else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
97 else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
98 else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
99@@ -995,12 +997,12 @@
100 else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
101 else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum;
102 else if (strcmp(name,"IntInput2")==0) return IntInput2Enum;
103- else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
104- else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
105 else stage=9;
106 }
107 if(stage==9){
108- if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
109+ if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
110+ else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
111+ else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
112 else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
113 else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
114 else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
115@@ -1118,12 +1120,12 @@
116 else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum;
117 else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum;
118 else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum;
119- else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
120- else if (strcmp(name,"Indexed")==0) return IndexedEnum;
121 else stage=10;
122 }
123 if(stage==10){
124- if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
125+ if (strcmp(name,"Incremental")==0) return IncrementalEnum;
126+ else if (strcmp(name,"Indexed")==0) return IndexedEnum;
127+ else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
128 else if (strcmp(name,"IntInput")==0) return IntInputEnum;
129 else if (strcmp(name,"ElementInput2")==0) return ElementInput2Enum;
130 else if (strcmp(name,"SegInput2")==0) return SegInput2Enum;
131@@ -1241,12 +1243,12 @@
132 else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
133 else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum;
134 else if (strcmp(name,"Regular")==0) return RegularEnum;
135- else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
136- else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
137 else stage=11;
138 }
139 if(stage==11){
140- if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
141+ if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
142+ else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
143+ else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
144 else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
145 else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
146 else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
147Index: ../trunk-jpl/src/c/shared/Enum/Enum.vim
148===================================================================
149--- ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24682)
150+++ ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24683)
151@@ -148,6 +148,7 @@
152 syn keyword cConstant FlowequationIsSSAEnum
153 syn keyword cConstant FrictionCouplingEnum
154 syn keyword cConstant FrictionDeltaEnum
155+syn keyword cConstant FrictionEffectivePressureLimitEnum
156 syn keyword cConstant FrictionFEnum
157 syn keyword cConstant FrictionGammaEnum
158 syn keyword cConstant FrictionLawEnum
159@@ -716,6 +717,7 @@
160 syn keyword cConstant SmbDziniEnum
161 syn keyword cConstant SmbEAirEnum
162 syn keyword cConstant SmbECEnum
163+syn keyword cConstant SmbECDtEnum
164 syn keyword cConstant SmbECiniEnum
165 syn keyword cConstant SmbElaEnum
166 syn keyword cConstant SmbEvaporationEnum
167@@ -773,6 +775,7 @@
168 syn keyword cConstant SmbVmeanEnum
169 syn keyword cConstant SmbVzEnum
170 syn keyword cConstant SmbWEnum
171+syn keyword cConstant SmbWAddEnum
172 syn keyword cConstant SmbWiniEnum
173 syn keyword cConstant SmbZMaxEnum
174 syn keyword cConstant SmbZMinEnum
175@@ -1334,7 +1337,6 @@
176 syn keyword cType Cfsurfacelogvel
177 syn keyword cType Cfsurfacesquare
178 syn keyword cType Channel
179-syn keyword cType classes
180 syn keyword cType Constraint
181 syn keyword cType Constraints
182 syn keyword cType Contour
183@@ -1341,8 +1343,8 @@
184 syn keyword cType Contours
185 syn keyword cType ControlInput2
186 syn keyword cType Covertree
187+syn keyword cType DataSetParam
188 syn keyword cType DatasetInput2
189-syn keyword cType DataSetParam
190 syn keyword cType Definition
191 syn keyword cType DependentObject
192 syn keyword cType DoubleMatArrayParam
193@@ -1354,8 +1356,8 @@
194 syn keyword cType ElementHook
195 syn keyword cType ElementInput2
196 syn keyword cType ElementMatrix
197+syn keyword cType ElementVector
198 syn keyword cType Elements
199-syn keyword cType ElementVector
200 syn keyword cType ExponentialVariogram
201 syn keyword cType ExternalResult
202 syn keyword cType FemModel
203@@ -1362,12 +1364,11 @@
204 syn keyword cType FileParam
205 syn keyword cType Friction
206 syn keyword cType Gauss
207-syn keyword cType GaussianVariogram
208-syn keyword cType gaussobjects
209 syn keyword cType GaussPenta
210 syn keyword cType GaussSeg
211 syn keyword cType GaussTetra
212 syn keyword cType GaussTria
213+syn keyword cType GaussianVariogram
214 syn keyword cType GenericExternalResult
215 syn keyword cType GenericOption
216 syn keyword cType GenericParam
217@@ -1382,7 +1383,6 @@
218 syn keyword cType IoModel
219 syn keyword cType IssmDirectApplicInterface
220 syn keyword cType IssmParallelDirectApplicInterface
221-syn keyword cType krigingobjects
222 syn keyword cType Load
223 syn keyword cType Loads
224 syn keyword cType Masscon
225@@ -1393,7 +1393,6 @@
226 syn keyword cType Matestar
227 syn keyword cType Matice
228 syn keyword cType Matlitho
229-syn keyword cType matrixobjects
230 syn keyword cType MatrixParam
231 syn keyword cType Misfit
232 syn keyword cType Moulin
233@@ -1406,8 +1405,8 @@
234 syn keyword cType Observation
235 syn keyword cType Observations
236 syn keyword cType Option
237+syn keyword cType OptionUtilities
238 syn keyword cType Options
239-syn keyword cType OptionUtilities
240 syn keyword cType Param
241 syn keyword cType Parameters
242 syn keyword cType Pengrid
243@@ -1421,12 +1420,12 @@
244 syn keyword cType Radar
245 syn keyword cType Regionaloutput
246 syn keyword cType Results
247+syn keyword cType RiftStruct
248 syn keyword cType Riftfront
249-syn keyword cType RiftStruct
250 syn keyword cType Seg
251 syn keyword cType SegInput2
252+syn keyword cType SegRef
253 syn keyword cType Segment
254-syn keyword cType SegRef
255 syn keyword cType SpcDynamic
256 syn keyword cType SpcStatic
257 syn keyword cType SpcTransient
258@@ -1445,6 +1444,10 @@
259 syn keyword cType VectorParam
260 syn keyword cType Vertex
261 syn keyword cType Vertices
262+syn keyword cType classes
263+syn keyword cType gaussobjects
264+syn keyword cType krigingobjects
265+syn keyword cType matrixobjects
266 syn keyword cType AdjointBalancethickness2Analysis
267 syn keyword cType AdjointBalancethicknessAnalysis
268 syn keyword cType AdjointHorizAnalysis
269@@ -1463,8 +1466,8 @@
270 syn keyword cType ExtrudeFromTopAnalysis
271 syn keyword cType FreeSurfaceBaseAnalysis
272 syn keyword cType FreeSurfaceTopAnalysis
273+syn keyword cType GLheightadvectionAnalysis
274 syn keyword cType GiaIvinsAnalysis
275-syn keyword cType GLheightadvectionAnalysis
276 syn keyword cType HydrologyDCEfficientAnalysis
277 syn keyword cType HydrologyDCInefficientAnalysis
278 syn keyword cType HydrologyGlaDSAnalysis
279Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
280===================================================================
281--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24682)
282+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24683)
283@@ -5,6 +5,8 @@
284 #include "./SurfaceMassBalancex.h"
285 #include "../../shared/shared.h"
286 #include "../../toolkits/toolkits.h"
287+#include "../modules.h"
288+#include "../../classes/Inputs2/TransientInput2.h"
289
290 const double Pi = 3.141592653589793;
291 const double CtoK = 273.15; // Kelvin to Celcius conversion/ice melt. point T in K
292@@ -31,9 +33,51 @@
293
294 void Gembx(FemModel* femmodel){ /*{{{*/
295
296- for(int i=0;i<femmodel->elements->Size();i++){
297- Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
298- element->SmbGemb();
299+ int count=0;
300+ IssmDouble time,dt,finaltime,starttime;
301+ IssmDouble timeclim=0.0;
302+ IssmDouble t,smb_dt;
303+ IssmDouble delta;
304+ bool isclimatology=false;
305+
306+ femmodel->parameters->FindParam(&time,TimeEnum); /*transient core time at which we run the smb core*/
307+ femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*transient core time step*/
308+ femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
309+ femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
310+ femmodel->parameters->FindParam(&smb_dt,SmbDtEnum); /*time period for the smb solution, usually smaller than the glaciological dt*/
311+ femmodel->parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
312+
313+ //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT.
314+ //go back to time - deltaT:
315+ time-=dt;
316+
317+ IssmDouble timeinputs = time;
318+
319+ /*Start loop: */
320+ count=1;
321+ for (t=time;t<time+dt;t=t+smb_dt){
322+
323+ for(int i=0;i<femmodel->elements->Size();i++){
324+ Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
325+
326+ timeclim=time;
327+ if (isclimatology){
328+ //If this is a climatology, we need to repeat the forcing after the final time
329+ TransientInput2* Ta_input_tr = element->inputs2->GetTransientInput(SmbTaEnum); _assert_(Ta_input_tr);
330+
331+ /*Get temperature climatology value*/
332+ int offsetend = Ta_input_tr->GetTimeInputOffset(finaltime);
333+ IssmDouble time0 = Ta_input_tr->GetTimeByOffset(-1);
334+ IssmDouble timeend = Ta_input_tr->GetTimeByOffset(offsetend);
335+ if (time>time0 & timeend>time0){
336+ delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
337+ timeclim=time0+delta;
338+ }
339+ }
340+ timeinputs = t-time+timeclim;
341+ element->SmbGemb(timeinputs,count);
342+ }
343+ count=count+1;
344 }
345
346 } /*}}}*/
347@@ -1044,7 +1088,7 @@
348
349 // SWs and SWss coefficients need to be better constranted. Greuell
350 // and Konzelmann 1994 used SWs = 0.36 and SWss = 0.64 as this the
351- // the // of SW radiation with wavelengths > and < 800 nm
352+ // the % of SW radiation with wavelengths > and < 800 nm
353 // respectively. This, however, may not account for the fact that
354 // the albedo of wavelengths > 800 nm has a much lower albedo.
355
356@@ -2106,9 +2150,9 @@
357 if(V< 0.01-Dtol)V=0.01;
358
359 // calculate the Bulk Richardson Number (Ri)
360- Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2));
361+ Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
362
363- // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
364+ // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
365
366 // do not allow Ri to exceed 0.19
367 Ri = min(Ri, 0.19);
368Index: ../trunk-jpl/src/c/classes/Elements/Element.h
369===================================================================
370--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 24682)
371+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 24683)
372@@ -171,7 +171,7 @@
373 void SetIntInput(Inputs2* inputs2,int enum_in,int value);
374 void SmbSemic();
375 int Sid();
376- void SmbGemb();
377+ void SmbGemb(IssmDouble timeinputs, int count);
378 void StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input);
379 void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input,Input2* vz_input);
380 void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input);
381Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
382===================================================================
383--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24682)
384+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24683)
385@@ -561,7 +561,7 @@
386 this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum);
387 this->parameters->FindParam(&isPrecipScaled,SmbIsprecipscaledEnum);
388
389- /*Recover present day temperature and precipitation*/
390+ /*Recover present day temperature and precipitation*/
391 DatasetInput2 *dinput3 = NULL;
392 DatasetInput2 *dinput4 = NULL;
393 int offset_t,offset_p,N;
394@@ -3587,7 +3587,7 @@
395
396 }
397 /*}}}*/
398-void Element::SmbGemb(){/*{{{*/
399+void Element::SmbGemb(IssmDouble timeinputs, int count){/*{{{*/
400
401 /*only compute SMB at the surface: */
402 if (!IsOnSurface()) return;
403@@ -3604,9 +3604,6 @@
404 IssmDouble Vmean=0.0;
405 IssmDouble C=0.0;
406 IssmDouble Tz,Vz=0.0;
407- IssmDouble time,dt,starttime,finaltime;
408- IssmDouble timeclim=0.0;
409- IssmDouble t,smb_dt;
410 IssmDouble yts;
411 IssmDouble Ta=0.0;
412 IssmDouble V=0.0;
413@@ -3617,6 +3614,7 @@
414 IssmDouble pAir=0.0;
415 IssmDouble teValue=1.0;
416 IssmDouble aValue=0.0;
417+ IssmDouble dt,time,smb_dt;
418 int aIdx=0;
419 int denIdx=0;
420 int dsnowIdx=0;
421@@ -3626,13 +3624,12 @@
422 IssmDouble shf=0.0;
423 IssmDouble dayEC=0.0;
424 IssmDouble initMass=0.0;
425- IssmDouble sumR=0.0;
426- IssmDouble sumM=0.0;
427- IssmDouble sumEC=0.0;
428- IssmDouble sumP=0.0;
429- IssmDouble sumW=0.0;
430- IssmDouble sumMassAdd=0.0;
431- IssmDouble sumdz_add=0.0;
432+ IssmDouble sumR=0.0;
433+ IssmDouble sumM=0.0;
434+ IssmDouble sumEC=0.0;
435+ IssmDouble sumP=0.0;
436+ IssmDouble sumW=0.0;
437+ IssmDouble sumMassAdd=0.0;
438 IssmDouble fac=0.0;
439 IssmDouble sumMass=0.0;
440 IssmDouble dMass=0.0;
441@@ -3641,8 +3638,6 @@
442 IssmDouble init_scaling=0.0;
443 IssmDouble thermo_scaling=1.0;
444 IssmDouble adThresh=1023.0;
445- int offsetend=-1;
446- IssmDouble time0, timeend, delta;
447 /*}}}*/
448 /*Output variables:{{{ */
449 IssmDouble* dz=NULL;
450@@ -3675,7 +3670,6 @@
451 IssmDouble* aini = NULL;
452 IssmDouble* Tini = NULL;
453 int m=0;
454- int count=0;
455 /*}}}*/
456
457 /*Retrieve material properties and parameters:{{{ */
458@@ -3686,8 +3680,6 @@
459 parameters->FindParam(&time,TimeEnum); /*transient core time at which we run the smb core*/
460 parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*transient core time step*/
461 parameters->FindParam(&yts,ConstantsYtsEnum);
462- parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
463- parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
464 parameters->FindParam(&smb_dt,SmbDtEnum); /*time period for the smb solution, usually smaller than the glaciological dt*/
465 parameters->FindParam(&aIdx,SmbAIdxEnum);
466 parameters->FindParam(&denIdx,SmbDenIdxEnum);
467@@ -3697,7 +3689,6 @@
468 parameters->FindParam(&t0wet,SmbT0wetEnum);
469 parameters->FindParam(&t0dry,SmbT0dryEnum);
470 parameters->FindParam(&K,SmbKEnum);
471- parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
472 parameters->FindParam(&isgraingrowth,SmbIsgraingrowthEnum);
473 parameters->FindParam(&isalbedo,SmbIsalbedoEnum);
474 parameters->FindParam(&isshortwave,SmbIsshortwaveEnum);
475@@ -3722,8 +3713,6 @@
476 Input2 *C_input = this->GetInput2(SmbCEnum); _assert_(C_input);
477 Input2 *Tz_input = this->GetInput2(SmbTzEnum); _assert_(Tz_input);
478 Input2 *Vz_input = this->GetInput2(SmbVzEnum); _assert_(Vz_input);
479- Input2 *teValue_input = this->GetInput2(SmbTeValueEnum); _assert_(teValue_input);
480- Input2 *aValue_input = this->GetInput2(SmbAValueEnum); _assert_(aValue_input);
481 Input2 *EC_input = NULL;
482
483 /*Retrieve input values:*/
484@@ -3741,8 +3730,6 @@
485 C_input->GetInputValue(&C,gauss);
486 Tz_input->GetInputValue(&Tz,gauss);
487 Vz_input->GetInputValue(&Vz,gauss);
488- teValue_input->GetInputValue(&teValue,gauss);
489- aValue_input->GetInputValue(&aValue,gauss);
490 /*}}}*/
491
492 /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/
493@@ -3762,7 +3749,7 @@
494 EC_input = this->GetInput2(SmbECiniEnum); _assert_(EC_input);
495 EC_input->GetInputAverage(&EC);
496
497- /*Retrive the correct value of m (without the zeroes at the end)*/
498+ /*Retrieve the correct value of m (without the zeroes at the end)*/
499 this->GetInput2Value(&m,SmbSizeiniEnum);
500
501 if(m==2){ //Snow properties are initialized with default values. Vertical grid has to be initialized too
502@@ -3815,9 +3802,8 @@
503 this->inputs2->GetArray(SmbWEnum,this->lid,&W,&m);
504 this->inputs2->GetArray(SmbAEnum,this->lid,&a,&m);
505 this->inputs2->GetArray(SmbTEnum,this->lid,&T,&m);
506- EC_input = this->GetInput2(SmbECEnum); _assert_(EC_input);
507+ EC_input = this->GetInput2(SmbECDtEnum); _assert_(EC_input);
508 EC_input->GetInputAverage(&EC);
509- EC=EC*dt*rho_ice;
510
511 //fixed lower temperature bounday condition - T is fixed
512 _assert_(m>0);
513@@ -3829,152 +3815,167 @@
514
515 // initialize cumulative variables
516 sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
517- sumdz_add=0;
518
519 //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT.
520- //go back to time - deltaT:
521- time-=dt;
522+ //go back to time - deltaT:
523+ time-=dt;
524
525- timeclim=time;
526- if (isclimatology){
527- //If this is a climatology, we need to repeat the forcing after the final time
528- TransientInput2* Ta_input_tr = this->inputs2->GetTransientInput(SmbTaEnum); _assert_(Ta_input_tr);
529- offsetend = Ta_input_tr->GetTimeInputOffset(finaltime);
530- time0 = Ta_input_tr->GetTimeByOffset(-1);
531- timeend = Ta_input_tr->GetTimeByOffset(offsetend);
532- if (time>time0 & timeend>time0){
533- delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
534- timeclim=time0+delta;
535- }
536- }
537+ if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << timeinputs/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n");
538
539- /*Start loop: */
540- count=1;
541- for (t=time;t<time+dt;t=t+smb_dt){
542+ /*Get daily accumulated inputs {{{*/
543+ if (count>1){
544+ Input2 *sumEC_input = this->GetInput2(SmbECEnum); _assert_(sumEC_input);
545+ Input2 *sumM_input = this->GetInput2(SmbMeltEnum); _assert_(sumM_input);
546+ Input2 *sumR_input = this->GetInput2(SmbRunoffEnum); _assert_(sumR_input);
547+ Input2 *sumP_input = this->GetInput2(SmbPrecipitationEnum); _assert_(sumP_input);
548+ Input2 *ULW_input = this->GetInput2(SmbMeanULWEnum); _assert_(ULW_input);
549+ Input2 *LW_input = this->GetInput2(SmbNetLWEnum); _assert_(LW_input);
550+ Input2 *SW_input = this->GetInput2(SmbNetSWEnum); _assert_(SW_input);
551+ Input2 *LHF_input = this->GetInput2(SmbMeanLHFEnum); _assert_(LHF_input);
552+ Input2 *SHF_input = this->GetInput2(SmbMeanSHFEnum); _assert_(SHF_input);
553+ Input2 *DzAdd_input = this->GetInput2(SmbDzAddEnum); _assert_(DzAdd_input);
554+ Input2 *MassAdd_input = this->GetInput2(SmbMAddEnum); _assert_(MassAdd_input);
555+ Input2 *InitMass_input = this->GetInput2(SmbMInitnum); _assert_(InitMass_input);
556
557- if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << t/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n");
558+ ULW_input->GetInputAverage(&meanULW);
559+ LW_input->GetInputAverage(&netLW);
560+ SW_input->GetInputAverage(&netSW);
561+ LHF_input->GetInputAverage(&meanLHF);
562+ SHF_input->GetInputAverage(&meanSHF);
563+ DzAdd_input->GetInputAverage(&dz_add);
564+ MassAdd_input->GetInputAverage(&sumMassAdd);
565+ sumMassAdd=sumMassAdd*dt;
566+ InitMass_input->GetInputAverage(&initMass);
567+ sumEC_input->GetInputAverage(&sumEC);
568+ sumEC=sumEC*dt*rho_ice;
569+ sumM_input->GetInputAverage(&sumM);
570+ sumM=sumM*dt*rho_ice;
571+ sumR_input->GetInputAverage(&sumR);
572+ sumR=sumR*dt*rho_ice;
573+ sumP_input->GetInputAverage(&sumP);
574+ sumP=sumP*dt*rho_ice;
575+ }
576+ /*}}}*/
577
578- Input2 *Ta_input = this->GetInput2(SmbTaEnum,t-time+timeclim); _assert_(Ta_input);
579- Input2 *V_input = this->GetInput2(SmbVEnum,t-time+timeclim); _assert_(V_input);
580- Input2 *Dlwr_input= this->GetInput2(SmbDlwrfEnum,t-time+timeclim); _assert_(Dlwr_input);
581- Input2 *Dswr_input= this->GetInput2(SmbDswrfEnum,t-time+timeclim); _assert_(Dswr_input);
582- Input2 *P_input = this->GetInput2(SmbPEnum,t-time+timeclim); _assert_(P_input);
583- Input2 *eAir_input= this->GetInput2(SmbEAirEnum,t-time+timeclim); _assert_(eAir_input);
584- Input2 *pAir_input= this->GetInput2(SmbPAirEnum,t-time+timeclim); _assert_(pAir_input);
585+ // Get time forcing inputs
586+ Input2 *Ta_input = this->GetInput2(SmbTaEnum,timeinputs); _assert_(Ta_input);
587+ Input2 *V_input = this->GetInput2(SmbVEnum,timeinputs); _assert_(V_input);
588+ Input2 *Dlwr_input= this->GetInput2(SmbDlwrfEnum,timeinputs); _assert_(Dlwr_input);
589+ Input2 *Dswr_input= this->GetInput2(SmbDswrfEnum,timeinputs); _assert_(Dswr_input);
590+ Input2 *P_input = this->GetInput2(SmbPEnum,timeinputs); _assert_(P_input);
591+ Input2 *eAir_input= this->GetInput2(SmbEAirEnum,timeinputs); _assert_(eAir_input);
592+ Input2 *pAir_input= this->GetInput2(SmbPAirEnum,timeinputs); _assert_(pAir_input);
593+ Input2 *teValue_input= this->GetInput2(SmbTeValueEnum,timeinputs); _assert_(teValue_input);
594+ Input2 *aValue_input= this->GetInput2(SmbAValueEnum,timeinputs); _assert_(aValue_input);
595
596- /*extract daily data:{{{*/
597- Ta_input->GetInputValue(&Ta,gauss);//screen level air temperature [K]
598- V_input->GetInputValue(&V,gauss); //wind speed [m s-1]
599- Dlwr_input->GetInputValue(&dlw,gauss); //downward longwave radiation flux [W m-2]
600- Dswr_input->GetInputValue(&dsw,gauss); //downward shortwave radiation flux [W m-2]
601- P_input->GetInputValue(&P,gauss); //precipitation [kg m-2]
602- eAir_input->GetInputValue(&eAir,gauss); //screen level vapor pressure [Pa]
603- pAir_input->GetInputValue(&pAir,gauss); // screen level air pressure [Pa]
604- teValue_input->GetInputValue(&teValue,gauss); // Emissivity [0-1]
605- aValue_input->GetInputValue(&aValue,gauss); // Albedo [0 1]
606- //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
607- /*}}}*/
608+ /*extract daily data:{{{*/
609+ Ta_input->GetInputValue(&Ta,gauss);//screen level air temperature [K]
610+ V_input->GetInputValue(&V,gauss); //wind speed [m s-1]
611+ Dlwr_input->GetInputValue(&dlw,gauss); //downward longwave radiation flux [W m-2]
612+ Dswr_input->GetInputValue(&dsw,gauss); //downward shortwave radiation flux [W m-2]
613+ P_input->GetInputValue(&P,gauss); //precipitation [kg m-2]
614+ eAir_input->GetInputValue(&eAir,gauss); //screen level vapor pressure [Pa]
615+ pAir_input->GetInputValue(&pAir,gauss); // screen level air pressure [Pa]
616+ teValue_input->GetInputValue(&teValue,gauss); // Emissivity [0-1]
617+ aValue_input->GetInputValue(&aValue,gauss); // Albedo [0 1]
618+ //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
619+ /*}}}*/
620
621- /*Snow grain metamorphism:*/
622- if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid());
623+ /*Snow grain metamorphism:*/
624+ if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid());
625
626- /*Snow, firn and ice albedo:*/
627- if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
628+ /*Snow, firn and ice albedo:*/
629+ if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
630
631- /*Distribution of absorbed short wave radation with depth:*/
632- if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
633+ /*Distribution of absorbed short wave radation with depth:*/
634+ if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
635
636- /*Calculate net shortwave [W m-2]*/
637- netSW = netSW + cellsum(swf,m)*smb_dt/dt;
638+ /*Calculate net shortwave [W m-2]*/
639+ netSW = netSW + cellsum(swf,m)*smb_dt/dt;
640
641- /*Thermal profile computation:*/
642- if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid());
643+ /*Thermal profile computation:*/
644+ if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid());
645
646- /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell.
647- * need to fix this in case all or more of cell evaporates */
648- dz[0] = dz[0] + EC / d[0];
649+ /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell.
650+ * need to fix this in case all or more of cell evaporates */
651+ dz[0] = dz[0] + EC / d[0];
652
653- /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
654- if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, dsnowIdx, Tmean, Ta, P, dzMin, aSnow, C, V, Vmean, rho_ice,this->Sid());
655+ /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
656+ if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, dsnowIdx, Tmean, Ta, P, dzMin, aSnow, C, V, Vmean, rho_ice,this->Sid());
657
658- /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
659- * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
660- if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid());
661+ /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
662+ * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
663+ if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid());
664
665- /*Allow non-melt densification and determine compaction [m]*/
666- if(isdensification)densification(&d,&dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid());
667+ /*Allow non-melt densification and determine compaction [m]*/
668+ if(isdensification)densification(&d,&dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid());
669
670- /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
671- * sub-time step in thermo equations*/
672- //ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here
673+ /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
674+ * sub-time step in thermo equations*/
675+ //ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here
676
677- /*Calculate net longwave [W m-2]*/
678- meanULW = meanULW + ulw*smb_dt/dt;
679- netLW = netLW + (dlw - ulw)*smb_dt/dt;
680+ /*Calculate net longwave [W m-2]*/
681+ meanULW = meanULW + ulw*smb_dt/dt;
682+ netLW = netLW + (dlw - ulw)*smb_dt/dt;
683
684- /*Calculate turbulent heat fluxes [W m-2]*/
685- if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());
686+ /*Calculate turbulent heat fluxes [W m-2]*/
687+ if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());
688
689- /*Verbose some results in debug mode: {{{*/
690- if(VerboseSmb() && 0){
691- _printf_("smb log: count[" << count << "] m[" << m << "] "
692- << setprecision(16) << "T[" << cellsum(T,m) << "] "
693- << "d[" << cellsum(d,m) << "] "
694- << "dz[" << cellsum(dz,m) << "] "
695- << "a[" << cellsum(a,m) << "] "
696- << "W[" << cellsum(W,m) << "] "
697- << "re[" << cellsum(re,m) << "] "
698- << "gdn[" << cellsum(gdn,m) << "] "
699- << "gsp[" << cellsum(gsp,m) << "] "
700- << "swf[" << netSW << "] "
701- << "lwf[" << netLW << "] "
702- << "a[" << a << "] "
703- << "te[" << teValue << "] "
704- << "\n");
705- } /*}}}*/
706+ /*Verbose some results in debug mode: {{{*/
707+ if(VerboseSmb() && 0){
708+ _printf_("smb log: count[" << count << "] m[" << m << "] "
709+ << setprecision(16) << "T[" << cellsum(T,m) << "] "
710+ << "d[" << cellsum(d,m) << "] "
711+ << "dz[" << cellsum(dz,m) << "] "
712+ << "a[" << cellsum(a,m) << "] "
713+ << "W[" << cellsum(W,m) << "] "
714+ << "re[" << cellsum(re,m) << "] "
715+ << "gdn[" << cellsum(gdn,m) << "] "
716+ << "gsp[" << cellsum(gsp,m) << "] "
717+ << "swf[" << netSW << "] "
718+ << "lwf[" << netLW << "] "
719+ << "a[" << a << "] "
720+ << "te[" << teValue << "] "
721+ << "\n");
722+ } /*}}}*/
723
724- meanLHF = meanLHF + lhf*smb_dt/dt;
725- meanSHF = meanSHF + shf*smb_dt/dt;
726+ meanLHF = meanLHF + lhf*smb_dt/dt;
727+ meanSHF = meanSHF + shf*smb_dt/dt;
728
729- /*Sum component mass changes [kg m-2]*/
730- sumMassAdd = mAdd + sumMassAdd;
731- sumM = M + sumM;
732- sumR = R + sumR;
733- sumW = cellsum(W,m);
734- sumP = P + sumP;
735- sumEC = sumEC + EC; // evap (-)/cond(+)
736+ /*Sum component mass changes [kg m-2]*/
737+ sumMassAdd = mAdd + sumMassAdd;
738+ sumM = M + sumM;
739+ sumR = R + sumR;
740+ sumW = cellsum(W,m);
741+ sumP = P + sumP;
742+ sumEC = sumEC + EC; // evap (-)/cond(+)
743
744- /*Calculate total system mass:*/
745- sumMass=0;
746- fac=0;
747- for(int i=0;i<m;i++){
748- sumMass += dz[i]*d[i];
749- fac += dz[i]*(rho_ice - fmin(d[i],rho_ice));
750- }
751+ /*Calculate total system mass:*/
752+ sumMass=0;
753+ fac=0;
754+ for(int i=0;i<m;i++){
755+ sumMass += dz[i]*d[i];
756+ fac += dz[i]*(rho_ice - fmin(d[i],rho_ice));
757+ }
758
759- #if defined(_HAVE_AD_)
760- /*we want to avoid the round operation at all cost. Not differentiable.*/
761- _error_("not implemented yet");
762- #else
763- dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
764- dMass = round(dMass * 100.0)/100.0;
765+ #if defined(_HAVE_AD_)
766+ /*we want to avoid the round operation at all cost. Not differentiable.*/
767+ _error_("not implemented yet");
768+ #else
769+ dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
770+ dMass = round(dMass * 100.0)/100.0;
771
772- /*Check mass conservation:*/
773- if (dMass != 0.0) _printf_("total system mass not conserved in MB function \n");
774- #endif
775+ /*Check mass conservation:*/
776+ if (dMass != 0.0){
777+ _printf_("total system mass not conserved in MB function \n");
778+ }
779+ #endif
780
781- /*Check bottom grid cell T is unchanged:*/
782- if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0){
783- if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
784- }
785+ /*Check bottom grid cell T is unchanged:*/
786+ if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0){
787+ if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
788+ }
789
790- /*Free ressources: */
791- xDelete<IssmDouble>(swf);
792-
793- /*increase counter:*/
794- count++;
795- } //for (t=time;t<time+dt;t=t+smb_dt)
796-
797 /*Save generated inputs: */
798 this->inputs2->SetArrayInput(SmbDzEnum,this->lid,dz,m);
799 this->inputs2->SetArrayInput(SmbDEnum,this->lid,d,m);
800@@ -3994,9 +3995,12 @@
801 this->SetElementInput(SmbNetSWEnum,netSW);
802 this->SetElementInput(SmbMeanLHFEnum,meanLHF);
803 this->SetElementInput(SmbMeanSHFEnum,meanSHF);
804- this->SetElementInput(SmbDzAddEnum,sumdz_add);
805+ this->SetElementInput(SmbDzAddEnum,dz_add);
806+ this->SetElementInput(SmbMInitnum,initMass);
807 this->SetElementInput(SmbMAddEnum,sumMassAdd/dt);
808+ this->SetElementInput(SmbWAddEnum,sumW/dt);
809 this->SetElementInput(SmbFACEnum,fac/1000.); // output in meters
810+ this->SetElementInput(SmbECDtEnum,EC);
811
812 /*Free allocations:{{{*/
813 if(dz) xDelete<IssmDouble>(dz);
814@@ -4015,6 +4019,7 @@
815 if(Wini) xDelete<IssmDouble>(Wini);
816 if(aini) xDelete<IssmDouble>(aini);
817 if(Tini) xDelete<IssmDouble>(Tini);
818+ if(swf) xDelete<IssmDouble>(swf);
819
820 delete gauss;
821 /*}}}*/
Note: See TracBrowser for help on using the repository browser.