source: issm/oecreview/Archive/22819-23185/ISSM-23065-23066.diff

Last change on this file was 23186, checked in by Mathieu Morlighem, 7 years ago

CHG: added Archive/22819-23185
File size: 181.5 KB
RevLine 
[23186]1Index: ../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 23065)
4+++ ../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (revision 23066)
5@@ -91,7 +91,7 @@
6
7 /* Intermediaries */
8 int numvertices = element->GetNumberOfVertices();
9-
10+
11 if(element->IsIceInElement()){
12 for(int i = 0;i<numvertices;i++){
13 vec_mask_ice->SetValue(element->vertices[i]->Sid(),1.,INS_VAL);
14Index: ../trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp
15===================================================================
16--- ../trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp (revision 23065)
17+++ ../trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp (revision 23066)
18@@ -21,7 +21,7 @@
19
20 /*This is the object that we have been chasing for. compute the response and return: */
21 IssmDouble return_value=definition->Response(femmodel);
22-
23+
24 /*cleanup: */
25 xDelete<char>(name);
26
27@@ -30,7 +30,7 @@
28 }
29 xDelete<char>(name);
30 }
31-
32+
33 /*If we are here, did not find the definition for this response, not good!: */
34 _error_("Could not find the response for output definition " << output_string << " because could not find the definition itself!");
35
36@@ -43,7 +43,7 @@
37
38 /*Now, go through the output definitions, and retrieve the object which corresponds to our requested response, output_enum: */
39 for(int i=0;i<output_definitions->Size();i++){
40-
41+
42 //Definition* definition=xDynamicCast<Definition*>(output_definitions->GetObjectByOffset(i));
43 Definition* definition=dynamic_cast<Definition*>(output_definitions->GetObjectByOffset(i));
44
45@@ -52,12 +52,12 @@
46
47 /*This is the object that we have been chasing for. compute the response and return: */
48 IssmDouble return_value=definition->Response(femmodel);
49-
50+
51 /*return:*/
52 return return_value;
53 }
54 }
55-
56+
57 /*If we are here, did not find the definition for this response, not good!: */
58 _error_("Could not find the response for output definition " << EnumToStringx(output_enum)
59 <<" ("<<output_enum<<")"
60Index: ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp
61===================================================================
62--- ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 23065)
63+++ ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp (revision 23066)
64@@ -14,13 +14,12 @@
65 IssmDouble *phi_ungrounding = NULL;
66 Element *element = NULL;
67
68-
69 /*retrieve parameters: */
70 parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
71 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
72
73 if(migration_style==NoneEnum) return;
74-
75+
76 if(VerboseModule()) _printf0_(" Migrating grounding line\n");
77
78 /*Set toolkit to default*/
79Index: ../trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp
80===================================================================
81--- ../trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp (revision 23065)
82+++ ../trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp (revision 23066)
83@@ -77,7 +77,7 @@
84 if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
85 if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
86 }
87-
88+
89 /*Free ressources: */
90 xDelete<char>(toolkittype);
91
92Index: ../trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp
93===================================================================
94--- ../trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp (revision 23065)
95+++ ../trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp (revision 23066)
96@@ -27,7 +27,6 @@
97 int size = 0;
98 for(int i=0;i<num_controls;i++) size+=M[i]*N[i];
99
100-
101 /*2. Allocate vector*/
102 vector=new Vector<IssmDouble>(size);
103
104@@ -59,7 +58,6 @@
105 parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
106 parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
107
108-
109 /*2. Allocate vector*/
110 vector=new Vector<IssmDouble>(num_controls*vertices->NumberOfVertices());
111
112Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp
113===================================================================
114--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp (revision 23065)
115+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp (revision 23066)
116@@ -165,10 +165,10 @@
117 edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
118
119 /*edge can be either IntersectEnum (one and only one intersection between the edge and the segment), ColinearEnum (edge and segment are collinear) and SeparateEnum (no intersection): */
120-
121+
122 /*if (el==318 && contouri==9){
123 _printf_(edge1 << " " << edge2 << " " << edge3 << " " << alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " " << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
124-
125+
126 _printf_("Bool" << (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum));
127 }*/
128
129@@ -198,7 +198,6 @@
130 }
131 segments_dataset->AddObject(new Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
132
133-
134 /*This case is impossible: not quite! */
135 //_printf_(alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " " << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
136 /* _error_("error: a line cannot go through 3 different vertices!");*/
137@@ -225,7 +224,7 @@
138 }
139 }
140 else if( (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum) ){
141-
142+
143 /*if (el==318 && contouri==9){
144 _printf_("hello" << " NodeInElement 0 " << (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) << " NodeInElement 1 " << (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])));
145 }*/
146@@ -250,7 +249,7 @@
147 if (IsIdenticalNode(xnodes[0],ynodes[0],xsegment[0],ysegment[0],tolerance) ||
148 IsIdenticalNode(xnodes[1],ynodes[1],xsegment[0],ysegment[0],tolerance) ||
149 IsIdenticalNode(xnodes[2],ynodes[2],xsegment[0],ysegment[0],tolerance)){
150-
151+
152 /*ok, segments[0] is common to one of our vertices: */
153 coord1=0;
154 if(edge1==IntersectEnum){coord2=alpha1;}
155Index: ../trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp
156===================================================================
157--- ../trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp (revision 23065)
158+++ ../trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp (revision 23066)
159@@ -38,7 +38,7 @@
160
161 IssmDouble d0=1.e+10;
162 IssmDouble xi,yi;
163-
164+
165 //recover vertex position:
166 xi=x[i]; yi=y[i];
167
168Index: ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp
169===================================================================
170--- ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 23065)
171+++ ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 23066)
172@@ -83,7 +83,7 @@
173 int interpolation_type;
174 /*this one is special: we don't specify the type, but let the nature of the inputs dictace.
175 * P0 -> ElementSIdEnum, P1 ->VertexSIdEnum: */
176-
177+
178 /*We go find the input of the first element, and query its interpolation type: */
179 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(0));
180 Input* input=element->GetInput(name);
181Index: ../trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp
182===================================================================
183--- ../trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp (revision 23065)
184+++ ../trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp (revision 23066)
185@@ -20,7 +20,7 @@
186 *element, figure out if they hold the vertex, and the data. If so, return it: */
187 cpu_found=-1;
188 found=0;
189-
190+
191 for(int i=0;i<elements->Size();i++){
192 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
193 found=element->NodalValue(&value,index,natureofdataenum);
194Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp
195===================================================================
196--- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 23065)
197+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 23066)
198@@ -8,7 +8,6 @@
199 #include "../../MeshPartitionx/MeshPartitionx.h"
200 #include "../ModelProcessorx.h"
201
202-
203 #if !defined(_HAVE_ADOLC_)
204 void UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters,Materials* materials, IoModel* iomodel){
205 /*Intermediary*/
206@@ -21,7 +20,6 @@
207 char **controls = NULL;
208 char **cost_functions = NULL;
209
210-
211 /*Fetch parameters: */
212 iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
213 if(control_analysis) iomodel->FindConstant(&num_controls,"md.inversion.num_control_parameters");
214@@ -28,7 +26,7 @@
215
216 /*Now, return if no control*/
217 if(!control_analysis) return;
218-
219+
220 /*Process controls and convert from string to enums*/
221 iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters");
222 if(num_controls<1) _error_("no controls found");
223@@ -47,7 +45,7 @@
224 }
225
226 iomodel->FetchData(3,"md.inversion.cost_functions_coefficients","md.inversion.min_parameters","md.inversion.max_parameters");
227-
228+
229 /*Fetch Observations */
230 iomodel->FindConstant(&domaintype,"md.mesh.domain_type");
231 for(int i=0;i<num_cost_functions;i++){
232@@ -153,7 +151,6 @@
233 _assert_(M==num_independent_objects);
234 iomodel->FetchData(&types,NULL,NULL,"md.autodiff.independent_object_types");
235
236-
237 M_all = xNew<int>(num_independent_objects);
238
239 /*create independent objects, and at the same time, fetch the corresponding independent variables,
240@@ -161,7 +158,7 @@
241 iomodel->FetchData(&independents_fullmin,&M_par,&N_par,"md.autodiff.independent_min_parameters");
242 iomodel->FetchData(&independents_fullmax,&M_par,&N_par,"md.autodiff.independent_max_parameters");
243 iomodel->FetchData(&control_sizes,NULL,NULL,"md.autodiff.independent_control_sizes");
244-
245+
246 int* start_point = NULL;
247 start_point = xNew<int>(num_independent_objects);
248 int counter = 0;
249@@ -179,7 +176,7 @@
250 int input_enum;
251 IssmDouble* independents_min = NULL;
252 IssmDouble* independents_max = NULL;
253-
254+
255 FieldAndEnumFromCode(&input_enum,&iofieldname,names[i]);
256
257 /*Fetch required data*/
258@@ -186,7 +183,7 @@
259 iomodel->FetchData(&independent,&M,&N,iofieldname);
260 _assert_(independent);
261 _assert_(N==control_sizes[i]);
262-
263+
264 independents_min = NULL; independents_min = xNew<IssmDouble>(M*N);
265 independents_max = NULL; independents_max = xNew<IssmDouble>(M*N);
266 for(int m=0;m<M;m++){
267@@ -205,7 +202,7 @@
268 xDelete<IssmDouble>(independent);
269 xDelete<IssmDouble>(independents_min);
270 xDelete<IssmDouble>(independents_max);
271-
272+
273 }
274 else{
275 _error_("Independent cannot be of size " << types[i]);
276Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp
277===================================================================
278--- ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 23065)
279+++ ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 23066)
280@@ -68,11 +68,11 @@
281 /*Intermediaries*/
282 int num_independent_objects,M;
283 char** names = NULL;
284-
285+
286 /*this is done somewhere else*/
287 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_independent_objects",InversionNumControlParametersEnum));
288 parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_dependent_objects",InversionNumCostFunctionsEnum));
289-
290+
291 /*Step1: create controls (independents)*/
292 iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects");
293 _assert_(num_independent_objects>0);
294@@ -95,10 +95,10 @@
295 }
296 xDelete<char*>(cm_responses);
297 parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc));
298-
299+
300 iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.autodiff.independent_scaling_factors");
301 parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_independent_objects));
302-
303+
304 /*cleanup*/
305 for(int i=0;i<num_independent_objects;i++){
306 xDelete<char>(names[i]);
307Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp
308===================================================================
309--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 23065)
310+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 23066)
311@@ -10,7 +10,7 @@
312 void CreateOutputDefinitions(Elements* elements, Parameters* parameters,IoModel* iomodel){
313
314 int i,j;
315-
316+
317 DataSet* output_definitions = NULL;
318 int* output_definition_enums = NULL;
319 int num_output_definitions;
320@@ -64,7 +64,7 @@
321 }
322 else if (output_definition_enums[i]==MisfitEnum){
323 /*Deal with misfits: {{{*/
324-
325+
326 /*misfit variables: */
327 int nummisfits;
328 char** misfit_name_s = NULL;
329@@ -114,7 +114,6 @@
330 else
331 _error_("misfit weight size not supported yet");
332
333-
334 /*First create a misfit object for that specific string (misfit_model_string_s[j]):*/
335 output_definitions->AddObject(new Misfit(misfit_name_s[j],StringToEnumx(misfit_definitionstring_s[j]),StringToEnumx(misfit_model_string_s[j]),StringToEnumx(misfit_observation_string_s[j]),misfit_timeinterpolation_s[j],misfit_local_s[j],StringToEnumx(misfit_weights_string_s[j])));
336
337@@ -158,7 +157,7 @@
338 }
339 else if (output_definition_enums[i]==CfsurfacesquareEnum){
340 /*Deal with cfsurfacesquare: {{{*/
341-
342+
343 /*cfsurfacesquare variables: */
344 int num_cfsurfacesquares;
345 char** cfsurfacesquare_name_s = NULL;
346@@ -213,7 +212,7 @@
347 for(int k=0;k<elements->Size();k++){
348
349 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
350-
351+
352 element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_observation_s[j], iomodel,cfsurfacesquare_observation_M_s[j],cfsurfacesquare_observation_N_s[j],obs_vector_type,StringToEnumx(cfsurfacesquare_observation_string_s[j]),7,SurfaceObservationEnum);
353 element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_weights_s[j], iomodel,cfsurfacesquare_weights_M_s[j],cfsurfacesquare_weights_N_s[j],weight_vector_type,StringToEnumx(cfsurfacesquare_weights_string_s[j]),7,WeightsSurfaceObservationEnum);
354
355@@ -250,7 +249,7 @@
356 }
357 else if (output_definition_enums[i]==CfdragcoeffabsgradEnum){
358 /*Deal with cfdragcoeffabsgrad: {{{*/
359-
360+
361 /*cfdragcoeffabsgrad variables: */
362 int num_cfdragcoeffabsgrads;
363 char** cfdragcoeffabsgrad_name_s = NULL;
364@@ -285,7 +284,7 @@
365 for(int k=0;k<elements->Size();k++){
366
367 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
368-
369+
370 element->DatasetInputAdd(StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),cfdragcoeffabsgrad_weights_s[j], iomodel,cfdragcoeffabsgrad_weights_M_s[j],cfdragcoeffabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),7,WeightsSurfaceObservationEnum);
371
372 }
373@@ -312,7 +311,7 @@
374 }
375 else if (output_definition_enums[i]==CfsurfacelogvelEnum){
376 /*Deal with cfsurfacelogvel: {{{*/
377-
378+
379 /*cfsurfacelogvel variables: */
380 int num_cfsurfacelogvels;
381 char** cfsurfacelogvel_name = NULL;
382@@ -368,7 +367,7 @@
383 for(int k=0;k<elements->Size();k++){
384
385 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
386-
387+
388 element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vxobs[j], iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vxobs_string[j]),7,VxObsEnum);
389 element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vyobs[j], iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vyobs_string[j]),7,VyObsEnum);
390 element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_weights[j], iomodel,cfsurfacelogvel_weights_M[j],cfsurfacelogvel_weights_N[j],weight_vector_type,StringToEnumx(cfsurfacelogvel_weightstring[j]),7,WeightsSurfaceObservationEnum);
391@@ -408,7 +407,7 @@
392 }
393 else if (output_definition_enums[i]==NodalvalueEnum){
394 /*Deal with nodal values: {{{*/
395-
396+
397 /*nodal value variables: */
398 int numnodalvalues;
399 char** nodalvalue_name_s = NULL;
400@@ -427,7 +426,7 @@
401 /*First create a nodalvalue object for that specific enum (nodalvalue_model_enum_s[j]):*/
402 output_definitions->AddObject(new Nodalvalue(nodalvalue_name_s[j],StringToEnumx(nodalvalue_definitionstrings[j]),StringToEnumx(nodalvalue_modelstrings[j]),nodalvalue_node_s[j]-1)); //-1 because matlab to c indexing.
403 }
404-
405+
406 /*Free ressources:*/
407 for(j=0;j<numnodalvalues;j++){
408 char* string=NULL;
409@@ -481,7 +480,7 @@
410 }
411 else if (output_definition_enums[i]==MassconaxpbyEnum){
412 /*Deal with masscon combinations: {{{*/
413-
414+
415 /*masscon variables: */
416 char** masscon_name_s = NULL;
417 char** masscon_definitionstring_s = NULL;
418@@ -616,7 +615,7 @@
419 }
420 output_definitions->AddObject(new Numberedcostfunction(ncf_name_s[j],StringToEnumx(ncf_definitionstring_s[j]),num_cost_functions,cost_function_enums));
421 }
422-
423+
424 /*Free data: */
425 iomodel->DeleteData(2,"md.numberedcostfunction.name","md.numberedcostfunction.definitionstring");
426 xDelete<int>(cost_function_enums);
427Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
428===================================================================
429--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23065)
430+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23066)
431@@ -22,7 +22,7 @@
432 char** requestedoutputs = NULL;
433 char* fieldname = NULL;
434 IssmDouble time;
435-
436+
437 /*parameters for mass flux:*/
438 int mass_flux_num_profiles = 0;
439 bool qmu_mass_flux_present = false;
440@@ -61,7 +61,6 @@
441 parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum));
442 parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum));
443 parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1));
444-
445
446 {/*This is specific to ice...*/
447 parameters->AddObject(iomodel->CopyConstantObject("md.mesh.elementtype",MeshElementtypeEnum));
448Index: ../trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp
449===================================================================
450--- ../trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 23065)
451+++ ../trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp (revision 23066)
452@@ -54,7 +54,6 @@
453 analysis->UpdateElements(elements,iomodel,i,analysis_enum);
454 delete analysis;
455
456-
457 /* Update counters, because we have created more nodes, loads and
458 * constraints, and ids for objects created in next call to CreateDataSets
459 * will need to start at the end of the updated counters: */
460Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
461===================================================================
462--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 23065)
463+++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 23066)
464@@ -34,7 +34,7 @@
465
466 /*Create elements*/
467 if(control_analysis && !adolc_analysis)iomodel->FetchData(2,"md.inversion.min_parameters","md.inversion.max_parameters");
468-
469+
470 switch(iomodel->meshelementtype){
471 case TriaEnum:
472 for(i=0;i<iomodel->numberofelements;i++){
473@@ -122,7 +122,7 @@
474 }
475 break;
476 case MaterialsEnum:
477-
478+
479 //we have several types of materials. Retrieve this info first:
480 iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature");
481
482@@ -234,7 +234,7 @@
483 if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->FetchData(3,"md.mesh.lat","md.mesh.long","md.mesh.r");
484 else iomodel->FetchDataToInput(elements,"md.mesh.scale_factor",MeshScaleFactorEnum,1.);
485 if (isoceancoupling) iomodel->FetchData(2,"md.mesh.lat","md.mesh.long");
486-
487+
488 CreateNumberNodeToElementConnectivity(iomodel,solution_type);
489
490 int lid = 0;
491Index: ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp
492===================================================================
493--- ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp (revision 23065)
494+++ ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp (revision 23066)
495@@ -31,7 +31,6 @@
496 offset += M[i]*N[i];
497 }
498
499-
500 xDelete<int>(control_type);
501 }
502 else{
503Index: ../trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp
504===================================================================
505--- ../trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 23065)
506+++ ../trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 23066)
507@@ -66,4 +66,3 @@
508 }
509 }
510 /*}}}*/
511-
512Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
513===================================================================
514--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23065)
515+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23066)
516@@ -6,7 +6,6 @@
517 #include "../../shared/shared.h"
518 #include "../../toolkits/toolkits.h"
519
520-
521 void SmbGradientsx(FemModel* femmodel){/*{{{*/
522
523 // void SurfaceMassBalancex(hd,agd,ni){
524@@ -317,7 +316,6 @@
525 correct for number of seconds in a year [s/yr]*/
526 smb = smb/yts + anomaly;
527
528-
529 /*Update array accordingly*/
530 smblist[v] = smb;
531
532Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
533===================================================================
534--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23065)
535+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23066)
536@@ -106,7 +106,7 @@
537 /*Free ressouces:*/
538 xDelete(dzT);
539 xDelete(dzB);
540-
541+
542 //---------NEED TO IMPLEMENT A PROPER GRID STRECHING ALGORITHM------------
543
544 /*assign ouput pointers: */
545@@ -201,7 +201,7 @@
546 re=*pre;
547 gdn=*pgdn;
548 gsp=*pgsp;
549-
550+
551 /*only when aIdx = 1 or 2 do we run grainGrowth: */
552 if(aIdx!=1 & aIdx!=2){
553 /*come out as we came in: */
554@@ -220,7 +220,6 @@
555 //set maximum water content by mass to 9 percent (Brun, 1980)
556 for(int i=0;i<m;i++)if(lwc[i]>9.0-Wtol) lwc[i]=9.0;
557
558-
559 /* Calculate temperature gradiant across grid cells
560 * Returns the avereage gradient across the upper and lower grid cell */
561
562@@ -243,7 +242,7 @@
563
564 // take absolute value of temperature gradient
565 for(int i=0;i<m;i++)dT[i]=fabs(dT[i]);
566-
567+
568 /*Snow metamorphism. Depends on value of snow dendricity and wetness of the snowpack: */
569 for(int i=0;i<m;i++){
570 if (gdn[i]>0.0+Gdntol){
571@@ -323,7 +322,7 @@
572 //convert grain size back to effective grain radius:
573 re[i]=gsz[i]/2.0;
574 }
575-
576+
577 /*Free ressources:*/
578 xDelete<IssmDouble>(gsz);
579 xDelete<IssmDouble>(dT);
580@@ -347,7 +346,7 @@
581
582 //// Inputs
583 // aIdx = albedo method to use
584-
585+
586 // Method 0
587 // aValue = direct input value for albedo, override all changes to albedo
588
589@@ -396,7 +395,7 @@
590
591 /*Recover pointers: */
592 a=*pa;
593-
594+
595 //some constants:
596 const IssmDouble dSnow = 300; // density of fresh snow [kg m-3]
597
598@@ -512,7 +511,7 @@
599 void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid) { /*{{{*/
600
601 /* ENGLACIAL THERMODYNAMICS*/
602-
603+
604 /* Description:
605 computes new temperature profile accounting for energy absorption and
606 thermal diffusion.*/
607@@ -536,7 +535,7 @@
608 // OUTPUTS
609 // T: grid cell temperature [k]
610 // EC: evaporation/condensation [kg]
611-
612+
613 /*intermediary: */
614 IssmDouble* K = NULL;
615 IssmDouble* KU = NULL;
616@@ -579,7 +578,7 @@
617 IssmDouble dT_ulw=0.0;
618 IssmDouble dlw=0.0;
619 IssmDouble dT_dlw=0.0;
620-
621+
622 /*outputs:*/
623 IssmDouble EC=0.0;
624 IssmDouble* T=*pT;
625@@ -596,7 +595,7 @@
626
627 //initialize Evaporation - Condenstation
628 EC = 0.0;
629-
630+
631 // check if all SW applied to surface or distributed throught subsurface
632 // swIdx = length(swf) > 1
633
634@@ -608,7 +607,7 @@
635
636 // if V = 0 goes to infinity therfore if V = 0 change
637 if(V<0.01-Dtol)V=0.01;
638-
639+
640 // Bulk-transfer coefficient for turbulent fluxes
641 An = pow(0.4,2) / pow(log(Tz/z0),2); // Bulk-transfer coefficient
642 C = An * dAir * V; // shf & lhf common coefficient
643@@ -626,16 +625,16 @@
644 for(int i=0;i<m;i++) if(d[i]>=dIce-Dtol) K[i] = 9.828 * exp(-5.7E-3*T[i]);
645
646 // THERMAL DIFFUSION COEFFICIENTS
647-
648+
649 // A discretization scheme which truncates the Taylor-Series expansion
650 // after the 3rd term is used. See Patankar 1980, Ch. 3&4
651-
652+
653 // discretized heat equation:
654-
655+
656 // Tp = (Au*Tu° + Ad*Td° + (Ap-Au-Ad)Tp° + S) / Ap
657-
658+
659 // where neighbor coefficients Au, Ap, & Ad are
660-
661+
662 // Au = [dz_u/2KP + dz_p/2KE]^-1
663 // Ad = [dz_d/2KP + dz_d/2KD]^-1
664 // Ap = d*CI*dz/Dt
665@@ -648,7 +647,7 @@
666 KU = xNew<IssmDouble>(m);
667 KD = xNew<IssmDouble>(m);
668 KP = xNew<IssmDouble>(m);
669-
670+
671 KU[0] = UNDEF;
672 KD[m-1] = UNDEF;
673 for(int i=1;i<m;i++) KU[i]= K[i-1];
674@@ -660,7 +659,7 @@
675 dzD = xNew<IssmDouble>(m);
676 dzU[0]=UNDEF;
677 dzD[m-1]=UNDEF;
678-
679+
680 for(int i=1;i<m;i++) dzU[i]= dz[i-1];
681 for(int i=0;i<m-1;i++) dzD[i] = dz[i+1];
682
683@@ -692,7 +691,7 @@
684 Ad[i] = pow((dzD[i]/2.0/KP[i] + dz[i]/2.0/KD[i]),-1.0);
685 Ap[i] = (d[i]*dz[i]*CI)/dt;
686 }
687-
688+
689 // create "neighbor" coefficient matrix
690 Nu = xNew<IssmDouble>(m);
691 Nd = xNew<IssmDouble>(m);
692@@ -702,24 +701,24 @@
693 Nd[i] = Ad[i] / Ap[i];
694 Np[i]= 1.0 - Nu[i] - Nd[i];
695 }
696-
697+
698 // specify boundary conditions: constant flux at bottom
699 Nu[m-1] = 0.0;
700 Np[m-1] = 1.0;
701-
702+
703 // zero flux at surface
704 Np[0] = 1.0 - Nd[0];
705-
706+
707 // Create neighbor arrays for diffusion calculations instead of a tridiagonal matrix
708 Nu[0] = 0.0;
709 Nd[m-1] = 0.0;
710-
711+
712 /* RADIATIVE FLUXES*/
713
714 // energy supplied by shortwave radiation [J]
715 sw = xNew<IssmDouble>(m);
716 for(int i=0;i<m;i++) sw[i]= swf[i] * dt;
717-
718+
719 // temperature change due to SW
720 dT_sw = xNew<IssmDouble>(m);
721 for(int i=0;i<m;i++) dT_sw[i]= sw[i] / (CI * d[i] * dz[i]);
722@@ -745,7 +744,7 @@
723 // PART OF ENERGY CONSERVATION CHECK
724 // store initial temperature
725 //T_init = T;
726-
727+
728 // calculate temperature of snow surface (Ts)
729 // when incoming SW radition is allowed to penetrate the surface,
730 // the modeled energy balance becomes very sensitive to how Ts is
731@@ -753,9 +752,9 @@
732 // less when Ts is taken as the mean of the x top grid cells.
733 Ts = (T[0] + T[1])/2.0;
734 Ts = fmin(CtoK,Ts); // don't allow Ts to exceed 273.15 K (0 degC)
735-
736+
737 //TURBULENT HEAT FLUX
738-
739+
740 // Monin-Obukhov Stability Correction
741 // Reference:
742 // Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
743@@ -763,9 +762,9 @@
744
745 // calculate the Bulk Richardson Number (Ri)
746 Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
747-
748+
749 // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
750-
751+
752 // do not allow Ri to exceed 0.19
753 Ri = fmin(Ri, 0.19);
754
755@@ -777,11 +776,11 @@
756 else {
757 coefM =pow (1.0-18*Ri,-0.25);
758 }
759-
760+
761 // calculate heat/wind 'coef_H' stability factor
762 if (Ri <= -0.03+Ttol) coefH = coefM/1.3;
763 else coefH = coefM;
764-
765+
766 //// Sensible Heat
767 // calculate the sensible heat flux [W m-2](Patterson, 1998)
768 shf = C * CA * (Ta - Ts);
769@@ -800,7 +799,7 @@
770 }
771 else{
772 L = LS; // latent heat of sublimation
773-
774+
775 // for an ice surface Murphy and Koop, 2005 [Equation 7]
776 eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
777 }
778@@ -822,19 +821,19 @@
779 // upward longwave contribution
780 ulw = - (SB * pow(Ts,4.0)* teValue) * dt ;
781 dT_ulw = ulw / TCs;
782-
783+
784 // new grid point temperature
785-
786+
787 //SW penetrates surface
788 for(int j=0;j<m;j++) T[j] = T[j] + dT_sw[j];
789 T[0] = T[0] + dT_dlw + dT_ulw + dT_turb;
790-
791+
792 // temperature diffusion
793 for(int j=0;j<m;j++) T0[1+j]=T[j];
794 for(int j=0;j<m;j++) Tu[j] = T0[j];
795 for(int j=0;j<m;j++) Td[j] = T0[2+j];
796 for(int j=0;j<m;j++) T[j] = (Np[j] * T[j]) + (Nu[j] * Tu[j]) + (Nd[j] * Td[j]);
797-
798+
799 // calculate cumulative evaporation (+)/condensation(-)
800 EC = EC + (EC_day/dts)*dt;
801
802@@ -872,7 +871,6 @@
803 xDelete<IssmDouble>(Tu);
804 xDelete<IssmDouble>(Td);
805
806-
807 /*Assign output pointers:*/
808 *pEC=EC;
809 *pT=T;
810@@ -900,7 +898,7 @@
811 // Outputs
812 // swf = absorbed shortwave radiation [W m-2]
813 //
814-
815+
816 /*outputs: */
817 IssmDouble* swf=NULL;
818
819@@ -909,10 +907,9 @@
820 /*Initialize and allocate: */
821 swf=xNewZeroInit<IssmDouble>(m);
822
823-
824 // SHORTWAVE FUNCTION
825 if (swIdx == 0) {// all sw radation is absorbed in by the top grid cell
826-
827+
828 // calculate surface shortwave radiation fluxes [W m-2]
829 swf[0] = (1.0 - as) * dsw;
830 }
831@@ -947,7 +944,6 @@
832 B2_cum[i]=1.0;
833 }
834
835-
836 // spectral albedos:
837 // 0.3 - 0.8um
838 IssmDouble a0 = fmin(0.98, 1.0 - 1.58 *pow(gsz[0],0.5));
839@@ -988,7 +984,6 @@
840 B2_cum[i+1]=cum2;
841 }
842
843-
844 // flux across grid cell boundaries
845 Qs1 = xNew<IssmDouble>(m+1);
846 Qs2 = xNew<IssmDouble>(m+1);
847@@ -1014,8 +1009,7 @@
848 xDelete<IssmDouble>(exp2);
849 xDelete<IssmDouble>(Qs1);
850 xDelete<IssmDouble>(Qs2);
851-
852-
853+
854 }
855 else{ //function of grid cell density
856
857@@ -1143,7 +1137,6 @@
858 for(int i=0;i<m;i++)massinit+=mInit[i];
859
860 if (P > 0.0+Dtol){
861-
862
863 if (T_air <= CtoK+Ttol){ // if snow
864
865@@ -1209,7 +1202,7 @@
866 mass=0; for(int i=0;i<m;i++)mass+=d[i]*dz[i];
867
868 mass_diff = mass - massinit - P;
869-
870+
871 #ifndef _HAVE_ADOLC_ //avoid round operation. only check in forward mode.
872 mass_diff = round(mass_diff * 100.0)/100.0;
873 if (mass_diff > 0) _error_("mass not conserved in accumulation function");
874@@ -1252,7 +1245,7 @@
875 IssmDouble* EW=NULL;
876 IssmDouble* M=NULL;
877 int* D=NULL;
878-
879+
880 IssmDouble sumM=0.0;
881 IssmDouble sumER=0.0;
882 IssmDouble addE=0.0;
883@@ -1264,7 +1257,7 @@
884 IssmDouble dm=0.0;
885 IssmDouble X=0.0;
886 IssmDouble Wi=0.0;
887-
888+
889 IssmDouble Ztot=0.0;
890 IssmDouble T_bot=0.0;
891 IssmDouble m_bot=0.0;
892@@ -1278,7 +1271,7 @@
893 IssmDouble EI_bot=0.0;
894 IssmDouble EW_bot=0.0;
895 bool top=false;
896-
897+
898 IssmDouble* Zcum=NULL;
899 IssmDouble* dzMin2=NULL;
900 IssmDouble zY2=1.025;
901@@ -1285,7 +1278,7 @@
902 bool lastCellFlag = false;
903 int X1=0;
904 int X2=0;
905-
906+
907 int D_size = 0;
908
909 /*outputs:*/
910@@ -1302,7 +1295,7 @@
911 IssmDouble* gsp=*pgsp;
912 int n=*pn;
913 IssmDouble* R=NULL;
914-
915+
916 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" melt module\n");
917
918 //// INITIALIZATION
919@@ -1334,7 +1327,7 @@
920 // new grid point center temperature, T [K]
921 // for(int i=0;i<n;i++) T[i]-=exsT[i];
922 for(int i=0;i<n;i++) T[i]=fmin(T[i],CtoK);
923-
924+
925 // specify irreducible water content saturation [fraction]
926 const IssmDouble Swi = 0.07; // assumed constant after Colbeck, 1974
927
928@@ -1492,7 +1485,6 @@
929 flxDn[i+1] = inM - F1 - dW[i]; // meltwater out
930 T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
931
932-
933 // check if an ice layer forms
934 if (fabs(d[i] - dIce) < Dtol){
935 // _printf_("ICE LAYER FORMS");
936@@ -1504,10 +1496,9 @@
937 }
938 }
939
940-
941 //// GRID CELL SPACING AND MODEL DEPTH
942 for(int i=0;i<n;i++)if (W[i] < 0.0 - Wtol) _error_("negative pore water generated in melt equations");
943-
944+
945 // delete all cells with zero mass
946 // adjust pore water
947 for(int i=0;i<n;i++)W[i] += dW[i];
948@@ -1532,7 +1523,7 @@
949 celldelete(&EW,n,D,D_size);
950 n=D_size;
951 xDelete<int>(D);
952-
953+
954 // calculate new grid lengths
955 for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
956
957@@ -1544,13 +1535,13 @@
958 //Merging of cells as they are burried under snow.
959 Zcum=xNew<IssmDouble>(n);
960 dzMin2=xNew<IssmDouble>(n);
961-
962+
963 Zcum[0]=dz[0]; // Compute a cumulative depth vector
964-
965+
966 for (int i=1;i<n;i++){
967 Zcum[i]=Zcum[i-1]+dz[i];
968 }
969-
970+
971 for (int i=0;i<n;i++){
972 if (Zcum[i]<=zTop+Dtol){
973 dzMin2[i]=dzMin;
974@@ -1568,7 +1559,7 @@
975 break;
976 }
977 }
978-
979+
980 //Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell)
981 if(X==n-1){
982 lastCellFlag = true;
983@@ -1588,13 +1579,13 @@
984 re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new;
985 gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
986 gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
987-
988+
989 // merge with underlying grid cell and delete old cell
990 dz[i+1] = dz[i] + dz[i+1]; // combine cell depths
991 d[i+1] = m_new / dz[i+1]; // combine top densities
992 W[i+1] = W[i+1] + W[i]; // combine liquid water
993 m[i+1] = m_new; // combine top masses
994-
995+
996 // set cell to 99999 for deletion
997 m[i] = -99999;
998 }
999@@ -1610,7 +1601,7 @@
1000 break;
1001 }
1002 }
1003-
1004+
1005 // adjust variables as a linearly weighted function of mass
1006 IssmDouble m_new = m[X2] + m[X1];
1007 T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
1008@@ -1618,13 +1609,13 @@
1009 re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
1010 gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
1011 gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
1012-
1013+
1014 // merge with underlying grid cell and delete old cell
1015 dz[X1] = dz[X2] + dz[X1]; // combine cell depths
1016 d[X1] = m_new / dz[X1]; // combine top densities
1017 W[X1] = W[X1] + W[X2]; // combine liquid water
1018 m[X1] = m_new; // combine top masses
1019-
1020+
1021 // set cell to 99999 for deletion
1022 m[X2] = -99999;
1023 }
1024@@ -1647,7 +1638,7 @@
1025 celldelete(&EW,n,D,D_size);
1026 n=D_size;
1027 xDelete<int>(D);
1028-
1029+
1030 // check if any of the top 10 cell depths are too large
1031 X=0;
1032 for(int i=9;i>=0;i--){
1033@@ -1656,7 +1647,7 @@
1034 break;
1035 }
1036 }
1037-
1038+
1039 int j=0;
1040 while(j<=X){
1041 if (dz[j] > dzMin*2.0+Dtol){
1042@@ -1685,13 +1676,13 @@
1043 // INTERPRETATION
1044 // calculate total model depth
1045 Ztot = cellsum(dz,n);
1046-
1047+
1048 if (Ztot < zMin-Dtol){
1049 // printf("Total depth < zMin %f \n", Ztot);
1050 // mass and energy to be added
1051 mAdd = m[n-1]+W[n-1];
1052 addE = T[n-1]*m[n-1]*CI + W[n-1]*(LF+CtoK*CI);
1053-
1054+
1055 // add a grid cell of the same size and temperature to the bottom
1056 dz_bot=dz[n-1];
1057 T_bot=T[n-1];
1058@@ -1704,9 +1695,9 @@
1059 gsp_bot=gsp[n-1];
1060 EI_bot=EI[n-1];
1061 EW_bot=EW[n-1];
1062-
1063+
1064 dz_add=dz_bot;
1065-
1066+
1067 newcell(&dz,dz_bot,top,n);
1068 newcell(&T,T_bot,top,n);
1069 newcell(&W,W_bot,top,n);
1070@@ -1726,11 +1717,11 @@
1071 mAdd = -(m[n-1]+W[n-1]);
1072 addE = -(T[n-1]*m[n-1]*CI) - (W[n-1]*(LF+CtoK*CI));
1073 dz_add=-(dz[n-1]);
1074-
1075+
1076 // remove a grid cell from the bottom
1077 D_size=n-1;
1078 D=xNew<int>(D_size);
1079-
1080+
1081 for(int i=0;i<n-1;i++) D[i]=i;
1082 celldelete(&dz,n,D,D_size);
1083 celldelete(&T,n,D,D_size);
1084@@ -1767,7 +1758,7 @@
1085 if (dm !=0 || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
1086 << "dm: " << dm << " dE: " << dE << "\n");
1087 #endif
1088-
1089+
1090 /*Free ressources:*/
1091 if(m)xDelete<IssmDouble>(m);
1092 if(EI)xDelete<IssmDouble>(EI);
1093@@ -1783,7 +1774,7 @@
1094 if(M)xDelete<IssmDouble>(M);
1095 xDelete<IssmDouble>(Zcum);
1096 xDelete<IssmDouble>(dzMin2);
1097-
1098+
1099 /*Assign output pointers:*/
1100 *pM=sumM;
1101 *pR=Rsum;
1102@@ -1861,7 +1852,7 @@
1103 /*output: */
1104 IssmDouble* dz=NULL;
1105 IssmDouble* d=NULL;
1106-
1107+
1108 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" densification module\n");
1109
1110 /*Recover pointers: */
1111@@ -1870,7 +1861,7 @@
1112
1113 // initial mass
1114 IssmDouble* mass_init = xNew<IssmDouble>(m);for(int i=0;i<m;i++) mass_init[i]=d[i] * dz[i];
1115-
1116+
1117 /*allocations and initialization of overburden pressure and factor H: */
1118 IssmDouble* cumdz = xNew<IssmDouble>(m-1);
1119 cumdz[0]=dz[0];
1120@@ -1879,12 +1870,11 @@
1121 IssmDouble* obp = xNew<IssmDouble>(m);
1122 obp[0]=0.0;
1123 for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];
1124-
1125+
1126 // calculate new snow/firn density for:
1127 // snow with densities <= 550 [kg m-3]
1128 // snow with densities > 550 [kg m-3]
1129-
1130-
1131+
1132 for(int i=0;i<m;i++){
1133 switch (denIdx){
1134 case 1: // Herron and Langway (1980)
1135@@ -2003,7 +1993,7 @@
1136 IssmDouble shf=0.0;
1137 IssmDouble lhf=0.0;
1138 IssmDouble EC=0.0;
1139-
1140+
1141 if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" turbulentFlux module\n");
1142
1143 // calculated air density [kg/m3]
1144@@ -2028,7 +2018,7 @@
1145 Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2));
1146
1147 // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
1148-
1149+
1150 // do not allow Ri to exceed 0.19
1151 Ri = fmin(Ri, 0.19);
1152
1153@@ -2072,7 +2062,6 @@
1154 eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
1155 }
1156
1157-
1158 // Latent heat flux [W m-2]
1159 lhf = C * L * (eAir - eS) * 0.622 / pAir;
1160
1161Index: ../trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp
1162===================================================================
1163--- ../trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp (revision 23065)
1164+++ ../trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp (revision 23066)
1165@@ -106,7 +106,7 @@
1166
1167 femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
1168 femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
1169-
1170+
1171 IssmDouble* boxareas=xNewZeroInit<IssmDouble>(num_basins*maxbox);
1172
1173 for(int i=0;i<femmodel->elements->Size();i++){
1174@@ -121,7 +121,7 @@
1175 IssmDouble* sumareas =xNew<IssmDouble>(num_basins*maxbox);
1176 ISSM_MPI_Allreduce(boxareas,sumareas,num_basins*maxbox,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
1177 //if(sumareas[0]==0){_error_("No elements in box 0, basal meltrates will be 0. Consider decreasing md.basalforcings.maxboxcount or refining your mesh!");}
1178-
1179+
1180 /*Update parameters to keep track of the new areas in future calculations*/
1181 femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoBoxAreaEnum,sumareas,maxbox*num_basins));
1182
1183Index: ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp
1184===================================================================
1185--- ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp (revision 23065)
1186+++ ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp (revision 23066)
1187@@ -51,7 +51,7 @@
1188 double x1,y1;
1189 double x2,y2;
1190 double mind;
1191-
1192+
1193 for(int i=i0;i<i1;i++){
1194
1195 /*Get current point*/
1196@@ -92,7 +92,7 @@
1197 /*Beyond the 'w' end of the segment*/
1198 return sqrt(pow(x2-x0,2)+pow(y2-y0,2));
1199 }
1200-
1201+
1202 /*Projection falls on segment*/
1203 double projx= x1 + t* (x2-x1);
1204 double projy= y1 + t* (y2-y1);
1205Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp
1206===================================================================
1207--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp (revision 23065)
1208+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp (revision 23066)
1209@@ -19,7 +19,7 @@
1210 Vector<IssmDouble>* ys = NULL;
1211 int configuration_type;
1212 IssmDouble solver_residue_threshold;
1213-
1214+
1215 /*solver convergence test: */
1216 IssmDouble nKUF;
1217 IssmDouble nF;
1218@@ -38,9 +38,9 @@
1219 femmodel->profiler->Start(SOLVER);
1220 Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters);
1221 femmodel->profiler->Stop(SOLVER);
1222-
1223+
1224 /*Check that the solver converged nicely: */
1225-
1226+
1227 //compute KUF = KU - F = K*U - F
1228 KU=uf->Duplicate(); Kff->MatMult(uf,KU);
1229 KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0);
1230@@ -52,7 +52,6 @@
1231
1232 if(!xIsNan(solver_residue_threshold) && solver_residue>solver_residue_threshold)_error_(" solver residue too high!: norm(KU-F)/norm(F)=" << solver_residue << "\n");
1233
1234-
1235 //clean up
1236 delete KU; delete KUF;
1237 delete Kff; delete pf; delete df;
1238@@ -62,7 +61,7 @@
1239 // ADOLC_DUMP_MACRO(uf->svector->vector[i]);
1240 // }
1241 //#endif
1242-
1243+
1244 Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
1245 InputUpdateFromSolutionx(femmodel,ug);
1246 delete ug;
1247Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
1248===================================================================
1249--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 23065)
1250+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 23066)
1251@@ -184,7 +184,7 @@
1252
1253 /*Initial guess*/
1254 VecZeroEntries(udot);
1255-
1256+
1257 /*Richardson's iterations*/
1258 for(int i=0;i<5;i++){
1259 MatMult(Mc,udot,temp1);
1260@@ -401,7 +401,7 @@
1261 /*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
1262 CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
1263 delete uf;
1264-
1265+
1266 /*Go solve lower order solution*/
1267 femmodel->profiler->Start(SOLVER);
1268 SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);
1269Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp
1270===================================================================
1271--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp (revision 23065)
1272+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp (revision 23066)
1273@@ -60,7 +60,7 @@
1274 }
1275
1276 count=1;
1277-
1278+
1279 for(;;){
1280 delete tf_old;tf_old=tf;
1281
1282Index: ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp
1283===================================================================
1284--- ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp (revision 23065)
1285+++ ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp (revision 23066)
1286@@ -177,7 +177,6 @@
1287 rhs=xNew<IssmDouble>(pf_M);
1288 #endif
1289
1290-
1291 recvcounts=xNew<int>(num_procs);
1292 displs=xNew<int>(num_procs);
1293
1294@@ -262,7 +261,6 @@
1295 rhs=xNew<IssmDouble>(pf_M);
1296 #endif
1297
1298-
1299 recvcounts=xNew<int>(num_procs);
1300 displs=xNew<int>(num_procs);
1301
1302Index: ../trunk-jpl/src/c/cores/controlad_core.cpp
1303===================================================================
1304--- ../trunk-jpl/src/c/cores/controlad_core.cpp (revision 23065)
1305+++ ../trunk-jpl/src/c/cores/controlad_core.cpp (revision 23066)
1306@@ -104,7 +104,7 @@
1307 case 7: _printf0_(" <g,d> > 0 or <y,s> <0\n"); break;
1308 default: _printf0_(" Unknown end condition\n");
1309 }
1310-
1311+
1312 /*Save results:*/
1313 femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,G,n,1,0,0.0));
1314 femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffXpEnum,X,intn,1,0,0.0));
1315@@ -129,7 +129,7 @@
1316 int intn;
1317 IssmPDouble* X=NULL;
1318 int i;
1319-
1320+
1321 /*Now things get complicated. The femmodel we recovered did not initialize an AD trace, so we can't compute gradients with it. We are going to recreate
1322 *a new femmodel, identical in all aspects to the first one, with trace on though, which will allow us to run the forward mode and get the gradient
1323 in one run of the solution core. So first recover the filenames required for the FemModel constructor, then call a new ad tailored constructor:*/
1324@@ -142,7 +142,6 @@
1325
1326 femmodel=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, restartfilename,IssmComm::GetComm(), femmodel->solution_type,NULL);
1327
1328-
1329 /*Get initial guess:*/
1330 femmodel->parameters->FindParam(&Xd,&intn,AutodiffXpEnum);
1331 X=xNew<IssmPDouble>(intn); for(i=0;i<intn;i++) X[i]=reCast<IssmPDouble>(Xd[i]);
1332@@ -174,7 +173,7 @@
1333 FemModel *femmodel = NULL;
1334 IssmDouble pfd;
1335 int i;
1336-
1337+
1338 /*Recover Femmodel*/
1339 femmodel = (FemModel*)dzs;
1340
1341@@ -194,7 +193,7 @@
1342 IssmDouble* Jlist = NULL;
1343 femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast<IssmPDouble>(pfd);
1344 _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<" | ");
1345-
1346+
1347 /*Compute gradient using AD. Gradient is in the results after the ad_core is called*/
1348 adgradient_core(femmodel);
1349
1350@@ -210,10 +209,10 @@
1351
1352 /*MPI broadcast results:*/
1353 ISSM_MPI_Bcast(G2,*n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
1354-
1355+
1356 /*Send gradient to m1qn3 core: */
1357 for(long i=0;i<*n;i++) G[i] = G2[i];
1358-
1359+
1360 /*Constrain X and G*/
1361 IssmDouble Gnorm = 0.;
1362 for(long i=0;i<*n;i++) Gnorm += G[i]*G[i];
1363@@ -227,7 +226,7 @@
1364 /*Clean-up and return*/
1365 xDelete<IssmDouble>(Jlist);
1366 xDelete<IssmPDouble>(G2);
1367-
1368+
1369 xDelete<char>(rootpath);
1370 xDelete<char>(inputfilename);
1371 xDelete<char>(outputfilename);
1372@@ -250,7 +249,7 @@
1373 FemModel *femmodelad = NULL;
1374 IssmDouble pfd;
1375 int i;
1376-
1377+
1378 /*Recover Femmodel*/
1379 femmodel = (FemModel*)dzs;
1380
1381@@ -270,7 +269,7 @@
1382
1383 femmodelad=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, restartfilename,IssmComm::GetComm(), femmodel->solution_type,X);
1384 femmodel=femmodelad; //We can do this, because femmodel is being called from outside, not by reference, so we won't erase it
1385-
1386+
1387 /*Recover some parameters*/
1388 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
1389
1390@@ -283,7 +282,7 @@
1391 IssmDouble* Jlist = NULL;
1392 femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast<IssmPDouble>(pfd);
1393 _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<" | ");
1394-
1395+
1396 /*Compute gradient using AD. Gradient is in the results after the ad_core is called*/
1397 adgradient_core(femmodel);
1398
1399@@ -299,10 +298,10 @@
1400
1401 /*MPI broadcast results:*/
1402 ISSM_MPI_Bcast(G2,*n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
1403-
1404+
1405 /*Send gradient to m1qn3 core: */
1406 for(long i=0;i<*n;i++) G[i] = G2[i];
1407-
1408+
1409 /*Recover Gnorm: */
1410 IssmDouble Gnorm = 0.;
1411 for(int i=0;i<*n;i++) Gnorm += G[i]*G[i];
1412@@ -316,7 +315,6 @@
1413 /*Clean-up and return*/
1414 xDelete<IssmDouble>(Jlist);
1415 xDelete<IssmPDouble>(G2);
1416-
1417 xDelete<char>(rootpath);
1418 xDelete<char>(inputfilename);
1419 xDelete<char>(outputfilename);
1420Index: ../trunk-jpl/src/c/cores/gia_core.cpp
1421===================================================================
1422--- ../trunk-jpl/src/c/cores/gia_core.cpp (revision 23065)
1423+++ ../trunk-jpl/src/c/cores/gia_core.cpp (revision 23066)
1424@@ -52,7 +52,7 @@
1425 int outputs[2] = {GiaWEnum,GiadWdtEnum};
1426 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
1427 }
1428-
1429+
1430 xDelete<IssmDouble>(x);
1431 xDelete<IssmDouble>(y);
1432 }
1433Index: ../trunk-jpl/src/c/cores/bmb_core.cpp
1434===================================================================
1435--- ../trunk-jpl/src/c/cores/bmb_core.cpp (revision 23065)
1436+++ ../trunk-jpl/src/c/cores/bmb_core.cpp (revision 23066)
1437@@ -11,7 +11,6 @@
1438
1439 void bmb_core(FemModel* femmodel){
1440
1441-
1442 /*First, get BMB model from parameters*/
1443 int basalforcing_model;
1444 bool isplume = false;
1445Index: ../trunk-jpl/src/c/cores/dakota_core.cpp
1446===================================================================
1447--- ../trunk-jpl/src/c/cores/dakota_core.cpp (revision 23065)
1448+++ ../trunk-jpl/src/c/cores/dakota_core.cpp (revision 23066)
1449@@ -234,7 +234,6 @@
1450 /*}}}*/
1451 void dakota_core(FemModel* femmodel){ /*{{{*/
1452
1453-
1454 int my_rank;
1455 char *dakota_input_file = NULL;
1456 char *dakota_output_file = NULL;
1457Index: ../trunk-jpl/src/c/cores/esa_core.cpp
1458===================================================================
1459--- ../trunk-jpl/src/c/cores/esa_core.cpp (revision 23065)
1460+++ ../trunk-jpl/src/c/cores/esa_core.cpp (revision 23066)
1461@@ -22,7 +22,7 @@
1462 int solution_type;
1463 int numoutputs = 0;
1464 char **requested_outputs = NULL;
1465-
1466+
1467 /*additional parameters: */
1468 int gsize;
1469 bool spherical=true;
1470@@ -79,7 +79,7 @@
1471 U_east = new Vector<IssmDouble>(gsize);
1472 U_x = new Vector<IssmDouble>(gsize);
1473 U_y = new Vector<IssmDouble>(gsize);
1474-
1475+
1476 /*call the geodetic main modlule:*/
1477 if(domaintype==Domain3DsurfaceEnum){
1478 femmodel->EsaGeodetic3D(U_radial,U_north,U_east,latitude,longitude,radius,xx,yy,zz);
1479@@ -94,13 +94,13 @@
1480 InputUpdateFromVectorx(femmodel,U_radial,EsaUmotionEnum,VertexSIdEnum); // radial displacement
1481 InputUpdateFromVectorx(femmodel,U_north,EsaNmotionEnum,VertexSIdEnum); // north motion
1482 InputUpdateFromVectorx(femmodel,U_east,EsaEmotionEnum,VertexSIdEnum); // east motion
1483-
1484+
1485 if(save_results){
1486 if(VerboseSolution()) _printf0_(" saving results\n");
1487 femmodel->parameters->FindParam(&requested_outputs,&numoutputs,EsaRequestedOutputsEnum);
1488 femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
1489 }
1490-
1491+
1492 if(solution_type==EsaSolutionEnum)femmodel->RequestedDependentsx();
1493
1494 /*Free ressources:*/
1495@@ -114,4 +114,3 @@
1496
1497 }
1498 /*}}}*/
1499-
1500Index: ../trunk-jpl/src/c/cores/masstransport_core.cpp
1501===================================================================
1502--- ../trunk-jpl/src/c/cores/masstransport_core.cpp (revision 23065)
1503+++ ../trunk-jpl/src/c/cores/masstransport_core.cpp (revision 23066)
1504@@ -30,7 +30,7 @@
1505 femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
1506 femmodel->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
1507 if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
1508-
1509+
1510 if(VerboseSolution()) _printf0_(" computing mass transport\n");
1511
1512 /*Transport mass or free surface*/
1513Index: ../trunk-jpl/src/c/cores/sealevelrise_core.cpp
1514===================================================================
1515--- ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 23065)
1516+++ ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 23066)
1517@@ -12,17 +12,16 @@
1518 /*cores:*/
1519 void sealevelrise_core(FemModel* femmodel){ /*{{{*/
1520
1521-
1522 /*Parameters, variables:*/
1523 bool save_results;
1524 bool isslr=0;
1525 int solution_type;
1526-
1527+
1528 /*Retrieve parameters:*/
1529 femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
1530 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
1531 femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
1532-
1533+
1534 /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
1535 if(solution_type==SealevelriseSolutionEnum)isslr=1;
1536
1537@@ -40,7 +39,7 @@
1538
1539 /*Run steric core for sure:*/
1540 steric_core(femmodel);
1541-
1542+
1543 /*Save results: */
1544 if(save_results){
1545 int numoutputs = 0;
1546@@ -58,7 +57,6 @@
1547 /*}}}*/
1548 void geodetic_core(FemModel* femmodel){ /*{{{*/
1549
1550-
1551 /*variables:*/
1552 Vector<IssmDouble> *RSLg = NULL;
1553 Vector<IssmDouble> *RSLg_rate = NULL;
1554@@ -87,7 +85,7 @@
1555
1556 /*Should we even be here?:*/
1557 femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return;
1558-
1559+
1560 /*Verbose: */
1561 if(VerboseSolution()) _printf0_(" computing geodetic sea level rise\n");
1562
1563@@ -154,7 +152,7 @@
1564
1565 /*recover N_esa = U_esa + RSLg:*/
1566 N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1);
1567-
1568+
1569 /*transform these values into rates (as we only run this once each frequency turn:*/
1570 N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency));
1571 U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency));
1572@@ -161,7 +159,7 @@
1573 N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency));
1574 U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency));
1575 RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency));
1576-
1577+
1578 /*get some results into elements:{{{*/
1579 InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum);
1580 InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum);
1581@@ -181,7 +179,6 @@
1582 } /*}}}*/
1583 }
1584
1585-
1586 if(iscoupler){
1587 /*transfer sea level back to ice caps:*/
1588 TransferSealevel(femmodel,SealevelNEsaRateEnum);
1589@@ -227,7 +224,7 @@
1590 Vector<IssmDouble> *N_esa_rate= NULL;
1591 Vector<IssmDouble> *U_gia_rate= NULL;
1592 Vector<IssmDouble> *N_gia_rate= NULL;
1593-
1594+
1595 /*parameters: */
1596 bool isslr=0;
1597 int solution_type;
1598@@ -242,7 +239,7 @@
1599
1600 /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
1601 if(solution_type==SealevelriseSolutionEnum)isslr=1;
1602-
1603+
1604 /*Should we be here?:*/
1605 if(!isslr)return;
1606
1607@@ -290,7 +287,7 @@
1608 }
1609 /*}}}*/
1610 Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/
1611-
1612+
1613 /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
1614
1615 Vector<IssmDouble> *RSLgi = NULL;
1616@@ -317,7 +314,7 @@
1617
1618 /*Figure out size of g-set deflection vector and allocate solution vector: */
1619 gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
1620-
1621+
1622 /*Initialize:*/
1623 RSLgi = new Vector<IssmDouble>(gsize);
1624
1625@@ -330,11 +327,10 @@
1626 /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the
1627 * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/
1628 RSLgi->Shift(-eustatic-RSLgi_oceanaverage);
1629-
1630+
1631 /*save eustatic value for results: */
1632 femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic));
1633
1634-
1635 /*clean up and return:*/
1636 xDelete<IssmDouble>(latitude);
1637 xDelete<IssmDouble>(longitude);
1638@@ -346,7 +342,6 @@
1639 /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
1640 non eustatic core of the SLR solution */
1641
1642-
1643 Vector<IssmDouble> *RSLg = NULL;
1644 Vector<IssmDouble> *RSLg_old = NULL;
1645
1646@@ -379,7 +374,7 @@
1647 femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum);
1648 femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
1649 femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
1650-
1651+
1652 /*computational flag: */
1653 femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum);
1654
1655@@ -388,7 +383,7 @@
1656
1657 /*Figure out size of g-set deflection vector and allocate solution vector: */
1658 gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
1659-
1660+
1661 /*Initialize:*/
1662 RSLg = new Vector<IssmDouble>(gsize);
1663 RSLg->Assemble();
1664@@ -399,7 +394,7 @@
1665
1666 count=1;
1667 converged=false;
1668-
1669+
1670 /*Start loop: */
1671 for(;;){
1672
1673@@ -412,7 +407,7 @@
1674
1675 /*call the non eustatic module: */
1676 femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop);
1677-
1678+
1679 /*assemble solution vector: */
1680 RSLgo->Assemble();
1681
1682@@ -422,7 +417,7 @@
1683 RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
1684 femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz,latitude,longitude,radius);
1685 RSLgo_rot->Assemble();
1686-
1687+
1688 /*save changes in inertia tensor as results: */
1689 femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz));
1690 femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz));
1691@@ -430,10 +425,10 @@
1692
1693 RSLgo->AXPY(RSLgo_rot,1);
1694 }
1695-
1696+
1697 /*we need to average RSLgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
1698 RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo);
1699-
1700+
1701 /*RSLg is the sum of the eustatic term, and the ocean terms: */
1702 RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1);
1703 RSLg->Shift(-RSLgo_oceanaverage);
1704@@ -454,10 +449,10 @@
1705 converged=true;
1706 break;
1707 }
1708-
1709+
1710 /*some minor verbosing adjustment:*/
1711 if(count>1)verboseconvolution=false;
1712-
1713+
1714 }
1715 if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n");
1716
1717@@ -473,7 +468,7 @@
1718 Vector<IssmDouble> *U_esa = NULL;
1719 Vector<IssmDouble> *U_north_esa = NULL;
1720 Vector<IssmDouble> *U_east_esa = NULL;
1721-
1722+
1723 /*parameters: */
1724 int configuration_type;
1725 int gsize;
1726@@ -506,7 +501,7 @@
1727 /*retrieve geometric information: */
1728 VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
1729 VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
1730-
1731+
1732 /*call the elastic main modlule:*/
1733 femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,latitude,longitude,radius,xx,yy,zz,loop,horiz);
1734
1735@@ -531,7 +526,7 @@
1736 /*variables:*/
1737 Vector<IssmDouble> *U_gia = NULL;
1738 Vector<IssmDouble> *N_gia = NULL;
1739-
1740+
1741 /*parameters:*/
1742 int frequency;
1743 IssmDouble dt;
1744@@ -569,7 +564,7 @@
1745 IssmDouble* forcing=NULL;
1746 Vector<IssmDouble>* forcingglobal=NULL;
1747 int* nvs=NULL;
1748-
1749+
1750 /*transition vectors:*/
1751 IssmDouble** transitions=NULL;
1752 int ntransitions;
1753@@ -576,7 +571,7 @@
1754 int* transitions_m=NULL;
1755 int* transitions_n=NULL;
1756 int nv;
1757-
1758+
1759 /*communicators:*/
1760 ISSM_MPI_Comm tocomm;
1761 ISSM_MPI_Comm* fromcomms=NULL;
1762@@ -590,7 +585,7 @@
1763 femmodel->parameters->FindParam(&earthid,EarthIdEnum);
1764 femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
1765 my_rank=IssmComm::GetRank();
1766-
1767+
1768 /*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */
1769 if(modelid==earthid){
1770 GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
1771@@ -619,7 +614,7 @@
1772 forcings[i]=xNew<IssmDouble>(nvs[i]);
1773 ISSM_MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status);
1774 }
1775-
1776+
1777 }
1778 else{
1779 ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, modelid, tocomm);
1780@@ -630,7 +625,7 @@
1781
1782 /*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/
1783 if(modelid==earthid){
1784-
1785+
1786 /*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions.
1787 *First, build the global delta thickness vector in the earth model: */
1788 nv=femmodel->vertices->NumberOfVertices();
1789@@ -652,7 +647,6 @@
1790 /*build index to plug values: */
1791 int* index=xNew<int>(M); for(int i=0;i<M;i++)index[i]=reCast<int>(transition[i])-1; //matlab indexing!
1792
1793-
1794 /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular
1795 * ice cap: */
1796 forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL);
1797@@ -662,7 +656,7 @@
1798
1799 /*Assemble vector:*/
1800 forcingglobal->Assemble();
1801-
1802+
1803 /*Plug into elements:*/
1804 InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum);
1805 }
1806@@ -695,7 +689,7 @@
1807 /*forcing being transferred from earth to ice caps: */
1808 IssmDouble* forcing=NULL;
1809 IssmDouble* forcingglobal=NULL;
1810-
1811+
1812 /*transition vectors:*/
1813 IssmDouble** transitions=NULL;
1814 int ntransitions;
1815@@ -702,7 +696,7 @@
1816 int* transitions_m=NULL;
1817 int* transitions_n=NULL;
1818 int nv;
1819-
1820+
1821 /*communicators:*/
1822 ISSM_MPI_Comm fromcomm;
1823 ISSM_MPI_Comm* tocomms=NULL;
1824@@ -717,7 +711,7 @@
1825 femmodel->parameters->FindParam(&earthid,EarthIdEnum);
1826 femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
1827 my_rank=IssmComm::GetRank();
1828-
1829+
1830 /*retrieve the inter communicators that will be used to send data from earth to ice caps:*/
1831 if(modelid==earthid){
1832 GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
1833@@ -732,7 +726,6 @@
1834 //femmodel->parameters->FindParam((int*)(&fromcomm), IcecapToEarthCommEnum);
1835 }
1836
1837-
1838 /*Retrieve sea-level on earth model: */
1839 if(modelid==earthid){
1840 nv=femmodel->vertices->NumberOfVertices();
1841@@ -741,12 +734,12 @@
1842
1843 /*Send the forcing to the ice caps:{{{*/
1844 if(my_rank==0){
1845-
1846+
1847 if(modelid==earthid){
1848-
1849+
1850 /*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */
1851 femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
1852-
1853+
1854 if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!");
1855
1856 for(int i=0;i<earthid;i++){
1857@@ -803,7 +796,7 @@
1858 Vector<IssmDouble> *deltathickness = NULL;
1859 Vector<IssmDouble> *cumdeltathickness = NULL;
1860 int nv;
1861-
1862+
1863 if(VerboseSolution()) _printf0_(" computing earth mass transport\n");
1864
1865 /*This mass transport module for the Earth is because we might have thickness variations as spcs
1866@@ -839,7 +832,7 @@
1867
1868 } /*}}}*/
1869 void slrconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
1870-
1871+
1872 bool converged=true;
1873 IssmDouble ndS,nS;
1874 Vector<IssmDouble> *dRSLg = NULL;
1875@@ -847,15 +840,14 @@
1876 //compute norm(du) and norm(u) if requested
1877 dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
1878 ndS=dRSLg->Norm(NORM_TWO);
1879-
1880+
1881 if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");
1882-
1883+
1884 if(!xIsNan<IssmDouble>(eps_rel)){
1885 nS=RSLg_old->Norm(NORM_TWO);
1886 if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
1887 }
1888
1889-
1890 //clean up
1891 delete dRSLg;
1892
1893Index: ../trunk-jpl/src/c/cores/stressbalance_core.cpp
1894===================================================================
1895--- ../trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 23065)
1896+++ ../trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 23066)
1897@@ -21,7 +21,6 @@
1898 int numoutputs = 0;
1899 char **requested_outputs = NULL;
1900 Analysis *analysis = NULL;
1901-
1902
1903 /* recover parameters:*/
1904 femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
1905@@ -35,7 +34,7 @@
1906 femmodel->parameters->FindParam(&numoutputs,StressbalanceNumRequestedOutputsEnum);
1907 femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum);
1908 if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,StressbalanceRequestedOutputsEnum);
1909-
1910+
1911 if(VerboseSolution()) _printf0_(" computing new velocity\n");
1912
1913 /*Compute slopes if necessary */
1914@@ -78,12 +77,11 @@
1915 delete analysis;
1916 }
1917
1918-
1919 if(save_results){
1920 if(VerboseSolution()) _printf0_(" saving results\n");
1921 femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
1922 }
1923-
1924+
1925 if(solution_type==StressbalanceSolutionEnum && !control_analysis)femmodel->RequestedDependentsx();
1926
1927 /*Free ressources:*/
1928Index: ../trunk-jpl/src/c/cores/adgradient_core.cpp
1929===================================================================
1930--- ../trunk-jpl/src/c/cores/adgradient_core.cpp (revision 23065)
1931+++ ../trunk-jpl/src/c/cores/adgradient_core.cpp (revision 23066)
1932@@ -38,7 +38,7 @@
1933 /*intermediary: */
1934 IssmPDouble *aWeightVector=NULL;
1935 IssmPDouble *weightVectorTimesJac=NULL;
1936-
1937+
1938 /*output: */
1939 IssmPDouble *totalgradient=NULL;
1940
1941@@ -52,7 +52,7 @@
1942
1943 /*First, stop tracing: */
1944 trace_off();
1945-
1946+
1947 if(VerboseAutodiff()){ /*{{{*/
1948 size_t tape_stats[15];
1949 tapestats(my_rank,tape_stats); //reading of tape statistics
1950@@ -95,11 +95,11 @@
1951 }
1952 delete [] sstats;
1953 } /*}}}*/
1954-
1955+
1956 /*retrieve parameters: */
1957 femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
1958 femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
1959-
1960+
1961 /*if no dependents, no point in running a driver: */
1962 if(!(num_dependents*num_independents)) return;
1963
1964@@ -108,10 +108,10 @@
1965 if (my_rank!=0){
1966 num_dependents=0; num_independents=0;
1967 }
1968-
1969+
1970 /*retrieve state variable: */
1971 femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
1972-
1973+
1974 /* driver argument */
1975 xp=xNew<double>(num_independents);
1976 for(i=0;i<num_independents;i++){
1977@@ -130,13 +130,13 @@
1978
1979 /*Initialize outputs: */
1980 totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
1981-
1982+
1983 for(aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
1984
1985 /*initialize direction index in the weights vector: */
1986 aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
1987 if (my_rank==0) aWeightVector[aDepIndex]=1.0;
1988-
1989+
1990 /*initialize output gradient: */
1991 weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
1992
1993@@ -161,12 +161,12 @@
1994 xDelete(weightVectorTimesJac);
1995 xDelete(aWeightVector);
1996 }
1997-
1998+
1999 /*add totalgradient to results*/
2000 femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,totalgradient,num_independents,1,1,0.0));
2001
2002 if(VerboseAutodiff())_printf0_(" end ad core\n");
2003-
2004+
2005 /* delete the allocated space for the parameters and free ressources:{{{*/
2006 xDelete(anEDF_for_solverx_p->dp_x);
2007 xDelete(anEDF_for_solverx_p->dp_X);
2008Index: ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp
2009===================================================================
2010--- ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 23065)
2011+++ ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 23066)
2012@@ -125,7 +125,7 @@
2013 ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
2014 SetControlInputsFromVectorx(femmodel,aX);
2015 xDelete<IssmDouble>(aX);
2016-
2017+
2018 /*Compute solution (forward)*/
2019 void (*solutioncore)(FemModel*)=NULL;
2020 CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
2021@@ -391,7 +391,7 @@
2022 /*Optimization criterions*/
2023 long niter = long(maxsteps); /*Maximum number of iterations*/
2024 long nsim = long(maxiter);/*Maximum number of function calls*/
2025-
2026+
2027 /*Get initial guess*/
2028 Vector<double> *Xpetsc = NULL;
2029
2030@@ -400,7 +400,7 @@
2031 Xpetsc->GetSize(&intn);
2032 delete Xpetsc;
2033 //_assert_(intn==numberofvertices*num_controls);
2034-
2035+
2036 /*Get problem dimension and initialize gradient and initial guess*/
2037 long n = long(intn);
2038 G = xNew<double>(n);
2039@@ -414,7 +414,7 @@
2040 }
2041 N_add+=N[c];
2042 }
2043-
2044+
2045 /*Allocate m1qn3 working arrays (see documentation)*/
2046 long m = 100;
2047 long ndz = 4*n+m*(2*n+1);
2048@@ -470,7 +470,7 @@
2049 }
2050 N_add +=N[c];
2051 }
2052-
2053+
2054 /*Set X as our new control*/
2055 IssmDouble* aX=xNew<IssmDouble>(intn);
2056 IssmDouble* aG=xNew<IssmDouble>(intn);
2057@@ -482,10 +482,9 @@
2058
2059 ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
2060 SetControlInputsFromVectorx(femmodel,aX);
2061-
2062+
2063 xDelete(aX);
2064
2065-
2066 if (solution_type == TransientSolutionEnum){
2067 int step = 1;
2068 femmodel->parameters->SetParam(step,StepEnum);
2069@@ -505,7 +504,7 @@
2070 /*Add to results*/
2071 femmodel->results->AddObject(G_output);
2072 femmodel->results->AddObject(X_output);
2073-
2074+
2075 offset += N[i]*numberofvertices;
2076 }
2077 }
2078Index: ../trunk-jpl/src/c/cores/love_core.cpp
2079===================================================================
2080--- ../trunk-jpl/src/c/cores/love_core.cpp (revision 23065)
2081+++ ../trunk-jpl/src/c/cores/love_core.cpp (revision 23066)
2082@@ -27,12 +27,12 @@
2083
2084 /*parameters: */
2085 bool save_results;
2086-
2087+
2088 if(VerboseSolution()) _printf0_(" computing LOVE numbers\n");
2089
2090 /*Recover some parameters: */
2091 femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
2092-
2093+
2094 /*recover love number parameters: */
2095 femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
2096 femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); _assert_(nfreq==dummy);
2097@@ -44,7 +44,7 @@
2098 femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
2099 femmodel->parameters->FindParam(&love_kernels,LoveKernelsEnum);
2100 femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
2101-
2102+
2103 /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
2104 Matlitho* matlitho=NULL;
2105 for (int i=femmodel->materials->Size()-1;i>=0;i--){
2106Index: ../trunk-jpl/src/c/cores/damage_core.cpp
2107===================================================================
2108--- ../trunk-jpl/src/c/cores/damage_core.cpp (revision 23065)
2109+++ ../trunk-jpl/src/c/cores/damage_core.cpp (revision 23066)
2110@@ -18,7 +18,6 @@
2111 int numoutputs = 0;
2112 char **requested_outputs = NULL;
2113
2114-
2115 //first recover parameters common to all solutions
2116 femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
2117 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
2118Index: ../trunk-jpl/src/c/cores/smb_core.cpp
2119===================================================================
2120--- ../trunk-jpl/src/c/cores/smb_core.cpp (revision 23065)
2121+++ ../trunk-jpl/src/c/cores/smb_core.cpp (revision 23066)
2122@@ -28,9 +28,9 @@
2123 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
2124 femmodel->parameters->FindParam(&numoutputs,SmbNumRequestedOutputsEnum);
2125 if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SmbRequestedOutputsEnum);
2126-
2127+
2128 if(VerboseSolution()) _printf0_(" computing smb \n");
2129-
2130+
2131 analysis = new SmbAnalysis();
2132 analysis->Core(femmodel);
2133 delete analysis;
2134Index: ../trunk-jpl/src/c/cores/ad_core.cpp
2135===================================================================
2136--- ../trunk-jpl/src/c/cores/ad_core.cpp (revision 23065)
2137+++ ../trunk-jpl/src/c/cores/ad_core.cpp (revision 23066)
2138@@ -45,7 +45,7 @@
2139
2140 /*First, stop tracing: */
2141 trace_off();
2142-
2143+
2144 /*Print tape statistics so that user can kill this run if something is off already:*/
2145 if(VerboseAutodiff()){ /*{{{*/
2146 tapestats(my_rank,tape_stats); //reading of tape statistics
2147@@ -92,7 +92,7 @@
2148 /*retrieve parameters: */
2149 femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
2150 femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
2151-
2152+
2153 /*if no dependents, no point in running a driver: */
2154 if(!(num_dependents*num_independents)) return;
2155
2156@@ -100,7 +100,7 @@
2157 if (my_rank!=0){
2158 num_dependents=0; num_independents=0;
2159 }
2160-
2161+
2162 /*retrieve state variable: */
2163 femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
2164
2165@@ -116,7 +116,6 @@
2166 /*Branch according to AD driver: */
2167 femmodel->parameters->FindParam(&driver,AutodiffDriverEnum);
2168
2169-
2170 if (strcmp(driver,"fos_forward")==0){ /*{{{*/
2171
2172 int anIndepIndex;
2173@@ -140,7 +139,6 @@
2174 anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx;
2175 #endif
2176
2177-
2178 /*call driver: */
2179 fos_forward(my_rank,num_dependents,num_independents, 0, xp, tangentDir, theOutput, jacTimesTangentDir );
2180
2181@@ -184,7 +182,6 @@
2182 #endif
2183 // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
2184
2185-
2186 /*seed matrix: */
2187 seed=xNewZeroInit<double>(num_independents,tangentDirNum);
2188
2189@@ -243,7 +240,6 @@
2190 anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
2191 #endif
2192
2193-
2194 /*call driver: */
2195 fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
2196
2197@@ -284,7 +280,6 @@
2198 anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
2199 #endif
2200
2201-
2202 /*seed matrix: */
2203 weights=xNewZeroInit<double>(weightNum,num_dependents);
2204
2205@@ -316,7 +311,6 @@
2206 } /*}}}*/
2207 else _error_("driver: " << driver << " not yet supported!");
2208
2209-
2210 if(VerboseAutodiff())_printf0_(" end AD core\n");
2211
2212 /*Free resources: */
2213Index: ../trunk-jpl/src/c/classes/Materials/Matestar.cpp
2214===================================================================
2215--- ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 23065)
2216+++ ../trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 23066)
2217@@ -108,7 +108,7 @@
2218 void Matestar::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
2219
2220 if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook();
2221-
2222+
2223 MARSHALLING_ENUM(MatestarEnum);
2224 MARSHALLING(mid);
2225 this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
2226Index: ../trunk-jpl/src/c/classes/Materials/Matlitho.cpp
2227===================================================================
2228--- ../trunk-jpl/src/c/classes/Materials/Matlitho.cpp (revision 23065)
2229+++ ../trunk-jpl/src/c/classes/Materials/Matlitho.cpp (revision 23066)
2230@@ -36,13 +36,13 @@
2231
2232 this->radius=xNew<IssmDouble>(this->numlayers+1);
2233 xMemCpy<IssmDouble>(this->radius, iomodel->Data("md.materials.radius"),this->numlayers+1);
2234-
2235+
2236 this->viscosity=xNew<IssmDouble>(this->numlayers);
2237 xMemCpy<IssmDouble>(this->viscosity, iomodel->Data("md.materials.viscosity"),this->numlayers);
2238-
2239+
2240 this->lame_lambda=xNew<IssmDouble>(this->numlayers);
2241 xMemCpy<IssmDouble>(this->lame_lambda, iomodel->Data("md.materials.lame_lambda"),this->numlayers);
2242-
2243+
2244 this->lame_mu=xNew<IssmDouble>(this->numlayers);
2245 xMemCpy<IssmDouble>(this->lame_mu, iomodel->Data("md.materials.lame_mu"),this->numlayers);
2246
2247@@ -60,17 +60,17 @@
2248
2249 this->issolid=xNew<IssmDouble>(this->numlayers);
2250 xMemCpy<IssmDouble>(this->issolid, iomodel->Data("md.materials.issolid"),this->numlayers);
2251-
2252+
2253 /*isburgersd= xNew<IssmDouble>(this->numlayers);
2254 this->isburgers=xNew<bool>(this->numlayers);
2255 xMemCpy<IssmDouble>(isburgersd, iomodel->Data("md.materials.isburgers"),this->numlayers);
2256 for (int i=0;i<this->numlayers;i++)this->isburgers[i]=reCast<bool,IssmDouble>(isburgersd[i]);
2257-
2258+
2259 issolidd= xNew<IssmDouble>(this->numlayers);
2260 this->issolid=xNew<bool>(this->numlayers);
2261 xMemCpy<IssmDouble>(issolidd, iomodel->Data("md.materials.issolid"),this->numlayers);
2262 for (int i=0;i<this->numlayers;i++)this->issolid[i]=reCast<bool,IssmDouble>(issolidd[i]);*/
2263-
2264+
2265 /*free ressources: */
2266 xDelete<IssmDouble>(isburgersd);
2267 xDelete<IssmDouble>(issolidd);
2268@@ -77,7 +77,7 @@
2269 }
2270 /*}}}*/
2271 Matlitho::~Matlitho(){/*{{{*/
2272-
2273+
2274 xDelete<IssmDouble>(radius);
2275 xDelete<IssmDouble>(viscosity);
2276 xDelete<IssmDouble>(lame_lambda);
2277Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp
2278===================================================================
2279--- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 23065)
2280+++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 23066)
2281@@ -134,7 +134,7 @@
2282 _printf_(" note: helement not printed to avoid recursion.\n");
2283 //if(helement) helement->DeepEcho();
2284 //else _printf_(" helement = NULL\n");
2285-
2286+
2287 _printf_(" element:\n");
2288 _printf_(" note: element not printed to avoid recursion.\n");
2289 //if(element) element->DeepEcho();
2290@@ -142,12 +142,12 @@
2291 }
2292 /*}}}*/
2293 void Matice::Echo(void){/*{{{*/
2294-
2295+
2296 _printf_("Matice:\n");
2297 _printf_(" mid: " << mid << "\n");
2298 _printf_(" isdamaged: " << isdamaged << "\n");
2299 _printf_(" isenhanced: " << isenhanced << "\n");
2300-
2301+
2302 /*helement and element Echo were commented to avoid recursion.*/
2303 /*Example: element->Echo calls matice->Echo which calls element->Echo etc*/
2304 _printf_(" helement:\n");
2305@@ -154,7 +154,7 @@
2306 _printf_(" note: helement not printed to avoid recursion.\n");
2307 //if(helement) helement->Echo();
2308 //else _printf_(" helement = NULL\n");
2309-
2310+
2311 _printf_(" element:\n");
2312 _printf_(" note: element not printed to avoid recursion.\n");
2313 //if(element) element->Echo();
2314@@ -166,7 +166,7 @@
2315 void Matice::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
2316
2317 if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook();
2318-
2319+
2320 MARSHALLING_ENUM(MaticeEnum);
2321 MARSHALLING(mid);
2322 MARSHALLING(isdamaged);
2323@@ -540,7 +540,6 @@
2324 /*input strain rate: */
2325 IssmDouble exx,eyy,exy,exz,eyz;
2326
2327-
2328 if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
2329 (epsilon[3]==0) && (epsilon[4]==0)){
2330 mu_prime=0.5*pow(10.,14);
2331Index: ../trunk-jpl/src/c/classes/Numberedcostfunction.cpp
2332===================================================================
2333--- ../trunk-jpl/src/c/classes/Numberedcostfunction.cpp (revision 23065)
2334+++ ../trunk-jpl/src/c/classes/Numberedcostfunction.cpp (revision 23066)
2335@@ -8,7 +8,6 @@
2336 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
2337 #endif
2338
2339-
2340 /*Headers:*/
2341 //#include "./Definition.h"
2342 //#include "../datastructures/datastructures.h"
2343@@ -29,7 +28,7 @@
2344 #include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
2345
2346 /*}}}*/
2347-
2348+
2349 /*Numberedcostfunction constructors, destructors :*/
2350 Numberedcostfunction::Numberedcostfunction(){/*{{{*/
2351
2352@@ -109,7 +108,7 @@
2353 }
2354 /*}}}*/
2355 IssmDouble Numberedcostfunction::Response(FemModel* femmodel){/*{{{*/
2356-
2357+
2358 _assert_(number_cost_functions>0 && number_cost_functions<1e3);
2359
2360 /*output:*/
2361@@ -156,4 +155,3 @@
2362 return value_sum;
2363 }
2364 /*}}}*/
2365-
2366Index: ../trunk-jpl/src/c/classes/Params/BoolParam.cpp
2367===================================================================
2368--- ../trunk-jpl/src/c/classes/Params/BoolParam.cpp (revision 23065)
2369+++ ../trunk-jpl/src/c/classes/Params/BoolParam.cpp (revision 23066)
2370@@ -62,4 +62,3 @@
2371
2372 }
2373 /*}}}*/
2374-
2375Index: ../trunk-jpl/src/c/classes/Params/StringParam.cpp
2376===================================================================
2377--- ../trunk-jpl/src/c/classes/Params/StringParam.cpp (revision 23065)
2378+++ ../trunk-jpl/src/c/classes/Params/StringParam.cpp (revision 23066)
2379@@ -54,7 +54,7 @@
2380 int size = 0;
2381
2382 if(marshall_direction==MARSHALLING_FORWARD || marshall_direction == MARSHALLING_SIZE)size=strlen(value)+1;
2383-
2384+
2385 MARSHALLING_ENUM(StringParamEnum);
2386 MARSHALLING(enum_type);
2387 MARSHALLING(size);
2388Index: ../trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp
2389===================================================================
2390--- ../trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp (revision 23065)
2391+++ ../trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp (revision 23066)
2392@@ -115,7 +115,7 @@
2393
2394 /*DoubleVecParam specific routines:*/
2395 void DoubleVecParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM){/*{{{*/
2396-
2397+
2398 /*Assign output pointers:*/
2399 if(pM) *pM=M;
2400 *pIssmDoublearray=values;
2401Index: ../trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp
2402===================================================================
2403--- ../trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp (revision 23065)
2404+++ ../trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp (revision 23066)
2405@@ -115,7 +115,7 @@
2406
2407 /*DoubleMatParam specific routines:*/
2408 void DoubleMatParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM,int* pN){/*{{{*/
2409-
2410+
2411 /*Assign output pointers:*/
2412 if(pM) *pM=M;
2413 if(pN) *pN=N;
2414Index: ../trunk-jpl/src/c/classes/Params/DataSetParam.cpp
2415===================================================================
2416--- ../trunk-jpl/src/c/classes/Params/DataSetParam.cpp (revision 23065)
2417+++ ../trunk-jpl/src/c/classes/Params/DataSetParam.cpp (revision 23066)
2418@@ -54,7 +54,7 @@
2419 void DataSetParam::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
2420
2421 if(marshall_direction==MARSHALLING_BACKWARD)value=new DataSet();
2422-
2423+
2424 MARSHALLING_ENUM(DataSetParamEnum);
2425 MARSHALLING(enum_type);
2426 value->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
2427Index: ../trunk-jpl/src/c/classes/Params/IntParam.cpp
2428===================================================================
2429--- ../trunk-jpl/src/c/classes/Params/IntParam.cpp (revision 23065)
2430+++ ../trunk-jpl/src/c/classes/Params/IntParam.cpp (revision 23066)
2431@@ -63,4 +63,3 @@
2432
2433 }
2434 /*}}}*/
2435-
2436Index: ../trunk-jpl/src/c/classes/Params/StringArrayParam.cpp
2437===================================================================
2438--- ../trunk-jpl/src/c/classes/Params/StringArrayParam.cpp (revision 23065)
2439+++ ../trunk-jpl/src/c/classes/Params/StringArrayParam.cpp (revision 23066)
2440@@ -82,7 +82,7 @@
2441 }
2442
2443 MARSHALLING_ENUM(StringArrayParamEnum);
2444-
2445+
2446 MARSHALLING(enum_type);
2447 MARSHALLING(numstrings);
2448
2449Index: ../trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp
2450===================================================================
2451--- ../trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp (revision 23065)
2452+++ ../trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp (revision 23066)
2453@@ -232,4 +232,3 @@
2454
2455 }
2456 /*}}}*/
2457-
2458Index: ../trunk-jpl/src/c/classes/AmrBamg.cpp
2459===================================================================
2460--- ../trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23065)
2461+++ ../trunk-jpl/src/c/classes/AmrBamg.cpp (revision 23066)
2462@@ -33,7 +33,7 @@
2463 this->deviatoricerror_threshold = -1;
2464 this->deviatoricerror_groupthreshold = -1;
2465 this->deviatoricerror_maximum = -1;
2466-
2467+
2468 /*Geometry and mesh as NULL*/
2469 this->geometry = NULL;
2470 this->fathermesh = NULL;
2471@@ -136,7 +136,7 @@
2472 this->options->hmaxVertices = NULL;
2473 this->options->hmaxVerticesSize[0] = 0;
2474 this->options->hmaxVerticesSize[1] = 0;
2475-
2476+
2477 /*verify if the metric will be reseted or not*/
2478 if(this->keepmetric==0){
2479 if(this->options->metric) xDelete<IssmDouble>(this->options->metric);
2480@@ -154,15 +154,15 @@
2481 int *datalist = xNew<int>(2);
2482 IssmDouble *xylist= xNew<IssmDouble>(nbv*2);
2483 int* elementslist = xNew<int>(nbt*3);
2484-
2485+
2486 datalist[0] = nbv;
2487 datalist[1] = nbt;
2488-
2489+
2490 for(int i=0;i<nbv;i++){
2491 xylist[2*i] = meshout->Vertices[i*3+0];
2492 xylist[2*i+1] = meshout->Vertices[i*3+1];
2493 }
2494-
2495+
2496 for(int i=0;i<nbt;i++){
2497 elementslist[3*i+0] = reCast<int>(meshout->Triangles[4*i+0]);
2498 elementslist[3*i+1] = reCast<int>(meshout->Triangles[4*i+1]);
2499@@ -173,7 +173,7 @@
2500 *pdatalist = datalist;
2501 *pxylist = xylist;
2502 *pelementslist = elementslist;
2503-
2504+
2505 /*Cleanup and return*/
2506 delete geomout;
2507 }/*}}}*/
2508@@ -180,7 +180,7 @@
2509 void AmrBamg::SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in){/*{{{*/
2510
2511 if(!this->options) _error_("AmrBamg->options is NULL!");
2512-
2513+
2514 this->options->hmin = hmin_in;
2515 this->options->hmax = hmax_in;
2516 this->options->err[0] = err_in;
2517Index: ../trunk-jpl/src/c/classes/kriging/Observations.cpp
2518===================================================================
2519--- ../trunk-jpl/src/c/classes/kriging/Observations.cpp (revision 23065)
2520+++ ../trunk-jpl/src/c/classes/kriging/Observations.cpp (revision 23066)
2521@@ -740,4 +740,3 @@
2522 /*Assign output pointer*/
2523 xDelete<IssmPDouble>(counter);
2524 }/*}}}*/
2525-
2526Index: ../trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp
2527===================================================================
2528--- ../trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp (revision 23065)
2529+++ ../trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp (revision 23066)
2530@@ -22,7 +22,7 @@
2531 #include "../classes/Inputs/Input.h"
2532 #include "../classes/gauss/Gauss.h"
2533 /*}}}*/
2534-
2535+
2536 /*Cfsurfacelogvel constructors, destructors :*/
2537 Cfsurfacelogvel::Cfsurfacelogvel(){/*{{{*/
2538
2539@@ -39,13 +39,13 @@
2540 Cfsurfacelogvel::Cfsurfacelogvel(char* in_name, int in_definitionenum, IssmDouble in_datatime, bool in_timepassedflag){/*{{{*/
2541
2542 this->definitionenum=in_definitionenum;
2543-
2544+
2545 this->name = xNew<char>(strlen(in_name)+1);
2546 xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
2547
2548 this->datatime=in_datatime;
2549 this->timepassedflag=in_timepassedflag;
2550-
2551+
2552 this->misfit=0;
2553 this->lock=0;
2554 }
2555@@ -99,10 +99,10 @@
2556 }
2557 /*}}}*/
2558 IssmDouble Cfsurfacelogvel::Response(FemModel* femmodel){/*{{{*/
2559-
2560+
2561 /*diverse: */
2562 IssmDouble time;
2563-
2564+
2565 /*recover time parameters: */
2566 femmodel->parameters->FindParam(&time,TimeEnum);
2567 if(time < last_time) timepassedflag = false;
2568@@ -111,7 +111,7 @@
2569 int i;
2570 IssmDouble J=0.;
2571 IssmDouble J_sum=0.;
2572-
2573+
2574 if(datatime<=time && !timepassedflag){
2575 for(i=0;i<femmodel->elements->Size();i++){
2576 Element* element=(Element*)femmodel->elements->GetObjectByOffset(i);
2577@@ -121,7 +121,7 @@
2578 ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
2579 ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
2580 J=J_sum;
2581-
2582+
2583 timepassedflag = true;
2584 return J;
2585 }
2586@@ -138,7 +138,7 @@
2587 IssmDouble misfit,Jdet;
2588 IssmDouble vx,vy,vxobs,vyobs,weight;
2589 IssmDouble* xyz_list = NULL;
2590-
2591+
2592 /*Get basal element*/
2593 if(!element->IsOnSurface()) return 0.;
2594
2595@@ -159,7 +159,7 @@
2596
2597 /* Get node coordinates*/
2598 topelement->GetVerticesCoordinates(&xyz_list);
2599-
2600+
2601 /*Get model values*/
2602 Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input);
2603 Input* vy_input =NULL;
2604@@ -174,7 +174,6 @@
2605 if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do");
2606 datasetinput = (DatasetInput*)tempinput;
2607
2608-
2609 /* Start looping on the number of gaussian points: */
2610 Gauss* gauss=topelement->NewGauss(2);
2611 for(int ig=gauss->begin();ig<gauss->end();ig++){
2612@@ -221,4 +220,3 @@
2613 delete gauss;
2614 return Jelem;
2615 }/*}}}*/
2616-
2617Index: ../trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp
2618===================================================================
2619--- ../trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp (revision 23065)
2620+++ ../trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp (revision 23066)
2621@@ -22,7 +22,7 @@
2622 #include "../classes/Inputs/Input.h"
2623 #include "../classes/gauss/Gauss.h"
2624 /*}}}*/
2625-
2626+
2627 /*Cfdragcoeffabsgrad constructors, destructors :*/
2628 Cfdragcoeffabsgrad::Cfdragcoeffabsgrad(){/*{{{*/
2629
2630@@ -39,13 +39,13 @@
2631 Cfdragcoeffabsgrad::Cfdragcoeffabsgrad(char* in_name, int in_definitionenum, int in_weights_enum, bool in_timepassedflag){/*{{{*/
2632
2633 this->definitionenum=in_definitionenum;
2634-
2635+
2636 this->name = xNew<char>(strlen(in_name)+1);
2637 xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
2638
2639 this->weights_enum=in_weights_enum;
2640 this->timepassedflag=in_timepassedflag;
2641-
2642+
2643 this->misfit=0;
2644 this->lock=0;
2645 }
2646@@ -101,7 +101,7 @@
2647 IssmDouble Cfdragcoeffabsgrad::Response(FemModel* femmodel){/*{{{*/
2648 /*diverse: */
2649 IssmDouble time;
2650-
2651+
2652 /*recover time parameters: */
2653 int i;
2654 IssmDouble J=0.;
2655@@ -115,7 +115,7 @@
2656 ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
2657 ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
2658 J=J_sum;
2659-
2660+
2661 timepassedflag = true;
2662 return J;
2663 }
2664@@ -182,4 +182,3 @@
2665 delete gauss;
2666 return Jelem;
2667 }/*}}}*/
2668-
2669Index: ../trunk-jpl/src/c/classes/Loads/Neumannflux.cpp
2670===================================================================
2671--- ../trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 23065)
2672+++ ../trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 23066)
2673@@ -30,7 +30,6 @@
2674 /*}}}*/
2675 Neumannflux::Neumannflux(int neumannflux_id,int i,IoModel* iomodel,int* segments,int in_analysis_type){/*{{{*/
2676
2677-
2678 /*Some sanity checks*/
2679 _assert_(segments);
2680
2681Index: ../trunk-jpl/src/c/classes/Loads/Moulin.cpp
2682===================================================================
2683--- ../trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23065)
2684+++ ../trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 23066)
2685@@ -177,7 +177,7 @@
2686 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
2687
2688 switch(analysis_type){
2689-
2690+
2691 case HydrologyShaktiAnalysisEnum:
2692 pe = this->CreatePVectorHydrologyShakti();
2693 break;
2694@@ -394,10 +394,10 @@
2695 if(!this->node->IsActive()) return NULL;
2696 IssmDouble moulin_load,dt;
2697 ElementVector* pe=new ElementVector(&node,1,this->parameters);
2698-
2699+
2700 this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
2701 parameters->FindParam(&dt,TimesteppingTimeStepEnum);
2702-
2703+
2704 pe->values[0]=moulin_load*dt;
2705 /*Clean up and return*/
2706 return pe;
2707Index: ../trunk-jpl/src/c/classes/Cfsurfacesquare.cpp
2708===================================================================
2709--- ../trunk-jpl/src/c/classes/Cfsurfacesquare.cpp (revision 23065)
2710+++ ../trunk-jpl/src/c/classes/Cfsurfacesquare.cpp (revision 23066)
2711@@ -22,7 +22,7 @@
2712 #include "../classes/Inputs/Input.h"
2713 #include "../classes/gauss/Gauss.h"
2714 /*}}}*/
2715-
2716+
2717 /*Cfsurfacesquare constructors, destructors :*/
2718 Cfsurfacesquare::Cfsurfacesquare(){/*{{{*/
2719
2720@@ -42,7 +42,7 @@
2721 Cfsurfacesquare::Cfsurfacesquare(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, int in_weights_enum, IssmDouble in_datatime, bool in_timepassedflag){/*{{{*/
2722
2723 this->definitionenum=in_definitionenum;
2724-
2725+
2726 this->name = xNew<char>(strlen(in_name)+1);
2727 xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
2728
2729@@ -51,7 +51,7 @@
2730 this->weights_enum=in_weights_enum;
2731 this->datatime=in_datatime;
2732 this->timepassedflag=in_timepassedflag;
2733-
2734+
2735 this->misfit=0;
2736 this->lock=0;
2737 }
2738@@ -110,7 +110,7 @@
2739 IssmDouble Cfsurfacesquare::Response(FemModel* femmodel){/*{{{*/
2740 /*diverse: */
2741 IssmDouble time;
2742-
2743+
2744 /*recover time parameters: */
2745 femmodel->parameters->FindParam(&time,TimeEnum);
2746 if(time < last_time) timepassedflag = false;
2747@@ -129,7 +129,7 @@
2748 ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
2749 ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
2750 J=J_sum;
2751-
2752+
2753 timepassedflag = true;
2754 return J;
2755 }
2756@@ -171,7 +171,7 @@
2757
2758 /*Get input if it already exists*/
2759 Input* tempinput = topelement->GetInput(definitionenum);
2760-
2761+
2762 /*Cast it to a Datasetinput*/
2763 if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do");
2764 datasetinput = (DatasetInput*)tempinput;
2765@@ -213,4 +213,3 @@
2766 delete gauss;
2767 return Jelem;
2768 }/*}}}*/
2769-
2770Index: ../trunk-jpl/src/c/classes/Profiler.cpp
2771===================================================================
2772--- ../trunk-jpl/src/c/classes/Profiler.cpp (revision 23065)
2773+++ ../trunk-jpl/src/c/classes/Profiler.cpp (revision 23066)
2774@@ -99,7 +99,6 @@
2775 _assert_(tag<MAXIMUMSIZE);
2776 if(this->running[tag]) _error_("Tag "<<tag<<" has not been stopped");
2777
2778-
2779 #ifdef _HAVE_MPI_
2780 return this->time[tag];
2781 #else
2782@@ -143,7 +142,6 @@
2783 _assert_(tag<MAXIMUMSIZE);
2784 if(this->running[tag]) _error_("Tag "<<tag<<" is already running");
2785
2786-
2787 /*If mpisync requested, make sure all the cpus are at the same point in the execution: */
2788 if(!dontmpisync){
2789 ISSM_MPI_Barrier(IssmComm::GetComm());
2790@@ -183,7 +181,6 @@
2791 _assert_(tag<MAXIMUMSIZE);
2792 if(!this->running[tag]) _error_("Tag "<<tag<<" is not running");
2793
2794-
2795 /*If mpisync requested, make sure all the cpus are at the same point in the execution: */
2796 if(!dontmpisync){
2797 ISSM_MPI_Barrier(IssmComm::GetComm());
2798Index: ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp
2799===================================================================
2800--- ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23065)
2801+++ ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp (revision 23066)
2802@@ -100,7 +100,7 @@
2803
2804 /*Mesh refinement methods*/
2805 void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxylist,int** pelementslist){/*{{{*/
2806-
2807+
2808 /*IMPORTANT! pelementslist are in Matlab indexing*/
2809 /*NEOPZ works only in C indexing*/
2810 if(!this->fathermesh || !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh is NULL!\n");
2811@@ -117,10 +117,10 @@
2812
2813 /*Intermediaries*/
2814 bool verbose=VerboseSolution();
2815-
2816+
2817 /*Execute refinement*/
2818 this->RefineMeshOneLevel(verbose,gl_distance,if_distance,deviatoricerror,thicknesserror);
2819-
2820+
2821 /*Get new geometric mesh in ISSM data structure*/
2822 this->GetMesh(this->previousmesh,pdatalist,pxylist,pelementslist);
2823
2824@@ -129,7 +129,7 @@
2825 }
2826 /*}}}*/
2827 void AdaptiveMeshRefinement::RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror){/*{{{*/
2828-
2829+
2830 /*Intermediaries*/
2831 int nelem =-1;
2832 int side2D = 6;
2833@@ -203,7 +203,7 @@
2834 if(criteria) index.push_back(gmesh->Element(this->specialelementsindex[i])->FatherIndex());
2835 }
2836 /*}}}*/
2837-
2838+
2839 /*Now, detele the special elements{{{*/
2840 if(this->refinement_type==1) this->DeleteSpecialElements(verbose,gmesh);
2841 else this->specialelementsindex.clear();
2842@@ -215,12 +215,12 @@
2843 gmesh=this->fathermesh;
2844 }
2845 /*}}}*/
2846-
2847+
2848 /*Unrefinement process: loop over level of refinements{{{*/
2849 if(verbose) _printf_("\tunrefinement process...\n");
2850 if(verbose) _printf_("\ttotal: ");
2851 count=0;
2852-
2853+
2854 nelem=gmesh->NElements();//must keep here
2855 for(int i=0;i<nelem;i++){
2856 if(!gmesh->Element(i)) continue;
2857@@ -265,7 +265,7 @@
2858 /*Adjust the connectivities before continue*/
2859 gmesh->BuildConnectivity();
2860 /*}}}*/
2861-
2862+
2863 /*Refinement process: loop over level of refinements{{{*/
2864 if(verbose) _printf_("\trefinement process (level max = "<<this->level_max<<")\n");
2865 if(verbose) _printf_("\ttotal: ");
2866@@ -319,14 +319,14 @@
2867 * 2 : uniform refinment
2868 * */
2869 if(!geoel) _error_("geoel is NULL!\n");
2870-
2871+
2872 /*Output*/
2873 int type=0;
2874 int count=0;
2875-
2876+
2877 /*Intermediaries*/
2878 TPZVec<TPZGeoEl*> sons;
2879-
2880+
2881 /*Loop over neighboors (sides 3, 4 and 5)*/
2882 for(int j=3;j<6;j++){
2883 sons.clear();
2884@@ -334,11 +334,11 @@
2885 if(sons.size()) count++; //if neighbour was refined
2886 if(sons.size()>4) count++; //if neighbour's level is > element level+1
2887 }
2888-
2889+
2890 /*Verify and return*/
2891 if(count>1) type=2;
2892 else type=count;
2893-
2894+
2895 return type;
2896 }
2897 /*}}}*/
2898@@ -390,14 +390,14 @@
2899 }
2900 /*}}}*/
2901 void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/
2902-
2903+
2904 /*Now, insert special elements to avoid hanging nodes*/
2905 if(verbose) _printf_("\trefining to avoid hanging nodes (total: ");
2906-
2907+
2908 /*Intermediaries*/
2909 int nelem=-1;
2910 int count= 1;
2911-
2912+
2913 while(count>0){
2914 nelem=gmesh->NElements();//must keep here
2915 count=0;
2916@@ -467,13 +467,13 @@
2917 long* vertex_index2sid = xNew<long>(ntotalvertices);
2918 this->index2sid.clear(); this->index2sid.resize(gmesh->NElements());
2919 this->sid2index.clear();
2920-
2921+
2922 /*Fill in the vertex_index2sid vector with non usual index value*/
2923 for(int i=0;i<gmesh->NNodes();i++) vertex_index2sid[i]=-1;
2924-
2925+
2926 /*Fill in the this->index2sid vector with non usual index value*/
2927 for(int i=0;i<gmesh->NElements();i++) this->index2sid[i]=-1;
2928-
2929+
2930 /*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */
2931 sid=0;
2932 for(int i=0;i<gmesh->NElements();i++){//over gmesh elements index
2933@@ -510,7 +510,7 @@
2934 newmeshXY[2*sid] = coords[0]; // X
2935 newmeshXY[2*sid+1] = coords[1]; // Y
2936 }
2937-
2938+
2939 for(int i=0;i<this->sid2index.size();i++){//over the sid (fill the ISSM elements)
2940 for(int j=0;j<this->GetNumberOfNodes();j++) {
2941 geoel = gmesh->ElementVec()[this->sid2index[i]];
2942@@ -532,7 +532,7 @@
2943 xc=newmeshXY[2*c]; yc=newmeshXY[2*c+1];
2944
2945 detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
2946-
2947+
2948 /*verify and swap, if necessary*/
2949 if(detJ<0) {
2950 newelements[i*this->GetNumberOfNodes()+0]=b+1;//a->b
2951@@ -544,7 +544,7 @@
2952 *pdata = newdata;
2953 *pxy = newmeshXY;
2954 *pelements = newelements;
2955-
2956+
2957 /*Cleanup*/
2958 xDelete<long>(vertex_index2sid);
2959 }
2960@@ -553,7 +553,7 @@
2961
2962 /* IMPORTANT! elements come in Matlab indexing
2963 NEOPZ works only in C indexing*/
2964-
2965+
2966 if(nvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n");
2967 if(nelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n");
2968 if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to create initial mesh: refinement type is not defined!\n");
2969@@ -560,7 +560,7 @@
2970
2971 /*Verify and creating initial mesh*/
2972 if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!");
2973-
2974+
2975 this->fathermesh = new TPZGeoMesh();
2976 this->fathermesh->NodeVec().Resize(nvertices);
2977
2978@@ -575,7 +575,7 @@
2979 this->fathermesh->NodeVec()[i].SetCoord(coord);
2980 this->fathermesh->NodeVec()[i].SetNodeId(i);
2981 }
2982-
2983+
2984 /*Generate the elements*/
2985 long index;
2986 const int mat = this->GetElemMaterialID();
2987@@ -603,10 +603,10 @@
2988 }
2989 /*}}}*/
2990 TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/
2991-
2992+
2993 TPZGeoMesh *newgmesh = new TPZGeoMesh();
2994 newgmesh->CleanUp();
2995-
2996+
2997 int nnodes = gmesh->NNodes();
2998 int nelem = gmesh->NElements();
2999 int mat = this->GetElemMaterialID();;
3000@@ -616,25 +616,25 @@
3001 //nodes
3002 newgmesh->NodeVec().Resize(nnodes);
3003 for(int i=0;i<nnodes;i++) newgmesh->NodeVec()[i] = gmesh->NodeVec()[i];
3004-
3005+
3006 //elements
3007 for(int i=0;i<nelem;i++){
3008 TPZGeoEl * geoel = gmesh->Element(i);
3009-
3010+
3011 if(!geoel){
3012 index=newgmesh->ElementVec().AllocateNewElement();
3013 newgmesh->ElementVec()[index] = NULL;
3014 continue;
3015 }
3016-
3017+
3018 TPZManVector<long> elem(3,0);
3019 for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j);
3020-
3021+
3022 newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);
3023 newgmesh->ElementVec()[index]->SetId(geoel->Id());
3024-
3025+
3026 TPZGeoElRefPattern<TPZGeoTriangle>* newgeoel = dynamic_cast<TPZGeoElRefPattern<TPZGeoTriangle>*>(newgmesh->ElementVec()[index]);
3027-
3028+
3029 //old neighbourhood
3030 const int nsides = TPZGeoTriangle::NSides;
3031 TPZVec< std::vector<TPZGeoElSide> > neighbourhood(nsides);
3032@@ -653,7 +653,7 @@
3033 neighS = neighS.Neighbour();
3034 }
3035 }
3036-
3037+
3038 //inserting in new element
3039 for(int s = 0; s < nsides; s++){
3040 TPZGeoEl * tempEl = newgeoel;
3041@@ -667,17 +667,17 @@
3042 }
3043 tempEl->SetNeighbour(byside, tempSide);
3044 }
3045-
3046+
3047 long fatherindex = geoel->FatherIndex();
3048 if(fatherindex>-1) newgeoel->SetFather(fatherindex);
3049-
3050+
3051 if(!geoel->HasSubElement()) continue;
3052-
3053+
3054 int nsons = geoel->NSubElements();
3055
3056 TPZAutoPointer<TPZRefPattern> ref = gRefDBase.GetUniformRefPattern(ETriangle);
3057 newgeoel->SetRefPattern(ref);
3058-
3059+
3060 for(int j=0;j<nsons;j++){
3061 TPZGeoEl* son = geoel->SubElement(j);
3062 if(!son){
3063@@ -689,7 +689,7 @@
3064
3065 /*Now, build connectivities*/
3066 newgmesh->BuildConnectivity();
3067-
3068+
3069 return newgmesh;
3070 }
3071 /*}}}*/
3072@@ -704,7 +704,7 @@
3073 }
3074 /*}}}*/
3075 void AdaptiveMeshRefinement::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/*{{{*/
3076-
3077+
3078 file.clear();
3079 long nelements = gmesh->NElements();
3080 TPZGeoEl *gel;
3081@@ -730,10 +730,10 @@
3082 }
3083 MElementType elt = gel->Type();
3084 int elNnodes = MElementType_NNodes(elt);
3085-
3086+
3087 size += (1+elNnodes);
3088 connectivity << elNnodes;
3089-
3090+
3091 for(int t = 0; t < elNnodes; t++)
3092 {
3093 for(int c = 0; c < 3; c++)
3094@@ -742,15 +742,15 @@
3095 node << coord << " ";
3096 }
3097 node << std::endl;
3098-
3099+
3100 actualNode++;
3101 connectivity << " " << actualNode;
3102 }
3103 connectivity << std::endl;
3104-
3105+
3106 int elType = this->GetVTK_ElType(gel);
3107 type << elType << std::endl;
3108-
3109+
3110 if(matColor == true)
3111 {
3112 material << gel->MaterialId() << std::endl;
3113@@ -759,7 +759,7 @@
3114 {
3115 material << gel->Index() << std::endl;
3116 }
3117-
3118+
3119 nVALIDelements++;
3120 }
3121 node << std::endl;
3122@@ -789,9 +789,9 @@
3123 }
3124 /*}}}*/
3125 int AdaptiveMeshRefinement::GetVTK_ElType(TPZGeoEl * gel){/*{{{*/
3126-
3127+
3128 MElementType pzElType = gel->Type();
3129-
3130+
3131 int elType = -1;
3132 switch (pzElType)
3133 {
3134@@ -848,8 +848,7 @@
3135 std::cout << "MIGHT BE CURVED ELEMENT (quadratic or quarter point)" << std::endl;
3136 DebugStop();
3137 }
3138-
3139+
3140 return elType;
3141 }
3142 /*}}}*/
3143-
3144Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
3145===================================================================
3146--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23065)
3147+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23066)
3148@@ -105,7 +105,7 @@
3149 for(i=0;i<num_nodes;i++) if(this->nodes[i]) tria->nodes[i]=this->nodes[i]; else tria->nodes[i] = NULL;
3150 }
3151 else tria->nodes = NULL;
3152-
3153+
3154 tria->vertices = (Vertex**)this->hvertices->deliverp();
3155 tria->material = (Material*)this->hmaterial->delivers();
3156 tria->matpar = (Matpar*)this->hmatpar->delivers();
3157@@ -114,9 +114,9 @@
3158 }
3159 /*}}}*/
3160 void Tria::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
3161-
3162+
3163 MARSHALLING_ENUM(TriaEnum);
3164-
3165+
3166 /*Call parent classes: */
3167 ElementHook::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
3168 Element::MarshallElement(pmarshalled_data,pmarshalled_data_size,marshall_direction,this->numanalyses);
3169@@ -134,7 +134,7 @@
3170
3171 /*Call inputs method*/
3172 _assert_(this->inputs);
3173-
3174+
3175 int domaintype;
3176 parameters->FindParam(&domaintype,DomainTypeEnum);
3177 switch(domaintype){
3178@@ -184,7 +184,7 @@
3179 if(N!=num_inputs) _error_("sizes are not consistent");
3180
3181 int tria_vertex_ids[3];
3182-
3183+
3184 for(int k=0;k<3;k++){
3185 tria_vertex_ids[k]=reCast<int>(iomodel->elements[3*this->Sid()+k]); //ids for vertices are in the elements array from Matlab
3186 }
3187@@ -327,7 +327,7 @@
3188 }
3189 /*}}}*/
3190 void Tria::CalvingCrevasseDepth(){/*{{{*/
3191-
3192+
3193 IssmDouble xyz_list[NUMVERTICES][3];
3194 IssmDouble calvingrate[NUMVERTICES];
3195 IssmDouble vx,vy,vel;
3196@@ -339,10 +339,10 @@
3197
3198 /* Get node coordinates and dof list: */
3199 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
3200-
3201+
3202 /*Get the critical fraction of thickness surface and basal crevasses penetrate for calving onset*/
3203 this->parameters->FindParam(&critical_fraction,CalvingCrevasseDepthEnum);
3204-
3205+
3206 IssmDouble rho_ice = this->GetMaterialParameter(MaterialsRhoIceEnum);
3207 IssmDouble rho_seawater = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
3208 IssmDouble rho_freshwater = this->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
3209@@ -360,12 +360,12 @@
3210 Input* s_xx_input = inputs->GetInput(DeviatoricStressxxEnum); _assert_(s_xx_input);
3211 Input* s_xy_input = inputs->GetInput(DeviatoricStressxyEnum); _assert_(s_xy_input);
3212 Input* s_yy_input = inputs->GetInput(DeviatoricStressyyEnum); _assert_(s_yy_input);
3213-
3214+
3215 /*Loop over all elements of this partition*/
3216 GaussTria* gauss=new GaussTria();
3217 for (int iv=0;iv<NUMVERTICES;iv++){
3218 gauss->GaussVertex(iv);
3219-
3220+
3221 H_input->GetInputValue(&thickness,gauss);
3222 bed_input->GetInputValue(&bed,gauss);
3223 surface_input->GetInputValue(&float_depth,gauss);
3224@@ -377,13 +377,13 @@
3225 s_xx_input->GetInputValue(&s_xx,gauss);
3226 s_xy_input->GetInputValue(&s_xy,gauss);
3227 s_yy_input->GetInputValue(&s_yy,gauss);
3228-
3229+
3230 vel=sqrt(vx*vx+vy*vy)+1.e-14;
3231
3232 s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
3233 s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
3234 if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
3235-
3236+
3237 Ho = thickness - (rho_seawater/rho_ice) * (-bed);
3238 if(Ho<0.) Ho=0.;
3239
3240@@ -398,19 +398,19 @@
3241 //if (surface_crevasse[iv]<water_height){
3242 // water_height = surface_crevasse[iv];
3243 //}
3244-
3245+
3246 /*basal crevasse*/
3247 //basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (rheology_B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho);
3248 basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
3249 if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
3250 if (bed>0.) basal_crevasse[iv] = 0.;
3251-
3252+
3253 H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth;
3254 H_surfbasal = (surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv])-(critical_fraction*thickness);
3255-
3256+
3257 crevasse_depth[iv] = max(H_surf,H_surfbasal);
3258 }
3259-
3260+
3261 this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum));
3262 this->inputs->AddInput(new TriaInput(BasalCrevasseEnum,&basal_crevasse[0],P1Enum));
3263 this->inputs->AddInput(new TriaInput(CrevasseDepthEnum,&crevasse_depth[0],P1Enum));
3264@@ -460,7 +460,7 @@
3265 }
3266 else
3267 calvingrate[iv]=0.;
3268-
3269+
3270 calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);
3271 calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14);
3272 }
3273@@ -561,7 +561,7 @@
3274 IssmDouble strain_xy[NUMVERTICES];
3275 IssmDouble vorticity_xy[NUMVERTICES];
3276 GaussTria* gauss=NULL;
3277-
3278+
3279 /* Get node coordinates and dof list: */
3280 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
3281
3282@@ -568,7 +568,7 @@
3283 /*Retrieve all inputs we will be needing: */
3284 Input* vx_input=this->GetInput(EsaXmotionEnum); _assert_(vx_input);
3285 Input* vy_input=this->GetInput(EsaYmotionEnum); _assert_(vy_input);
3286-
3287+
3288 /* Start looping on the number of vertices: */
3289 gauss=new GaussTria();
3290 for (int iv=0;iv<NUMVERTICES;iv++){
3291@@ -772,7 +772,6 @@
3292 }
3293 }
3294
3295-
3296 }/*}}}*/
3297 void Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
3298
3299@@ -1399,7 +1398,7 @@
3300 }
3301 /*}}}*/
3302 void Tria::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
3303-
3304+
3305 /* Intermediaries */
3306 int i, dir,nrfrontnodes;
3307 IssmDouble levelset[NUMVERTICES];
3308@@ -1488,7 +1487,7 @@
3309
3310 }/*}}}*/
3311 void Tria::GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/
3312-
3313+
3314 /* GetLevelsetIntersection computes:
3315 * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion,
3316 * 2. fraction of intersected triangle edges intersected by levelset, lying below level*/
3317@@ -1564,7 +1563,7 @@
3318 }
3319 /*}}}*/
3320 void Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
3321-
3322+
3323 /*Computeportion of the element that has a positive levelset*/
3324
3325 bool negative=true;
3326@@ -1678,11 +1677,11 @@
3327 Input* input=(Input*)this->inputs->GetInput(control_enum); _assert_(input);
3328
3329 parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
3330-
3331+
3332 /*Cast to Controlinput*/
3333 if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
3334 ControlInput* controlinput = xDynamicCast<ControlInput*>(input);
3335-
3336+
3337 if(strcmp(data,"value")==0){
3338 input = controlinput->values;
3339 }
3340@@ -1699,7 +1698,7 @@
3341 _error_("Data " << data << " not supported yet");
3342 }
3343 /*Check what input we are dealing with*/
3344-
3345+
3346 switch(input->ObjectEnum()){
3347 case TriaInputEnum:
3348 {
3349@@ -1716,7 +1715,7 @@
3350 vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
3351 break;
3352 }
3353-
3354+
3355 case TransientInputEnum:
3356 {
3357 TransientInput* transientinput = xDynamicCast<TransientInput*>(input);
3358@@ -1981,7 +1980,7 @@
3359 surface_input->GetInputAverage(&surface);
3360 base_input->GetInputAverage(&bed);
3361 bed_input->GetInputAverage(&bathymetry);
3362-
3363+
3364 /*Return: */
3365 return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.));
3366 }
3367@@ -2026,8 +2025,6 @@
3368 if(control_analysis && ad_analysis) iomodel->FindConstant(&num_control_type,"md.autodiff.num_independent_objects");
3369 if(control_analysis && ad_analysis) iomodel->FindConstant(&num_responses,"md.autodiff.num_dependent_objects");
3370
3371-
3372-
3373 /*Recover vertices ids needed to initialize inputs*/
3374 for(i=0;i<3;i++){
3375 tria_vertex_ids[i]=reCast<int>(iomodel->elements[3*index+i]); //ids for vertices are in the elements array from Matlab
3376@@ -2393,7 +2390,6 @@
3377 /*}}}*/
3378 IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/
3379
3380-
3381 /*intermediary: */
3382 IssmDouble* values=NULL;
3383 Input* thickness_input=NULL;
3384@@ -2406,7 +2402,7 @@
3385 int point1;
3386 IssmDouble fraction1,fraction2;
3387 bool mainlynegative=true;
3388-
3389+
3390 /*Output:*/
3391 volume=0;
3392
3393@@ -2415,7 +2411,7 @@
3394
3395 /*Retrieve inputs required:*/
3396 thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
3397-
3398+
3399 /*Retrieve material parameters: */
3400 rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
3401
3402@@ -2424,7 +2420,7 @@
3403 for(int i=0;i<NUMVERTICES;i++){
3404 values[i]=levelset[this->vertices[i]->Sid()];
3405 }
3406-
3407+
3408 /*Ok, use the level set values to figure out where we put our gaussian points:*/
3409 this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
3410 Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
3411@@ -2876,7 +2872,7 @@
3412 dist_gl_input->GetInputValue(&dist_gl,gauss);
3413 dist_cf_input->GetInputValue(&dist_cf,gauss);
3414 delete gauss;
3415-
3416+
3417 /*Ensure values are positive for floating ice*/
3418 dist_gl = fabs(dist_gl);
3419 dist_cf = fabs(dist_cf);
3420@@ -3196,7 +3192,6 @@
3421 }
3422 }
3423
3424-
3425 }
3426 /*}}}*/
3427 void Tria::ResetFSBasalBoundaryCondition(void){/*{{{*/
3428@@ -3277,7 +3272,6 @@
3429 IssmDouble values[NUMVERTICES];
3430 int idlist[NUMVERTICES],control_init;
3431
3432-
3433 /*Get Domain type*/
3434 int domaintype;
3435 parameters->FindParam(&domaintype,DomainTypeEnum);
3436@@ -3297,7 +3291,7 @@
3437
3438 /*Get out if this is not an element input*/
3439 if(!IsInputEnum(control_enum)) return;
3440-
3441+
3442 Input* input = (Input*)this->inputs->GetInput(control_enum); _assert_(input);
3443 if(input->ObjectEnum()!=ControlInputEnum){
3444 _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
3445@@ -3330,7 +3324,6 @@
3446 IssmDouble values[NUMVERTICES];
3447 int vertexpidlist[NUMVERTICES],control_init;
3448
3449-
3450 /*Get Domain type*/
3451 int domaintype;
3452 parameters->FindParam(&domaintype,DomainTypeEnum);
3453@@ -4166,7 +4159,7 @@
3454 IssmDouble* U_elastic_precomputed = NULL;
3455 IssmDouble* H_elastic_precomputed = NULL;
3456 int M, hemi;
3457-
3458+
3459 /*computation of Green functions:*/
3460 IssmDouble* U_elastic= NULL;
3461 IssmDouble* N_elastic= NULL;
3462@@ -4173,7 +4166,7 @@
3463 IssmDouble* E_elastic= NULL;
3464 IssmDouble* X_elastic= NULL;
3465 IssmDouble* Y_elastic= NULL;
3466-
3467+
3468 /*optimization:*/
3469 bool store_green_functions=false;
3470
3471@@ -4181,7 +4174,7 @@
3472 Input* deltathickness_input=inputs->GetInput(EsaDeltathicknessEnum);
3473 if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
3474 deltathickness_input->GetInputAverage(&I);
3475-
3476+
3477 /*early return if we are not on the (ice) loading point: */
3478 if(I==0) return;
3479
3480@@ -4194,7 +4187,7 @@
3481
3482 /*which hemisphere? for north-south, east-west components*/
3483 this->parameters->FindParam(&hemi,EsaHemisphereEnum);
3484-
3485+
3486 /*compute area of element:*/
3487 area=GetArea();
3488
3489@@ -4253,7 +4246,7 @@
3490 U_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*U_elastic[i];
3491 Y_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*Y_elastic[i];
3492 X_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*X_elastic[i];
3493-
3494+
3495 /*North-south, East-west components */
3496 if (hemi == -1) {
3497 ang2 = PI/2 - atan2(yy[i],xx[i]);
3498@@ -4276,7 +4269,7 @@
3499 pEast->SetValues(gsize,indices,E_values,ADD_VAL);
3500 pX->SetValues(gsize,indices,X_values,ADD_VAL);
3501 pY->SetValues(gsize,indices,Y_values,ADD_VAL);
3502-
3503+
3504 /*free ressources:*/
3505 xDelete<int>(indices);
3506 xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
3507@@ -4306,12 +4299,12 @@
3508 IssmDouble* U_elastic_precomputed = NULL;
3509 IssmDouble* H_elastic_precomputed = NULL;
3510 int M;
3511-
3512+
3513 /*computation of Green functions:*/
3514 IssmDouble* U_elastic= NULL;
3515 IssmDouble* N_elastic= NULL;
3516 IssmDouble* E_elastic= NULL;
3517-
3518+
3519 /*optimization:*/
3520 bool store_green_functions=false;
3521
3522@@ -4319,7 +4312,7 @@
3523 Input* deltathickness_input=inputs->GetInput(EsaDeltathicknessEnum);
3524 if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
3525 deltathickness_input->GetInputAverage(&I);
3526-
3527+
3528 /*early return if we are not on the (ice) loading point: */
3529 if(I==0) return;
3530
3531@@ -4418,7 +4411,7 @@
3532 dx = x_element-x; dy = y_element-y; dz = z_element-z;
3533 N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
3534 E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
3535-
3536+
3537 /*Elastic component (from Eq 17 in Adhikari et al, GMD 2015): */
3538 int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
3539 U_elastic[i] += U_elastic_precomputed[index];
3540@@ -4433,7 +4426,7 @@
3541 pUp->SetValues(gsize,indices,U_values,ADD_VAL);
3542 pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
3543 pEast->SetValues(gsize,indices,E_values,ADD_VAL);
3544-
3545+
3546 /*free ressources:*/
3547 xDelete<int>(indices);
3548 xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
3549@@ -4454,7 +4447,7 @@
3550 IssmDouble Tria::OceanAverage(IssmDouble* Sg){ /*{{{*/
3551
3552 if(IsWaterInElement()){
3553-
3554+
3555 IssmDouble area;
3556
3557 /*Compute area of element:*/
3558@@ -4527,13 +4520,13 @@
3559
3560 if(IsWaterInElement()){
3561 IssmDouble rho_water, S;
3562-
3563+
3564 /*recover material parameters: */
3565 rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
3566-
3567+
3568 /*From Sg_old, recover water sea level rise:*/
3569 S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
3570-
3571+
3572 /* Perturbation terms for moment of inertia (moi_list):
3573 * computed analytically (see Wu & Peltier, eqs 10 & 32)
3574 * also consistent with my GMD formulation!
3575@@ -4545,20 +4538,20 @@
3576 }
3577 else if(this->inputs->Max(MaskIceLevelsetEnum)<0){
3578 IssmDouble rho_ice, I;
3579-
3580+
3581 /*recover material parameters: */
3582 rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
3583-
3584+
3585 /*Compute ice thickness change: */
3586 Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
3587 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
3588 deltathickness_input->GetInputAverage(&I);
3589-
3590+
3591 dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea;
3592 dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea;
3593 dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea;
3594 }
3595-
3596+
3597 return;
3598 }/*}}}*/
3599 void Tria::SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
3600@@ -4591,7 +4584,7 @@
3601 bool computerigid = true;
3602 bool computeelastic= true;
3603 bool scaleoceanarea= false;
3604-
3605+
3606 /*early return if we are not on an ice cap:*/
3607 if(!(this->inputs->Max(MaskIceLevelsetEnum)<=0)){
3608 constant=0; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum));
3609@@ -4605,12 +4598,12 @@
3610 *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
3611 return;
3612 }
3613-
3614+
3615 /*If we are here, we are on ice that is fully grounded or half-way to floating: */
3616 if ((this->inputs->Min(MaskGroundediceLevelsetEnum))<0){
3617 notfullygrounded=true; //used later on.
3618 }
3619-
3620+
3621 /*Inform mask: */
3622 constant=1; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum));
3623
3624@@ -4635,10 +4628,10 @@
3625 this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
3626
3627 /* Where is the centroid of this element?:{{{*/
3628-
3629+
3630 /*retrieve coordinates: */
3631 ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
3632-
3633+
3634 IssmDouble minlong=400;
3635 IssmDouble maxlong=-20;
3636 for (int i=0;i<NUMVERTICES;i++){
3637@@ -4652,12 +4645,12 @@
3638 if (llr_list[1][1]==0)llr_list[1][1]=360;
3639 if (llr_list[2][1]==0)llr_list[2][1]=360;
3640 }
3641-
3642+
3643 // correction at the north pole
3644 if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
3645 if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
3646 if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
3647-
3648+
3649 //correction at the south pole
3650 if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
3651 if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
3652@@ -4683,7 +4676,7 @@
3653 phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
3654 area*=phi;
3655 }
3656-
3657+
3658 /*Compute ice thickness change: */
3659 Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
3660 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
3661@@ -4695,7 +4688,7 @@
3662 bool mainlyfloating = true;
3663 int point1;
3664 IssmDouble fraction1,fraction2;
3665-
3666+
3667 /*Recover portion of element that is grounded*/
3668 this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
3669 Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
3670@@ -4749,12 +4742,12 @@
3671 values[i]=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
3672 }
3673 pSgi->SetValues(gsize,indices,values,ADD_VAL);
3674-
3675+
3676 /*free ressources:*/
3677 xDelete<IssmDouble>(values);
3678 xDelete<int>(indices);
3679 }
3680-
3681+
3682 /*Assign output pointer:*/
3683 _assert_(!xIsNan<IssmDouble>(eustatic));
3684 _assert_(!xIsInf<IssmDouble>(eustatic));
3685@@ -4779,7 +4772,7 @@
3686 /*precomputed elastic green functions:*/
3687 IssmDouble* G_elastic_precomputed = NULL;
3688 int M;
3689-
3690+
3691 /*computation of Green functions:*/
3692 IssmDouble* G_elastic= NULL;
3693 IssmDouble* G_rigid= NULL;
3694@@ -4856,9 +4849,9 @@
3695 late=late/180*PI;
3696 longe=longe/180*PI;
3697 /*}}}*/
3698-
3699+
3700 if(computeelastic){
3701-
3702+
3703 /*recover elastic green function:*/
3704 DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
3705 parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
3706@@ -4896,7 +4889,7 @@
3707 values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
3708 }
3709 }
3710-
3711+
3712 pSgo->SetValues(gsize,indices,values,ADD_VAL);
3713
3714 /*free ressources:*/
3715@@ -4929,7 +4922,7 @@
3716 IssmDouble* U_elastic_precomputed = NULL;
3717 IssmDouble* H_elastic_precomputed = NULL;
3718 int M;
3719-
3720+
3721 /*computation of Green functions:*/
3722 IssmDouble* U_elastic= NULL;
3723 IssmDouble* N_elastic= NULL;
3724@@ -4956,7 +4949,7 @@
3725 /*recover computational flags: */
3726 this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
3727 this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
3728-
3729+
3730 /*early return if elastic not requested:*/
3731 if(!computeelastic) return;
3732
3733@@ -5025,12 +5018,12 @@
3734
3735 /*From Sg, recover water sea level rise:*/
3736 S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
3737-
3738+
3739 /*Compute ice thickness change: */
3740 Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
3741 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
3742 deltathickness_input->GetInputAverage(&I);
3743-
3744+
3745 /*initialize: */
3746 U_elastic=xNewZeroInit<IssmDouble>(gsize);
3747 if(horiz){
3748@@ -5072,7 +5065,7 @@
3749 N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
3750 E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
3751 }
3752-
3753+
3754 /*Elastic component (from Eq 17 in Adhikari et al, GMD 2015): */
3755 int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
3756 U_elastic[i] += U_elastic_precomputed[index];
3757Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
3758===================================================================
3759--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23065)
3760+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 23066)
3761@@ -286,7 +286,6 @@
3762 IssmDouble calvingratey[NUMVERTICES];
3763 IssmDouble calvingrate[NUMVERTICES];
3764
3765-
3766 /* Get node coordinates and dof list: */
3767 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
3768
3769@@ -1047,7 +1046,7 @@
3770 }
3771 /*}}}*/
3772 void Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
3773-
3774+
3775 /* Intermediaries */
3776 const int dim=3;
3777 int i, dir,nrfrontnodes;
3778@@ -1181,7 +1180,6 @@
3779 if(!IsInputEnum(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");
3780 Input* input=(Input*)this->inputs->GetInput(control_enum); _assert_(input);
3781
3782-
3783 /*Cast to Controlinput*/
3784 if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
3785 ControlInput* controlinput = xDynamicCast<ControlInput*>(input);
3786@@ -2594,7 +2592,7 @@
3787
3788 Penta* penta=this;
3789 for(;;){
3790-
3791+
3792 IssmDouble xyz_list[NUMVERTICES][3];
3793 /* Get node coordinates and dof list: */
3794 ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
3795@@ -2603,7 +2601,7 @@
3796 Jdet[0]=(xyz_list[3][2]-xyz_list[0][2])*0.5;
3797 Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5;
3798 Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5;
3799-
3800+
3801 /*Retrieve all inputs we will need*/
3802 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3803 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3804@@ -2614,7 +2612,7 @@
3805 Input* deviayy_input=inputs->GetInput(DeviatoricStressyyEnum); _assert_(deviayy_input);
3806 Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
3807 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3808-
3809+
3810 /* Start looping on the number of 2D vertices: */
3811 for(int ig=0;ig<3;ig++){
3812 GaussPenta* gauss=new GaussPenta(ig,3+ig,11);
3813@@ -2643,10 +2641,10 @@
3814 }
3815 delete gauss;
3816 }
3817-
3818+
3819 /*Stop if we have reached the surface/base*/
3820 if(penta->IsOnSurface()) break;
3821-
3822+
3823 /*get upper Penta*/
3824 penta=penta->GetUpperPenta();
3825 _assert_(penta->Id()!=this->id);
3826Index: ../trunk-jpl/src/c/classes/Elements/ElementHook.cpp
3827===================================================================
3828--- ../trunk-jpl/src/c/classes/Elements/ElementHook.cpp (revision 23065)
3829+++ ../trunk-jpl/src/c/classes/Elements/ElementHook.cpp (revision 23066)
3830@@ -103,7 +103,7 @@
3831 MARSHALLING_DYNAMIC(hnodesi_null,bool,numanalyses);
3832
3833 if(marshall_direction==MARSHALLING_BACKWARD){
3834-
3835+
3836 if (!hnodes_null)this->hnodes = new Hook*[numanalyses];
3837 else this->hnodes=NULL;
3838 this->hvertices = new Hook();
3839@@ -139,7 +139,7 @@
3840
3841 _printf_(" ElementHook DeepEcho:\n");
3842 _printf_(" numanalyses : "<< this->numanalyses <<"\n");
3843-
3844+
3845 _printf_(" hnodes:\n");
3846 if(hnodes){
3847 for(int i=0;i<this->numanalyses;i++) {
3848@@ -148,11 +148,11 @@
3849 }
3850 }
3851 else _printf_(" hnodes = NULL\n");
3852-
3853+
3854 _printf_(" hvertices:\n");
3855 if(hvertices) hvertices->DeepEcho();
3856 else _printf_(" hvertices = NULL\n");
3857-
3858+
3859 _printf_(" hmaterial:\n");
3860 if(hmaterial) hmaterial->DeepEcho();
3861 else _printf_(" hmaterial = NULL\n");
3862@@ -169,10 +169,10 @@
3863 }
3864 /*}}}*/
3865 void ElementHook::Echo(){/*{{{*/
3866-
3867+
3868 _printf_(" ElementHook Echo:\n");
3869 _printf_(" numanalyses : "<< this->numanalyses <<"\n");
3870-
3871+
3872 _printf_(" hnodes:\n");
3873 if(hnodes){
3874 for(int i=0;i<this->numanalyses;i++) {
3875@@ -180,11 +180,11 @@
3876 }
3877 }
3878 else _printf_(" hnodes = NULL\n");
3879-
3880+
3881 _printf_(" hvertices:\n");
3882 if(hvertices) hvertices->Echo();
3883 else _printf_(" hvertices = NULL\n");
3884-
3885+
3886 _printf_(" hmaterial:\n");
3887 if(hmaterial) hmaterial->Echo();
3888 else _printf_(" hmaterial = NULL\n");
3889Index: ../trunk-jpl/src/c/classes/Elements/Seg.cpp
3890===================================================================
3891--- ../trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 23065)
3892+++ ../trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 23066)
3893@@ -103,7 +103,6 @@
3894
3895 return seg;
3896
3897-
3898 }
3899 /*}}}*/
3900 void Seg::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
3901@@ -140,7 +139,7 @@
3902 }
3903 /*}}}*/
3904 void Seg::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
3905-
3906+
3907 /* Intermediaries */
3908 int nrfrontnodes,index;
3909 IssmDouble levelset[NUMVERTICES];
3910Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
3911===================================================================
3912--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23065)
3913+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23066)
3914@@ -1104,7 +1104,6 @@
3915 /*}}}*/
3916 void Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug){/*{{{*/
3917
3918-
3919 /*Get number of nodes for this element*/
3920 int numnodes = this->GetNumberOfNodes();
3921
3922@@ -1121,7 +1120,6 @@
3923 if(ug[nodes[i]->GetDof(0,GsetEnum)] > input_max) input_max = ug[nodes[i]->GetDof(0,GsetEnum)];
3924 }
3925
3926-
3927 /*Second loop to reassign min and max with local extrema*/
3928 for(int i=0;i<numnodes;i++){
3929 if(min[nodes[i]->GetDof(0,GsetEnum)]>input_min) min[nodes[i]->GetDof(0,GsetEnum)] = input_min;
3930@@ -1388,7 +1386,7 @@
3931 default:
3932 _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
3933 }
3934-
3935+
3936 /*Clean up*/
3937 xDelete<int>(doflist);
3938 xDelete<IssmDouble>(values);
3939@@ -1476,7 +1474,6 @@
3940 parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
3941 parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
3942
3943-
3944 /*Get number of vertices*/
3945 int numvertices = this->GetNumberOfVertices();
3946 if(isautodiff){
3947@@ -1807,7 +1804,6 @@
3948 this->inputs->AddInput(datasetinput);
3949 }
3950
3951-
3952 /*Branch on type of vector: nodal or elementary: */
3953 if(vector_type==1){ //nodal vector
3954
3955@@ -2066,7 +2062,7 @@
3956 this->GetInputListOnVertices(deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum);
3957 this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum);
3958 this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum);
3959-
3960+
3961 for(int i=0;i<numvertices;i++){
3962 if(base[i]>upperwaterel[i]) values[i]=0;
3963 else if (base[i]<deepwaterel[i]) values[i]=deepwatermelt[i];
3964@@ -2869,7 +2865,6 @@
3965 /*only compute SMB at the surface: */
3966 if (!IsOnSurface()) return;
3967
3968-
3969 /*Retrieve material properties and parameters:{{{ */
3970 rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
3971 rho_water = matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
3972@@ -3071,7 +3066,6 @@
3973 /*Snow, firn and ice albedo:*/
3974 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());
3975
3976-
3977 /*Distribution of absorbed short wave radation with depth:*/
3978 if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
3979
3980Index: ../trunk-jpl/src/c/classes/gauss/GaussSeg.cpp
3981===================================================================
3982--- ../trunk-jpl/src/c/classes/gauss/GaussSeg.cpp (revision 23065)
3983+++ ../trunk-jpl/src/c/classes/gauss/GaussSeg.cpp (revision 23066)
3984@@ -30,7 +30,7 @@
3985 /*Get gauss points*/
3986 this->numgauss = order;
3987 GaussLegendreLinear(&pcoords1,&pweights,order);
3988-
3989+
3990 this->coords1=xNew<IssmDouble>(numgauss);
3991 this->weights=xNew<IssmDouble>(numgauss);
3992
3993Index: ../trunk-jpl/src/c/classes/Regionaloutput.cpp
3994===================================================================
3995--- ../trunk-jpl/src/c/classes/Regionaloutput.cpp (revision 23065)
3996+++ ../trunk-jpl/src/c/classes/Regionaloutput.cpp (revision 23066)
3997@@ -154,4 +154,3 @@
3998 return val_t;
3999 }
4000 /*}}}*/
4001-
4002Index: ../trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp
4003===================================================================
4004--- ../trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 23065)
4005+++ ../trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp (revision 23066)
4006@@ -19,7 +19,7 @@
4007
4008 int world_rank;
4009 ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&world_rank);
4010-
4011+
4012 /*Build an femmodel if you are a slave, using the corresponding communicator:*/
4013 if(world_rank!=0){
4014 femmodel_init= new FemModel(argc,argv,evaluation_comm);
4015@@ -32,7 +32,7 @@
4016
4017 int world_rank;
4018 ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&world_rank);
4019-
4020+
4021 if(world_rank!=0){
4022
4023 /*Wrap up: */
4024Index: ../trunk-jpl/src/c/classes/Node.cpp
4025===================================================================
4026--- ../trunk-jpl/src/c/classes/Node.cpp (revision 23065)
4027+++ ../trunk-jpl/src/c/classes/Node.cpp (revision 23066)
4028@@ -497,7 +497,7 @@
4029 }
4030 /*}}}*/
4031 void Node::SetApproximation(int in_approximation){/*{{{*/
4032-
4033+
4034 this->approximation = in_approximation;
4035 }
4036 /*}}}*/
4037Index: ../trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp
4038===================================================================
4039--- ../trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp (revision 23065)
4040+++ ../trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp (revision 23066)
4041@@ -36,7 +36,7 @@
4042
4043 /*Object virtual functions definitions:*/
4044 Object* SpcStatic::copy() {/*{{{*/
4045-
4046+
4047 SpcStatic* spcstat = new SpcStatic(*this);
4048
4049 spcstat->id=this->id;
4050Index: ../trunk-jpl/src/c/classes/IoModel.cpp
4051===================================================================
4052--- ../trunk-jpl/src/c/classes/IoModel.cpp (revision 23065)
4053+++ ../trunk-jpl/src/c/classes/IoModel.cpp (revision 23066)
4054@@ -156,7 +156,7 @@
4055
4056 bool autodiff=false;
4057 bool iscontrol=false;
4058-
4059+
4060 /*First, keep track of the file handle: */
4061 this->fid=iomodel_handle;
4062
4063@@ -421,7 +421,7 @@
4064 /*Initialize array detecting whether data[i] is an independent AD mode variable: */
4065 this->FetchData(&autodiff,"md.autodiff.isautodiff");
4066 this->FetchData(&iscontrol,"md.inversion.iscontrol");
4067-
4068+
4069 if(trace || (autodiff && !iscontrol)){
4070
4071 #ifdef _HAVE_ADOLC_
4072@@ -505,7 +505,7 @@
4073 void IoModel::DeleteData(char*** pstringarray, int numstrings,const char* data_name){/*{{{*/
4074
4075 char** stringarray=*pstringarray;
4076-
4077+
4078 if(numstrings){
4079 for(int i=0;i<numstrings;i++){
4080 char* string=stringarray[i];
4081@@ -596,7 +596,7 @@
4082 case 2:
4083 /*Read the integer and broadcast it to other cpus:*/
4084 if(fread(&integer,sizeof(int),1,this->fid)!=1) _error_("could not read integer ");
4085-
4086+
4087 /*Convert codes to Enums if needed*/
4088 if(strcmp(record_name,"md.smb.model")==0) integer = IoCodeToEnumSMB(integer);
4089 if(strcmp(record_name,"md.basalforcings.model")==0) integer = IoCodeToEnumBasal(integer);
4090@@ -962,7 +962,7 @@
4091
4092 /*recover my_rank:*/
4093 int my_rank=IssmComm::GetRank();
4094-
4095+
4096 /*Set file pointer to beginning of the data: */
4097 fid=this->SetFilePointerToData(&code,NULL,data_name);
4098
4099@@ -1133,7 +1133,7 @@
4100 return;
4101 }
4102 }
4103-
4104+
4105 /*output: */
4106 int M,N;
4107 IssmPDouble *matrix = NULL;
4108@@ -1824,7 +1824,7 @@
4109 if(my_rank==0){
4110 /*check we are indeed finding a string, not something else: */
4111 if(codes[i]!=4)_error_("expecting a string for \""<<data_name<<"\"");
4112-
4113+
4114 /*We have to read a string from disk. First read the dimensions of the string, then the string: */
4115 fsetpos(fid,file_positions+i);
4116 if(fread(&string_size,sizeof(int),1,fid)!=1) _error_("could not read length of string ");
4117@@ -1853,7 +1853,7 @@
4118 /*Free ressources:*/
4119 xDelete<int>(codes);
4120 xDelete<fpos_t>(file_positions);
4121-
4122+
4123 /*Assign output pointers: */
4124 *pstrings=strings;
4125 *pnumstrings=num_instances;
4126@@ -1874,7 +1874,7 @@
4127
4128 /*recover my_rank:*/
4129 int my_rank=IssmComm::GetRank();
4130-
4131+
4132 /*Get file pointers to beginning of the data (multiple instances of it): */
4133 file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
4134
4135@@ -1889,7 +1889,7 @@
4136 code=codes[i];
4137
4138 if(code!=2)_error_("expecting an integer for \""<<data_name<<"\"");
4139-
4140+
4141 /*We have to read a integer from disk. First read the dimensions of the integer, then the integer: */
4142 fsetpos(fid,file_positions+i);
4143 if(my_rank==0){
4144@@ -1902,7 +1902,7 @@
4145 vector[i]=integer;
4146 }
4147 }
4148-
4149+
4150 /*Free ressources:*/
4151 xDelete<fpos_t>(file_positions);
4152 xDelete<int>(codes);
4153@@ -1927,7 +1927,7 @@
4154
4155 /*recover my_rank:*/
4156 int my_rank=IssmComm::GetRank();
4157-
4158+
4159 /*Get file pointers to beginning of the data (multiple instances of it): */
4160 file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
4161
4162@@ -1942,7 +1942,7 @@
4163 code=codes[i];
4164
4165 if(code!=3)_error_("expecting a double for \""<<data_name<<"\"");
4166-
4167+
4168 /*We have to read a double from disk: */
4169 fsetpos(fid,file_positions+i);
4170 if(my_rank==0){
4171@@ -1955,7 +1955,7 @@
4172 vector[i]=scalar;
4173 }
4174 }
4175-
4176+
4177 /*Free ressources:*/
4178 xDelete<fpos_t>(file_positions);
4179 xDelete<int>(codes);
4180@@ -1984,7 +1984,7 @@
4181
4182 /*recover my_rank:*/
4183 int my_rank=IssmComm::GetRank();
4184-
4185+
4186 /*Get file pointers to beginning of the data (multiple instances of it): */
4187 file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
4188
4189@@ -2014,7 +2014,6 @@
4190 }
4191 ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
4192
4193-
4194 /*Now allocate matrix: */
4195 if(M*N){
4196 pmatrix=xNew<IssmPDouble>(M*N);
4197@@ -2038,8 +2037,7 @@
4198 }
4199 else
4200 matrix=NULL;
4201-
4202-
4203+
4204 /*Assign: */
4205 mdims[i]=M;
4206 matrices[i]=matrix;
4207@@ -2046,7 +2044,7 @@
4208 ndims[i]=N;
4209 }
4210 }
4211-
4212+
4213 /*Free ressources:*/
4214 xDelete<fpos_t>(file_positions);
4215 xDelete<int>(codes);
4216@@ -2088,7 +2086,7 @@
4217
4218 /*recover my_rank:*/
4219 int my_rank=IssmComm::GetRank();
4220-
4221+
4222 /*Get file pointers to beginning of the data (multiple instances of it): */
4223 file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
4224
4225@@ -2118,7 +2116,6 @@
4226 }
4227 ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
4228
4229-
4230 /*Now allocate matrix: */
4231 if(M*N){
4232 pmatrix=xNew<IssmPDouble>(M*N);
4233@@ -2143,8 +2140,7 @@
4234 }
4235 else
4236 integer_matrix=NULL;
4237-
4238-
4239+
4240 /*Assign: */
4241 mdims[i]=M;
4242 matrices[i]=integer_matrix;
4243@@ -2151,7 +2147,7 @@
4244 ndims[i]=N;
4245 }
4246 }
4247-
4248+
4249 /*Free ressources:*/
4250 xDelete<fpos_t>(file_positions);
4251 xDelete<int>(codes);
4252@@ -2375,7 +2371,7 @@
4253 codes = xNew<int>(num_instances);
4254 vector_types = xNew<int>(num_instances);
4255 }
4256-
4257+
4258 /*Reset FILE* position to the beginning of the file, and start again, this time saving the data information
4259 * as we find it: */
4260 counter=0;
4261@@ -2418,12 +2414,12 @@
4262 codes[counter] = record_code;
4263 vector_types[counter] = vector_type;
4264 fgetpos(fid,file_positions+counter);
4265-
4266+
4267 /*backup and skip over the record, as we have more work to do: */
4268 if(5<=record_code && record_code<=7) fseek(fid,-sizeof(int),SEEK_CUR);
4269 fseek(fid,-sizeof(int),SEEK_CUR);
4270 fseek(fid,-sizeof(int),SEEK_CUR);
4271-
4272+
4273 /*increment counter: */
4274 counter++;
4275 }
4276Index: ../trunk-jpl/src/c/classes/Nodalvalue.cpp
4277===================================================================
4278--- ../trunk-jpl/src/c/classes/Nodalvalue.cpp (revision 23065)
4279+++ ../trunk-jpl/src/c/classes/Nodalvalue.cpp (revision 23066)
4280@@ -86,7 +86,7 @@
4281 }
4282 /*}}}*/
4283 IssmDouble Nodalvalue::Response(FemModel* femmodel){/*{{{*/
4284-
4285+
4286 /*output:*/
4287 IssmDouble value;
4288
4289Index: ../trunk-jpl/src/c/classes/Misfit.cpp
4290===================================================================
4291--- ../trunk-jpl/src/c/classes/Misfit.cpp (revision 23065)
4292+++ ../trunk-jpl/src/c/classes/Misfit.cpp (revision 23066)
4293@@ -22,7 +22,7 @@
4294 #include "../classes/Inputs/Input.h"
4295 #include "../classes/gauss/Gauss.h"
4296 /*}}}*/
4297-
4298+
4299 /*Misfit constructors, destructors :*/
4300 Misfit::Misfit(){/*{{{*/
4301
4302@@ -41,18 +41,18 @@
4303 Misfit::Misfit(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, int in_local, int in_weights_enum){/*{{{*/
4304
4305 this->definitionenum=in_definitionenum;
4306-
4307+
4308 this->name = xNew<char>(strlen(in_name)+1);
4309 xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
4310
4311 this->timeinterpolation = xNew<char>(strlen(in_timeinterpolation)+1);
4312 xMemCpy<char>(this->timeinterpolation,in_timeinterpolation,strlen(in_timeinterpolation)+1);
4313-
4314+
4315 this->model_enum=in_model_enum;
4316 this->observation_enum=in_observation_enum;
4317 this->weights_enum=in_weights_enum;
4318 this->local=in_local;
4319-
4320+
4321 this->misfit=0;
4322 this->lock=0;
4323 }
4324@@ -110,11 +110,11 @@
4325 }
4326 /*}}}*/
4327 IssmDouble Misfit::Response(FemModel* femmodel){/*{{{*/
4328-
4329+
4330 /*diverse: */
4331 IssmDouble time,starttime,finaltime;
4332 IssmDouble dt;
4333-
4334+
4335 /*recover time parameters: */
4336 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
4337 femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
4338@@ -129,7 +129,6 @@
4339 IssmDouble area_t=0.;
4340 IssmDouble all_area_t;
4341
4342-
4343 /*If we are locked, return time average: */
4344 if(this->lock)return misfit/(time-starttime);
4345
4346@@ -143,7 +142,7 @@
4347 ISSM_MPI_Allreduce ( (void*)&area_t,(void*)&all_area_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
4348 area_t=all_area_t;
4349 misfit_t=all_misfit_t;
4350-
4351+
4352 /*Divide by surface area if not nill!: */
4353 if (area_t!=0) misfit_t=misfit_t/area_t;
4354
4355@@ -163,11 +162,11 @@
4356 IssmDouble* observation= NULL;
4357 IssmDouble* weights= NULL;
4358 int msize,osize,wsize;
4359-
4360+
4361 /*Are we transient?:*/
4362 if (time==0){
4363 IssmDouble misfit_t=0.;
4364-
4365+
4366 /*get global vectors: */
4367 GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
4368 GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
4369@@ -189,7 +188,7 @@
4370 return misfit;
4371 }
4372 else{
4373-
4374+
4375 IssmDouble misfit_t=0.;
4376 IssmDouble all_misfit_t=0.;
4377
4378@@ -200,7 +199,7 @@
4379 GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
4380 GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
4381 GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize);
4382-
4383+
4384 int count=0;
4385 for (int i=0;i<msize;i++){
4386 misfit_t += pow(model[i]-observation[i],2)*weights[i];
4387@@ -214,7 +213,7 @@
4388
4389 /*Do we lock? i.e. are we at final_time? :*/
4390 if(time==finaltime)this->lock=1;
4391-
4392+
4393 /*Free ressources:*/
4394 xDelete<IssmDouble>(model);
4395 xDelete<IssmDouble>(observation);
4396@@ -226,9 +225,9 @@
4397
4398 } /*}}}*/
4399 else{ /*global computation: {{{ */
4400-
4401+
4402 IssmDouble model, observation;
4403-
4404+
4405 /*If we are locked, return time average: */
4406 if(this->lock)return misfit/(time-starttime);
4407
4408@@ -238,13 +237,13 @@
4409 Element* element=(Element*)femmodel->elements->GetObjectByOffset(0); _assert_(element);
4410 Input* input = element->GetInput(observation_enum); _assert_(input);
4411 input->GetInputAverage(&observation);
4412-
4413+
4414 /*Add this time's contribution to curent misfit: */
4415 misfit+=dt*(model-observation);
4416-
4417+
4418 /*Do we lock? i.e. are we at final_time? :*/
4419 if(time==finaltime)this->lock=1;
4420-
4421+
4422 /*What we return is the value of misfit / time: */
4423 return misfit/(time-starttime);
4424 } /*}}}*/
4425Index: ../trunk-jpl/src/c/datastructures/DataSet.cpp
4426===================================================================
4427--- ../trunk-jpl/src/c/datastructures/DataSet.cpp (revision 23065)
4428+++ ../trunk-jpl/src/c/datastructures/DataSet.cpp (revision 23066)
4429@@ -91,7 +91,7 @@
4430
4431 /*Specific methods*/
4432 void DataSet::Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
4433-
4434+
4435 vector<Object*>::iterator obj;
4436 int obj_size=0;
4437 int obj_enum=0;
4438Index: ../trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp
4439===================================================================
4440--- ../trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 23065)
4441+++ ../trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp (revision 23066)
4442@@ -593,7 +593,6 @@
4443 xDelete<int>(responses);
4444 delete gauss;
4445
4446-
4447 }/*}}}*/
4448 void BalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
4449
4450Index: ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp
4451===================================================================
4452--- ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 23065)
4453+++ ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 23066)
4454@@ -260,7 +260,6 @@
4455 basis,1,numnodes,0,
4456 &Ke->values[0],1);
4457
4458-
4459 /*Transfer EPL part*/
4460 transfer=GetHydrologyKMatrixTransfer(basalelement);
4461 D_scalar=dt*transfer*gauss->weight*Jdet;
4462Index: ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
4463===================================================================
4464--- ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 23065)
4465+++ ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 23066)
4466@@ -1152,7 +1152,7 @@
4467
4468 /*Clean up and return*/
4469 xDelete<int>(responses);
4470-
4471+
4472 }/*}}}*/
4473 void AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
4474 /*WARNING: We use SSA as an estimate for now*/
4475@@ -2111,7 +2111,6 @@
4476 IssmDouble vx,vy,lambda,mu;
4477 IssmDouble *xyz_list= NULL;
4478
4479-
4480 /*Fetch number of vertices for this finite element*/
4481 int numvertices = basalelement->GetNumberOfVertices();
4482
4483@@ -2131,8 +2130,6 @@
4484 Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input);
4485 Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input);
4486
4487-
4488-
4489 IssmDouble q_exp;
4490 IssmDouble C_param;
4491 IssmDouble As;
4492@@ -2181,11 +2178,11 @@
4493 vmag = sqrt(vx*vx + vy*vy);
4494 Chi = vmag/(pow(C_param,n)*pow(Neff,n)*As);
4495 Gamma = (Chi/(1.+alpha*pow(Chi,q_exp)));
4496-
4497+
4498 Uder =Neff*C_param/(vmag*vmag*n) *
4499 (Gamma-alpha*q_exp*pow(Chi,q_exp-1.)*Gamma*Gamma* pow(Gamma,(1.-n)/n) -
4500 n* pow(Gamma,1./n));
4501-
4502+
4503 /*Build gradient vector (actually -dJ/dD): */
4504 for(int i=0;i<numvertices;i++){
4505 ge[i]+=-dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
4506@@ -2318,7 +2315,6 @@
4507 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
4508 }
4509
4510-
4511 /*Fetch number of nodes and dof for this finite element*/
4512 int vnumnodes = element->NumberofNodesVelocity();
4513 int pnumnodes = element->NumberofNodesPressure();
4514Index: ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
4515===================================================================
4516--- ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23065)
4517+++ ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 23066)
4518@@ -42,7 +42,7 @@
4519 counter++;
4520 }
4521 }
4522-
4523+
4524 iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
4525 iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
4526 iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
4527@@ -176,7 +176,7 @@
4528 }
4529
4530 /*Calving threshold*/
4531-
4532+
4533 /*Fetch number of nodes and dof for this finite element*/
4534 int numnodes = basalelement->GetNumberOfNodes();
4535
4536@@ -325,7 +325,7 @@
4537 c[i]=0.; m[i]=0.;
4538 }
4539 break;
4540-
4541+
4542 case CalvingLevermannEnum:
4543 calvingratex_input->GetInputValue(&c[0],gauss);
4544 if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
4545@@ -356,7 +356,7 @@
4546 m[i]=0.;
4547 }
4548 break;
4549-
4550+
4551 case CalvingHabEnum:
4552 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
4553 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
4554@@ -377,7 +377,7 @@
4555 m[i]=0.;
4556 }
4557 break;
4558-
4559+
4560 case CalvingCrevasseDepthEnum:
4561 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
4562 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
4563@@ -384,7 +384,7 @@
4564 meltingrate_input->GetInputValue(&meltingrate,gauss);
4565
4566 if(groundedice<0) meltingrate = 0.;
4567-
4568+
4569 norm_dlsf=0.;
4570 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
4571 norm_dlsf=sqrt(norm_dlsf);
4572@@ -515,7 +515,7 @@
4573 return Ke;
4574 }/*}}}*/
4575 ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
4576-
4577+
4578 if(!element->IsOnBase()) return NULL;
4579 Element* basalelement = element->SpawnBasalElement();
4580
4581@@ -524,7 +524,7 @@
4582 IssmDouble Jdet,dt;
4583 IssmDouble lsf;
4584 IssmDouble* xyz_list = NULL;
4585-
4586+
4587 /*Fetch number of nodes and dof for this finite element*/
4588 int numnodes = basalelement->GetNumberOfNodes();
4589
4590@@ -531,7 +531,7 @@
4591 /*Initialize Element vector*/
4592 ElementVector* pe = basalelement->NewElementVector();
4593 basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
4594-
4595+
4596 if(dt!=0.){
4597 /*Initialize basis vector*/
4598 IssmDouble* basis = xNew<IssmDouble>(numnodes);
4599@@ -622,7 +622,7 @@
4600 // returns distance d of point q to straight going through points s0, s1
4601 // d=|a x b|/|b|
4602 // with a=q-s0, b=s1-s0
4603-
4604+
4605 /* Intermediaries */
4606 const int dim=2;
4607 int i;
4608@@ -633,7 +633,7 @@
4609 a[i]=q[i]-s0[i];
4610 b[i]=s1[i]-s0[i];
4611 }
4612-
4613+
4614 norm_b=0.;
4615 for(i=0;i<dim;i++)
4616 norm_b+=b[i]*b[i];
4617@@ -670,7 +670,7 @@
4618 IssmDouble rho_ice,rho_water;
4619 IssmDouble bed,water_depth;
4620 IssmDouble levelset;
4621-
4622+
4623 femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
4624
4625 if(calvinglaw==CalvingMinthicknessEnum){
4626@@ -702,17 +702,17 @@
4627 delete gauss;
4628 }
4629 }
4630-
4631+
4632 if(calvinglaw==CalvingHabEnum){
4633
4634 /*Get the fraction of the flotation thickness at the terminus*/
4635 femmodel->elements->InputDuplicate(MaskIceLevelsetEnum, DistanceToCalvingfrontEnum);
4636 femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
4637-
4638+
4639 /*Loop over all elements of this partition*/
4640 for(int i=0;i<femmodel->elements->Size();i++){
4641 Element* element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
4642-
4643+
4644 rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
4645 rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
4646
4647@@ -743,13 +743,13 @@
4648 delete gauss;
4649 }
4650 }
4651-
4652+
4653 if(calvinglaw==CalvingCrevasseDepthEnum){
4654-
4655+
4656 int nflipped,local_nflipped;
4657 Vector<IssmDouble>* vec_constraint_nodes = NULL;
4658 IssmDouble* constraint_nodes = NULL;
4659-
4660+
4661 /*Get the DistanceToCalvingfront*/
4662 femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
4663 femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
4664@@ -841,7 +841,7 @@
4665 for(int in=0;in<numnodes;in++){
4666 gauss->GaussNode(element->GetElementType(),in);
4667 Node* node=element->GetNode(in);
4668-
4669+
4670 if(constraint_nodes[node->Sid()]>0.){
4671 node->ApplyConstraint(0,+1.);
4672 }
4673Index: ../trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
4674===================================================================
4675--- ../trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 23065)
4676+++ ../trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 23066)
4677@@ -29,7 +29,6 @@
4678 else
4679 iscoupling = false;
4680
4681-
4682 /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
4683 if(!iscoupling){
4684 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvz",StressbalanceVerticalAnalysisEnum,P1Enum,0);
4685@@ -158,7 +157,6 @@
4686 }/*}}}*/
4687 ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixBase(Element* element){/*{{{*/
4688
4689-
4690 if(!element->IsOnBase()) return NULL;
4691
4692 /*Intermediaries*/
4693@@ -199,7 +197,6 @@
4694 }/*}}}*/
4695 ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/
4696
4697-
4698 if(!element->IsOnSurface()) return NULL;
4699
4700 /*Intermediaries*/
4701Index: ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp
4702===================================================================
4703--- ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 23065)
4704+++ ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 23066)
4705@@ -294,7 +294,7 @@
4706 B_input->GetInputValue(&B,gauss);
4707 n_input->GetInputValue(&n,gauss);
4708 A=pow(B,-n);
4709-
4710+
4711 /*Compute beta term*/
4712 if(gap<br)
4713 beta = (br-gap)/lr;
4714@@ -323,7 +323,7 @@
4715
4716 meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
4717 _assert_(meltrate>0.);
4718-
4719+
4720 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
4721 (
4722 meltrate*(1/rho_water-1/rho_ice)
4723@@ -397,7 +397,7 @@
4724
4725 /*Calculate effective pressure*/
4726 eff_pressure[i] = rho_ice*g*thickness[i] - rho_water*g*(values[i]-bed[i]);
4727-
4728+
4729 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
4730 if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
4731 }
4732@@ -443,7 +443,7 @@
4733 Input* gap_input = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input);
4734 reynolds_input->GetInputAverage(&reynolds);
4735 gap_input->GetInputAverage(&gap);
4736-
4737+
4738 /*Compute conductivity*/
4739 IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds));
4740 _assert_(conductivity>0);
4741@@ -453,7 +453,6 @@
4742 }/*}}}*/
4743 void HydrologyShaktiAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/
4744
4745-
4746 for(int j=0;j<femmodel->elements->Size();j++){
4747 Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
4748 UpdateGapHeight(element);
4749@@ -546,7 +545,7 @@
4750 IssmDouble pressure_ice = rho_ice*g*thickness; _assert_(pressure_ice>0.);
4751 IssmDouble pressure_water = rho_water*g*(head-bed);
4752 if(pressure_water>pressure_ice) pressure_water = pressure_ice;
4753-
4754+
4755 /* Compute change in sensible heat due to changes in pressure melting point*/
4756 dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
4757 dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
4758@@ -580,7 +579,7 @@
4759
4760 /*Add new gap as an input*/
4761 element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum);
4762-
4763+
4764 /*Divide by connectivity, add basal flux as an input*/
4765 q = q/totalweights;
4766 element->AddInput(HydrologyBasalFluxEnum,&q,P0Enum);
4767Index: ../trunk-jpl/src/c/analyses/EsaAnalysis.cpp
4768===================================================================
4769--- ../trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 23065)
4770+++ ../trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 23066)
4771@@ -39,7 +39,7 @@
4772 int nl;
4773 IssmDouble* love_h=NULL;
4774 IssmDouble* love_l=NULL;
4775-
4776+
4777 IssmDouble* U_elastic = NULL;
4778 IssmDouble* U_elastic_local = NULL;
4779 IssmDouble* H_elastic = NULL;
4780@@ -68,7 +68,7 @@
4781 M=reCast<int,IssmDouble>(180./degacc+1.);
4782 U_elastic=xNew<IssmDouble>(M);
4783 H_elastic=xNew<IssmDouble>(M);
4784-
4785+
4786 /*compute combined legendre + love number (elastic green function:*/
4787 m=DetermineLocalSize(M,IssmComm::GetComm());
4788 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
4789@@ -95,7 +95,7 @@
4790 IssmDouble deltalove_U;
4791
4792 deltalove_U = (love_h[n]-love_h[nl-1]);
4793-
4794+
4795 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
4796 if(n==0){
4797 Pn=1;
4798@@ -198,7 +198,7 @@
4799 void EsaAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
4800 /*Default, do nothing*/
4801 return;
4802-
4803+
4804 }/*}}}*/
4805 void EsaAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
4806 /*Default, do nothing*/
4807Index: ../trunk-jpl/src/c/analyses/SmoothAnalysis.cpp
4808===================================================================
4809--- ../trunk-jpl/src/c/analyses/SmoothAnalysis.cpp (revision 23065)
4810+++ ../trunk-jpl/src/c/analyses/SmoothAnalysis.cpp (revision 23066)
4811@@ -152,7 +152,6 @@
4812 default: input = element->GetInput(input_enum);
4813 }
4814
4815-
4816 /* Start looping on the number of gaussian points: */
4817 Gauss* gauss=element->NewGauss(2);
4818 for(int ig=gauss->begin();ig<gauss->end();ig++){
4819@@ -161,7 +160,6 @@
4820 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
4821 element->NodalFunctions(basis,gauss);
4822
4823-
4824 switch(input_enum){
4825 case DrivingStressXEnum:
4826 case DrivingStressYEnum:{
4827Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
4828===================================================================
4829--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 23065)
4830+++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 23066)
4831@@ -164,7 +164,7 @@
4832 }
4833 iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
4834 iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
4835-
4836+
4837 /*Initialize cumdeltalthickness input*/
4838 InputUpdateFromConstantx(elements,0.,SealevelriseCumDeltathicknessEnum);
4839
4840@@ -202,7 +202,7 @@
4841 iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
4842 iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
4843 }
4844-
4845+
4846 }/*}}}*/
4847 void MasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
4848
4849@@ -220,8 +220,6 @@
4850 parameters->AddObject(new IntParam(MasstransportNumRequestedOutputsEnum,numoutputs));
4851 if(numoutputs)parameters->AddObject(new StringArrayParam(MasstransportRequestedOutputsEnum,requestedoutputs,numoutputs));
4852 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.masstransport.requested_outputs");
4853-
4854-
4855
4856 }/*}}}*/
4857
4858@@ -804,7 +802,7 @@
4859 xDelete<IssmDouble>(phi);
4860 xDelete<IssmDouble>(sealevel);
4861 xDelete<IssmDouble>(bed);
4862-
4863+
4864 xDelete<int>(doflist);
4865 if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
4866 }/*}}}*/
4867Index: ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp
4868===================================================================
4869--- ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 23065)
4870+++ ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 23066)
4871@@ -68,7 +68,7 @@
4872 IssmDouble* love_h=NULL;
4873 IssmDouble* love_k=NULL;
4874 IssmDouble* love_l=NULL;
4875-
4876+
4877 bool elastic=false;
4878 IssmDouble* G_elastic = NULL;
4879 IssmDouble* G_elastic_local = NULL;
4880@@ -131,7 +131,6 @@
4881 H_elastic=xNew<IssmDouble>(M);
4882 #endif
4883
4884-
4885 /*compute combined legendre + love number (elastic green function:*/
4886 m=DetermineLocalSize(M,IssmComm::GetComm());
4887 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
4888@@ -167,7 +166,7 @@
4889
4890 deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
4891 deltalove_U = (love_h[n]-love_h[nl-1]);
4892-
4893+
4894 /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
4895 if(n==0){
4896 Pn=1;
4897@@ -229,7 +228,7 @@
4898 xDelete<IssmDouble>(H_elastic);
4899 xDelete<IssmDouble>(H_elastic_local);
4900 }
4901-
4902+
4903 /*Transitions: */
4904 iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions");
4905 if(transitions){
4906@@ -249,7 +248,6 @@
4907 if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
4908 iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs");
4909
4910-
4911 }/*}}}*/
4912
4913 /*Finite Element Analysis*/
4914Index: ../trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp
4915===================================================================
4916--- ../trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (revision 23065)
4917+++ ../trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (revision 23066)
4918@@ -198,8 +198,6 @@
4919 xDelete<IssmDouble>(lset);
4920 }
4921
4922-
4923-
4924 return;
4925 }/*}}}*/
4926
4927Index: ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp
4928===================================================================
4929--- ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 23065)
4930+++ ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp (revision 23066)
4931@@ -378,7 +378,6 @@
4932 return;
4933 }/*}}}*/
4934
4935-
4936 /*Needed changes to switch to the Johnson formulation*//*{{{*/
4937 /*All the changes are to be done in the velocity computation.
4938 The new velocity needs some new parameter that should be introduce in the hydrologyshreve class:
4939Index: ../trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp
4940===================================================================
4941--- ../trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp (revision 23065)
4942+++ ../trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp (revision 23066)
4943@@ -199,7 +199,6 @@
4944 Ny[i] = -H[i]*Ny[i]/norm;
4945 }
4946
4947-
4948 /* Start looping on the number of gaussian points: */
4949 Gauss* gauss=basalelement->NewGauss(2);
4950 for(int ig=gauss->begin();ig<gauss->end();ig++){
4951Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp
4952===================================================================
4953--- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 23065)
4954+++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 23066)
4955@@ -213,7 +213,6 @@
4956 return NULL;
4957 }
4958
4959-
4960 /*Intermediaries */
4961 bool active_element,isefficientlayer;
4962 IssmDouble D_scalar,Jdet,dt;
4963@@ -503,7 +502,6 @@
4964 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
4965 }
4966
4967-
4968 /*Fetch number of nodes for this finite element*/
4969 int numnodes = basalelement->GetNumberOfNodes();
4970
4971Index: ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp
4972===================================================================
4973--- ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23065)
4974+++ ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23066)
4975@@ -18,10 +18,10 @@
4976 return 1;
4977 }/*}}}*/
4978 void SmbAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
4979-
4980+
4981 int smb_model;
4982 bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled;
4983-
4984+
4985 /*Update elements: */
4986 int counter=0;
4987 for(int i=0;i<iomodel->numberofelements;i++){
4988@@ -31,10 +31,10 @@
4989 counter++;
4990 }
4991 }
4992-
4993+
4994 /*Figure out smb model: */
4995 iomodel->FindConstant(&smb_model,"md.smb.model");
4996-
4997+
4998 switch(smb_model){
4999 case SMBforcingEnum:
5000 iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
5001@@ -141,8 +141,6 @@
5002 default:
5003 _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
5004 }
5005-
5006-
5007
5008 }/*}}}*/
5009 void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
5010@@ -153,12 +151,12 @@
5011 int smb_model;
5012 IssmDouble *temp = NULL;
5013 int N,M;
5014-
5015+
5016 parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
5017-
5018+
5019 iomodel->FindConstant(&smb_model,"md.smb.model");
5020 iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
5021-
5022+
5023 switch(smb_model){
5024 case SMBforcingEnum:
5025 /*Nothing to add to parameters*/
5026@@ -197,7 +195,7 @@
5027 iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
5028 parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M));
5029 iomodel->DeleteData(temp,"md.smb.Pfac");
5030-
5031+
5032 iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
5033 parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M));
5034 iomodel->DeleteData(temp,"md.smb.Tdiff");
5035@@ -265,7 +263,7 @@
5036
5037 /*Figure out smb model: */
5038 femmodel->parameters->FindParam(&smb_model,SmbEnum);
5039-
5040+
5041 /*branch to correct module*/
5042 switch(smb_model){
5043 case SMBforcingEnum:
5044Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
5045===================================================================
5046--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 23065)
5047+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 23066)
5048@@ -125,7 +125,7 @@
5049
5050 /*Add input*/
5051 element->AddInput(DamageFEnum,f,element->GetElementType());
5052-
5053+
5054 /*Clean up and return*/
5055 xDelete<IssmDouble>(f);
5056 }/*}}}*/
5057@@ -169,7 +169,7 @@
5058 Gauss* gauss=element->NewGauss();
5059 for (int i=0;i<numnodes;i++){
5060 gauss->GaussNode(element->GetElementType(),i);
5061-
5062+
5063 eps_xx_input->GetInputValue(&eps_xx,gauss);
5064 eps_xy_input->GetInputValue(&eps_xy,gauss);
5065 eps_yy_input->GetInputValue(&eps_yy,gauss);
5066@@ -176,7 +176,7 @@
5067 B_input->GetInputValue(&B,gauss);
5068 n_input->GetInputValue(&n,gauss);
5069 damage_input->GetInputValue(&damage,gauss);
5070-
5071+
5072 /*Calculate principal effective strain rates*/
5073 eps1=(eps_xx+eps_yy)/2.+sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
5074 eps2=(eps_xx+eps_yy)/2.-sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
5075@@ -194,7 +194,7 @@
5076
5077 /*Add input*/
5078 element->AddInput(DamageFEnum,f,element->GetElementType());
5079-
5080+
5081 /*Clean up and return*/
5082 xDelete<IssmDouble>(f);
5083 delete gauss;
5084@@ -261,7 +261,7 @@
5085 Gauss* gauss=element->NewGauss();
5086 for (int i=0;i<numnodes;i++){
5087 gauss->GaussNode(element->GetElementType(),i);
5088-
5089+
5090 damage_input->GetInputValue(&damage,gauss);
5091 tau_xx_input->GetInputValue(&tau_xx,gauss);
5092 tau_xy_input->GetInputValue(&tau_xy,gauss);
5093@@ -313,7 +313,7 @@
5094 }
5095 /*Add input*/
5096 element->AddInput(DamageFEnum,f,element->GetElementType());
5097-
5098+
5099 /*Clean up and return*/
5100 xDelete<IssmDouble>(f);
5101 delete gauss;
5102@@ -374,7 +374,7 @@
5103
5104 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
5105 element->NodalFunctions(basis,gauss);
5106-
5107+
5108 vx_input->GetInputValue(&vx,gauss);
5109 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
5110 vy_input->GetInputValue(&vy,gauss);
5111@@ -537,7 +537,6 @@
5112 damaged_input = element->GetInput(DamageDEnum); _assert_(damaged_input);
5113 }
5114
5115-
5116 /* Start looping on the number of gaussian points: */
5117 Gauss* gauss=element->NewGauss(2);
5118 for(int ig=gauss->begin();ig<gauss->end();ig++){
5119Index: ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp
5120===================================================================
5121--- ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 23065)
5122+++ ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp (revision 23066)
5123@@ -7,7 +7,6 @@
5124 /*Model processing*/
5125 void Balancethickness2Analysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
5126
5127-
5128 int finiteelement = P1Enum;
5129 IoModelToConstraintsx(constraints,iomodel,"md.balancethickness.spcthickness",Balancethickness2AnalysisEnum,finiteelement);
5130
5131Index: ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
5132===================================================================
5133--- ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 23065)
5134+++ ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 23066)
5135@@ -102,7 +102,7 @@
5136 bool dakota_analysis,ismovingfront,isenthalpy;
5137 int frictionlaw,basalforcing_model,materialstype;
5138 int FrictionCoupling;
5139-
5140+
5141 /*Now, is the model 3d? otherwise, do nothing: */
5142 if(iomodel->domaintype==Domain2DhorizontalEnum)return;
5143
5144@@ -188,7 +188,7 @@
5145 default:
5146 _error_("not supported");
5147 }
5148-
5149+
5150 /*Friction law variables*/
5151 switch(frictionlaw){
5152 case 1:
5153@@ -804,7 +804,7 @@
5154
5155 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
5156 element->NodalFunctions(basis,gauss);
5157-
5158+
5159 /*viscous dissipation*/
5160 element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
5161
5162@@ -822,7 +822,7 @@
5163 enthalpypicard_input->GetInputDerivativeValue(&d1enthalpypicard[0],xyz_list,gauss);
5164 pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
5165 d2pressure=0.; // for linear elements, 2nd derivative is zero
5166-
5167+
5168 d1H_d1P=0.;
5169 for(i=0;i<3;i++) d1H_d1P+=d1enthalpypicard[i]*d1pressure[i];
5170 d1P2=0.;
5171@@ -1055,7 +1055,7 @@
5172 IssmDouble dt;
5173 int* basalnodeindices=NULL;
5174 Element* element= NULL;
5175-
5176+
5177 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
5178
5179 /*depth-integrate the drained water fraction */
5180@@ -1387,7 +1387,7 @@
5181 for(int i=0;i<numindices;i++){
5182 gauss->GaussNode(element->GetElementType(),indices[i]);
5183 gaussup->GaussNode(element->GetElementType(),indicesup[i]);
5184-
5185+
5186 enthalpy_input->GetInputValue(&enthalpy,gauss);
5187 enthalpy_input->GetInputValue(&enthalpyup,gaussup);
5188 pressure_input->GetInputValue(&pressure,gauss);
5189Index: ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
5190===================================================================
5191--- ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 23065)
5192+++ ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp (revision 23066)
5193@@ -144,7 +144,7 @@
5194 workelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
5195 GetB(B,workelement,xyz_list,gauss,dim);
5196 GetBprime(Bprime,workelement,xyz_list,gauss,dim);
5197-
5198+
5199 D_scalar=gauss->weight*Jdet;
5200
5201 if(extrapolatebydiffusion){
5202@@ -173,7 +173,7 @@
5203 for(i=0;i<dim;i++) normal[i]=dlsf[i]/norm_dlsf;
5204 else
5205 for(i=0;i<dim;i++) normal[i]=0.;
5206-
5207+
5208 for(row=0;row<dim;row++)
5209 for(col=0;col<dim;col++)
5210 if(row==col)
5211@@ -336,7 +336,7 @@
5212
5213 /* Get parameters */
5214 element->FindParam(&extvar_enum, ExtrapolationVariableEnum);
5215-
5216+
5217 Input* active_input=element->GetInput(IceMaskNodeActivationEnum); _assert_(active_input);
5218 Input* extvar_input=element->GetInput(extvar_enum); _assert_(extvar_input);
5219
5220Index: ../trunk-jpl/src/c/main/issm_slr.cpp
5221===================================================================
5222--- ../trunk-jpl/src/c/main/issm_slr.cpp (revision 23065)
5223+++ ../trunk-jpl/src/c/main/issm_slr.cpp (revision 23066)
5224@@ -27,7 +27,7 @@
5225
5226 /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/
5227 worldcomm=EnvironmentInit(argc,argv);
5228-
5229+
5230 /*What is my rank?:*/
5231 ISSM_MPI_Comm_rank(worldcomm,&my_rank);
5232
5233@@ -39,11 +39,11 @@
5234 rankzeros=xNew<int>(nummodels);
5235 for(int i=0;i<nummodels;i++){
5236 char* string=NULL;
5237-
5238+
5239 string=xNew<char>(strlen(argv[5+3*i])+1);
5240 xMemCpy<char>(string,argv[5+3*i],strlen(argv[5+3*i])+1);
5241 dirnames[i]=string;
5242-
5243+
5244 string=xNew<char>(strlen(argv[5+3*i+1])+1);
5245 xMemCpy<char>(string,argv[5+3*i+1],strlen(argv[5+3*i+1])+1);
5246 modelnames[i]=string;
5247@@ -93,7 +93,7 @@
5248
5249 /*Initialize femmodel from arguments provided command line: */
5250 FemModel *femmodel = new FemModel(4,arguments,modelcomm);
5251-
5252+
5253 /*Now that the models are initialized, keep communicator information in the parameters datasets of each model: */
5254 femmodel->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(worldcomm,WorldCommEnum));
5255 femmodel->parameters->AddObject(new IntParam(NumModelsEnum,nummodels));
5256Index: ../trunk-jpl/src/c/main/issm_dakota.cpp
5257===================================================================
5258--- ../trunk-jpl/src/c/main/issm_dakota.cpp (revision 23065)
5259+++ ../trunk-jpl/src/c/main/issm_dakota.cpp (revision 23066)
5260@@ -15,7 +15,6 @@
5261
5262 int main(int argc,char **argv){
5263
5264-
5265 #if defined(_HAVE_DAKOTA_) && _DAKOTA_MAJOR_ >= 6
5266
5267 bool parallel=true;
5268@@ -28,11 +27,11 @@
5269
5270 /*Initialize MPI: */
5271 ISSM_MPI_Init(&argc, &argv); // initialize MPI
5272-
5273+
5274 /*Recover file name for dakota input file:*/
5275 dakota_input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.in")+2));
5276 sprintf(dakota_input_file,"%s/%s%s",argv[2],argv[3],".qmu.in");
5277-
5278+
5279 dakota_output_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.out")+2));
5280 sprintf(dakota_output_file,"%s/%s%s",argv[2],argv[3],".qmu.out");
5281
5282@@ -56,7 +55,7 @@
5283 << std::endl;
5284 Dakota::abort_handler(-1);
5285 }
5286-
5287+
5288 Dakota::ProblemDescDB& problem_db = env.problem_description_db();
5289 Dakota::ModelLIter ml_iter;
5290 size_t model_index = problem_db.get_db_model_node(); // for restoration
5291Index: ../trunk-jpl/src/c/main/issm_ocean.cpp
5292===================================================================
5293--- ../trunk-jpl/src/c/main/issm_ocean.cpp (revision 23065)
5294+++ ../trunk-jpl/src/c/main/issm_ocean.cpp (revision 23066)
5295@@ -21,7 +21,7 @@
5296
5297 /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/
5298 worldcomm=EnvironmentInit(argc,argv);
5299-
5300+
5301 /*What is my rank?:*/
5302 ISSM_MPI_Comm_rank(worldcomm,&my_rank);
5303 ISSM_MPI_Comm_size(worldcomm,&my_size);
5304@@ -38,7 +38,7 @@
5305 ISSM_MPI_Intercomm_create( modelcomm, 0, worldcomm, my_local_size, 0, &tomitgcmcomm);
5306
5307 FemModel *femmodel = new FemModel(argc,argv,modelcomm);
5308-
5309+
5310 /*Now that the models are initialized, keep communicator information in the parameters datasets of each model: */
5311 femmodel->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(worldcomm,WorldCommEnum));
5312 femmodel->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(tomitgcmcomm,ToMITgcmCommEnum));
5313Index: ../trunk-jpl/src/c/main/esmfbinders.cpp
5314===================================================================
5315--- ../trunk-jpl/src/c/main/esmfbinders.cpp (revision 23065)
5316+++ ../trunk-jpl/src/c/main/esmfbinders.cpp (revision 23066)
5317@@ -2,7 +2,6 @@
5318 * \brief: ESMF binders for ISSM. Binders developed initially for the GEOS-5 framework.
5319 */
5320
5321-
5322 #include "./issm.h"
5323
5324 /*GCM specific declarations:*/
5325@@ -35,7 +34,7 @@
5326
5327 /*Some specific code here for the binding: */
5328 femmodel->parameters->SetParam(SMBgcmEnum,SmbEnum); //bypass SMB model, will be provided by GCM!
5329-
5330+
5331 /*Restart file: */
5332 femmodel->Restart();
5333
5334@@ -112,11 +111,11 @@
5335 {
5336
5337 IssmDouble surface;
5338-
5339+
5340 /*Recover surface from the ISSM element: */
5341 Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
5342 surface_input->GetInputAverage(&surface);
5343-
5344+
5345 *(issmoutputs+f*numberofelements+i) = surface;
5346
5347 }
5348@@ -146,7 +145,7 @@
5349
5350 /*Output results: */
5351 OutputResultsx(femmodel);
5352-
5353+
5354 /*Check point: */
5355 femmodel->CheckPoint();
5356
5357Index: ../trunk-jpl/src/c/bamg/BamgQuadtree.cpp
5358===================================================================
5359--- ../trunk-jpl/src/c/bamg/BamgQuadtree.cpp (revision 23065)
5360+++ ../trunk-jpl/src/c/bamg/BamgQuadtree.cpp (revision 23066)
5361@@ -284,7 +284,6 @@
5362 int xiv = b->v[k]->i.x;
5363 int yiv = b->v[k]->i.y;
5364
5365-
5366 int h0 = Norm(xi2,xiv,yi2,yiv);
5367
5368 /*is it smaller than previous value*/
5369Index: ../trunk-jpl/src/c/bamg/Mesh.cpp
5370===================================================================
5371--- ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 23065)
5372+++ ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 23066)
5373@@ -500,7 +500,7 @@
5374 int* connectivitysize_1=NULL;
5375 int connectivitymax_1=0;
5376 int verbose=0;
5377-
5378+
5379 /*Check needed pointer*/
5380 _assert_(bamgmesh);
5381
5382@@ -1033,7 +1033,7 @@
5383 // ---------------- !
5384 // s0 tt2 s1
5385 //-------------------------------
5386-
5387+
5388 /*Intermediaries*/
5389 Triangle* tt[3]; //the three triangles
5390 long long det3local[3]; //three determinants (integer)
5391@@ -1157,10 +1157,10 @@
5392
5393 int verbose=0;
5394 double lminaniso = 1./ (Max(hminaniso*hminaniso,1e-100));
5395-
5396+
5397 //Get options
5398 if(bamgopts) verbose=bamgopts->verbose;
5399-
5400+
5401 //display info
5402 if (verbose > 1) _printf_(" BoundAnisotropy by " << anisomax << "\n");
5403
5404@@ -1853,7 +1853,7 @@
5405
5406 /*Check pointer*/
5407 _assert_(bamgopts);
5408-
5409+
5410 /*Recover options*/
5411 verbose=bamgopts->verbose;
5412 NbJacobi=bamgopts->nbjacobi;
5413@@ -2727,7 +2727,7 @@
5414 long k0=this->RandomNumber(nbv);
5415 if (verbose>4) _printf_(" Prime Number = "<<PrimeNumber<<"\n");
5416 if (verbose>4) _printf_(" k0 = "<<k0<<"\n");
5417-
5418+
5419 //Build orderedvertices
5420 for (i=0; i<nbv; i++){
5421 orderedvertices[i]=&vertices[k0=(k0+PrimeNumber)%nbv];
5422@@ -2952,7 +2952,7 @@
5423 /*}}}*/
5424 void Mesh::MaxSubDivision(BamgOpts* bamgopts,double maxsubdiv) {/*{{{*/
5425 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)*/
5426-
5427+
5428 /*Intermediaries*/
5429 int verbose=0;
5430
5431@@ -3739,10 +3739,10 @@
5432 /*}}}*/
5433 void Mesh::SmoothingVertex(BamgOpts* bamgopts,int nbiter,double omega ) { /*{{{*/
5434 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/
5435-
5436+
5437 /*Intermediaries*/
5438 int verbose=0;
5439-
5440+
5441 /*Get options*/
5442 if(bamgopts) verbose=bamgopts->verbose;
5443
5444@@ -3784,7 +3784,7 @@
5445 /*}}}*/
5446 void Mesh::SmoothMetric(BamgOpts* bamgopts,double raisonmax) { /*{{{*/
5447 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)*/
5448-
5449+
5450 /*Intermediaries*/
5451 int verbose=0;
5452
5453@@ -4086,7 +4086,7 @@
5454 int verbose=0;
5455 int i;
5456 Metric M1(1);
5457-
5458+
5459 /*Get options*/
5460 if(bamgopts) verbose=bamgopts->verbose;
5461
5462Index: ../trunk-jpl/src/c/shared/Elements/Paterson.cpp
5463===================================================================
5464--- ../trunk-jpl/src/c/shared/Elements/Paterson.cpp (revision 23065)
5465+++ ../trunk-jpl/src/c/shared/Elements/Paterson.cpp (revision 23066)
5466@@ -15,7 +15,6 @@
5467 /*Switch to celsius from Kelvin: */
5468 T=temperature-273.15;
5469
5470-
5471 if(T<=-45.0){
5472 B=1.e+8*(-0.000292866376675*pow(T+50.,3)+ 0.011672640664130*pow(T+50.,2) -0.325004442485481*(T+50.)+ 6.524779401948101);
5473 }
5474Index: ../trunk-jpl/src/c/shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp
5475===================================================================
5476--- ../trunk-jpl/src/c/shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp (revision 23065)
5477+++ ../trunk-jpl/src/c/shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp (revision 23066)
5478@@ -15,7 +15,6 @@
5479 IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
5480 IssmDouble tdiffh;
5481
5482-
5483 for (int imonth = 0; imonth<12; imonth++){
5484 tdiffh = TdiffTime*( TemperaturesLgm[imonth] - TemperaturesPresentday[imonth] );
5485 monthlytemperaturestmp[imonth] = tdiffh + TemperaturesPresentday[imonth] ;
5486@@ -28,4 +27,3 @@
5487 }
5488 // printf(" tempera %f\n",monthlytemperaturestmp[1]);
5489 }
5490-
5491Index: ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp
5492===================================================================
5493--- ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 23065)
5494+++ ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp (revision 23066)
5495@@ -65,7 +65,7 @@
5496
5497 /*Get norm of epsprime*/
5498 epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]);
5499-
5500+
5501 /*Assign output pointers*/
5502 *pepsprime_norm=epsprime_norm;
5503 }/*}}}*/
5504Index: ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp
5505===================================================================
5506--- ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp (revision 23065)
5507+++ ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp (revision 23066)
5508@@ -10,18 +10,18 @@
5509 bool isPrecipScaled, IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
5510 IssmDouble* PrecipitationReconstructed,IssmDouble* TemperatureReconstructed, IssmDouble* monthlytemperaturesout,
5511 IssmDouble* monthlyprecout){
5512-
5513+
5514 IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
5515 IssmDouble deltaTemp;
5516-
5517+
5518 /* Constants */
5519 // dpermil = 2.4;/*degrees C per mil*/
5520-
5521+
5522 /*Create Delta Temp to be applied to monthly temps and used in precip scaling*/
5523 deltaTemp = dpermil * (d018+34.83);
5524-
5525+
5526 for (int imonth = 0; imonth<12; imonth++){
5527-
5528+
5529 if (isTemperatureScaled)monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth] + deltaTemp;
5530 else{
5531 monthlytemperaturestmp[imonth] = TemperatureReconstructed[imonth];
5532@@ -30,7 +30,7 @@
5533
5534 if (isPrecipScaled)monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp);
5535 else monthlyprectmp[imonth] = PrecipitationReconstructed[imonth];
5536-
5537+
5538 /*Assign output pointer*/
5539 *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth];
5540 *(monthlyprecout+imonth) = monthlyprectmp[imonth];
5541Index: ../trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp
5542===================================================================
5543--- ../trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp (revision 23065)
5544+++ ../trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp (revision 23066)
5545@@ -18,7 +18,7 @@
5546 /*get drainage function value*/
5547 if((w0==w1)||(w1==w2)||(w0==w2))
5548 _error_("Error: equal ordinates in DrainageFunctionWaterfraction -> division by zero. Abort");
5549-
5550+
5551 if(waterfraction<=w0)
5552 Dret=D0;
5553 else if((waterfraction>w0) && (waterfraction<=w1))
5554Index: ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp
5555===================================================================
5556--- ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 23065)
5557+++ ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 23066)
5558@@ -82,11 +82,11 @@
5559 // seasonal loop
5560 for (iqj = 0; iqj < 12; iqj++){
5561 imonth = ismon[iqj];
5562-
5563+
5564 /*********compute lapse rate ****************/
5565 st=(s-s0t)/1000.;
5566 rtlaps=TdiffTime*rlapslgm + (1.-TdiffTime)*rlaps; // lapse rate
5567-
5568+
5569 /*********compute Surface temperature *******/
5570 monthlytemperatures[imonth]=monthlytemperatures[imonth] - rtlaps *max(st,sealevTime*0.001);
5571 tstar = monthlytemperatures[imonth];
5572@@ -98,12 +98,12 @@
5573 if (tstar >= -siglimc){ pd = pds[reCast<int,IssmDouble>(tstar/DT + siglim0c)];}}
5574 else {
5575 pd = 0.;}
5576-
5577+
5578 /******exp des/elev precip reduction*******/
5579 sp=(s-s0p)/1000.-deselcut; // deselev effect is wrt chng in topo
5580 if (sp>0.0){q = exp(-desfac*sp);}
5581 else {q = 1.0;}
5582-
5583+
5584 qmt= qmt + monthlyprec[imonth]*sconv; //*sconv to convert in m of ice equivalent per month
5585 qmpt= q*monthlyprec[imonth]*sconv;
5586 qmp= qmp + qmpt;
5587@@ -113,7 +113,7 @@
5588 // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
5589 // gaussian=T_m, so ndd=-(Tsurf-pdd)
5590 if (iqj>5 && iqj<9){ Tsum=Tsum+tstar;}
5591-
5592+
5593 if (tstar >= siglim) {pdd = pdd + tstar*deltm;}
5594 else if (tstar> -siglim){
5595 pddsig=pdds[reCast<int,IssmDouble>(tstar/DT + siglim0)];
5596@@ -120,7 +120,7 @@
5597 pdd = pdd + pddsig*deltm;
5598 frzndd = frzndd - (tstar-pddsig)*deltm;}
5599 else{frzndd = frzndd - tstar*deltm; }
5600-
5601+
5602 /*Assign output pointer*/
5603 *(monthlytemperatures+imonth) = monthlytemperatures[imonth];
5604 *(monthlyprec+imonth) = monthlyprec[imonth];
5605@@ -149,7 +149,7 @@
5606 icefac=pddicefac;
5607 }
5608 }
5609-
5610+
5611 /***** determine PDD factors *****/
5612 if(Tsum< -1.) {
5613 snwmf=(2.65+snowfac-pddsnowfac0)*0.001; // ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
5614@@ -165,12 +165,12 @@
5615 }
5616 snwmf=0.95*snwmf;
5617 smf=0.95*smf;
5618-
5619+
5620 /***** compute PDD ablation and refreezing *****/
5621 pddt = pdd *365;
5622 snwm = snwmf*pddt; // snow that could have been melted in a year
5623 hmx2 = min(h,dfrz); // refreeze active layer max depth: dfrz
5624-
5625+
5626 if(snwm < saccu) {
5627 water=prect-saccu + snwm; //water=rain + snowmelt
5628 // l 2.2= capillary factor
5629Index: ../trunk-jpl/src/c/shared/Numerics/legendre.cpp
5630===================================================================
5631--- ../trunk-jpl/src/c/shared/Numerics/legendre.cpp (revision 23065)
5632+++ ../trunk-jpl/src/c/shared/Numerics/legendre.cpp (revision 23066)
5633@@ -120,11 +120,11 @@
5634 v(1:m,2) = x;
5635
5636 for i = 2 : n
5637-
5638+
5639 v(:,i+1) = ( ( 2 * i - 1 ) * x .* v(:,i) ...
5640 - ( i - 1 ) * v(:,i-1) ) ...
5641 / ( i );
5642-
5643+
5644 end
5645 }}} */
5646
5647@@ -239,7 +239,7 @@
5648 Output, double P_POLYNOMIAL_VALUE[M*(N+1)], the values of the Legendre
5649 polynomials of order 0 through N.
5650 }}}*/
5651-
5652+
5653 int i,j;
5654
5655 if(n<0) return NULL;
5656@@ -253,7 +253,7 @@
5657
5658 for ( i = 0; i < m; i++ ) {
5659 for ( j = 2; j <= n; j++ ) {
5660-
5661+
5662 v[j+i*(n+1)] = ( ( IssmDouble ) ( 2 * j - 1 ) * x[i] * v[(j-1)+i*(n+1)]
5663 - ( IssmDouble ) ( j - 1 ) * v[(j-2)+i*(n+1)] )
5664 / ( IssmDouble ) ( j );
5665Index: ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp
5666===================================================================
5667--- ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp (revision 23065)
5668+++ ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp (revision 23066)
5669@@ -67,7 +67,7 @@
5670 /*Get current Gradient at xmin=0*/
5671 _printf0_(" x = "<<setw(9)<<xmin<<" | ");
5672 fxmin = (*g)(&G,X0,usr); if(xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
5673-
5674+
5675 /*Get f(xmax)*/
5676 _printf0_(" x = "<<setw(9)<<xmax<<" | ");
5677 for(int i=0;i<nsize;i++) X[i]=X0[i]+xmax*G[i];
5678@@ -243,7 +243,7 @@
5679 xDelete<IssmDouble>(G);
5680 J[n]=fxbest;
5681 }
5682-
5683+
5684 /*return*/
5685 xDelete<IssmDouble>(X);
5686 *pJ=J;
5687Index: ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
5688===================================================================
5689--- ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 23065)
5690+++ ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 23066)
5691@@ -62,7 +62,7 @@
5692 else{
5693 dtemp = &dtemp_static[0];
5694 }
5695-
5696+
5697 /* perform the matrix triple product in the order that minimizes the
5698 number of multiplies and the temporary space used, noting that
5699 (a*b)*c requires ac(b+d) multiplies and ac IssmDoubles, and a*(b*c)
5700@@ -334,7 +334,7 @@
5701 /*Compute determinant*/
5702 Matrix2x2Determinant(&det,A);
5703 if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
5704-
5705+
5706 /*Multiplication is faster than divsion, so we multiply by the reciprocal*/
5707 det_reciprocal = 1./det;
5708
5709@@ -426,7 +426,6 @@
5710 *pvx = vx;
5711 *pvy = vy;
5712
5713-
5714 }/*}}}*/
5715
5716 void Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A){/*{{{*/
5717@@ -576,16 +575,16 @@
5718 }/*}}}*/
5719
5720 void newcell(IssmDouble** pcell, IssmDouble newvalue, bool top, int m){ /*{{{*/
5721-
5722+
5723 IssmDouble* cell=NULL;
5724 IssmDouble* dummy=NULL;
5725-
5726+
5727 /*recover pointer: */
5728 cell=*pcell;
5729-
5730+
5731 /*reallocate:*/
5732 dummy=xNew<IssmDouble>(m+1);
5733-
5734+
5735 /*copy data:*/
5736 if(top){
5737 dummy[0]=newvalue;
5738@@ -595,10 +594,10 @@
5739 dummy[m]=newvalue;
5740 for(int i=0;i<m;i++)dummy[i]=cell[i];
5741 }
5742-
5743+
5744 /*delete and reassign: */
5745 xDelete<IssmDouble>(cell); cell=dummy;
5746-
5747+
5748 /*assign output pointer:*/
5749 *pcell=cell;
5750 } /*}}}*/
5751@@ -614,7 +613,7 @@
5752
5753 /*input: */
5754 IssmDouble* cell=*pcell;
5755-
5756+
5757 /*output: */
5758 IssmDouble* newcell=xNew<IssmDouble>(nind);
5759
5760@@ -621,7 +620,7 @@
5761 for(int i=0;i<nind;i++){
5762 newcell[i]=cell[indices[i]];
5763 }
5764-
5765+
5766 /*free allocation:*/
5767 xDelete<IssmDouble>(cell);
5768
5769@@ -632,7 +631,7 @@
5770
5771 /*input: */
5772 IssmDouble* cell=*pcell;
5773-
5774+
5775 /*output: */
5776 IssmDouble* newcell=xNew<IssmDouble>(m+1);
5777
5778@@ -640,7 +639,7 @@
5779 newcell[i]=scale*cell[i];
5780 newcell[i+1]=scale* cell[i];
5781 for(int j=i+2;j<m+1;j++)newcell[j]=cell[j-1];
5782-
5783+
5784 /*free allocation:*/
5785 xDelete<IssmDouble>(cell);
5786
5787@@ -661,7 +660,7 @@
5788 celllist[x]= va_arg ( arguments, IssmDouble*);
5789 }
5790 va_end ( arguments );
5791-
5792+
5793 _printf_("Echo of cell: \n");
5794 for(int i=0;i<m;i++){
5795 _printf_(i << ": ");
5796Index: ../trunk-jpl/src/c/classes/Materials/Matpar.cpp
5797===================================================================
5798--- ../trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23065)
5799+++ ../trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23066)
5800@@ -58,7 +58,7 @@
5801 lithosphere_density = 0;
5802 mantle_shear_modulus = 0;
5803 mantle_density = 0;
5804-
5805+
5806 earth_density = 0;
5807
5808 int nnat,dummy;
5809@@ -337,7 +337,7 @@
5810 matpar->lithosphere_density=this->lithosphere_density;
5811 matpar->mantle_shear_modulus=this->mantle_shear_modulus;
5812 matpar->mantle_density=this->mantle_density;
5813-
5814+
5815 matpar->earth_density=this->earth_density;
5816
5817 return matpar;
5818@@ -438,7 +438,7 @@
5819 MARSHALLING(lithosphere_density);
5820 MARSHALLING(mantle_shear_modulus);
5821 MARSHALLING(mantle_density);
5822-
5823+
5824 //slr:
5825 MARSHALLING(earth_density);
5826 }
5827Index: ../trunk-jpl/src/c/classes/Loads/Friction.cpp
5828===================================================================
5829--- ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23065)
5830+++ ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23066)
5831@@ -195,7 +195,7 @@
5832 /*Check to prevent dividing by zero if vmag==0*/
5833 if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
5834 else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
5835-
5836+
5837 _assert_(!xIsNan<IssmDouble>(alpha_complement));
5838 _assert_(!xIsInf<IssmDouble>(alpha_complement));
5839
5840@@ -270,7 +270,6 @@
5841 r=drag_q/drag_p;
5842 s=1./drag_p;
5843
5844-
5845 /*Get effective pressure*/
5846 IssmDouble Neff = EffectivePressure(gauss);
5847
5848Index: ../trunk-jpl/src/c/classes/FemModel.cpp
5849===================================================================
5850--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23065)
5851+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23066)
5852@@ -242,7 +242,6 @@
5853 char *lockfilename = NULL;
5854 bool waitonlock = false;
5855
5856-
5857 /*Write lock file if requested: */
5858 this->parameters->FindParam(&waitonlock,SettingsWaitonlockEnum);
5859 this->parameters->FindParam(&lockfilename,LockFileNameEnum);
5860@@ -4392,7 +4391,6 @@
5861 /*Free ressources:*/
5862 xDelete<IssmDouble>(RSLg_old);
5863
5864-
5865 }
5866 /*}}}*/
5867 void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
5868@@ -4806,7 +4804,6 @@
5869 char** poutputbuffer;
5870 size_t* poutputbuffersize;
5871
5872-
5873 /*Before we delete the profiler, report statistics for this run: */
5874 profiler->Stop(TOTAL); //final tagging
5875 _printf0_("\n");
5876@@ -4868,7 +4865,6 @@
5877 }/*}}}*/
5878 #endif
5879
5880-
5881 #if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_)
5882 void FemModel::ReMeshBamg(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/
5883
5884@@ -5245,7 +5241,7 @@
5885 newx[i] = newxylist[2*i];
5886 newy[i] = newxylist[2*i+1];
5887 }
5888-
5889+
5890 /*Assign the pointers*/
5891 (*pnewelementslist) = newelementslist; //Matlab indexing
5892 (*pnewx) = newx;
5893Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
5894===================================================================
5895--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23065)
5896+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23066)
5897@@ -2735,7 +2735,6 @@
5898 }
5899 }
5900
5901-
5902 /*Clean-up*/
5903 xDelete<IssmDouble>(dbasis);
5904 }/*}}}*/
Note: See TracBrowser for help on using the repository browser.