Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23807) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23808) @@ -214,4 +214,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdensification",SmbIsdensificationEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum)); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23807) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23808) @@ -3215,5 +3215,6 @@ IssmDouble Tz,Vz=0.0; IssmDouble rho_ice, rho_water,aSnow,aIce; - IssmDouble time,dt; + IssmDouble time,dt,starttime,finaltime; + IssmDouble timeclim=0.0; IssmDouble t,smb_dt; IssmDouble yts; @@ -3250,7 +3251,10 @@ IssmDouble dMass=0.0; bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux; + bool isclimatology=false; IssmDouble init_scaling=0.0; IssmDouble thermo_scaling=1.0; IssmDouble adThresh=1023.0; + int offset0, offsetend; + IssmDouble time0, timeend; /*}}}*/ @@ -3296,4 +3300,6 @@ parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*transient core time step*/ parameters->FindParam(&yts,ConstantsYtsEnum); + parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); + parameters->FindParam(&starttime,TimesteppingStartTimeEnum); parameters->FindParam(&smb_dt,SmbDtEnum); /*time period for the smb solution, usually smaller than the glaciological dt*/ parameters->FindParam(&aIdx,SmbAIdxEnum); @@ -3305,4 +3311,5 @@ parameters->FindParam(&t0dry,SmbT0dryEnum); parameters->FindParam(&K,SmbKEnum); + parameters->FindParam(&isclimatology,SmbIsclimatologyEnum); parameters->FindParam(&isgraingrowth,SmbIsgraingrowthEnum); parameters->FindParam(&isalbedo,SmbIsalbedoEnum); @@ -3466,4 +3473,15 @@ //go back to time - deltaT: time-=dt; + if (isclimatology){ + //If this is a climatology, we need to repeat the forcing after the final time + offset0=dynamic_cast(Ta_input)->GetTimeInputOffset(starttime); + offsetend=dynamic_cast(Ta_input)->GetTimeInputOffset(finaltime); + time0=dynamic_cast(Ta_input)->GetTimeByOffset(offset0); + timeend=dynamic_cast(Ta_input)->GetTimeByOffset(offsetend); + timeclim=time0+fmod(time-time0,timeend-time0); + } + else{ + timeclim=time; + } /*Start loop: */ @@ -3474,11 +3492,11 @@ /*extract daily data:{{{*/ - Ta_input->GetInputValue(&Ta,gauss,t);//screen level air temperature [K] - V_input->GetInputValue(&V,gauss,t); //wind speed [m s-1] - Dlwr_input->GetInputValue(&dlw,gauss,t); //downward longwave radiation flux [W m-2] - Dswr_input->GetInputValue(&dsw,gauss,t); //downward shortwave radiation flux [W m-2] - P_input->GetInputValue(&P,gauss,t); //precipitation [kg m-2] - eAir_input->GetInputValue(&eAir,gauss,t); //screen level vapor pressure [Pa] - pAir_input->GetInputValue(&pAir,gauss,t); // screen level air pressure [Pa] + Ta_input->GetInputValue(&Ta,gauss,t-time+timeclim);//screen level air temperature [K] + V_input->GetInputValue(&V,gauss,t-time+timeclim); //wind speed [m s-1] + Dlwr_input->GetInputValue(&dlw,gauss,t-time+timeclim); //downward longwave radiation flux [W m-2] + Dswr_input->GetInputValue(&dsw,gauss,t-time+timeclim); //downward shortwave radiation flux [W m-2] + P_input->GetInputValue(&P,gauss,t-time+timeclim); //precipitation [kg m-2] + eAir_input->GetInputValue(&eAir,gauss,t-time+timeclim); //screen level vapor pressure [Pa] + pAir_input->GetInputValue(&pAir,gauss,t-time+timeclim); // screen level air pressure [Pa] teValue_input->GetInputValue(&teValue,gauss); // screen level air pressure [Pa] aValue_input->GetInputValue(&aValue,gauss); // screen level air pressure [Pa] Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23807) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23808) @@ -335,4 +335,5 @@ SmbIsaccumulationEnum, SmbIsalbedoEnum, + SmbIsclimatologyEnum, SmbIsd18opdEnum, SmbIsdelta18oEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23807) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23808) @@ -343,4 +343,5 @@ case SmbIsaccumulationEnum : return "SmbIsaccumulation"; case SmbIsalbedoEnum : return "SmbIsalbedo"; + case SmbIsclimatologyEnum : return "SmbIsclimatology"; case SmbIsd18opdEnum : return "SmbIsd18opd"; case SmbIsdelta18oEnum : return "SmbIsdelta18o"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23807) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23808) @@ -349,4 +349,5 @@ else if (strcmp(name,"SmbIsaccumulation")==0) return SmbIsaccumulationEnum; else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum; + else if (strcmp(name,"SmbIsclimatology")==0) return SmbIsclimatologyEnum; else if (strcmp(name,"SmbIsd18opd")==0) return SmbIsd18opdEnum; else if (strcmp(name,"SmbIsdelta18o")==0) return SmbIsdelta18oEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum; else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum; - else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; + if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum; + else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; else if (strcmp(name,"Step")==0) return StepEnum; else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"EplHeadHydrostep")==0) return EplHeadHydrostepEnum; else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum; - else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum; + if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum; + else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum; else if (strcmp(name,"EsaDeltathickness")==0) return EsaDeltathicknessEnum; else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum; else if (strcmp(name,"SmbC")==0) return SmbCEnum; - else if (strcmp(name,"SmbDailysnowfall")==0) return SmbDailysnowfallEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SmbDailyrainfall")==0) return SmbDailyrainfallEnum; + if (strcmp(name,"SmbDailysnowfall")==0) return SmbDailysnowfallEnum; + else if (strcmp(name,"SmbDailyrainfall")==0) return SmbDailyrainfallEnum; else if (strcmp(name,"SmbDailydsradiation")==0) return SmbDailydsradiationEnum; else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum; @@ -751,9 +752,9 @@ 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 stage=7; } if(stage==7){ - if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum; + 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; @@ -874,9 +875,9 @@ 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 stage=8; } if(stage==8){ - if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; + 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; @@ -997,9 +998,9 @@ else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; else if (strcmp(name,"None")==0) return NoneEnum; - else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; + if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum; + else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"P2xP1")==0) return P2xP1Enum; else if (strcmp(name,"P2xP4")==0) return P2xP4Enum; - else if (strcmp(name,"Paterson")==0) return PatersonEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Pengrid")==0) return PengridEnum; + if (strcmp(name,"Paterson")==0) return PatersonEnum; + else if (strcmp(name,"Pengrid")==0) return PengridEnum; else if (strcmp(name,"Penpair")==0) return PenpairEnum; else if (strcmp(name,"Penta")==0) return PentaEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum; else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum; - else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"StrainRate")==0) return StrainRateEnum; + if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum; + else if (strcmp(name,"StrainRate")==0) return StrainRateEnum; else if (strcmp(name,"StressbalanceViscosityOvershoot")==0) return StressbalanceViscosityOvershootEnum; else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 23807) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 23808) @@ -22,4 +22,5 @@ isdensification; isturbulentflux; + isclimatology; %inputs: @@ -159,4 +160,5 @@ self.isdensification=1; self.isturbulentflux=1; + self.isclimatology=0; self.aIdx = 1; @@ -213,4 +215,5 @@ md = checkfield(md,'fieldname','smb.isdensification','values',[0 1]); md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]); + md = checkfield(md,'fieldname','smb.isclimatology','values',[0 1]); md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value @@ -277,4 +280,5 @@ fielddisplay(self,'isdensification','run densification module (default true)'); fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)'); + fielddisplay(self,'isclimatology','repeat all forcings when past last forcing time (default false)'); fielddisplay(self,'Ta','2 m air temperature, in Kelvin'); fielddisplay(self,'V','wind speed (m s-1)'); @@ -370,4 +374,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','isdensification','format','Boolean'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','isturbulentflux','format','Boolean'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','isclimatology','format','Boolean'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 23807) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 23808) @@ -132,4 +132,5 @@ string = "%s\n%s"%(string,fielddisplay(self,'isdensification','run densification module (default true)')) string = "%s\n%s"%(string,fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)')) + string = "%s\n%s"%(string,fielddisplay(self,'isclimatology','repeat all forcings when past last forcing time (default false)')) string = "%s\n%s"%(string,fielddisplay(self,'Ta','2 m air temperature, in Kelvin')) string = "%s\n%s"%(string,fielddisplay(self,'V','wind speed (m s-1)')) @@ -237,4 +238,5 @@ self.isdensification = 1 self.isturbulentflux = 1 + self.isclimatology = 0 self.aIdx = 1 @@ -291,4 +293,5 @@ md = checkfield(md,'fieldname','smb.isdensification','values',[0,1]) md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0,1]) + md = checkfield(md,'fieldname','smb.isclimatology','values',[0,1]) md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100) #-100/100 celsius min/max value @@ -357,4 +360,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','isdensification','format','Boolean') WriteData(fid,prefix,'object',self,'class','smb','fieldname','isturbulentflux','format','Boolean') + WriteData(fid,prefix,'object',self,'class','smb','fieldname','isclimatology','format','Boolean') WriteData(fid,prefix,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)