source: issm/oecreview/Archive/22819-23185/ISSM-22851-22852.diff

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

CHG: added Archive/22819-23185
File size: 27.0 KB
RevLine 
[23186]1Index: ../trunk-jpl/src/m/classes/SMBd18opdd.m
2===================================================================
3--- ../trunk-jpl/src/m/classes/SMBd18opdd.m (revision 22851)
4+++ ../trunk-jpl/src/m/classes/SMBd18opdd.m (revision 22852)
5@@ -17,11 +17,14 @@
6 ismungsm = 0;
7 isd18opd = 0;
8 issetpddfac = 0;
9- istemperaturescaled = 0;
10+ istemperaturescaled = 1;
11+ isprecipscaled = 1;
12 delta18o = NaN;
13 delta18o_surface = NaN;
14 temperatures_presentday = NaN;
15 precipitations_presentday = NaN;
16+ temperatures_reconstructed = NaN;
17+ precipitations_reconstructed = NaN;
18 pddfac_snow = NaN;
19 pddfac_ice = NaN;
20 requested_outputs = {};
21@@ -38,6 +41,8 @@
22 function self = extrude(self,md) % {{{
23 if(self.isd18opd),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end
24 if(self.isd18opd),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end
25+ if(self.istemperaturescaled==0),self.temperatures_reconstructed=project3d(md,'vector',self.temperatures_reconstructed,'type','node');end
26+ if(self.isprecipscaled),self.precipitations_reconstructed=project3d(md,'vector',self.precipitations_reconstructed,'type','node');end
27 if(self.issetpddfac), self.pddfac_snow=project3d(md,'vector',self.pddfac_snow,'type','node');end
28 if(self.issetpddfac), self.pddfac_ice=project3d(md,'vector',self.pddfac_ice,'type','node');end
29 self.s0p=project3d(md,'vector',self.s0p,'type','node');
30@@ -66,6 +71,7 @@
31 self.ismungsm = 0;
32 self.isd18opd = 1;
33 self.istemperaturescaled = 1;
34+ self.isprecipscaled = 1;
35 self.desfac = 0.5;
36 self.rlaps = 6.5;
37 self.rlapslgm = 6.5;
38@@ -92,6 +98,16 @@
39 md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1);
40 md = checkfield(md,'fieldname','smb.dpermil','>=',0,'numel',1);
41 md = checkfield(md,'fieldname','smb.f','>=',0,'numel',1);
42+ if(self.istemperaturescaled==0)
43+ lent=size(self.temperatures_reconstructed,2);
44+ multt=ceil(lent/12)*12;
45+ md = checkfield(md,'fieldname','smb.temperatures_reconstructed','size',[md.mesh.numberofvertices+1 multt],'NaN',1,'Inf',1,'timeseries',1);
46+ end
47+ if(self.isprecipscaled==0)
48+ lenp=size(self.temperatures_reconstructed,2);
49+ multp=ceil(lent/12)*12;
50+ md = checkfield(md,'fieldname','smb.precipitations_reconstructed','size',[md.mesh.numberofvertices+1 multt],'NaN',1,'Inf',1,'timeseries',1);
51+ end
52 end
53 if(self.issetpddfac==1)
54 md = checkfield(md,'fieldname','smb.pddfac_snow','>=',0,'NaN',1,'Inf',1);
55@@ -113,7 +129,14 @@
56 if(self.isd18opd==1)
57 fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm/d18opd is activated');
58 fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated');
59- fielddisplay(self,'istemperaturescaled','if delta18o parametrisation from present day temperature and precipitation is activated, is temperature scaled to delta18o value? (0 or 1, default is 0)');
60+ fielddisplay(self,'istemperaturescaled','if delta18o parametrisation from present day temperature and precipitation is activated, is temperature scaled to delta18o value? (0 or 1, default is 1)');
61+ fielddisplay(self,'isprecipscaled','if delta18o parametrisation from present day temperature and precipitation is activated, is precip scaled to delta18o value? (0 or 1, default is 1)');
62+ if(self.istemperaturescaled==0)
63+ fielddisplay(self,'temperatures_reconstructed','monthly historical surface temperatures [K], required if delta18o/mungsm/d18opd is activated and istemperaturescaled is not activated');
64+ end
65+ if(self.isprecipscaled==0)
66+ fielddisplay(self,'precipitations_reconstructed','monthly historical precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated and isprecipscaled is not activated');
67+ end
68 fielddisplay(self,'delta18o','delta18o [per mil], required if pdd is activated and d18opd activated');
69 fielddisplay(self,'dpermil','degree per mil, required if d18opd is activated');
70 fielddisplay(self,'f','precip/temperature scaling factor, required if d18opd is activated');
71@@ -148,9 +171,16 @@
72 WriteData(fid,prefix,'object',self,'class','smb','fieldname','temperatures_presentday','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
73 WriteData(fid,prefix,'object',self,'class','smb','fieldname','precipitations_presentday','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
74 WriteData(fid,prefix,'object',self,'class','smb','fieldname','istemperaturescaled','format','Boolean');
75+ WriteData(fid,prefix,'object',self,'class','smb','fieldname','isprecipscaled','format','Boolean');
76 WriteData(fid,prefix,'object',self,'class','smb','fieldname','delta18o','format','DoubleMat','mattype',1,'timeserieslength',2,'yts',md.constants.yts);
77 WriteData(fid,prefix,'object',self,'class','smb','fieldname','dpermil','format','Double');
78 WriteData(fid,prefix,'object',self,'class','smb','fieldname','f','format','Double');
79+ if self.istemperaturescaled==0
80+ WriteData(fid,prefix,'object',self,'class','smb','fieldname','temperatures_reconstructed','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
81+ end
82+ if self.isprecipscaled==0
83+ WriteData(fid,prefix,'object',self,'class','smb','fieldname','precipitations_reconstructed','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
84+ end
85 end
86 if self.issetpddfac==1
87 WriteData(fid,prefix,'object',self,'class','smb','fieldname','pddfac_snow','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
88Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
89===================================================================
90--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22851)
91+++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22852)
92@@ -278,6 +278,7 @@
93 SmbIssetpddfacEnum,
94 SmbIsshortwaveEnum,
95 SmbIstemperaturescaledEnum,
96+ SmbIsprecipscaledEnum,
97 SmbIsthermalEnum,
98 SmbIsturbulentfluxEnum,
99 SmbKEnum,
100@@ -532,6 +533,7 @@
101 SmbPrecipitationEnum,
102 SmbPrecipitationsLgmEnum,
103 SmbPrecipitationsPresentdayEnum,
104+ SmbPrecipitationsReconstructedEnum,
105 SmbReEnum,
106 SmbRefreezeEnum,
107 SmbReiniEnum,
108@@ -543,6 +545,7 @@
109 SmbTaEnum,
110 SmbTemperaturesLgmEnum,
111 SmbTemperaturesPresentdayEnum,
112+ SmbTemperaturesReconstructedEnum,
113 SmbTEnum,
114 SmbTeValueEnum,
115 SmbTiniEnum,
116Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
117===================================================================
118--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22851)
119+++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22852)
120@@ -286,6 +286,7 @@
121 case SmbIssetpddfacEnum : return "SmbIssetpddfac";
122 case SmbIsshortwaveEnum : return "SmbIsshortwave";
123 case SmbIstemperaturescaledEnum : return "SmbIstemperaturescaled";
124+ case SmbIsprecipscaledEnum : return "SmbIsprecipscaled";
125 case SmbIsthermalEnum : return "SmbIsthermal";
126 case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux";
127 case SmbKEnum : return "SmbK";
128@@ -538,6 +539,7 @@
129 case SmbPrecipitationEnum : return "SmbPrecipitation";
130 case SmbPrecipitationsLgmEnum : return "SmbPrecipitationsLgm";
131 case SmbPrecipitationsPresentdayEnum : return "SmbPrecipitationsPresentday";
132+ case SmbPrecipitationsReconstructedEnum : return "SmbPrecipitationsReconstructed";
133 case SmbReEnum : return "SmbRe";
134 case SmbRefreezeEnum : return "SmbRefreeze";
135 case SmbReiniEnum : return "SmbReini";
136@@ -549,6 +551,7 @@
137 case SmbTaEnum : return "SmbTa";
138 case SmbTemperaturesLgmEnum : return "SmbTemperaturesLgm";
139 case SmbTemperaturesPresentdayEnum : return "SmbTemperaturesPresentday";
140+ case SmbTemperaturesReconstructedEnum : return "SmbTemperaturesReconstructed";
141 case SmbTEnum : return "SmbT";
142 case SmbTeValueEnum : return "SmbTeValue";
143 case SmbTiniEnum : return "SmbTini";
144Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
145===================================================================
146--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22851)
147+++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22852)
148@@ -292,6 +292,7 @@
149 else if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum;
150 else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum;
151 else if (strcmp(name,"SmbIstemperaturescaled")==0) return SmbIstemperaturescaledEnum;
152+ else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum;
153 else if (strcmp(name,"SmbIsthermal")==0) return SmbIsthermalEnum;
154 else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum;
155 else if (strcmp(name,"SmbK")==0) return SmbKEnum;
156@@ -381,11 +382,11 @@
157 else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum;
158 else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
159 else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum;
160- else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum;
161 else stage=4;
162 }
163 if(stage==4){
164- if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
165+ if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum;
166+ else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
167 else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
168 else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
169 else if (strcmp(name,"Base")==0) return BaseEnum;
170@@ -504,11 +505,11 @@
171 else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
172 else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
173 else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
174- else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
175 else stage=5;
176 }
177 if(stage==5){
178- if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
179+ if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
180+ else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
181 else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
182 else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum;
183 else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum;
184@@ -550,6 +551,7 @@
185 else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum;
186 else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum;
187 else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum;
188+ else if (strcmp(name,"SmbPrecipitationsReconstructed")==0) return SmbPrecipitationsReconstructedEnum;
189 else if (strcmp(name,"SmbRe")==0) return SmbReEnum;
190 else if (strcmp(name,"SmbRefreeze")==0) return SmbRefreezeEnum;
191 else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
192@@ -561,6 +563,7 @@
193 else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
194 else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
195 else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
196+ else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum;
197 else if (strcmp(name,"SmbT")==0) return SmbTEnum;
198 else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
199 else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
200@@ -625,13 +628,13 @@
201 else if (strcmp(name,"VzSSA")==0) return VzSSAEnum;
202 else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum;
203 else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
204- else if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum;
205- else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum;
206- else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
207 else stage=6;
208 }
209 if(stage==6){
210- if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
211+ if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum;
212+ else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum;
213+ else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
214+ else if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
215 else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
216 else if (strcmp(name,"InputsEND")==0) return InputsENDEnum;
217 else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
218@@ -748,13 +751,13 @@
219 else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
220 else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
221 else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
222- else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
223- else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
224- else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
225 else stage=7;
226 }
227 if(stage==7){
228- if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
229+ if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
230+ else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
231+ else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
232+ else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
233 else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
234 else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum;
235 else if (strcmp(name,"GiaIvinsAnalysis")==0) return GiaIvinsAnalysisEnum;
236@@ -871,13 +874,13 @@
237 else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
238 else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
239 else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
240- else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
241- else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
242- else if (strcmp(name,"Melange")==0) return MelangeEnum;
243 else stage=8;
244 }
245 if(stage==8){
246- if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
247+ if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
248+ else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
249+ else if (strcmp(name,"Melange")==0) return MelangeEnum;
250+ else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
251 else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
252 else if (strcmp(name,"MeshX")==0) return MeshXEnum;
253 else if (strcmp(name,"MeshY")==0) return MeshYEnum;
254@@ -994,13 +997,13 @@
255 else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum;
256 else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum;
257 else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum;
258- else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
259- else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
260- else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
261 else stage=9;
262 }
263 if(stage==9){
264- if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
265+ if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
266+ else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
267+ else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
268+ else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
269 else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum;
270 else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum;
271 else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
272@@ -1117,13 +1120,13 @@
273 else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
274 else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
275 else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
276- else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
277- else if (strcmp(name,"Tria")==0) return TriaEnum;
278- else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
279 else stage=10;
280 }
281 if(stage==10){
282- if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;
283+ if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
284+ else if (strcmp(name,"Tria")==0) return TriaEnum;
285+ else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
286+ else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;
287 else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
288 else if (strcmp(name,"Vertex")==0) return VertexEnum;
289 else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;
290Index: ../trunk-jpl/src/c/shared/Elements/elements.h
291===================================================================
292--- ../trunk-jpl/src/c/shared/Elements/elements.h (revision 22851)
293+++ ../trunk-jpl/src/c/shared/Elements/elements.h (revision 22852)
294@@ -33,8 +33,9 @@
295 IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
296 IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
297 void ComputeD18OTemperaturePrecipitationFromPD(IssmDouble d018,IssmDouble dpermil,bool isTemperatureScaled,
298- IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
299- IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
300+ bool isPrecipScaled, IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
301+ IssmDouble* PrecipitationReconstructed,IssmDouble* TemperatureReconstructed, IssmDouble* monthlytemperaturesout,
302+ IssmDouble* monthlyprecout);
303 IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.);
304 IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness);
305
306Index: ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp
307===================================================================
308--- ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp (revision 22851)
309+++ ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp (revision 22852)
310@@ -7,8 +7,9 @@
311 #include "../Numerics/numerics.h"
312
313 void ComputeD18OTemperaturePrecipitationFromPD(IssmDouble d018,IssmDouble dpermil,bool isTemperatureScaled,
314- IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
315- IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout){
316+ bool isPrecipScaled, IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
317+ IssmDouble* PrecipitationReconstructed,IssmDouble* TemperatureReconstructed, IssmDouble* monthlytemperaturesout,
318+ IssmDouble* monthlyprecout){
319
320 IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
321 IssmDouble deltaTemp;
322@@ -22,10 +23,14 @@
323 for (int imonth = 0; imonth<12; imonth++){
324
325 if (isTemperatureScaled)monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth] + deltaTemp;
326- else monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth];
327+ else{
328+ monthlytemperaturestmp[imonth] = TemperatureReconstructed[imonth];
329+ deltaTemp=TemperatureReconstructed[imonth]-TemperaturePresentday[imonth];
330+ }
331
332- monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp);
333-
334+ if (isPrecipScaled)monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp);
335+ else monthlyprectmp[imonth] = PrecipitationReconstructed[imonth];
336+
337 /*Assign output pointer*/
338 *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth];
339 *(monthlyprecout+imonth) = monthlyprectmp[imonth];
340Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
341===================================================================
342--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22851)
343+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22852)
344@@ -534,11 +534,14 @@
345 IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
346 IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
347 IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
348+ IssmDouble* TemperaturesReconstructed=xNew<IssmDouble>(12*numvertices);
349+ IssmDouble* PrecipitationsReconstructed=xNew<IssmDouble>(12*numvertices);
350 IssmDouble* tmp=xNew<IssmDouble>(numvertices);
351 IssmDouble Delta18oTime;
352 IssmDouble dpermil,f;
353 IssmDouble time,yts,time_yr,month,time_clim;
354 bool isTemperatureScaled=true;
355+ bool isPrecipScaled=true;
356 this->parameters->FindParam(&time,TimeEnum);
357 this->parameters->FindParam(&yts,ConstantsYtsEnum);
358 this->parameters->FindParam(&f,SmbFEnum);
359@@ -548,9 +551,20 @@
360 /*Get some pdd parameters*/
361 dpermil=this->matpar->GetMaterialParameter(SmbDpermilEnum);
362
363+ this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum);
364+ this->parameters->FindParam(&isPrecipScaled,SmbIsprecipscaledEnum);
365+
366 /*Recover present day temperature and precipitation*/
367 Input* input=this->inputs->GetInput(SmbTemperaturesPresentdayEnum); _assert_(input);
368 Input* input2=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input2);
369+ Input* input3=NULL;
370+ if(!isTemperatureScaled){
371+ input3=this->inputs->GetInput(SmbTemperaturesPresentdayEnum); _assert_(input3);
372+ }
373+ Input* input4=NULL;
374+ if(!isPrecipScaled){
375+ input4=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input4);
376+ }
377 int offset;
378
379 offset=dynamic_cast<TransientInput*>(input)->GetTimeInputOffset(time_yr);
380@@ -565,19 +579,27 @@
381 gauss->GaussVertex(iv);
382 input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,time_clim+month/12.*yts);
383 input2->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,time_clim+month/12.*yts);
384+ PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts;
385
386- PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts;
387+ if(!isTemperatureScaled){
388+ input3->GetInputValue(&TemperaturesReconstructed[iv*12+month],gauss,time_clim+month/12.*yts);
389+ }
390+ if(!isPrecipScaled){
391+ input4->GetInputValue(&PrecipitationsReconstructed[iv*12+month],gauss,time_clim+month/12.*yts);
392+ PrecipitationsReconstructed[iv*12+month]=PrecipitationsReconstructed[iv*12+month]*yts;
393+ }
394+
395 }
396 }
397
398 /*Recover interpolation parameters at time t*/
399 this->parameters->FindParam(&Delta18oTime,SmbDelta18oEnum,time);
400- this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum);
401
402 /*Compute the temperature and precipitation*/
403 for(int iv=0;iv<numvertices;iv++){
404- ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,isTemperatureScaled,f,
405- &PrecipitationsPresentday[iv*12], &TemperaturesPresentday[iv*12],
406+ ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,isTemperatureScaled,isPrecipScaled,
407+ f,&PrecipitationsPresentday[iv*12], &TemperaturesPresentday[iv*12],
408+ &PrecipitationsReconstructed[iv*12], &TemperaturesReconstructed[iv*12],
409 &monthlytemperatures[iv*12], &monthlyprec[iv*12]);
410 }
411
412@@ -622,6 +644,8 @@
413 xDelete<IssmDouble>(monthlyprec);
414 xDelete<IssmDouble>(TemperaturesPresentday);
415 xDelete<IssmDouble>(PrecipitationsPresentday);
416+ xDelete<IssmDouble>(TemperaturesReconstructed);
417+ xDelete<IssmDouble>(PrecipitationsReconstructed);
418 xDelete<IssmDouble>(tmp);
419
420 }
421Index: ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp
422===================================================================
423--- ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22851)
424+++ ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22852)
425@@ -20,7 +20,7 @@
426 void SmbAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
427
428 int smb_model;
429- bool isdelta18o,ismungsm,isd18opd,issetpddfac;
430+ bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled;
431
432 /*Update elements: */
433 int counter=0;
434@@ -82,13 +82,14 @@
435 iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
436 iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
437 iomodel->FetchDataToInput(elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum);
438- }
439- else{
440+ }else{
441 iomodel->FetchDataToInput(elements,"md.smb.precipitation",SmbPrecipitationEnum);
442 iomodel->FetchDataToInput(elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
443 }
444 break;
445 case SMBd18opddEnum:
446+ iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
447+ iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
448 iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
449 iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
450 iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
451@@ -98,6 +99,12 @@
452 if(isd18opd){
453 iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
454 iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
455+ if(!istemperaturescaled){
456+ iomodel->FetchDataToInput(elements,"md.smb.temperatures_reconstructed",SmbTemperaturesReconstructedEnum);
457+ }
458+ if(!isprecipscaled){
459+ iomodel->FetchDataToInput(elements,"md.smb.precipitations_reconstructed",SmbPrecipitationsReconstructedEnum);
460+ }
461 }
462 if(issetpddfac){
463 iomodel->FetchDataToInput(elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
464@@ -219,6 +226,7 @@
465 if(isd18opd){
466 parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum));
467 parameters->AddObject(iomodel->CopyConstantObject("md.smb.istemperaturescaled",SmbIstemperaturescaledEnum));
468+ parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipscaled",SmbIsprecipscaledEnum));
469 iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
470 parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
471 iomodel->DeleteData(temp,"md.smb.delta18o");
Note: See TracBrowser for help on using the repository browser.