source: issm/oecreview/Archive/25834-26739/ISSM-26614-26615.diff@ 26740

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

CHG: added 25834-26739
File size: 65.3 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26614)
4+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26615)
5@@ -82,12 +82,18 @@
6 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
7
8 /*Get moving front parameters*/
9- int calvinglaw;
10- iomodel->FindConstant(&calvinglaw,"md.calving.law");
11- switch(calvinglaw){
12- case DefaultCalvingEnum:
13- iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",CalvingCalvingrateEnum);
14- break;
15+ bool isstochastic;
16+ int calvinglaw;
17+ iomodel->FindConstant(&calvinglaw,"md.calving.law");
18+ iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing");
19+ switch(calvinglaw){
20+ case DefaultCalvingEnum:
21+ iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",CalvingCalvingrateEnum);
22+ if(isstochastic){
23+ iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
24+ iomodel->FetchDataToInput(inputs,elements,"md.calving.calvingrate",BaselineCalvingCalvingrateEnum);
25+ }
26+ break;
27 case CalvingLevermannEnum:
28 iomodel->FetchDataToInput(inputs,elements,"md.calving.coeff",CalvinglevermannCoeffEnum);
29 break;
30Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
31===================================================================
32--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26614)
33+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26615)
34@@ -98,26 +98,49 @@
35 xDelete<IssmDouble>(varspin);
36 xDelete<IssmDouble>(phi_basin);
37 }/*}}}*/
38-void Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type){/*{{{*/
39-
40- const int numvertices = this->GetNumberOfVertices();
41- int basinid,M,N;
42- IssmDouble beta0_basin,beta1_basin,noise_basin;
43+void Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type){/*{{{*/
44+
45+ const int numvertices = this->GetNumberOfVertices();
46+ int basinid,M,N,arenum_type,basinenum_type,noiseenum_type,outenum_type;
47+ IssmDouble beta0_basin,beta1_basin,noiseterm;
48 IssmDouble* phi_basin = xNew<IssmDouble>(arorder);
49 IssmDouble* varlist = xNew<IssmDouble>(numvertices);
50 IssmDouble* valuesautoregression = NULL;
51+ Input* noiseterm_input = NULL;
52
53- /*Get Basin ID and Basin coefficients*/
54- if(enum_type==SMBautoregressionEnum) this->GetInputValue(&basinid,SmbBasinsIdEnum);
55- if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum);
56+ /*Get field-specific enums*/
57+ switch(enum_type){
58+ case(SMBautoregressionEnum):
59+ arenum_type = SmbValuesAutoregressionEnum;
60+ basinenum_type = SmbBasinsIdEnum;
61+ noiseenum_type = SmbAutoregressionNoiseEnum;
62+ outenum_type = SmbMassBalanceEnum;
63+ break;
64+ case(FrontalForcingsRignotAutoregressionEnum):
65+ arenum_type = ThermalforcingValuesAutoregressionEnum;
66+ basinenum_type = FrontalForcingsBasinIdEnum;
67+ noiseenum_type = ThermalforcingAutoregressionNoiseEnum;
68+ outenum_type = FrontalForcingsThermalForcingEnum;
69+ break;
70+ }
71+
72+ /*Get noise and autoregressive terms*/
73+ this->GetInputValue(&basinid,basinenum_type);
74+ if(isfieldstochastic){
75+ noiseterm_input = this->GetInput(noiseenum_type);
76+ Gauss* gauss = this->NewGauss();
77+ noiseterm_input->GetInputValue(&noiseterm,gauss);
78+ delete gauss;
79+ }
80+ else noiseterm = 0.0;
81+ this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M);
82+
83+ /*Get basin coefficients*/
84 for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii];
85 beta0_basin = beta0[basinid];
86 beta1_basin = beta1[basinid];
87- noise_basin = noiseterms[basinid];
88- if(enum_type==SMBautoregressionEnum) this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
89- if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->GetArray(ThermalforcingValuesAutoregressionEnum,this->lid,&valuesautoregression,&M);
90
91- /*If not AR model timestep: take the old values of variable*/
92+ /*If not AR model timestep: take the old values of variable*/
93 if(isstepforar==false){
94 for(int i=0;i<numvertices;i++) varlist[i]=valuesautoregression[i];
95 }
96@@ -130,21 +153,19 @@
97 for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*valuesautoregression[v+ii*numvertices];
98
99 /*Stochastic variable value*/
100- varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin;
101+ varlist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noiseterm;
102 }
103- /*Update autoregression TF values*/
104+ /*Update autoregression values*/
105 IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
106 /*Assign newest values and shift older values*/
107 for(int i=0;i<numvertices;i++) temparray[i] = varlist[i];
108 for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = valuesautoregression[i];
109- if(enum_type==SMBautoregressionEnum) this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder);
110- if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder);
111+ this->inputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder);
112 xDelete<IssmDouble>(temparray);
113 }
114
115 /*Add input to element*/
116- if(enum_type==SMBautoregressionEnum) this->AddInput(SmbMassBalanceEnum,varlist,P1Enum);
117- if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->AddInput(FrontalForcingsThermalForcingEnum,varlist,P1Enum);
118+ this->AddInput(outenum_type,varlist,P1Enum);
119
120 /*Cleanup*/
121 xDelete<IssmDouble>(phi_basin);
122Index: ../trunk-jpl/src/c/classes/Elements/Element.h
123===================================================================
124--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 26614)
125+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 26615)
126@@ -68,7 +68,7 @@
127 /*bool AnyActive(void);*/
128 bool AnyFSet(void);
129 void AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type);
130- void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms,int enum_type);
131+ void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type);
132 void ComputeLambdaS(void);
133 void ComputeNewDamage();
134 void ComputeStrainRate();
135Index: ../trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp
136===================================================================
137--- ../trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26614)
138+++ ../trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 26615)
139@@ -79,7 +79,8 @@
140
141 /*Load parameters*/
142 bool isstochastic;
143- int M,N,Nphi,arorder,numbasins,my_rank;
144+ bool istfstochastic = false;
145+ int M,N,Nphi,arorder,numbasins,my_rank;
146 femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
147 femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
148 IssmDouble tinit_ar;
149@@ -86,7 +87,6 @@
150 IssmDouble* beta0 = NULL;
151 IssmDouble* beta1 = NULL;
152 IssmDouble* phi = NULL;
153- IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
154
155 femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
156 femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum); _assert_(M==numbasins);
157@@ -93,7 +93,6 @@
158 femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum); _assert_(M==numbasins);
159 femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder);
160
161- /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/
162 femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
163 if(isstochastic){
164 int numstochasticfields;
165@@ -101,10 +100,9 @@
166 femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
167 femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
168 for(int i=0;i<numstochasticfields;i++){
169- if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum){
170- femmodel->parameters->FindParam(&noiseterms,&M,ThermalforcingAutoregressionNoiseEnum); _assert_(M==numbasins);
171- }
172+ if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum) istfstochastic = true;
173 }
174+ xDelete<int>(stochasticfields);
175 }
176 /*Time elapsed with respect to AR model initial time*/
177 IssmDouble telapsed_ar = time-tinit_ar;
178@@ -112,7 +110,7 @@
179 /*Loop over each element to compute Thermal Forcing at vertices*/
180 for(Object* &object:femmodel->elements->objects){
181 Element* element = xDynamicCast<Element*>(object);
182- element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,FrontalForcingsRignotAutoregressionEnum);
183+ element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum);
184 }
185
186 /*Cleanup*/
187@@ -119,5 +117,4 @@
188 xDelete<IssmDouble>(beta0);
189 xDelete<IssmDouble>(beta1);
190 xDelete<IssmDouble>(phi);
191- xDelete<IssmDouble>(noiseterms);
192 }/*}}}*/
193Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
194===================================================================
195--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26614)
196+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26615)
197@@ -499,30 +499,30 @@
198 }
199
200 bool isstochasticforcing;
201- parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
202- if(isstochasticforcing){
203- int num_fields;
204- char** fields;
205- parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum));
206- iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields");
207- if(num_fields<1) _error_("no stochasticforcing fields found");
208- int* stochasticforcing_enums = xNew<int>(num_fields);
209- for(int i=0;i<num_fields;i++){
210- stochasticforcing_enums[i] = StringToEnumx(fields[i]);
211- xDelete<char>(fields[i]);
212- }
213- xDelete<char*>(fields);
214- parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields));
215- xDelete<int>(stochasticforcing_enums);
216-
217+ parameters->FindParam(&isstochasticforcing,StochasticForcingIsStochasticForcingEnum);
218+ if(isstochasticforcing){
219+ int num_fields,stochastic_dim;
220+ char** fields;
221+ parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.num_fields",StochasticForcingNumFieldsEnum));
222+ parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.defaultdimension",StochasticForcingDefaultDimensionEnum));
223+ iomodel->FindConstant(&fields,&num_fields,"md.stochasticforcing.fields");
224+ if(num_fields<1) _error_("no stochasticforcing fields found");
225+ int* stochasticforcing_enums = xNew<int>(num_fields);
226+ for(int i=0;i<num_fields;i++){
227+ stochasticforcing_enums[i] = StringToEnumx(fields[i]);
228+ xDelete<char>(fields[i]);
229+ }
230+ xDelete<char*>(fields);
231+ parameters->AddObject(new IntVecParam(StochasticForcingFieldsEnum,stochasticforcing_enums,num_fields));
232+ xDelete<int>(stochasticforcing_enums);
233 parameters->AddObject(iomodel->CopyConstantObject("md.stochasticforcing.randomflag",StochasticForcingRandomflagEnum));
234- iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.dimensions");
235- parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N));
236- xDelete<IssmDouble>(transparam);
237- iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance");
238- parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N));
239- xDelete<IssmDouble>(transparam);
240- }
241+ iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.dimensions");
242+ parameters->AddObject(new IntVecParam(StochasticForcingDimensionsEnum,transparam,N));
243+ xDelete<IssmDouble>(transparam);
244+ iomodel->FetchData(&transparam,&M,&N,"md.stochasticforcing.covariance");
245+ parameters->AddObject(new DoubleMatParam(StochasticForcingCovarianceEnum,transparam,M,N));
246+ xDelete<IssmDouble>(transparam);
247+ }
248
249 /*Deal with mass flux segments: {{{*/
250 iomodel->FetchData(&qmu_mass_flux_present,"md.qmu.mass_flux_segments_present");
251Index: ../trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp
252===================================================================
253--- ../trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 26614)
254+++ ../trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp (revision 26615)
255@@ -35,7 +35,7 @@
256 /*Determine whether random seed is fixed to time step (randomflag==false) or random seed truly random (randomflag==true)*/
257 if(randomflag) fixedseed=-1;
258 else fixedseed = reCast<int,IssmDouble>((time-starttime)/dt);
259- /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
260+ /*multivariateNormal needs to be passed a NULL pointer to avoid memory leak issues*/
261 IssmDouble* temparray = NULL;
262 multivariateNormal(&temparray,dimtot,0.0,covariance,fixedseed);
263 for(int i=0;i<dimtot;i++) noiseterms[i]=temparray[i];
264@@ -45,17 +45,68 @@
265
266 int i=0;
267 for(int j=0;j<numfields;j++){
268- int enum_type;
269+ int dimenum_type,noiseenum_type;
270 IssmDouble* noisefield = xNew<IssmDouble>(dimensions[j]);
271 for(int k=0;k<dimensions[j];k++){
272 noisefield[k]=noiseterms[i+k];
273 }
274-
275- if(fields[j]==SMBautoregressionEnum) enum_type = SmbAutoregressionNoiseEnum;
276- else if(fields[j]==FrontalForcingsRignotAutoregressionEnum) enum_type = ThermalforcingAutoregressionNoiseEnum;
277- else _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.\n");
278- femmodel->parameters->SetParam(noisefield,dimensions[j],enum_type);
279- i=i+dimensions[j];
280+
281+ int dimensionid;
282+
283+ /*Deal with the autoregressive models*/
284+ if(fields[j]==SMBautoregressionEnum || fields[j]==FrontalForcingsRignotAutoregressionEnum){
285+ switch(fields[j]){
286+ case SMBautoregressionEnum:
287+ dimenum_type = SmbBasinsIdEnum;
288+ noiseenum_type = SmbAutoregressionNoiseEnum;
289+ break;
290+ case FrontalForcingsRignotAutoregressionEnum:
291+ dimenum_type = FrontalForcingsBasinIdEnum;
292+ noiseenum_type = ThermalforcingAutoregressionNoiseEnum;
293+ break;
294+ }
295+ for(Object* &object:femmodel->elements->objects){
296+ Element* element = xDynamicCast<Element*>(object);
297+ int numvertices = element->GetNumberOfVertices();
298+ IssmDouble* noise_element = xNew<IssmDouble>(numvertices);
299+ element->GetInputValue(&dimensionid,dimenum_type);
300+ for(int i=0;i<numvertices;i++) noise_element[i] = noisefield[dimensionid];
301+ element->AddInput(noiseenum_type,noise_element,P0Enum);
302+ xDelete<IssmDouble>(noise_element);
303+ }
304+ }
305+ else{
306+ switch(fields[j]){
307+ case SMBautoregressionEnum:
308+ case FrontalForcingsRignotAutoregressionEnum:
309+ /*Already done above*/
310+ break;
311+ case DefaultCalvingEnum:
312+ /*Delete CalvingCalvingrateEnum at previous time step (required if it is transient)*/
313+ femmodel->inputs->DeleteInput(CalvingCalvingrateEnum);
314+ for(Object* &object:femmodel->elements->objects){
315+ Element* element = xDynamicCast<Element*>(object);
316+ int numvertices = element->GetNumberOfVertices();
317+ IssmDouble baselinecalvingrate;
318+ IssmDouble calvingrate_tot[numvertices];
319+ Input* baselinecalvingrate_input = NULL;
320+ baselinecalvingrate_input = element->GetInput(BaselineCalvingCalvingrateEnum); _assert_(baselinecalvingrate_input);
321+ element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
322+ Gauss* gauss = element->NewGauss();
323+ for(int i=0;i<numvertices;i++){
324+ gauss->GaussVertex(i);
325+ baselinecalvingrate_input->GetInputValue(&baselinecalvingrate,gauss);
326+ calvingrate_tot[i] = max(0.0,baselinecalvingrate+noisefield[dimensionid]);
327+ }
328+ element->AddInput(CalvingCalvingrateEnum,&calvingrate_tot[0],P1DGEnum);
329+ delete gauss;
330+ }
331+ break;
332+ default:
333+ _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");
334+ }
335+ }
336+ i=i+dimensions[j];
337 xDelete<IssmDouble>(noisefield);
338 }
339
340Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
341===================================================================
342--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26614)
343+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26615)
344@@ -176,68 +176,64 @@
345 }/*}}}*/
346 void Smbautoregressionx(FemModel* femmodel){/*{{{*/
347
348- /*Get time parameters*/
349- IssmDouble time,dt,starttime,tstep_ar;
350- femmodel->parameters->FindParam(&time,TimeEnum);
351- femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
352- femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
353+ /*Get time parameters*/
354+ IssmDouble time,dt,starttime,tstep_ar;
355+ femmodel->parameters->FindParam(&time,TimeEnum);
356+ femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
357+ femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
358 femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
359
360- /*Initialize module at first time step*/
361- if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
362- /*Determine if this is a time step for the AR model*/
363- bool isstepforar = false;
364+ /*Initialize module at first time step*/
365+ if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
366+ /*Determine if this is a time step for the AR model*/
367+ bool isstepforar = false;
368
369- #ifndef _HAVE_AD_
370+ #ifndef _HAVE_AD_
371 if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt) isstepforar = true;
372- #else
373- _error_("not implemented yet");
374- #endif
375+ #else
376+ _error_("not implemented yet");
377+ #endif
378
379- /*Load parameters*/
380- bool isstochastic;
381- int M,N,Nphi,arorder,numbasins,my_rank;
382- femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
383- femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
384- IssmDouble tinit_ar;
385- IssmDouble* beta0 = NULL;
386- IssmDouble* beta1 = NULL;
387- IssmDouble* phi = NULL;
388+ /*Load parameters*/
389+ bool isstochastic;
390+ bool issmbstochastic = false;
391+ int M,N,Nphi,arorder,numbasins,my_rank;
392+ femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
393+ femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
394+ IssmDouble tinit_ar;
395+ IssmDouble* beta0 = NULL;
396+ IssmDouble* beta1 = NULL;
397+ IssmDouble* phi = NULL;
398
399- femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
400- femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins);
401- femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins);
402- femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder);
403+ femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
404+ femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins);
405+ femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins);
406+ femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder);
407
408- /*Retrieve noise terms if stochasticity, otherwise leave noiseterms as 0*/
409- IssmDouble* noiseterms = xNewZeroInit<IssmDouble>(numbasins);
410- femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
411+ femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
412 if(isstochastic){
413- int numstochasticfields;
414+ int numstochasticfields;
415 int* stochasticfields;
416 femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
417 femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
418 for(int i=0;i<numstochasticfields;i++){
419- if(stochasticfields[i]==SMBautoregressionEnum){
420- femmodel->parameters->FindParam(&noiseterms,&M,SmbAutoregressionNoiseEnum); _assert_(M==numbasins);
421- }
422- }
423- xDelete<int>(stochasticfields);
424- }
425- /*Time elapsed with respect to AR model initial time*/
426- IssmDouble telapsed_ar = time-tinit_ar;
427+ if(stochasticfields[i]==SMBautoregressionEnum) issmbstochastic = true;
428+ }
429+ xDelete<int>(stochasticfields);
430+ }
431+ /*Time elapsed with respect to AR model initial time*/
432+ IssmDouble telapsed_ar = time-tinit_ar;
433
434- /*Loop over each element to compute SMB at vertices*/
435- for(Object* &object:femmodel->elements->objects){
436- Element* element = xDynamicCast<Element*>(object);
437- element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,noiseterms,SMBautoregressionEnum);
438- }
439+ /*Loop over each element to compute SMB at vertices*/
440+ for(Object* &object:femmodel->elements->objects){
441+ Element* element = xDynamicCast<Element*>(object);
442+ element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum);
443+ }
444
445- /*Cleanup*/
446- xDelete<IssmDouble>(beta0);
447- xDelete<IssmDouble>(beta1);
448- xDelete<IssmDouble>(phi);
449- xDelete<IssmDouble>(noiseterms);
450+ /*Cleanup*/
451+ xDelete<IssmDouble>(beta0);
452+ xDelete<IssmDouble>(beta1);
453+ xDelete<IssmDouble>(phi);
454 }/*}}}*/
455 void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
456
457Index: ../trunk-jpl/src/c/shared/Enum/Enum.vim
458===================================================================
459--- ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26614)
460+++ ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26615)
461@@ -397,6 +397,7 @@
462 syn keyword cConstant SolidearthSettingsGrdOceanEnum
463 syn keyword cConstant SolidearthSettingsOceanAreaScalingEnum
464 syn keyword cConstant StochasticForcingCovarianceEnum
465+syn keyword cConstant StochasticForcingDefaultDimensionEnum
466 syn keyword cConstant StochasticForcingDimensionsEnum
467 syn keyword cConstant StochasticForcingFieldsEnum
468 syn keyword cConstant StochasticForcingIsStochasticForcingEnum
469@@ -427,7 +428,6 @@
470 syn keyword cConstant SmbAccurefEnum
471 syn keyword cConstant SmbAdThreshEnum
472 syn keyword cConstant SmbAutoregressionInitialTimeEnum
473-syn keyword cConstant SmbAutoregressionNoiseEnum
474 syn keyword cConstant SmbAutoregressionTimestepEnum
475 syn keyword cConstant SmbAutoregressiveOrderEnum
476 syn keyword cConstant SmbAveragingEnum
477@@ -503,7 +503,6 @@
478 syn keyword cConstant StressbalanceRestolEnum
479 syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum
480 syn keyword cConstant StressbalanceShelfDampeningEnum
481-syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
482 syn keyword cConstant ThermalIsdrainicecolumnEnum
483 syn keyword cConstant ThermalIsdynamicbasalspcEnum
484 syn keyword cConstant ThermalIsenthalpyEnum
485@@ -601,6 +600,7 @@
486 syn keyword cConstant BaseOldEnum
487 syn keyword cConstant BaseSlopeXEnum
488 syn keyword cConstant BaseSlopeYEnum
489+syn keyword cConstant BaselineCalvingCalvingrateEnum
490 syn keyword cConstant BedEnum
491 syn keyword cConstant BedGRDEnum
492 syn keyword cConstant BedEastEnum
493@@ -894,6 +894,7 @@
494 syn keyword cConstant SmbAccumulationEnum
495 syn keyword cConstant SmbAdiffiniEnum
496 syn keyword cConstant SmbAiniEnum
497+syn keyword cConstant SmbAutoregressionNoiseEnum
498 syn keyword cConstant SmbBasinsIdEnum
499 syn keyword cConstant SmbBMaxEnum
500 syn keyword cConstant SmbBMinEnum
501@@ -996,6 +997,7 @@
502 syn keyword cConstant SolidearthExternalDisplacementNorthRateEnum
503 syn keyword cConstant SolidearthExternalDisplacementUpRateEnum
504 syn keyword cConstant SolidearthExternalGeoidRateEnum
505+syn keyword cConstant StochasticForcingDefaultIdEnum
506 syn keyword cConstant StrainRateeffectiveEnum
507 syn keyword cConstant StrainRateparallelEnum
508 syn keyword cConstant StrainRateperpendicularEnum
509@@ -1031,6 +1033,7 @@
510 syn keyword cConstant TemperaturePDDEnum
511 syn keyword cConstant TemperaturePicardEnum
512 syn keyword cConstant TemperatureSEMICEnum
513+syn keyword cConstant ThermalforcingAutoregressionNoiseEnum
514 syn keyword cConstant ThermalforcingValuesAutoregressionEnum
515 syn keyword cConstant ThermalSpctemperatureEnum
516 syn keyword cConstant ThicknessAbsGradientEnum
517Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
518===================================================================
519--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26614)
520+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26615)
521@@ -391,6 +391,7 @@
522 SolidearthSettingsGrdOceanEnum,
523 SolidearthSettingsOceanAreaScalingEnum,
524 StochasticForcingCovarianceEnum,
525+ StochasticForcingDefaultDimensionEnum,
526 StochasticForcingDimensionsEnum,
527 StochasticForcingFieldsEnum,
528 StochasticForcingIsStochasticForcingEnum,
529@@ -421,7 +422,6 @@
530 SmbAccurefEnum,
531 SmbAdThreshEnum,
532 SmbAutoregressionInitialTimeEnum,
533- SmbAutoregressionNoiseEnum,
534 SmbAutoregressionTimestepEnum,
535 SmbAutoregressiveOrderEnum,
536 SmbAveragingEnum,
537@@ -497,7 +497,6 @@
538 StressbalanceRestolEnum,
539 StressbalanceRiftPenaltyThresholdEnum,
540 StressbalanceShelfDampeningEnum,
541- ThermalforcingAutoregressionNoiseEnum,
542 ThermalIsdrainicecolumnEnum,
543 ThermalIsdynamicbasalspcEnum,
544 ThermalIsenthalpyEnum,
545@@ -597,6 +596,7 @@
546 BaseOldEnum,
547 BaseSlopeXEnum,
548 BaseSlopeYEnum,
549+ BaselineCalvingCalvingrateEnum,
550 BedEnum,
551 BedGRDEnum,
552 BedEastEnum,
553@@ -890,6 +890,7 @@
554 SmbAccumulationEnum,
555 SmbAdiffiniEnum,
556 SmbAiniEnum,
557+ SmbAutoregressionNoiseEnum,
558 SmbBasinsIdEnum,
559 SmbBMaxEnum,
560 SmbBMinEnum,
561@@ -993,6 +994,7 @@
562 SolidearthExternalDisplacementNorthRateEnum,
563 SolidearthExternalDisplacementUpRateEnum,
564 SolidearthExternalGeoidRateEnum,
565+ StochasticForcingDefaultIdEnum,
566 StrainRateeffectiveEnum,
567 StrainRateparallelEnum,
568 StrainRateperpendicularEnum,
569@@ -1028,6 +1030,7 @@
570 TemperaturePDDEnum,
571 TemperaturePicardEnum,
572 TemperatureSEMICEnum,
573+ ThermalforcingAutoregressionNoiseEnum,
574 ThermalforcingValuesAutoregressionEnum,
575 ThermalSpctemperatureEnum,
576 ThicknessAbsGradientEnum,
577Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
578===================================================================
579--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26614)
580+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26615)
581@@ -399,6 +399,7 @@
582 case SolidearthSettingsGrdOceanEnum : return "SolidearthSettingsGrdOcean";
583 case SolidearthSettingsOceanAreaScalingEnum : return "SolidearthSettingsOceanAreaScaling";
584 case StochasticForcingCovarianceEnum : return "StochasticForcingCovariance";
585+ case StochasticForcingDefaultDimensionEnum : return "StochasticForcingDefaultDimension";
586 case StochasticForcingDimensionsEnum : return "StochasticForcingDimensions";
587 case StochasticForcingFieldsEnum : return "StochasticForcingFields";
588 case StochasticForcingIsStochasticForcingEnum : return "StochasticForcingIsStochasticForcing";
589@@ -429,7 +430,6 @@
590 case SmbAccurefEnum : return "SmbAccuref";
591 case SmbAdThreshEnum : return "SmbAdThresh";
592 case SmbAutoregressionInitialTimeEnum : return "SmbAutoregressionInitialTime";
593- case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
594 case SmbAutoregressionTimestepEnum : return "SmbAutoregressionTimestep";
595 case SmbAutoregressiveOrderEnum : return "SmbAutoregressiveOrder";
596 case SmbAveragingEnum : return "SmbAveraging";
597@@ -505,7 +505,6 @@
598 case StressbalanceRestolEnum : return "StressbalanceRestol";
599 case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold";
600 case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening";
601- case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
602 case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn";
603 case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc";
604 case ThermalIsenthalpyEnum : return "ThermalIsenthalpy";
605@@ -603,6 +602,7 @@
606 case BaseOldEnum : return "BaseOld";
607 case BaseSlopeXEnum : return "BaseSlopeX";
608 case BaseSlopeYEnum : return "BaseSlopeY";
609+ case BaselineCalvingCalvingrateEnum : return "BaselineCalvingCalvingrate";
610 case BedEnum : return "Bed";
611 case BedGRDEnum : return "BedGRD";
612 case BedEastEnum : return "BedEast";
613@@ -896,6 +896,7 @@
614 case SmbAccumulationEnum : return "SmbAccumulation";
615 case SmbAdiffiniEnum : return "SmbAdiffini";
616 case SmbAiniEnum : return "SmbAini";
617+ case SmbAutoregressionNoiseEnum : return "SmbAutoregressionNoise";
618 case SmbBasinsIdEnum : return "SmbBasinsId";
619 case SmbBMaxEnum : return "SmbBMax";
620 case SmbBMinEnum : return "SmbBMin";
621@@ -998,6 +999,7 @@
622 case SolidearthExternalDisplacementNorthRateEnum : return "SolidearthExternalDisplacementNorthRate";
623 case SolidearthExternalDisplacementUpRateEnum : return "SolidearthExternalDisplacementUpRate";
624 case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate";
625+ case StochasticForcingDefaultIdEnum : return "StochasticForcingDefaultId";
626 case StrainRateeffectiveEnum : return "StrainRateeffective";
627 case StrainRateparallelEnum : return "StrainRateparallel";
628 case StrainRateperpendicularEnum : return "StrainRateperpendicular";
629@@ -1033,6 +1035,7 @@
630 case TemperaturePDDEnum : return "TemperaturePDD";
631 case TemperaturePicardEnum : return "TemperaturePicard";
632 case TemperatureSEMICEnum : return "TemperatureSEMIC";
633+ case ThermalforcingAutoregressionNoiseEnum : return "ThermalforcingAutoregressionNoise";
634 case ThermalforcingValuesAutoregressionEnum : return "ThermalforcingValuesAutoregression";
635 case ThermalSpctemperatureEnum : return "ThermalSpctemperature";
636 case ThicknessAbsGradientEnum : return "ThicknessAbsGradient";
637Index: ../trunk-jpl/src/c/shared/Enum/Enumjl.vim
638===================================================================
639--- ../trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26614)
640+++ ../trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26615)
641@@ -169,8 +169,14 @@
642 syn keyword juliaConstC FrictionThresholdSpeedEnum
643 syn keyword juliaConstC FrictionVoidRatioEnum
644 syn keyword juliaConstC FrontalForcingsBasinIcefrontAreaEnum
645+syn keyword juliaConstC FrontalForcingsAutoregressionInitialTimeEnum
646+syn keyword juliaConstC FrontalForcingsAutoregressionTimestepEnum
647+syn keyword juliaConstC FrontalForcingsAutoregressiveOrderEnum
648+syn keyword juliaConstC FrontalForcingsBeta0Enum
649+syn keyword juliaConstC FrontalForcingsBeta1Enum
650 syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum
651 syn keyword juliaConstC FrontalForcingsParamEnum
652+syn keyword juliaConstC FrontalForcingsPhiEnum
653 syn keyword juliaConstC GrdModelEnum
654 syn keyword juliaConstC GroundinglineFrictionInterpolationEnum
655 syn keyword juliaConstC GroundinglineMeltInterpolationEnum
656@@ -383,6 +389,13 @@
657 syn keyword juliaConstC SolidearthSettingsMaxiterEnum
658 syn keyword juliaConstC SolidearthSettingsGrdOceanEnum
659 syn keyword juliaConstC SolidearthSettingsOceanAreaScalingEnum
660+syn keyword juliaConstC StochasticForcingCovarianceEnum
661+syn keyword juliaConstC StochasticForcingDefaultDimensionEnum
662+syn keyword juliaConstC StochasticForcingDimensionsEnum
663+syn keyword juliaConstC StochasticForcingFieldsEnum
664+syn keyword juliaConstC StochasticForcingIsStochasticForcingEnum
665+syn keyword juliaConstC StochasticForcingNumFieldsEnum
666+syn keyword juliaConstC StochasticForcingRandomflagEnum
667 syn keyword juliaConstC RotationalPolarMoiEnum
668 syn keyword juliaConstC SolidearthSettingsReltolEnum
669 syn keyword juliaConstC SealevelchangeRequestedOutputsEnum
670@@ -413,7 +426,6 @@
671 syn keyword juliaConstC SmbAveragingEnum
672 syn keyword juliaConstC SmbBeta0Enum
673 syn keyword juliaConstC SmbBeta1Enum
674-syn keyword juliaConstC SmbCovmatEnum
675 syn keyword juliaConstC SmbDesfacEnum
676 syn keyword juliaConstC SmbDpermilEnum
677 syn keyword juliaConstC SmbDsnowIdxEnum
678@@ -423,6 +435,7 @@
679 syn keyword juliaConstC SmbDenIdxEnum
680 syn keyword juliaConstC SmbDtEnum
681 syn keyword juliaConstC SmbEnum
682+syn keyword juliaConstC SmbEIdxEnum
683 syn keyword juliaConstC SmbFEnum
684 syn keyword juliaConstC SmbInitDensityScalingEnum
685 syn keyword juliaConstC SmbIsaccumulationEnum
686@@ -431,6 +444,7 @@
687 syn keyword juliaConstC SmbIsd18opdEnum
688 syn keyword juliaConstC SmbIsdelta18oEnum
689 syn keyword juliaConstC SmbIsdensificationEnum
690+syn keyword juliaConstC SmbIsdeltaLWupEnum
691 syn keyword juliaConstC SmbIsfirnwarmingEnum
692 syn keyword juliaConstC SmbIsgraingrowthEnum
693 syn keyword juliaConstC SmbIsmeltEnum
694@@ -446,7 +460,6 @@
695 syn keyword juliaConstC SmbNumRequestedOutputsEnum
696 syn keyword juliaConstC SmbPfacEnum
697 syn keyword juliaConstC SmbPhiEnum
698-syn keyword juliaConstC SmbRandomflagEnum
699 syn keyword juliaConstC SmbRdlEnum
700 syn keyword juliaConstC SmbRequestedOutputsEnum
701 syn keyword juliaConstC SmbRlapsEnum
702@@ -459,6 +472,7 @@
703 syn keyword juliaConstC SmbSwIdxEnum
704 syn keyword juliaConstC SmbT0dryEnum
705 syn keyword juliaConstC SmbT0wetEnum
706+syn keyword juliaConstC SmbTeThreshEnum
707 syn keyword juliaConstC SmbTdiffEnum
708 syn keyword juliaConstC SmbThermoDeltaTScalingEnum
709 syn keyword juliaConstC SmbTemperaturesReconstructedYearsEnum
710@@ -579,6 +593,7 @@
711 syn keyword juliaConstC BaseOldEnum
712 syn keyword juliaConstC BaseSlopeXEnum
713 syn keyword juliaConstC BaseSlopeYEnum
714+syn keyword juliaConstC BaselineCalvingCalvingrateEnum
715 syn keyword juliaConstC BedEnum
716 syn keyword juliaConstC BedGRDEnum
717 syn keyword juliaConstC BedEastEnum
718@@ -872,6 +887,7 @@
719 syn keyword juliaConstC SmbAccumulationEnum
720 syn keyword juliaConstC SmbAdiffiniEnum
721 syn keyword juliaConstC SmbAiniEnum
722+syn keyword juliaConstC SmbAutoregressionNoiseEnum
723 syn keyword juliaConstC SmbBasinsIdEnum
724 syn keyword juliaConstC SmbBMaxEnum
725 syn keyword juliaConstC SmbBMinEnum
726@@ -893,6 +909,7 @@
727 syn keyword juliaConstC SmbDailywindspeedEnum
728 syn keyword juliaConstC SmbDiniEnum
729 syn keyword juliaConstC SmbDlwrfEnum
730+syn keyword juliaConstC SmbDulwrfValueEnum
731 syn keyword juliaConstC SmbDswrfEnum
732 syn keyword juliaConstC SmbDswdiffrfEnum
733 syn keyword juliaConstC SmbDzAddEnum
734@@ -973,6 +990,7 @@
735 syn keyword juliaConstC SolidearthExternalDisplacementNorthRateEnum
736 syn keyword juliaConstC SolidearthExternalDisplacementUpRateEnum
737 syn keyword juliaConstC SolidearthExternalGeoidRateEnum
738+syn keyword juliaConstC StochasticForcingDefaultIdEnum
739 syn keyword juliaConstC StrainRateeffectiveEnum
740 syn keyword juliaConstC StrainRateparallelEnum
741 syn keyword juliaConstC StrainRateperpendicularEnum
742@@ -1008,6 +1026,8 @@
743 syn keyword juliaConstC TemperaturePDDEnum
744 syn keyword juliaConstC TemperaturePicardEnum
745 syn keyword juliaConstC TemperatureSEMICEnum
746+syn keyword juliaConstC ThermalforcingAutoregressionNoiseEnum
747+syn keyword juliaConstC ThermalforcingValuesAutoregressionEnum
748 syn keyword juliaConstC ThermalSpctemperatureEnum
749 syn keyword juliaConstC ThicknessAbsGradientEnum
750 syn keyword juliaConstC ThicknessAbsMisfitEnum
751@@ -1254,6 +1274,7 @@
752 syn keyword juliaConstC FreeSurfaceTopAnalysisEnum
753 syn keyword juliaConstC FrontalForcingsDefaultEnum
754 syn keyword juliaConstC FrontalForcingsRignotEnum
755+syn keyword juliaConstC FrontalForcingsRignotAutoregressionEnum
756 syn keyword juliaConstC FsetEnum
757 syn keyword juliaConstC FullMeltOnPartiallyFloatingEnum
758 syn keyword juliaConstC GLheightadvectionAnalysisEnum
759Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
760===================================================================
761--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26614)
762+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26615)
763@@ -408,6 +408,7 @@
764 else if (strcmp(name,"SolidearthSettingsGrdOcean")==0) return SolidearthSettingsGrdOceanEnum;
765 else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum;
766 else if (strcmp(name,"StochasticForcingCovariance")==0) return StochasticForcingCovarianceEnum;
767+ else if (strcmp(name,"StochasticForcingDefaultDimension")==0) return StochasticForcingDefaultDimensionEnum;
768 else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
769 else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum;
770 else if (strcmp(name,"StochasticForcingIsStochasticForcing")==0) return StochasticForcingIsStochasticForcingEnum;
771@@ -438,7 +439,6 @@
772 else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum;
773 else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
774 else if (strcmp(name,"SmbAutoregressionInitialTime")==0) return SmbAutoregressionInitialTimeEnum;
775- else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum;
776 else if (strcmp(name,"SmbAutoregressionTimestep")==0) return SmbAutoregressionTimestepEnum;
777 else if (strcmp(name,"SmbAutoregressiveOrder")==0) return SmbAutoregressiveOrderEnum;
778 else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum;
779@@ -517,7 +517,6 @@
780 else if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum;
781 else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
782 else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
783- else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum;
784 else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum;
785 else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
786 else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
787@@ -615,6 +614,7 @@
788 else if (strcmp(name,"BaseOld")==0) return BaseOldEnum;
789 else if (strcmp(name,"BaseSlopeX")==0) return BaseSlopeXEnum;
790 else if (strcmp(name,"BaseSlopeY")==0) return BaseSlopeYEnum;
791+ else if (strcmp(name,"BaselineCalvingCalvingrate")==0) return BaselineCalvingCalvingrateEnum;
792 else if (strcmp(name,"Bed")==0) return BedEnum;
793 else if (strcmp(name,"BedGRD")==0) return BedGRDEnum;
794 else if (strcmp(name,"BedEast")==0) return BedEastEnum;
795@@ -917,6 +917,7 @@
796 else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
797 else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
798 else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
799+ else if (strcmp(name,"SmbAutoregressionNoise")==0) return SmbAutoregressionNoiseEnum;
800 else if (strcmp(name,"SmbBasinsId")==0) return SmbBasinsIdEnum;
801 else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
802 else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
803@@ -996,11 +997,11 @@
804 else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
805 else if (strcmp(name,"SmbT")==0) return SmbTEnum;
806 else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
807- else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
808 else stage=9;
809 }
810 if(stage==9){
811- if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
812+ if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
813+ else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
814 else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
815 else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
816 else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum;
817@@ -1022,6 +1023,7 @@
818 else if (strcmp(name,"SolidearthExternalDisplacementNorthRate")==0) return SolidearthExternalDisplacementNorthRateEnum;
819 else if (strcmp(name,"SolidearthExternalDisplacementUpRate")==0) return SolidearthExternalDisplacementUpRateEnum;
820 else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
821+ else if (strcmp(name,"StochasticForcingDefaultId")==0) return StochasticForcingDefaultIdEnum;
822 else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
823 else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
824 else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
825@@ -1057,6 +1059,7 @@
826 else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
827 else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
828 else if (strcmp(name,"TemperatureSEMIC")==0) return TemperatureSEMICEnum;
829+ else if (strcmp(name,"ThermalforcingAutoregressionNoise")==0) return ThermalforcingAutoregressionNoiseEnum;
830 else if (strcmp(name,"ThermalforcingValuesAutoregression")==0) return ThermalforcingValuesAutoregressionEnum;
831 else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
832 else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
833@@ -1117,13 +1120,13 @@
834 else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
835 else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
836 else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
837- else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
838- else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
839- else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
840 else stage=10;
841 }
842 if(stage==10){
843- if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
844+ if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
845+ else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
846+ else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
847+ else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
848 else if (strcmp(name,"Outputdefinition30")==0) return Outputdefinition30Enum;
849 else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
850 else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
851@@ -1240,13 +1243,13 @@
852 else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
853 else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
854 else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
855- else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
856- else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
857- else if (strcmp(name,"Channel")==0) return ChannelEnum;
858 else stage=11;
859 }
860 if(stage==11){
861- if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
862+ if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
863+ else if (strcmp(name,"Cflevelsetmisfit")==0) return CflevelsetmisfitEnum;
864+ else if (strcmp(name,"Channel")==0) return ChannelEnum;
865+ else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
866 else if (strcmp(name,"ChannelAreaOld")==0) return ChannelAreaOldEnum;
867 else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum;
868 else if (strcmp(name,"Closed")==0) return ClosedEnum;
869@@ -1363,13 +1366,13 @@
870 else if (strcmp(name,"IntParam")==0) return IntParamEnum;
871 else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
872 else if (strcmp(name,"Inputs")==0) return InputsEnum;
873- else if (strcmp(name,"Internal")==0) return InternalEnum;
874- else if (strcmp(name,"Intersect")==0) return IntersectEnum;
875- else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
876 else stage=12;
877 }
878 if(stage==12){
879- if (strcmp(name,"J")==0) return JEnum;
880+ if (strcmp(name,"Internal")==0) return InternalEnum;
881+ else if (strcmp(name,"Intersect")==0) return IntersectEnum;
882+ else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
883+ else if (strcmp(name,"J")==0) return JEnum;
884 else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
885 else if (strcmp(name,"MLHOApproximation")==0) return MLHOApproximationEnum;
886 else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
887@@ -1486,13 +1489,13 @@
888 else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
889 else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
890 else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
891- else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
892- else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
893- else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
894 else stage=13;
895 }
896 if(stage==13){
897- if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
898+ if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
899+ else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
900+ else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
901+ else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
902 else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
903 else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
904 else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum;
905Index: ../trunk-jpl/test/NightlyRun/test257.m
906===================================================================
907--- ../trunk-jpl/test/NightlyRun/test257.m (revision 26614)
908+++ ../trunk-jpl/test/NightlyRun/test257.m (revision 26615)
909@@ -48,7 +48,6 @@
910 %Stochastic forcing
911 md.stochasticforcing.isstochasticforcing = 1;
912 md.stochasticforcing.fields = [{'SMBautoregression'}];
913-md.stochasticforcing.dimensions = [md.smb.num_basins]; %dimension of each field
914 md.stochasticforcing.covariance = [[0.15 0.08 -0.02];[0.08 0.12 -0.05];[-0.02 -0.05 0.1]]; %global covariance among- and between-fields
915 md.stochasticforcing.randomflag = 0; %fixed random seeds
916
917Index: ../trunk-jpl/test/NightlyRun/test257.py
918===================================================================
919--- ../trunk-jpl/test/NightlyRun/test257.py (revision 26614)
920+++ ../trunk-jpl/test/NightlyRun/test257.py (revision 26615)
921@@ -55,7 +55,6 @@
922 # Stochastic forcing
923 md.stochasticforcing.isstochasticforcing = 1
924 md.stochasticforcing.fields = ['SMBautoregression']
925-md.stochasticforcing.dimensions = [md.smb.num_basins] # dimension of each field
926 md.stochasticforcing.covariance = np.array([[0.15, 0.08, -0.02], [0.08, 0.12, -0.05], [-0.02, -0.05, 0.1]]) # global covariance among- and between-fields
927 md.stochasticforcing.randomflag = 0 # fixed random seeds
928
929Index: ../trunk-jpl/test/NightlyRun/test543.m
930===================================================================
931--- ../trunk-jpl/test/NightlyRun/test543.m (revision 26614)
932+++ ../trunk-jpl/test/NightlyRun/test543.m (revision 26615)
933@@ -38,7 +38,6 @@
934 %Stochastic forcing
935 md.stochasticforcing.isstochasticforcing = 1;
936 md.stochasticforcing.fields = [{'FrontalForcingsRignotAutoregression'}];
937-md.stochasticforcing.dimensions = [md.frontalforcings.num_basins]; %dimension of each field
938 md.stochasticforcing.covariance = 1e-4*[[1.5,0.5];[0.5,0.4]]; %global covariance among- and between-fields
939 md.stochasticforcing.randomflag = 0; %determines true/false randomness
940
941Index: ../trunk-jpl/test/NightlyRun/test543.py
942===================================================================
943--- ../trunk-jpl/test/NightlyRun/test543.py (revision 26614)
944+++ ../trunk-jpl/test/NightlyRun/test543.py (revision 26615)
945@@ -48,7 +48,6 @@
946 # Stochastic forcing
947 md.stochasticforcing.isstochasticforcing = 1
948 md.stochasticforcing.fields = ['FrontalForcingsRignotAutoregression']
949-md.stochasticforcing.dimensions = [md.frontalforcings.num_basins] # dimension of each field
950 md.stochasticforcing.covariance = 1e-4 * np.array([[1.5, 0.5], [0.5, 0.4]]) # global covariance among- and between-fields
951 md.stochasticforcing.randomflag = 0 # determines true/false randomness
952
953Index: ../trunk-jpl/src/m/classes/stochasticforcing.m
954===================================================================
955--- ../trunk-jpl/src/m/classes/stochasticforcing.m (revision 26614)
956+++ ../trunk-jpl/src/m/classes/stochasticforcing.m (revision 26615)
957@@ -7,7 +7,8 @@
958 properties (SetAccess=public)
959 isstochasticforcing = 0;
960 fields = NaN;
961- dimensions = NaN;
962+ defaultdimension = 0;
963+ default_id = NaN;
964 covariance = NaN;
965 randomflag = 1;
966 end
967@@ -33,33 +34,68 @@
968 if ~self.isstochasticforcing, return; end
969
970 num_fields = numel(self.fields);
971- size_tot = sum(self.dimensions);
972+
973+ %Check that covariance matrix is positive definite
974+ try
975+ chol(self.covariance);
976+ catch
977+ error('md.stochasticforcing.covariance is not positive definite');
978+ end
979
980+ %Check that all fields agree with the corresponding md class and if any field needs the default params
981+ checkdefaults = false; %need to check defaults only if one of the field does not have its own dimensionality
982+ for field=self.fields
983+ %Checking agreement of classes
984+ if(contains(field,'SMB'))
985+ if~(isequal(class(md.smb),char(field)))
986+ error('md.smb does not agree with stochasticforcing field %s', char(field));
987+ end
988+ end
989+ if(contains(field,'frontalforcings'))
990+ if~(isequal(class(md.frontalforcings),char(field)))
991+ error('md.frontalforcings does not agree with stochasticforcing field %s', char(field));
992+ end
993+ end
994+ %Checking for specific dimensions
995+ if ~(strcmp(field,'SMBautoregression') || strcmp(field,'FrontalForcingsRignotAutoregression'))
996+ checkdefaults = true; %field with non-specific dimensionality
997+ end
998+ end
999+
1000+ %Retrieve sum of all the field dimensionalities
1001+ size_tot = self.defaultdimension*num_fields;
1002+ indSMBar = -1; %about to check for index of SMBautoregression
1003+ indTFar = -1; %about to check for index of FrontalForcingsRignotAutoregression
1004+ if any(contains(self.fields,'SMBautoregression'))
1005+ size_tot = size_tot-self.defaultdimension+md.smb.num_basins;
1006+ indSMBar = find(contains(self.fields,'SMBautoregression')); %index of SMBar, now check for consistency with TFar timestep (08Nov2021)
1007+ end
1008+ if any(contains(self.fields,'FrontalForcingsRignotAutoregression'))
1009+ size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins;
1010+ indTFar = find(contains(self.fields,'FrontalForcingsRignotAutoregression')); %index of TFar, now check for consistency with SMBar timestep (08Nov2021)
1011+ end
1012+
1013+ if(indSMBar~=-1 && indTFar~=-1) %both autoregressive models are used: check autoregressive time step consistency
1014+ if((md.smb.ar_timestep~=md.frontalforcings.ar_timestep) && any(self.covariance(1+sum(self.dimensions(1:indSMBar-1)):sum(self.dimensions(1:indSMBar)),1+sum(self.dimensions(1:indTFar-1)):sum(self.dimensions(1:indTFar))))~=0)
1015+ error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance');
1016+ end
1017+ end
1018+
1019 md = checkfield(md,'fieldname','stochasticforcing.isstochasticforcing','values',[0 1]);
1020- md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings()); %VV check here 'cell' (19Oct2021)
1021- md = checkfield(md,'fieldname','stochasticforcing.dimensions','NaN',1,'Inf',1,'>',0,'size',[1,num_fields]); %specific dimension for each field
1022+ md = checkfield(md,'fieldname','stochasticforcing.fields','numel',num_fields,'cell',1,'values',supportedstochforcings());
1023 md = checkfield(md,'fieldname','stochasticforcing.covariance','NaN',1,'Inf',1,'size',[size_tot,size_tot]); %global covariance matrix
1024 md = checkfield(md,'fieldname','stochasticforcing.randomflag','numel',[1],'values',[0 1]);
1025-
1026- %Check that all fields agree with the corresponding md class
1027- for field=self.fields
1028- if(contains(field,'SMB'))
1029- if~(isequal(class(md.smb),char(field)))
1030- error('md.smb does not agree with stochasticforcing field %s', char(field));
1031- end
1032- end
1033- if(contains(field,'frontalforcings'))
1034- if~(isequal(class(md.frontalforcings),char(field)))
1035- error('md.frontalforcings does not agree with stochasticforcing field %s', char(field));
1036- end
1037- end
1038- end
1039+ if(checkdefaults) %need to check the defaults
1040+ md = checkfield(md,'fieldname','stochasticforcing.defaultdimension','numel',1,'NaN',1,'Inf',1,'>',0);
1041+ md = checkfield(md,'fieldname','stochasticforcing.default_id','Inf',1,'>=',0,'<=',self.defaultdimension,'size',[md.mesh.numberofelements,1]);
1042+ end
1043 end % }}}
1044 function disp(self) % {{{
1045 disp(sprintf(' stochasticforcing parameters:'));
1046 fielddisplay(self,'isstochasticforcing','is stochasticity activated?');
1047 fielddisplay(self,'fields','fields with stochasticity applied, ex: {''SMBautoregression''}, or {''FrontalForcingsRignotAutoregression''}');
1048- fielddisplay(self,'dimensions','dimensionality of each field');
1049+ fielddisplay(self,'defaultdimension','dimensionality of the noise terms (does not apply to fields with their specific dimension)');
1050+ fielddisplay(self,'default_id','id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)');
1051 fielddisplay(self,'covariance','covariance matrix for within- and between-fields covariance (units must be squared field units)');
1052 fielddisplay(self,'randomflag','whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)');
1053 disp('Available fields:');
1054@@ -76,12 +112,27 @@
1055 if ~self.isstochasticforcing
1056 return
1057 else
1058+
1059+ %Retrieve dimensionality of each field
1060+ dimensions = self.defaultdimension*ones(1,num_fields);
1061+ ind = 1;
1062+ for field=self.fields
1063+ %Checking for specific dimensions
1064+ if(strcmp(field,'SMBautoregression'))
1065+ dimensions(ind) = md.smb.num_basins;
1066+ end
1067+ if(strcmp(field,'FrontalForcingsRignotAutoregression'))
1068+ dimensions(ind) = md.frontalforcings.num_basins;
1069+ end
1070+ ind = ind+1;
1071+ end
1072+
1073 %Scaling covariance matrix (scale column-by-column and row-by-row)
1074- scaledfields = {'SMBautoregression'}; %list of fields that need scaling *1/yts
1075+ scaledfields = {'DefaultCalving','SMBautoregression'}; %list of fields that need scaling *1/yts
1076 tempcovariance = self.covariance; %copy of covariance to avoid writing back in member variable
1077 for i=1:num_fields
1078 if any(strcmp(scaledfields,self.fields(i)))
1079- inds = [1+sum(self.dimensions(1:i-1)):1:sum(self.dimensions(1:i))];
1080+ inds = [1+sum(dimensions(1:i-1)):1:sum(dimensions(1:i))];
1081 for row=inds %scale rows corresponding to scaled field
1082 tempcovariance(row,:) = 1./yts*tempcovariance(row,:);
1083 end
1084@@ -92,7 +143,9 @@
1085 end
1086 WriteData(fid,prefix,'data',num_fields,'name','md.stochasticforcing.num_fields','format','Integer');
1087 WriteData(fid,prefix,'object',self,'fieldname','fields','format','StringArray');
1088- WriteData(fid,prefix,'object',self,'fieldname','dimensions','format','IntMat');
1089+ WriteData(fid,prefix,'data',dimensions,'name','md.stochasticforcing.dimensions','format','IntMat');
1090+ WriteData(fid,prefix,'object',self,'fieldname','default_id','data',self.default_id-1,'format','IntMat','mattype',2); %0-indexed
1091+ WriteData(fid,prefix,'object',self,'fieldname','defaultdimension','format','Integer');
1092 WriteData(fid,prefix,'data',tempcovariance,'name','md.stochasticforcing.covariance','format','DoubleMat');
1093 WriteData(fid,prefix,'object',self,'fieldname','randomflag','format','Boolean');
1094 end
1095@@ -104,7 +157,8 @@
1096 % by the class md.stochasticforcing
1097
1098 list = {...
1099- 'SMBautoregression',...
1100- 'FrontalForcingsRignotAutoregression'
1101+ 'DefaultCalving',...
1102+ 'FrontalForcingsRignotAutoregression',...
1103+ 'SMBautoregression'
1104 };
1105 end % }}}
1106Index: ../trunk-jpl/src/m/classes/stochasticforcing.py
1107===================================================================
1108--- ../trunk-jpl/src/m/classes/stochasticforcing.py (revision 26614)
1109+++ ../trunk-jpl/src/m/classes/stochasticforcing.py (revision 26615)
1110@@ -14,10 +14,11 @@
1111
1112 def __init__(self, *args): # {{{
1113 self.isstochasticforcing = 0
1114- self.fields = np.nan
1115- self.dimensions = np.nan
1116- self.covariance = np.nan
1117- self.randomflag = 1
1118+ self.fields = np.nan
1119+ self.defaultdimension = 0
1120+ self.default_id = np.nan
1121+ self.covariance = np.nan
1122+ self.randomflag = 1
1123
1124 if len(args) == 0:
1125 self.setdefaultparameters()
1126@@ -28,9 +29,12 @@
1127 s = ' stochasticforcing parameters:\n'
1128 s += '{}\n'.format(fielddisplay(self, 'isstochasticforcing', 'is stochasticity activated?'))
1129 s += '{}\n'.format(fielddisplay(self, 'fields', 'fields with stochasticity applied, ex: [\'SMBautoregression\'], or [\'FrontalForcingsRignotAutoregression\']'))
1130+ s += '{}\n'.format(fielddisplay(self, 'defaultdimension', 'dimensionality of the noise terms (does not apply to fields with their specific dimension)'))
1131+ s += '{}\n'.format(fielddisplay(self, 'default_id', 'id of each element for partitioning of the noise terms (does not apply to fields with their specific partition)'))
1132 s += '{}\n'.format(fielddisplay(self, 'covariance', 'covariance matrix for within- and between-fields covariance (units must be squared field units)'))
1133 s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'))
1134 s += 'Available fields:\n'
1135+ s += ' DefaultCalving\n'
1136 s += ' SMBautoregression\n'
1137 s += ' FrontalForcingsRignotAutoregression (thermal forcing)\n'
1138 return s
1139@@ -50,15 +54,15 @@
1140 return md
1141
1142 num_fields = numel(self.fields)
1143- size_tot = np.sum(self.dimensions)
1144
1145- md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
1146- md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', supportedstochforcings()) # VV check here 'cell' (19Oct2021)
1147- md = checkfield(md, 'fieldname', 'stochasticforcing.dimensions', 'NaN', 1, 'Inf', 1, '>', 0, 'size', [num_fields]) # specific dimension for each field; NOTE: As opposed to MATLAB implementation, pass list
1148- md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix
1149- md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1])
1150+ #Check that covariance matrix is positive definite
1151+ try:
1152+ np.linalg.cholesky(self.covariance);
1153+ except:
1154+ error('md.stochasticforcing.covariance is not positive definite');
1155
1156 # Check that all fields agree with the corresponding md class
1157+ checkdefaults = False
1158 for field in self.fields:
1159 if (contains(field, 'SMB')):
1160 if not (type(md.smb) == field):
1161@@ -66,6 +70,32 @@
1162 if (contains(field, 'frontalforcings')):
1163 if not (type(md.frontalforcings) == field):
1164 error('md.frontalforcings does not agree with stochasticforcing field {}'.format(field))
1165+ #Checking for specific dimensions
1166+ if (field!='SMBautoregression' or field!='FrontalForcingsRignotAutoregression'):
1167+ checkdefaults = True #field with non-specific dimensionality
1168+
1169+ #Retrieve sum of all the field dimensionalities
1170+ size_tot = self.defaultdimension*num_fields;
1171+ indSMBar = -1; #about to check for index of SMBautoregression
1172+ indTFar = -1; #about to check for index of FrontalForcingsRignotAutoregression
1173+ if ('SMBautoregression' in self.fields):
1174+ size_tot = size_tot-self.defaultdimension+md.smb.num_basins
1175+ indSMBar = np.where(self.fields=='SMBautoregression')[0][0]
1176+ if ('FrontalForcingsRignotAutoregression' in self.fields):
1177+ size_tot = size_tot-self.defaultdimension+md.frontalforcings.num_basins
1178+ indSMBar = np.where(self.fields=='FrontalForcingsRignotAutoregression')[0][0]
1179+ if (indSMBar!=-1 and indTFar!=-1):
1180+ if((md.smb.ar_timestep!=md.frontalforcings.ar_timestep) and np.any(self.covariance[np.sum(self.dimensions[0:indSMBar]).astype(int):np.sum(self.dimensions[0:indSMBar+1]).astype(int),np.sum(self.dimensions[0:indTFar]).astype(int):np.sum(self.dimensions[0:indTFar+1]).astype(int)]!=0)):
1181+ error('SMBautoregression and FrontalForcingsRignotAutoregression have different ar_timestep and non-zero covariance')
1182+
1183+ md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
1184+ md = checkfield(md, 'fieldname', 'stochasticforcing.fields', 'numel', num_fields, 'cell', 1, 'values', supportedstochforcings())
1185+ #md = checkfield(md, 'fieldname', 'stochasticforcing.dimensions', 'NaN', 1, 'Inf', 1, '>', 0, 'size', [num_fields]) # specific dimension for each field; NOTE: As opposed to MATLAB implementation, pass list
1186+ md = checkfield(md, 'fieldname', 'stochasticforcing.covariance', 'NaN', 1, 'Inf', 1, 'size', [size_tot, size_tot]) # global covariance matrix
1187+ md = checkfield(md, 'fieldname', 'stochasticforcing.randomflag', 'numel', [1], 'values', [0, 1])
1188+ if (checkdefaults):
1189+ md = checkfield(md, 'fieldname', 'stochasticforcing.defaultdimension', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0)
1190+ md = checkfield(md, 'fieldname', 'stochasticforcing.default_id', 'Inf', 1, '>=', 0, '<=', self.defaultdimension, 'size', [md.mesh.numberofelements,1])
1191 return md
1192 # }}}
1193
1194@@ -82,19 +112,30 @@
1195 if not self.isstochasticforcing:
1196 return md
1197 else:
1198+ #Retrieve dimensionality of each field
1199+ dimensions = self.defaultdimension*np.ones(num_fields)
1200+ for ind,field in enumerate(self.fields):
1201+ #Checking for specific dimensions
1202+ if (field=='SMBautoregression'):
1203+ dimensions[ind] = md.smb.num_basins
1204+ if (field=='FrontalForcingsRignotAutoregression'):
1205+ dimensions[ind] = md.frontalforcings.num_basins
1206+
1207 # Scaling covariance matrix (scale column-by-column and row-by-row)
1208- scaledfields = ['SMBautoregression'] # list of fields that need scaling * 1/yts
1209- tempcovariance = self.covariance
1210+ scaledfields = ['DefaultCalving','SMBautoregression'] # list of fields that need scaling * 1/yts
1211+ tempcovariance = np.copy(self.covariance)
1212 for i in range(num_fields):
1213 if self.fields[i] in scaledfields:
1214 inds = range(int(np.sum(self.dimensions[0:i])), int(np.sum(self.dimensions[0:i + 1])))
1215 for row in inds: # scale rows corresponding to scaled field
1216- tempcovariance[row, :] = 1 / yts * self.covariance[row, :]
1217+ tempcovariance[row, :] = 1 / yts * tempcovariance[row, :]
1218 for col in inds: # scale columns corresponding to scaled field
1219- tempcovariance[:, col] = 1 / yts * self.covariance[:, col]
1220+ tempcovariance[:, col] = 1 / yts * tempcovariance[:, col]
1221 WriteData(fid, prefix, 'data', num_fields, 'name', 'md.stochasticforcing.num_fields', 'format', 'Integer')
1222 WriteData(fid, prefix, 'object', self, 'fieldname', 'fields', 'format', 'StringArray')
1223- WriteData(fid, prefix, 'object', self, 'fieldname','dimensions', 'format', 'IntMat')
1224+ WriteData(fid, prefix, 'data', 'dimensions', 'name', 'md.stochasticforcing.dimensions', 'format', 'IntMat')
1225+ WriteData(fid, prefix, 'object', self, 'fieldname', 'default_id', 'format', 'IntMat') #12Nov2021 make sure this is zero-indexed!
1226+ WriteData(fid, prefix, 'object', self, 'fieldname', 'defaultdimension', 'format', 'Integer')
1227 WriteData(fid, prefix, 'data', tempcovariance, 'name', 'md.stochasticforcing.covariance', 'format', 'DoubleMat')
1228 WriteData(fid, prefix, 'object', self, 'fieldname', 'randomflag', 'format', 'Boolean')
1229 # }}}
1230@@ -103,6 +144,7 @@
1231 """ Defines list of fields supported by the class stochasticforcings
1232 """
1233 return [
1234- 'SMBautoregression',
1235- 'FrontalForcingsRignotAutoregression'
1236+ 'DefaultCalving',
1237+ 'FrontalForcingsRignotAutoregression',
1238+ 'SMBautoregression'
1239 ]
Note: See TracBrowser for help on using the repository browser.