source: issm/oecreview/Archive/24684-25833/ISSM-24939-24940.diff@ 25834

Last change on this file since 25834 was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 52.0 KB
RevLine 
[25834]1Index: ../trunk-jpl/src/m/classes/slr.m
2===================================================================
3--- ../trunk-jpl/src/m/classes/slr.m (revision 24939)
4+++ ../trunk-jpl/src/m/classes/slr.m (revision 24940)
5@@ -32,6 +32,7 @@
6 horiz = 0;
7 Ngia = NaN;
8 Ugia = NaN;
9+ eartharea = 4*pi*planetradius('earth')^2;
10
11 requested_outputs = {};
12 transitions = {};
13@@ -101,7 +102,8 @@
14 return;
15 end
16
17- md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
18+ %md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
19+ md = checkfield(md,'fieldname','slr.deltathickness','timeseries',1,'NaN',1,'Inf',1);
20 md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
21 md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1);
22 md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
23@@ -131,12 +133,12 @@
24 end
25
26 %cross check that whereever we have an ice load, the mask is <0 on each vertex:
27- pos=find(self.deltathickness);
28- maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:));
29- [els,vertices]=find(maskpos>0);
30- if length(els),
31- warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
32- end
33+ %pos=find(self.deltathickness);
34+ %maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:));
35+ %[els,vertices]=find(maskpos>0);
36+ %if length(els),
37+ % warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
38+ %end
39
40 %check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not,
41 %a coupler to a planet model is provided.
42@@ -220,6 +222,7 @@
43 WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer');
44 WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer');
45 WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer');
46+ WriteData(fid,prefix,'object',self,'fieldname','eartharea','format','Double');
47
48 %process requested outputs
49 outputs = self.requested_outputs;
50Index: ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp
51===================================================================
52--- ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24939)
53+++ ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24940)
54@@ -320,6 +320,7 @@
55 parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
56 parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
57 parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum));
58+ parameters->AddObject(iomodel->CopyConstantObject("md.slr.eartharea",SealevelEarthAreaEnum));
59
60 /*Deal with dsl multi-model ensembles: {{{*/
61 iomodel->FetchData(&dslmodel,"md.dsl.model");
62Index: ../trunk-jpl/src/c/shared/Enum/Enum.vim
63===================================================================
64--- ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24939)
65+++ ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24940)
66@@ -301,6 +301,7 @@
67 syn keyword cConstant ResultsEnum
68 syn keyword cConstant RootPathEnum
69 syn keyword cConstant SaveResultsEnum
70+syn keyword cConstant SealevelEarthAreaEnum
71 syn keyword cConstant SealevelEustaticEnum
72 syn keyword cConstant SealevelriseAbstolEnum
73 syn keyword cConstant SealevelriseAngularVelocityEnum
74@@ -1430,6 +1431,7 @@
75 syn keyword cType Results
76 syn keyword cType RiftStruct
77 syn keyword cType Riftfront
78+syn keyword cType SealevelMasks
79 syn keyword cType Seg
80 syn keyword cType SegInput2
81 syn keyword cType SegRef
82Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
83===================================================================
84--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24939)
85+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24940)
86@@ -295,6 +295,7 @@
87 ResultsEnum,
88 RootPathEnum,
89 SaveResultsEnum,
90+ SealevelEarthAreaEnum,
91 SealevelEustaticEnum,
92 SealevelriseAbstolEnum,
93 SealevelriseAngularVelocityEnum,
94Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
95===================================================================
96--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24939)
97+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24940)
98@@ -303,6 +303,7 @@
99 case ResultsEnum : return "Results";
100 case RootPathEnum : return "RootPath";
101 case SaveResultsEnum : return "SaveResults";
102+ case SealevelEarthAreaEnum : return "SealevelEarthArea";
103 case SealevelEustaticEnum : return "SealevelEustatic";
104 case SealevelriseAbstolEnum : return "SealevelriseAbstol";
105 case SealevelriseAngularVelocityEnum : return "SealevelriseAngularVelocity";
106Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
107===================================================================
108--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24939)
109+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24940)
110@@ -309,6 +309,7 @@
111 else if (strcmp(name,"Results")==0) return ResultsEnum;
112 else if (strcmp(name,"RootPath")==0) return RootPathEnum;
113 else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
114+ else if (strcmp(name,"SealevelEarthArea")==0) return SealevelEarthAreaEnum;
115 else if (strcmp(name,"SealevelEustatic")==0) return SealevelEustaticEnum;
116 else if (strcmp(name,"SealevelriseAbstol")==0) return SealevelriseAbstolEnum;
117 else if (strcmp(name,"SealevelriseAngularVelocity")==0) return SealevelriseAngularVelocityEnum;
118@@ -381,11 +382,11 @@
119 else if (strcmp(name,"SmbK")==0) return SmbKEnum;
120 else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
121 else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum;
122- else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum;
123 else stage=4;
124 }
125 if(stage==4){
126- if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
127+ if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum;
128+ else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
129 else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum;
130 else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum;
131 else if (strcmp(name,"SmbRunoffalti")==0) return SmbRunoffaltiEnum;
132@@ -504,11 +505,11 @@
133 else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
134 else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
135 else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
136- else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
137 else stage=5;
138 }
139 if(stage==5){
140- if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
141+ if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
142+ else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
143 else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
144 else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
145 else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
146@@ -627,11 +628,11 @@
147 else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
148 else if (strcmp(name,"Ice")==0) return IceEnum;
149 else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
150- else if (strcmp(name,"Input")==0) return InputEnum;
151 else stage=6;
152 }
153 if(stage==6){
154- if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
155+ if (strcmp(name,"Input")==0) return InputEnum;
156+ else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
157 else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum;
158 else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum;
159 else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
160@@ -750,11 +751,11 @@
161 else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
162 else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
163 else if (strcmp(name,"SmbMassBalanceClimate")==0) return SmbMassBalanceClimateEnum;
164- else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
165 else stage=7;
166 }
167 if(stage==7){
168- if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum;
169+ if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
170+ else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum;
171 else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum;
172 else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum;
173 else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
174@@ -873,11 +874,11 @@
175 else if (strcmp(name,"Outputdefinition10")==0) return Outputdefinition10Enum;
176 else if (strcmp(name,"Outputdefinition11")==0) return Outputdefinition11Enum;
177 else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum;
178- else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
179 else stage=8;
180 }
181 if(stage==8){
182- if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
183+ if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
184+ else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
185 else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
186 else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
187 else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
188@@ -996,11 +997,11 @@
189 else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
190 else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum;
191 else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum;
192- else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
193 else stage=9;
194 }
195 if(stage==9){
196- if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
197+ if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
198+ else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
199 else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
200 else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
201 else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum;
202@@ -1119,11 +1120,11 @@
203 else if (strcmp(name,"Hydrologyshakti")==0) return HydrologyshaktiEnum;
204 else if (strcmp(name,"Hydrologyshreve")==0) return HydrologyshreveEnum;
205 else if (strcmp(name,"IceMass")==0) return IceMassEnum;
206- else if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum;
207 else stage=10;
208 }
209 if(stage==10){
210- if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
211+ if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum;
212+ else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
213 else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum;
214 else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
215 else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum;
216@@ -1242,11 +1243,11 @@
217 else if (strcmp(name,"Paterson")==0) return PatersonEnum;
218 else if (strcmp(name,"Pengrid")==0) return PengridEnum;
219 else if (strcmp(name,"Penpair")==0) return PenpairEnum;
220- else if (strcmp(name,"Penta")==0) return PentaEnum;
221 else stage=11;
222 }
223 if(stage==11){
224- if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
225+ if (strcmp(name,"Penta")==0) return PentaEnum;
226+ else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
227 else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
228 else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
229 else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
230Index: ../trunk-jpl/src/c/cores/cores.h
231===================================================================
232--- ../trunk-jpl/src/c/cores/cores.h (revision 24939)
233+++ ../trunk-jpl/src/c/cores/cores.h (revision 24940)
234@@ -57,9 +57,9 @@
235 void steric_core(FemModel* femmodel);
236 void sealevelrise_core_geometry(FemModel* femmodel);
237 SealevelMasks* sealevelrise_core_masks(FemModel* femmodel);
238-Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* peartharea,IssmDouble* poceanarea);
239-Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea);
240-void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg,IssmDouble eartharea, SealevelMasks* masks);
241+Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea);
242+Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble oceanarea);
243+void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks);
244 void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia,Vector<IssmDouble>** pN_gia,FemModel* femmodel,Vector<IssmDouble>* RSLg);
245 void sealevelrise_diagnostics(FemModel* femmodel,Vector<IssmDouble>* RSLg);
246 IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);
247Index: ../trunk-jpl/src/c/cores/sealevelrise_core.cpp
248===================================================================
249--- ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24939)
250+++ ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24940)
251@@ -95,7 +95,7 @@
252 int horiz;
253 int geodetic=0;
254 IssmDouble dt;
255- IssmDouble eartharea,oceanarea;
256+ IssmDouble oceanarea;
257
258 /*Should we even be here?:*/
259 femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return;
260@@ -152,13 +152,13 @@
261 masks=sealevelrise_core_masks(femmodel);
262
263 /*call eustatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
264- RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&eartharea,&oceanarea);
265+ RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&oceanarea);
266
267 /*call non-eustatic core (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */
268- RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,eartharea,oceanarea);
269+ RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,oceanarea);
270
271 /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */
272- sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,eartharea,masks);
273+ sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,masks);
274
275 /*compute viscosus (GIA) geodetic signatures:*/
276 sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg);
277@@ -353,7 +353,7 @@
278
279
280 }/*}}}*/
281-Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/
282+Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/
283
284 /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
285
286@@ -367,7 +367,7 @@
287 IssmDouble *longitude = NULL;
288 IssmDouble *radius = NULL;
289 int loop;
290- IssmDouble eartharea,oceanarea;
291+ IssmDouble oceanarea;
292
293 /*outputs:*/
294 IssmDouble eustatic;
295@@ -387,7 +387,7 @@
296 RSLgi = new Vector<IssmDouble>(gsize);
297
298 /*call the eustatic main module: */
299- femmodel->SealevelriseEustatic(RSLgi,&eartharea, &oceanarea,&eustatic, masks, latitude, longitude, radius,loop); //this computes
300+ femmodel->SealevelriseEustatic(RSLgi,&oceanarea,&eustatic, masks, latitude, longitude, radius,loop); //this computes
301
302 /*we need to average RSLgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
303 RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,masks, oceanarea);
304@@ -405,11 +405,10 @@
305 xDelete<IssmDouble>(radius);
306
307 /*Assign output pointers and return: */
308- *peartharea=eartharea;
309 *poceanarea=oceanarea;
310 return RSLgi;
311 }/*}}}*/
312-Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea){ /*{{{*/
313+Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble oceanarea){ /*{{{*/
314
315 /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
316 non eustatic core of the SLR solution */
317@@ -478,7 +477,7 @@
318 RSLgo = new Vector<IssmDouble>(gsize); RSLgo->Assemble();
319
320 /*call the non eustatic module: */
321- femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, eartharea, masks, latitude,longitude, radius,verboseconvolution,loop);
322+ femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, masks, latitude,longitude, radius,verboseconvolution,loop);
323
324 /*assemble solution vector: */
325 RSLgo->Assemble();
326@@ -487,7 +486,7 @@
327
328 /*call rotational feedback module: */
329 RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
330- femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, eartharea, masks, latitude,longitude,radius);
331+ femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks, latitude,longitude,radius);
332 RSLgo_rot->Assemble();
333
334 /*save changes in inertia tensor as results: */
335@@ -536,7 +535,7 @@
336
337 return RSLg;
338 } /*}}}*/
339-void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg,IssmDouble eartharea, SealevelMasks* masks){ /*{{{*/
340+void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks){ /*{{{*/
341
342 Vector<IssmDouble> *U_esa = NULL;
343 Vector<IssmDouble> *U_north_esa = NULL;
344@@ -576,7 +575,7 @@
345 VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
346
347 /*call the elastic main modlule:*/
348- femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,eartharea, masks, latitude,longitude,radius,xx,yy,zz,loop,horiz);
349+ femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg, masks, latitude,longitude,radius,xx,yy,zz,loop,horiz);
350
351 /*Assign output pointers:*/
352 *pU_esa=U_esa;
353Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
354===================================================================
355--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 24939)
356+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 24940)
357@@ -208,16 +208,16 @@
358 #endif
359 #ifdef _HAVE_ESA_
360 void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
361- void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
362+ void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
363 #endif
364 #ifdef _HAVE_SEALEVELRISE_
365 IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
366 void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
367- void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");};
368+ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
369 void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
370- void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
371- void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
372- void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
373+ void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");};
374+ void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
375+ void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");};
376 #endif
377
378 /*}}}*/
379Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
380===================================================================
381--- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 24939)
382+++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 24940)
383@@ -168,15 +168,15 @@
384 #endif
385 #ifdef _HAVE_ESA_
386 void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
387- void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
388+ void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
389 #endif
390 #ifdef _HAVE_SEALEVELRISE_
391- void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");};
392+ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
393 void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
394 void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
395- void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
396- void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
397- void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
398+ void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");};
399+ void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
400+ void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");};
401 IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
402 #endif
403
404Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h
405===================================================================
406--- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24939)
407+++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24940)
408@@ -174,15 +174,15 @@
409 #endif
410 #ifdef _HAVE_ESA_
411 void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
412- void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
413+ void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
414 #endif
415 #ifdef _HAVE_SEALEVELRISE_
416 void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
417- void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");};
418+ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
419 void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
420- void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
421- void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
422- void SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
423+ void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");};
424+ void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
425+ void SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");};
426 IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
427 #endif
428
429Index: ../trunk-jpl/src/c/classes/Elements/Element.h
430===================================================================
431--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 24939)
432+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 24940)
433@@ -370,7 +370,7 @@
434 #endif
435 #ifdef _HAVE_ESA_
436 virtual void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy)=0;
437- virtual void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0;
438+ virtual void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz)=0;
439 #endif
440 #ifdef _HAVE_SEALEVELRISE_
441 virtual void SetSealevelMasks(SealevelMasks* masks)=0;
442@@ -377,11 +377,11 @@
443 virtual IssmDouble GetArea3D(void)=0;
444 virtual IssmDouble GetAreaSpherical(void)=0;
445 virtual IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks)=0;
446- virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea)=0;
447- virtual void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
448+ virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0;
449+ virtual void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea)=0;
450 virtual void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0;
451- virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0;
452- virtual void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0;
453+ virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0;
454+ virtual void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz)=0;
455 #endif
456
457 };
458Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
459===================================================================
460--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 24939)
461+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 24940)
462@@ -159,18 +159,18 @@
463 #endif
464 #ifdef _HAVE_ESA_
465 void EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy);
466- void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea);
467+ void EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz);
468 #endif
469 #ifdef _HAVE_SEALEVELRISE_
470 IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks);
471- void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble eartharea);
472+ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
473 void SetSealevelMasks(SealevelMasks* masks);
474 void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius);
475- void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
476- void SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
477- void SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
478- void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea);
479- void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz);
480+ void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea);
481+ void SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea);
482+ void SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea);
483+ void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius);
484+ void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz);
485 #endif
486 /*}}}*/
487 /*Tria specific routines:{{{*/
488Index: ../trunk-jpl/src/c/classes/FemModel.h
489===================================================================
490--- ../trunk-jpl/src/c/classes/FemModel.h (revision 24939)
491+++ ../trunk-jpl/src/c/classes/FemModel.h (revision 24940)
492@@ -162,11 +162,11 @@
493 void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
494 #endif
495 #ifdef _HAVE_SEALEVELRISE_
496- void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
497+ void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
498 void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
499- void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
500- void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea,SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
501- void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble eartharea,SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
502+ void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
503+ void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
504+ void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
505 IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg,SealevelMasks* masks, IssmDouble oceanarea);
506 #endif
507 void HydrologyEPLupdateDomainx(IssmDouble* pEplcount);
508Index: ../trunk-jpl/src/c/classes/FemModel.cpp
509===================================================================
510--- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 24939)
511+++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 24940)
512@@ -4634,9 +4634,6 @@
513 /*}}}*/
514 void FemModel::EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
515
516- IssmDouble eartharea=0;
517- IssmDouble eartharea_cpu=0;
518-
519 int ns,nsmax;
520
521 /*Go through elements, and add contribution from each element to the deflection vector wg:*/
522@@ -4645,10 +4642,7 @@
523 /*First, figure out the surface area of Earth: */
524 for(int i=0;i<ns;i++){
525 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
526- eartharea_cpu += element->GetAreaSpherical();
527 }
528- ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
529- ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
530
531 /*Figure out max of ns: */
532 ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
533@@ -4658,7 +4652,7 @@
534 for(int i=0;i<nsmax;i++){
535 if(i<ns){
536 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
537- element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea);
538+ element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz);
539 }
540 if(i%100==0){
541 pUp->Assemble();
542@@ -4693,7 +4687,7 @@
543
544 }
545 /*}}}*/
546-void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
547+void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
548
549 /*serialized vectors:*/
550 IssmDouble eustatic = 0.;
551@@ -4701,8 +4695,6 @@
552 IssmDouble eustatic_cpu_e = 0.;
553 IssmDouble oceanarea = 0.;
554 IssmDouble oceanarea_cpu = 0.;
555- IssmDouble eartharea = 0.;
556- IssmDouble eartharea_cpu = 0.;
557
558 /*Initialize temporary vector that will be used to sum eustatic components
559 * on all local elements, prior to assembly:*/
560@@ -4715,7 +4707,6 @@
561 for(int i=0;i<elements->Size();i++){
562 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
563 IssmDouble area=element->GetAreaSpherical();
564- eartharea_cpu += area;
565 if (masks->isoceanin[i]) oceanarea_cpu += area;
566 }
567 ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
568@@ -4722,13 +4713,10 @@
569 ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
570 _assert_(oceanarea>0.);
571
572- ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
573- ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
574-
575 /*Call the sea level rise core: */
576 for(int i=0;i<elements->Size();i++){
577 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
578- element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, latitude,longitude,radius,oceanarea,eartharea);
579+ element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, latitude,longitude,radius,oceanarea);
580 eustatic_cpu+=eustatic_cpu_e;
581 }
582
583@@ -4746,13 +4734,12 @@
584 xDelete<IssmDouble>(RSLgi);
585
586 /*Assign output pointers:*/
587- *peartharea = eartharea;
588 *poceanarea = oceanarea;
589 *peustatic = eustatic;
590
591 }
592 /*}}}*/
593-void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
594+void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
595
596 /*serialized vectors:*/
597 IssmDouble* RSLg_old=NULL;
598@@ -4773,7 +4760,7 @@
599 /*Call the sea level rise core: */
600 for(int i=0;i<elements->Size();i++){
601 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
602- element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius,eartharea);
603+ element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius);
604 }
605 pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL);
606 pRSLgo->Assemble();
607@@ -4785,7 +4772,7 @@
608
609 }
610 /*}}}*/
611-void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
612+void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
613
614 /*serialized vectors:*/
615 IssmDouble* RSLg_old=NULL;
616@@ -4801,7 +4788,7 @@
617 IssmDouble moi_list_cpu[3]={0,0,0};
618 for(int i=0;i<elements->Size();i++){
619 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
620- element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks, eartharea);
621+ element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks );
622 moi_list_cpu[0] += moi_list[0];
623 moi_list_cpu[1] += moi_list[1];
624 moi_list_cpu[2] += moi_list[2];
625@@ -4861,7 +4848,7 @@
626
627 }
628 /*}}}*/
629-void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
630+void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
631
632 /*serialized vectors:*/
633 IssmDouble* RSLg=NULL;
634@@ -4887,7 +4874,7 @@
635 /*Call the sea level rise core: */
636 for(int i=0;i<elements->Size();i++){
637 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
638- element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
639+ element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,horiz);
640 }
641
642 pUp->SetValues(gsize,indices,Up,ADD_VAL);
643Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
644===================================================================
645--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24939)
646+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24940)
647@@ -5288,7 +5288,7 @@
648 return;
649 }
650 /*}}}*/
651-void Tria::EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/
652+void Tria::EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){ /*{{{*/
653
654 /*diverse:*/
655 int gsize;
656@@ -5295,7 +5295,7 @@
657 bool spherical=true;
658 IssmDouble llr_list[NUMVERTICES][3];
659 IssmDouble xyz_list[NUMVERTICES][3];
660- IssmDouble area;
661+ IssmDouble area,eartharea;
662 IssmDouble I; //ice/water loading
663 IssmDouble late,longe,re;
664 IssmDouble lati,longi,ri;
665@@ -5327,6 +5327,9 @@
666 /*recover material parameters: */
667 rho_ice=FindParam(MaterialsRhoIceEnum);
668 rho_earth=FindParam(MaterialsEarthDensityEnum);
669+
670+ /*recover earth area: */
671+ this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
672
673 /*how many dofs are we working with here? */
674 this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
675@@ -5464,7 +5467,7 @@
676
677 }
678 /*}}}*/
679-void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks, IssmDouble eartharea){/*{{{*/
680+void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
681 /*early return if we are not on an ice cap OR ocean:*/
682 if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){
683 dI_list[0] = 0.0; // this is important!!!
684@@ -5474,9 +5477,12 @@
685 }
686
687 /*Compute area of element:*/
688- IssmDouble area;
689+ IssmDouble area,eartharea;
690 area=GetAreaSpherical();
691
692+ /*recover earth area: */
693+ this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
694+
695 /*Compute lat,long,radius of elemental centroid: */
696 bool spherical=true;
697 IssmDouble llr_list[NUMVERTICES][3];
698@@ -5660,7 +5666,7 @@
699 return;
700 }
701 /*}}}*/
702-void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
703+void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){ /*{{{*/
704
705 /*Computational flags:*/
706 int bp_compute_fingerprints= 0;
707@@ -5671,22 +5677,22 @@
708 if(!masks->isoceanin[this->lid]){
709 /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested
710 *bottom pressure fingerprints:*/
711- if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
712+ if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea);
713 }
714 //if(!IsIceInElement()){
715 /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/
716- this->SealevelriseEustaticIce(Sgi,peustatic,masks, latitude,longitude,radius,oceanarea,eartharea);
717+ this->SealevelriseEustaticIce(Sgi,peustatic,masks, latitude,longitude,radius,oceanarea);
718 //}
719
720 }
721 /*}}}*/
722-void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
723+void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){ /*{{{*/
724
725 /*diverse:*/
726 int gsize,dummy;
727 bool spherical=true;
728 IssmDouble llr_list[NUMVERTICES][3];
729- IssmDouble area;
730+ IssmDouble area,eartharea;
731 IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
732 IssmDouble rho;
733 IssmDouble late,longe,re;
734@@ -5744,6 +5750,9 @@
735 this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
736 this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum);
737
738+ /*recover earth area: */
739+ this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
740+
741 /*recover precomputed green function kernels:*/
742 DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
743 parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M);
744@@ -5835,13 +5844,13 @@
745 return;
746 }
747 /*}}}*/
748-void Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
749+void Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){ /*{{{*/
750
751 /*diverse:*/
752 int gsize;
753 bool spherical=true;
754 IssmDouble llr_list[NUMVERTICES][3];
755- IssmDouble area;
756+ IssmDouble area,eartharea;
757 IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
758 IssmDouble rho;
759 IssmDouble late,longe,re;
760@@ -5875,6 +5884,9 @@
761 rho_water=FindParam(MaterialsRhoSeawaterEnum);
762 rho_earth=FindParam(MaterialsEarthDensityEnum);
763
764+ /*recover earth area: */
765+ this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
766+
767 /*recover love numbers and computational flags: */
768 this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
769 this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
770@@ -5981,13 +5993,13 @@
771 return;
772 }
773 /*}}}*/
774-void Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/
775+void Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){ /*{{{*/
776
777 /*diverse:*/
778 int gsize,dummy;
779 bool spherical=true;
780 IssmDouble llr_list[NUMVERTICES][3];
781- IssmDouble area;
782+ IssmDouble area,eartharea;
783 IssmDouble S; //change in water water level(Farrel and Clarke, Equ. 4)
784 IssmDouble late,longe;
785 IssmDouble lati,longi,ri;
786@@ -6025,6 +6037,9 @@
787 this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
788 this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
789
790+ /*recover earth area: */
791+ this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
792+
793 /*early return if rigid or elastic not requested:*/
794 if(!computerigid && !computeelastic) return;
795
796@@ -6080,7 +6095,7 @@
797 return;
798 }
799 /*}}}*/
800-void Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/
801+void Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){ /*{{{*/
802
803 /*diverse:*/
804 int gsize,dummy;
805@@ -6087,7 +6102,7 @@
806 bool spherical=true;
807 IssmDouble llr_list[NUMVERTICES][3];
808 IssmDouble xyz_list[NUMVERTICES][3];
809- IssmDouble area;
810+ IssmDouble area,eartharea;
811 IssmDouble I, S; //change in relative ice thickness and sea level
812 IssmDouble late,longe,re;
813 IssmDouble lati,longi,ri;
814@@ -6140,6 +6155,9 @@
815 /*how many dofs are we working with here? */
816 this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
817
818+ /*recover earth area: */
819+ this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
820+
821 /*retrieve indices:*/
822 if(computerigid){this->inputs2->GetArrayPtr(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);}
823
Note: See TracBrowser for help on using the repository browser.