Index: ../trunk-jpl/src/m/classes/SMBd18opdd.m =================================================================== --- ../trunk-jpl/src/m/classes/SMBd18opdd.m (revision 22851) +++ ../trunk-jpl/src/m/classes/SMBd18opdd.m (revision 22852) @@ -17,11 +17,14 @@ ismungsm = 0; isd18opd = 0; issetpddfac = 0; - istemperaturescaled = 0; + istemperaturescaled = 1; + isprecipscaled = 1; delta18o = NaN; delta18o_surface = NaN; temperatures_presentday = NaN; precipitations_presentday = NaN; + temperatures_reconstructed = NaN; + precipitations_reconstructed = NaN; pddfac_snow = NaN; pddfac_ice = NaN; requested_outputs = {}; @@ -38,6 +41,8 @@ function self = extrude(self,md) % {{{ if(self.isd18opd),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end if(self.isd18opd),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end + if(self.istemperaturescaled==0),self.temperatures_reconstructed=project3d(md,'vector',self.temperatures_reconstructed,'type','node');end + if(self.isprecipscaled),self.precipitations_reconstructed=project3d(md,'vector',self.precipitations_reconstructed,'type','node');end if(self.issetpddfac), self.pddfac_snow=project3d(md,'vector',self.pddfac_snow,'type','node');end if(self.issetpddfac), self.pddfac_ice=project3d(md,'vector',self.pddfac_ice,'type','node');end self.s0p=project3d(md,'vector',self.s0p,'type','node'); @@ -66,6 +71,7 @@ self.ismungsm = 0; self.isd18opd = 1; self.istemperaturescaled = 1; + self.isprecipscaled = 1; self.desfac = 0.5; self.rlaps = 6.5; self.rlapslgm = 6.5; @@ -92,6 +98,16 @@ md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1); md = checkfield(md,'fieldname','smb.dpermil','>=',0,'numel',1); md = checkfield(md,'fieldname','smb.f','>=',0,'numel',1); + if(self.istemperaturescaled==0) + lent=size(self.temperatures_reconstructed,2); + multt=ceil(lent/12)*12; + md = checkfield(md,'fieldname','smb.temperatures_reconstructed','size',[md.mesh.numberofvertices+1 multt],'NaN',1,'Inf',1,'timeseries',1); + end + if(self.isprecipscaled==0) + lenp=size(self.temperatures_reconstructed,2); + multp=ceil(lent/12)*12; + md = checkfield(md,'fieldname','smb.precipitations_reconstructed','size',[md.mesh.numberofvertices+1 multt],'NaN',1,'Inf',1,'timeseries',1); + end end if(self.issetpddfac==1) md = checkfield(md,'fieldname','smb.pddfac_snow','>=',0,'NaN',1,'Inf',1); @@ -113,7 +129,14 @@ if(self.isd18opd==1) fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm/d18opd is activated'); fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated'); - 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)'); + 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)'); + 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)'); + if(self.istemperaturescaled==0) + fielddisplay(self,'temperatures_reconstructed','monthly historical surface temperatures [K], required if delta18o/mungsm/d18opd is activated and istemperaturescaled is not activated'); + end + if(self.isprecipscaled==0) + fielddisplay(self,'precipitations_reconstructed','monthly historical precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated and isprecipscaled is not activated'); + end fielddisplay(self,'delta18o','delta18o [per mil], required if pdd is activated and d18opd activated'); fielddisplay(self,'dpermil','degree per mil, required if d18opd is activated'); fielddisplay(self,'f','precip/temperature scaling factor, required if d18opd is activated'); @@ -148,9 +171,16 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','temperatures_presentday','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); 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); WriteData(fid,prefix,'object',self,'class','smb','fieldname','istemperaturescaled','format','Boolean'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','isprecipscaled','format','Boolean'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','delta18o','format','DoubleMat','mattype',1,'timeserieslength',2,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','dpermil','format','Double'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','f','format','Double'); + if self.istemperaturescaled==0 + WriteData(fid,prefix,'object',self,'class','smb','fieldname','temperatures_reconstructed','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + end + if self.isprecipscaled==0 + 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); + end end if self.issetpddfac==1 WriteData(fid,prefix,'object',self,'class','smb','fieldname','pddfac_snow','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22851) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22852) @@ -278,6 +278,7 @@ SmbIssetpddfacEnum, SmbIsshortwaveEnum, SmbIstemperaturescaledEnum, + SmbIsprecipscaledEnum, SmbIsthermalEnum, SmbIsturbulentfluxEnum, SmbKEnum, @@ -532,6 +533,7 @@ SmbPrecipitationEnum, SmbPrecipitationsLgmEnum, SmbPrecipitationsPresentdayEnum, + SmbPrecipitationsReconstructedEnum, SmbReEnum, SmbRefreezeEnum, SmbReiniEnum, @@ -543,6 +545,7 @@ SmbTaEnum, SmbTemperaturesLgmEnum, SmbTemperaturesPresentdayEnum, + SmbTemperaturesReconstructedEnum, SmbTEnum, SmbTeValueEnum, SmbTiniEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22851) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22852) @@ -286,6 +286,7 @@ case SmbIssetpddfacEnum : return "SmbIssetpddfac"; case SmbIsshortwaveEnum : return "SmbIsshortwave"; case SmbIstemperaturescaledEnum : return "SmbIstemperaturescaled"; + case SmbIsprecipscaledEnum : return "SmbIsprecipscaled"; case SmbIsthermalEnum : return "SmbIsthermal"; case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux"; case SmbKEnum : return "SmbK"; @@ -538,6 +539,7 @@ case SmbPrecipitationEnum : return "SmbPrecipitation"; case SmbPrecipitationsLgmEnum : return "SmbPrecipitationsLgm"; case SmbPrecipitationsPresentdayEnum : return "SmbPrecipitationsPresentday"; + case SmbPrecipitationsReconstructedEnum : return "SmbPrecipitationsReconstructed"; case SmbReEnum : return "SmbRe"; case SmbRefreezeEnum : return "SmbRefreeze"; case SmbReiniEnum : return "SmbReini"; @@ -549,6 +551,7 @@ case SmbTaEnum : return "SmbTa"; case SmbTemperaturesLgmEnum : return "SmbTemperaturesLgm"; case SmbTemperaturesPresentdayEnum : return "SmbTemperaturesPresentday"; + case SmbTemperaturesReconstructedEnum : return "SmbTemperaturesReconstructed"; case SmbTEnum : return "SmbT"; case SmbTeValueEnum : return "SmbTeValue"; case SmbTiniEnum : return "SmbTini"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22851) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22852) @@ -292,6 +292,7 @@ else if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum; else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum; else if (strcmp(name,"SmbIstemperaturescaled")==0) return SmbIstemperaturescaledEnum; + else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum; else if (strcmp(name,"SmbIsthermal")==0) return SmbIsthermalEnum; else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum; else if (strcmp(name,"SmbK")==0) return SmbKEnum; @@ -381,11 +382,11 @@ else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum; else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum; else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum; - else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum; + if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum; + else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum; else if (strcmp(name,"Base")==0) return BaseEnum; @@ -504,11 +505,11 @@ else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum; else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; - else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; + if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum; + else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum; else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum; @@ -550,6 +551,7 @@ else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum; else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum; else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum; + else if (strcmp(name,"SmbPrecipitationsReconstructed")==0) return SmbPrecipitationsReconstructedEnum; else if (strcmp(name,"SmbRe")==0) return SmbReEnum; else if (strcmp(name,"SmbRefreeze")==0) return SmbRefreezeEnum; else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum; @@ -561,6 +563,7 @@ else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum; else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum; + else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum; else if (strcmp(name,"SmbT")==0) return SmbTEnum; else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum; else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; @@ -625,13 +628,13 @@ else if (strcmp(name,"VzSSA")==0) return VzSSAEnum; else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum; else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum; - else if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum; - else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum; - else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"Waterheight")==0) return WaterheightEnum; + if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum; + else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum; + else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum; + else if (strcmp(name,"Waterheight")==0) return WaterheightEnum; else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum; else if (strcmp(name,"InputsEND")==0) return InputsENDEnum; else if (strcmp(name,"Absolute")==0) return AbsoluteEnum; @@ -748,13 +751,13 @@ else if (strcmp(name,"FSSolver")==0) return FSSolverEnum; else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; - else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; - else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; - else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; + if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; + else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; + else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; + else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; else if (strcmp(name,"GenericParam")==0) return GenericParamEnum; else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum; else if (strcmp(name,"GiaIvinsAnalysis")==0) return GiaIvinsAnalysisEnum; @@ -871,13 +874,13 @@ else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum; else if (strcmp(name,"MaxVel")==0) return MaxVelEnum; else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; - else if (strcmp(name,"MaxVy")==0) return MaxVyEnum; - else if (strcmp(name,"MaxVz")==0) return MaxVzEnum; - else if (strcmp(name,"Melange")==0) return MelangeEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; + if (strcmp(name,"MaxVy")==0) return MaxVyEnum; + else if (strcmp(name,"MaxVz")==0) return MaxVzEnum; + else if (strcmp(name,"Melange")==0) return MelangeEnum; + else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum; else if (strcmp(name,"MeshX")==0) return MeshXEnum; else if (strcmp(name,"MeshY")==0) return MeshYEnum; @@ -994,13 +997,13 @@ else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; - else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; - else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; - else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum; + if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; + else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; + else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum; + else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum; else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum; else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum; else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; @@ -1117,13 +1120,13 @@ else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; - else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; - else if (strcmp(name,"Tria")==0) return TriaEnum; - else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum; + if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; + else if (strcmp(name,"Tria")==0) return TriaEnum; + else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; + else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum; else if (strcmp(name,"VectorParam")==0) return VectorParamEnum; else if (strcmp(name,"Vertex")==0) return VertexEnum; else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum; Index: ../trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- ../trunk-jpl/src/c/shared/Elements/elements.h (revision 22851) +++ ../trunk-jpl/src/c/shared/Elements/elements.h (revision 22852) @@ -33,8 +33,9 @@ IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); void ComputeD18OTemperaturePrecipitationFromPD(IssmDouble d018,IssmDouble dpermil,bool isTemperatureScaled, - IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday, - IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); + bool isPrecipScaled, IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday, + IssmDouble* PrecipitationReconstructed,IssmDouble* TemperatureReconstructed, IssmDouble* monthlytemperaturesout, + IssmDouble* monthlyprecout); IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.); IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness); Index: ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp (revision 22851) +++ ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp (revision 22852) @@ -7,8 +7,9 @@ #include "../Numerics/numerics.h" void ComputeD18OTemperaturePrecipitationFromPD(IssmDouble d018,IssmDouble dpermil,bool isTemperatureScaled, - IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday, - IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout){ + bool isPrecipScaled, IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday, + IssmDouble* PrecipitationReconstructed,IssmDouble* TemperatureReconstructed, IssmDouble* monthlytemperaturesout, + IssmDouble* monthlyprecout){ IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12]; IssmDouble deltaTemp; @@ -22,10 +23,14 @@ for (int imonth = 0; imonth<12; imonth++){ if (isTemperatureScaled)monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth] + deltaTemp; - else monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth]; + else{ + monthlytemperaturestmp[imonth] = TemperatureReconstructed[imonth]; + deltaTemp=TemperatureReconstructed[imonth]-TemperaturePresentday[imonth]; + } - monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp); - + if (isPrecipScaled)monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp); + else monthlyprectmp[imonth] = PrecipitationReconstructed[imonth]; + /*Assign output pointer*/ *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth]; *(monthlyprecout+imonth) = monthlyprectmp[imonth]; Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22851) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22852) @@ -534,11 +534,14 @@ IssmDouble* monthlyprec=xNew(12*numvertices); IssmDouble* TemperaturesPresentday=xNew(12*numvertices); IssmDouble* PrecipitationsPresentday=xNew(12*numvertices); + IssmDouble* TemperaturesReconstructed=xNew(12*numvertices); + IssmDouble* PrecipitationsReconstructed=xNew(12*numvertices); IssmDouble* tmp=xNew(numvertices); IssmDouble Delta18oTime; IssmDouble dpermil,f; IssmDouble time,yts,time_yr,month,time_clim; bool isTemperatureScaled=true; + bool isPrecipScaled=true; this->parameters->FindParam(&time,TimeEnum); this->parameters->FindParam(&yts,ConstantsYtsEnum); this->parameters->FindParam(&f,SmbFEnum); @@ -548,9 +551,20 @@ /*Get some pdd parameters*/ dpermil=this->matpar->GetMaterialParameter(SmbDpermilEnum); + this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum); + this->parameters->FindParam(&isPrecipScaled,SmbIsprecipscaledEnum); + /*Recover present day temperature and precipitation*/ Input* input=this->inputs->GetInput(SmbTemperaturesPresentdayEnum); _assert_(input); Input* input2=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input2); + Input* input3=NULL; + if(!isTemperatureScaled){ + input3=this->inputs->GetInput(SmbTemperaturesPresentdayEnum); _assert_(input3); + } + Input* input4=NULL; + if(!isPrecipScaled){ + input4=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input4); + } int offset; offset=dynamic_cast(input)->GetTimeInputOffset(time_yr); @@ -565,19 +579,27 @@ gauss->GaussVertex(iv); input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,time_clim+month/12.*yts); input2->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,time_clim+month/12.*yts); + PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts; - PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts; + if(!isTemperatureScaled){ + input3->GetInputValue(&TemperaturesReconstructed[iv*12+month],gauss,time_clim+month/12.*yts); + } + if(!isPrecipScaled){ + input4->GetInputValue(&PrecipitationsReconstructed[iv*12+month],gauss,time_clim+month/12.*yts); + PrecipitationsReconstructed[iv*12+month]=PrecipitationsReconstructed[iv*12+month]*yts; + } + } } /*Recover interpolation parameters at time t*/ this->parameters->FindParam(&Delta18oTime,SmbDelta18oEnum,time); - this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum); /*Compute the temperature and precipitation*/ for(int iv=0;iv(monthlyprec); xDelete(TemperaturesPresentday); xDelete(PrecipitationsPresentday); + xDelete(TemperaturesReconstructed); + xDelete(PrecipitationsReconstructed); xDelete(tmp); } Index: ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22851) +++ ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 22852) @@ -20,7 +20,7 @@ void SmbAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ int smb_model; - bool isdelta18o,ismungsm,isd18opd,issetpddfac; + bool isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled; /*Update elements: */ int counter=0; @@ -82,13 +82,14 @@ iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum); iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum); iomodel->FetchDataToInput(elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum); - } - else{ + }else{ iomodel->FetchDataToInput(elements,"md.smb.precipitation",SmbPrecipitationEnum); iomodel->FetchDataToInput(elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum); } break; case SMBd18opddEnum: + iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled"); + iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled"); iomodel->FindConstant(&ismungsm,"md.smb.ismungsm"); iomodel->FindConstant(&isd18opd,"md.smb.isd18opd"); iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac"); @@ -98,6 +99,12 @@ if(isd18opd){ iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum); iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum); + if(!istemperaturescaled){ + iomodel->FetchDataToInput(elements,"md.smb.temperatures_reconstructed",SmbTemperaturesReconstructedEnum); + } + if(!isprecipscaled){ + iomodel->FetchDataToInput(elements,"md.smb.precipitations_reconstructed",SmbPrecipitationsReconstructedEnum); + } } if(issetpddfac){ iomodel->FetchDataToInput(elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.); @@ -219,6 +226,7 @@ if(isd18opd){ parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.istemperaturescaled",SmbIstemperaturescaledEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipscaled",SmbIsprecipscaledEnum)); iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2); parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M)); iomodel->DeleteData(temp,"md.smb.delta18o");