Changeset 23814

Timestamp:

03/29/19 15:50:05 (6 years ago)

Author:

schlegel

Message:

NEW: Add option to repeat smb forcing as a climatology for GEMB and smb components forcing

Location:

issm/trunk-jpl/src

Files:

• c/analyses/SmbAnalysis.cpp (modified) (3 diffs)
• c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (5 diffs)
• c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (6 diffs)
• m/classes/SMBcomponents.m (modified) (4 diffs)
• m/classes/SMBcomponents.py (modified) (4 diffs)
• m/classes/SMBforcing.m (modified) (5 diffs)
• m/classes/SMBforcing.py (modified) (3 diffs)
• m/classes/SMBgemb.m (modified) (1 diff)
• m/classes/SMBgemb.py (modified) (1 diff)
• m/classes/SMBmeltcomponents.m (modified) (4 diffs)
• m/classes/SMBmeltcomponents.py (modified) (4 diffs)
• m/classes/model.m (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp ¶

@@ -192 +192 @@
         switch(smb_model){
                 case SMBforcingEnum:
-                        /*Nothing to add to parameters*/
+                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
                         break;
                 case SMBgembEnum:
@@ -278 +278 @@
                         break;
                 case SMBcomponentsEnum:
-                        /*Nothing to add to parameters*/
+                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
                         break;
                 case SMBmeltcomponentsEnum:
-                        /*Nothing to add to parameters*/
+                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
                         break;
                 case SMBgradientscomponentsEnum:
@@ -322 +322 @@
         switch(smb_model){
                 case SMBforcingEnum:
-                        /*Nothing to be done*/
+                        SmbForcingx(femmodel);
                         break;
                 case SMBgembEnum:

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp ¶

@@ -6 +6 @@
 #include "../../shared/shared.h"
 #include "../../toolkits/toolkits.h"
-
+#include "../modules.h"
+
+void SmbForcingx(FemModel* femmodel){/*{{{*/
+
+        // void SmbForcingx(smb,ni){
+        //    INPUT parameters: ni: working size of arrays
+        //    OUTPUT: mass-balance (m/yr ice): agd(NA)
+        bool isclimatology=false;
+        femmodel->parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
+
+        if (isclimatology){
+                int v;
+                IssmDouble time,dt,delta,starttime,finaltime;
+                int offsetend;
+                IssmDouble time0, timeend, timeclim;
+
+                femmodel->parameters->FindParam(&time,TimeEnum); 
+                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+                femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
+                femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+                
+                if (time<=starttime+dt){
+                        InputDuplicatex(femmodel,SmbMassBalanceEnum,SmbMassBalanceClimateEnum);
+                }
+
+                /*Loop over all the elements of this partition*/
+                for(int i=0;i<femmodel->elements->Size();i++){
+                        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+
+                        int         numvertices = element->GetNumberOfVertices();
+                        IssmDouble* smb         = xNew<IssmDouble>(numvertices);
+
+                        /*Recover Smb*/
+                        //If this is a climatology, we need to repeat the forcing after the final time
+                        Input* smb_input=element->GetInput(SmbMassBalanceClimateEnum); _assert_(smb_input);
+
+                        //Get accumulation climatology value
+                        offsetend=dynamic_cast<TransientInput*>(smb_input)->GetTimeInputOffset(finaltime);
+                        time0=dynamic_cast<TransientInput*>(smb_input)->GetTimeByOffset(0);
+                        timeend=dynamic_cast<TransientInput*>(smb_input)->GetTimeByOffset(offsetend);
+                        if (time>time0 & timeend>time0){
+                                delta=(time-time0) - (timeend-time0)*((int)((time-time0)/(timeend-time0)));
+                                if (delta==0){ timeclim=timeend;}
+                                else{ timeclim=time0+delta;}
+                        }
+
+                        element->GetInputListOnVerticesAtTime(smb,SmbMassBalanceClimateEnum,timeclim);
+
+                        /*Add input to element and Free memory*/
+                        element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
+                        xDelete<IssmDouble>(smb);
+                }
+        }
+
+}/*}}}*/
 void SmbGradientsx(FemModel* femmodel){/*{{{*/
 
@@ -350 +404 @@
         //    OUTPUT: mass-balance (m/yr ice): agd(NA)
         int v;
+        bool isclimatology=false;
+        IssmDouble time,delta,starttime,finaltime;
+        int offsetend;
+        IssmDouble time0, timeend, timeclim;
+
+        femmodel->parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
+        femmodel->parameters->FindParam(&time,TimeEnum); 
+        femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
+        femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
 
         /*Loop over all the elements of this partition*/
@@ -363 +426 @@
 
                 /*Recover Smb Components*/
-                element->GetInputListOnVertices(acc,SmbAccumulationEnum);
-                element->GetInputListOnVertices(evap,SmbEvaporationEnum);
-                element->GetInputListOnVertices(runoff,SmbRunoffEnum);
+                if (isclimatology){
+
+                        //If this is a climatology, we need to repeat the forcing after the final time
+                        Input* acc_input=element->GetInput(SmbAccumulationEnum); _assert_(acc_input);
+                        Input* evap_input=element->GetInput(SmbEvaporationEnum); _assert_(evap_input);
+                        Input* runoff_input=element->GetInput(SmbRunoffEnum); _assert_(runoff_input);
+
+                        //Get accumulation climatology value
+                        offsetend=dynamic_cast<TransientInput*>(acc_input)->GetTimeInputOffset(finaltime);
+                        time0=dynamic_cast<TransientInput*>(acc_input)->GetTimeByOffset(0);
+                        timeend=dynamic_cast<TransientInput*>(acc_input)->GetTimeByOffset(offsetend);
+                        if (time>time0 & timeend>time0){
+                                delta=(time-time0) - (timeend-time0)*((int)((time-time0)/(timeend-time0)));
+                                if (delta==0){ timeclim=timeend;}
+                                else{ timeclim=time0+delta;}
+                        }
+
+                        element->GetInputListOnVerticesAtTime(acc,SmbAccumulationEnum,timeclim);
+
+                        //Get evaporation climatology value
+                        offsetend=dynamic_cast<TransientInput*>(evap_input)->GetTimeInputOffset(finaltime);
+                        time0=dynamic_cast<TransientInput*>(evap_input)->GetTimeByOffset(0);
+                        timeend=dynamic_cast<TransientInput*>(evap_input)->GetTimeByOffset(offsetend);
+                        if (time>time0 & timeend>time0){
+                                delta=(time-time0) - (timeend-time0)*((int)((time-time0)/(timeend-time0)));
+                                if (delta==0){ timeclim=timeend;}
+                                else{ timeclim=time0+delta;}
+                        }
+
+                        element->GetInputListOnVerticesAtTime(evap,SmbEvaporationEnum,timeclim);
+
+                        //Get runoff climatology value
+                        offsetend=dynamic_cast<TransientInput*>(runoff_input)->GetTimeInputOffset(finaltime);
+                        time0=dynamic_cast<TransientInput*>(runoff_input)->GetTimeByOffset(0);
+                        timeend=dynamic_cast<TransientInput*>(runoff_input)->GetTimeByOffset(offsetend);
+                        if (time>time0 & timeend>time0){
+                                delta=(time-time0) - (timeend-time0)*((int)((time-time0)/(timeend-time0)));
+                                if (delta==0){ timeclim=timeend;}
+                                else{ timeclim=time0+delta;}
+                        }
+
+                        element->GetInputListOnVerticesAtTime(runoff,SmbRunoffEnum,timeclim);
+                }
+                else{
+                        element->GetInputListOnVertices(acc,SmbAccumulationEnum);
+                        element->GetInputListOnVertices(evap,SmbEvaporationEnum);
+                        element->GetInputListOnVertices(runoff,SmbRunoffEnum);
+                }
 
                 // loop over all vertices
@@ -391 +499 @@
         //    OUTPUT: mass-balance (m/yr ice): agd(NA)
         int v;
+        bool isclimatology=false;
+        IssmDouble time,delta,starttime,finaltime;
+        int offsetend;
+        IssmDouble time0, timeend, timeclim;
+
+        femmodel->parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
+        femmodel->parameters->FindParam(&time,TimeEnum);
+        femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
+        femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
 
         /*Loop over all the elements of this partition*/
@@ -405 +522 @@
 
                 /*Recover Smb Components*/
-                element->GetInputListOnVertices(acc,SmbAccumulationEnum);
-                element->GetInputListOnVertices(evap,SmbEvaporationEnum);
-                element->GetInputListOnVertices(melt,SmbMeltEnum);
-                element->GetInputListOnVertices(refreeze,SmbRefreezeEnum);
+                if (isclimatology){
+
+                        //If this is a climatology, we need to repeat the forcing after the final time
+                        Input* acc_input=element->GetInput(SmbAccumulationEnum); _assert_(acc_input);
+                        Input* evap_input=element->GetInput(SmbEvaporationEnum); _assert_(evap_input);
+                        Input* melt_input=element->GetInput(SmbMeltEnum); _assert_(melt_input);
+                        Input* refreeze_input=element->GetInput(SmbRefreezeEnum); _assert_(refreeze_input);
+
+                        //Get accumulation climatology value
+                        offsetend=dynamic_cast<TransientInput*>(acc_input)->GetTimeInputOffset(finaltime);
+                        time0=dynamic_cast<TransientInput*>(acc_input)->GetTimeByOffset(0);
+                        timeend=dynamic_cast<TransientInput*>(acc_input)->GetTimeByOffset(offsetend);
+                        if (time>time0 & timeend>time0){
+                                delta=(time-time0) - (timeend-time0)*((int)((time-time0)/(timeend-time0)));
+                                if (delta==0){ timeclim=timeend;}
+                                else{ timeclim=time0+delta;}
+                        }
+
+                        element->GetInputListOnVerticesAtTime(acc,SmbAccumulationEnum,timeclim);
+
+                        //Get evaporation climatology value
+                        offsetend=dynamic_cast<TransientInput*>(evap_input)->GetTimeInputOffset(finaltime);
+                        time0=dynamic_cast<TransientInput*>(evap_input)->GetTimeByOffset(0);
+                        timeend=dynamic_cast<TransientInput*>(evap_input)->GetTimeByOffset(offsetend);
+                        if (time>time0){
+                                delta=(time-time0) - (timeend-time0)*((int)((time-time0)/(timeend-time0)));
+                                if (delta==0){ timeclim=timeend;}
+                                else{ timeclim=time0+delta;}
+                        }
+
+                        element->GetInputListOnVerticesAtTime(evap,SmbEvaporationEnum,timeclim);
+
+                        //Get melt climatology value
+                        offsetend=dynamic_cast<TransientInput*>(melt_input)->GetTimeInputOffset(finaltime);
+                        time0=dynamic_cast<TransientInput*>(melt_input)->GetTimeByOffset(0);
+                        timeend=dynamic_cast<TransientInput*>(melt_input)->GetTimeByOffset(offsetend);
+                        if (time>time0){
+                                delta=(time-time0) - (timeend-time0)*((int)((time-time0)/(timeend-time0)));
+                                if (delta==0){ timeclim=timeend;}
+                                else{ timeclim=time0+delta;}
+                        }
+
+                        element->GetInputListOnVerticesAtTime(melt,SmbMeltEnum,timeclim);
+
+                        //Get refreeze climatology value
+                        offsetend=dynamic_cast<TransientInput*>(refreeze_input)->GetTimeInputOffset(finaltime);
+                        time0=dynamic_cast<TransientInput*>(refreeze_input)->GetTimeByOffset(0);
+                        timeend=dynamic_cast<TransientInput*>(refreeze_input)->GetTimeByOffset(offsetend);
+                        if (time>time0){
+                                delta=(time-time0) - (timeend-time0)*((int)((time-time0)/(timeend-time0)));
+                                if (delta==0){ timeclim=timeend;}
+                                else{ timeclim=time0+delta;}
+                        }
+
+                        element->GetInputListOnVerticesAtTime(refreeze,SmbRefreezeEnum,timeclim);
+                }
+                else{
+                        element->GetInputListOnVertices(acc,SmbAccumulationEnum);
+                        element->GetInputListOnVertices(evap,SmbEvaporationEnum);
+                        element->GetInputListOnVertices(melt,SmbMeltEnum);
+                        element->GetInputListOnVertices(refreeze,SmbRefreezeEnum);
+                }
 
                 // loop over all vertices

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h ¶

@@ -10 +10 @@
 /* local prototypes: */
 void SurfaceMassBalancex(FemModel* femmodel);
+void SmbForcingx(FemModel* femmodel);
 void SmbGradientsx(FemModel* femmodel);
 void SmbGradientsElax(FemModel* femmodel);

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -640 +640 @@
         SmbIsInitializedEnum,
         SmbMassBalanceEnum,
+        SmbMassBalanceClimateEnum,
         SmbMeltEnum,
         SmbMonthlytemperaturesEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp ¶

@@ -646 +646 @@
                 case SmbIsInitializedEnum : return "SmbIsInitialized";
                 case SmbMassBalanceEnum : return "SmbMassBalance";
+                case SmbMassBalanceClimateEnum : return "SmbMassBalanceClimate";
                 case SmbMeltEnum : return "SmbMelt";
                 case SmbMonthlytemperaturesEnum : return "SmbMonthlytemperatures";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp ¶

@@ -661 +661 @@
               else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
               else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
+              else if (strcmp(name,"SmbMassBalanceClimate")==0) return SmbMassBalanceClimateEnum;
               else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
               else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
@@ -751 +752 @@
               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 stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum;
+              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;
@@ -874 +875 @@
               else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
               else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
-              else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
+              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;
@@ -997 +998 @@
               else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum;
               else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
-              else if (strcmp(name,"None")==0) return NoneEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
+              if (strcmp(name,"None")==0) return NoneEnum;
+              else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
               else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
               else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
@@ -1120 +1121 @@
               else if (strcmp(name,"P2")==0) return P2Enum;
               else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
-              else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"Paterson")==0) return PatersonEnum;
+              if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
+              else if (strcmp(name,"Paterson")==0) return PatersonEnum;
               else if (strcmp(name,"Pengrid")==0) return PengridEnum;
               else if (strcmp(name,"Penpair")==0) return PenpairEnum;
@@ -1243 +1244 @@
               else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
               else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
-              else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum;
+              if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
+              else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum;
               else if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
               else if (strcmp(name,"StressbalanceViscosityOvershoot")==0) return StressbalanceViscosityOvershootEnum;

issm/trunk-jpl/src/m/classes/SMBcomponents.m ¶

@@ -6 +6 @@
 classdef SMBcomponents
         properties (SetAccess=public) 
+                isclimatology = 0;
                 accumulation = NaN;
                 runoff = NaN;
@@ -68 +69 @@
                         end
                         md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
+                        md = checkfield(md,'fieldname','smb.isclimatology','values',[0 1]);
+                        if (self.isclimatology)
+                                md = checkfield(md,'fieldname', 'smb.accumulation', 'size',[md.mesh.numberofvertices+1],...
+                                        'message',['accumulation must have md.mesh.numberofvertices+1 rows in order to force a climatology']);
+                                md = checkfield(md,'fieldname', 'smb.runoff', 'size',[md.mesh.numberofvertices+1],...
+                                        'message',['runoff must have md.mesh.numberofvertices+1 rows in order to force a climatology']);
+                                md = checkfield(md,'fieldname', 'smb.evaporation', 'size',[md.mesh.numberofvertices+1],...
+                                        'message',['evaporation must have md.mesh.numberofvertices+1 rows in order to force a climatology']);
+                        end
                 end % }}}
                 function disp(self) % {{{
@@ -74 +84 @@
                         fielddisplay(self,'runoff','amount of ice melt lost from the ice column [m/yr ice eq]');
                         fielddisplay(self,'evaporation','amount of ice lost to evaporative processes [m/yr ice eq]');
+                        fielddisplay(self,'isclimatology','repeat all forcings when past last forcing time (default false)');
                         fielddisplay(self,'requested_outputs','additional outputs requested');
                 end % }}}
@@ -93 +104 @@
                         end
                         WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray');
+                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','isclimatology','format','Boolean');
 
                 end % }}}

issm/trunk-jpl/src/m/classes/SMBcomponents.py ¶

@@ -16 +16 @@
                 self.runoff = float('NaN')
                 self.evaporation = float('NaN')
+                self.isclimatology = 0
                 self.requested_outputs      = []
                 #}}}
@@ -23 +24 @@
                 string="%s\n%s"%(string,fielddisplay(self,'runoff','amount of ice melt lost from the ice column [m/yr ice eq]'))
                 string="%s\n%s"%(string,fielddisplay(self,'evaporation','mount of ice lost to evaporative processes [m/yr ice eq]'))
+                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,'requested_outputs','additional outputs requested'))
                 return string
@@ -73 +75 @@
                 
                 md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1)
+                md = checkfield(md,'fieldname','smb.isclimatology','values',[0,1])
 
                 return md
@@ -92 +95 @@
                         outputs    =outputscopy
                 WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray')
+                WriteData(fid,prefix,'object',self,'class','smb','fieldname','isclimatology','format','Boolean')
+                if (self.isclimatology>0):
+                        md = checkfield(md,'fieldname', 'smb.accumulation', 'size',[md.mesh.numberofvertices+1],'message','accumulation must have md.mesh.numberofvertices+1 rows in order to force a climatology')
+                        md = checkfield(md,'fieldname', 'smb.runoff', 'size',[md.mesh.numberofvertices+1],'message','runoff must have md.mesh.numberofvertices+1 rows in order to force a climatology')
+                        md = checkfield(md,'fieldname', 'smb.evaporation', 'size',[md.mesh.numberofvertices+1],'message','evaporation must have md.mesh.numberofvertices+1 rows in order to force a climatology')
 
         # }}}

issm/trunk-jpl/src/m/classes/SMBforcing.m ¶

@@ -6 +6 @@
 classdef SMBforcing
         properties (SetAccess=public) 
+                isclimatology = 0;
                 mass_balance = NaN;
                 requested_outputs      = {};
@@ -49 +50 @@
 
                         if (strcmp(solution,'TransientSolution') & md.transient.issmb == 0), return; end
-                        
+
                         if ismember('MasstransportAnalysis',analyses),
                                 md = checkfield(md,'fieldname','smb.mass_balance','timeseries',1,'NaN',1,'Inf',1);
@@ -57 +58 @@
                         end
                         md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
+                        md = checkfield(md,'fieldname','smb.isclimatology','values',[0 1]);
+                        if (self.isclimatology)
+                                md = checkfield(md,'fieldname', 'smb.mass_balance', 'size',[md.mesh.numberofvertices+1],...
+                                        'message',['mass_balance must have md.mesh.numberofvertices+1 rows in order to force a climatology']);
+                        end
                 end % }}}
                 function disp(self) % {{{
                         disp(sprintf('   surface forcings parameters:'));
                         fielddisplay(self,'mass_balance','surface mass balance [m/yr ice eq]');
+                        fielddisplay(self,'isclimatology','repeat all forcings when past last forcing time (default false)');
                         fielddisplay(self,'requested_outputs','additional outputs requested');
                 end % }}}
@@ -79 +86 @@
                         end
                         WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray');
+                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','isclimatology','format','Boolean');
 
                 end % }}}
@@ -85 +93 @@
                         writejs1Darray(fid,[modelname '.smb.mass_balance'],self.mass_balance);
                         writejscellstring(fid,[modelname '.smb.requested_outputs'],self.requested_outputs);
+                        writejs1Darray(fid,[modelname '.smb.isclimatology'],self.isclimatology);
 
                 end % }}}

issm/trunk-jpl/src/m/classes/SMBforcing.py ¶

@@ -16 +16 @@
                 self.mass_balance = float('NaN')
                 self.requested_outputs      = []
+                self.isclimatology = 0
                 #}}}
         def __repr__(self): # {{{
                 string="   surface forcings parameters:"
                 string="%s\n%s"%(string,fielddisplay(self,'mass_balance','surface mass balance [m/yr ice eq]'))
+                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,'requested_outputs','additional outputs requested'))
                 return string
@@ -48 +50 @@
 
                 md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1)
+                md = checkfield(md,'fieldname','smb.isclimatology','values',[0,1])
+                if (self.isclimatology>0):
+                    md = checkfield(md,'fieldname', 'smb.mass_balance', 'size',[md.mesh.numberofvertices+1],'message','mass_balance must have md.mesh.numberofvertices+1 rows in order to force a climatology')
+
                 return md
         # }}}
@@ -65 +71 @@
                         outputs    =outputscopy
                 WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray')
+                WriteData(fid,prefix,'object',self,'class','smb','fieldname','isclimatology','format','Boolean')
 
         # }}}

issm/trunk-jpl/src/m/classes/SMBgemb.m ¶

@@ -218 +218 @@
 
                         md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value
-                        md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<',45); %max 500 km/h
-                        md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400);
-                        md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'>=',0);
-                        md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',100);
-                        md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1);
+                        md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<',45,'size',size(self.Ta)); %max 500 km/h
+                        md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400,'size',size(self.Ta));
+                        md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'size',size(self.Ta));
+                        md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',100,'size',size(self.Ta));
+                        md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1,'size',size(self.Ta));
+
+                        if (self.isclimatology)
+                                md = checkfield(md,'fieldname', 'smb.Ta', 'size',[md.mesh.numberofelements+1],...
+                                        'message',['Ta must have md.mesh.numberofelements+1 rows in order to force a climatology']);
+                                md = checkfield(md,'fieldname', 'smb.V', 'size',[md.mesh.numberofelements+1],...
+                                        'message',['V must have md.mesh.numberofelements+1 rows in order to force a climatology']);
+                                md = checkfield(md,'fieldname', 'smb.dswrf', 'size',[md.mesh.numberofelements+1],...
+                                        'message',['dswrf must have md.mesh.numberofelements+1 rows in order to force a climatology']);
+                                md = checkfield(md,'fieldname', 'smb.dlwrf', 'size',[md.mesh.numberofelements+1],...
+                                        'message',['dlwrf must have md.mesh.numberofelements+1 rows in order to force a climatology']);
+                                md = checkfield(md,'fieldname', 'smb.P', 'size',[md.mesh.numberofelements+1],...
+                                        'message',['P must have md.mesh.numberofelements+1 rows in order to force a climatology']);
+                                md = checkfield(md,'fieldname', 'smb.eAir', 'size',[md.mesh.numberofelements+1],...
+                                        'message',['eAir must have md.mesh.numberofelements+1 rows in order to force a climatology']);
+                        end
 
                         md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value

issm/trunk-jpl/src/m/classes/SMBgemb.py ¶

@@ -296 +296 @@
 
                 md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100) #-100/100 celsius min/max value
-                md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'<',45) #max 500 km/h
-                md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'< = ',1400)
-                md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'> = ',0)
-                md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'< = ',100)
-                md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1)
+                md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'<',45,'size',np.shape(self.Ta)) #max 500 km/h
+                md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'< = ',1400,'size',np.shape(self.Ta))
+                md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'size',np.shape(self.Ta))
+                md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'> = ',0,'< = ',100,'size',np.shape(self.Ta))
+                md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1,'size',np.shape(self.Ta))
+
+                if (self.isclimatology>0):
+                        md = checkfield(md,'fieldname', 'smb.Ta', 'size',[md.mesh.numberofelements+1],'message','Ta must have md.mesh.numberofelements+1 rows in order to force a climatology')
+                        md = checkfield(md,'fieldname', 'smb.V', 'size',[md.mesh.numberofelements+1],'message','V must have md.mesh.numberofelements+1 rows in order to force a climatology')
+                        md = checkfield(md,'fieldname', 'smb.dswrf', 'size',[md.mesh.numberofelements+1],'message','dswrf must have md.mesh.numberofelements+1 rows in order to force a climatology')
+                        md = checkfield(md,'fieldname', 'smb.dlwrf', 'size',[md.mesh.numberofelements+1],'message','dlwrf must have md.mesh.numberofelements+1 rows in order to force a climatology')
+                        md = checkfield(md,'fieldname', 'smb.P', 'size',[md.mesh.numberofelements+1],'message','P must have md.mesh.numberofelements+1 rows in order to force a climatology')
+                        md = checkfield(md,'fieldname', 'smb.eAir', 'size',[md.mesh.numberofelements+1],'message','eAir must have md.mesh.numberofelements+1 rows in order to force a climatology')
 
                 md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements],'NaN',1,'Inf',1,'>',273-100,'<',273+100) #-100/100 celsius min/max value

issm/trunk-jpl/src/m/classes/SMBmeltcomponents.m ¶

@@ -6 +6 @@
 classdef SMBmeltcomponents
         properties (SetAccess=public) 
+                isclimatology = 0;
                 accumulation = NaN;
                 evaporation = NaN;
@@ -78 +79 @@
                         end
                         md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
+                        md = checkfield(md,'fieldname','smb.isclimatology','values',[0 1]);
+                        if (self.isclimatology)
+                                md = checkfield(md,'fieldname', 'smb.accumulation', 'size',[md.mesh.numberofvertices+1],...
+                                        'message',['accumulation must have md.mesh.numberofvertices+1 rows in order to force a climatology']);
+                                md = checkfield(md,'fieldname', 'smb.melt', 'size',[md.mesh.numberofvertices+1],...
+                                        'message',['melt must have md.mesh.numberofvertices+1 rows in order to force a climatology']);
+                                md = checkfield(md,'fieldname', 'smb.refreeze', 'size',[md.mesh.numberofvertices+1],...
+                                        'message',['refreeze must have md.mesh.numberofvertices+1 rows in order to force a climatology']);
+                                md = checkfield(md,'fieldname', 'smb.evaporation', 'size',[md.mesh.numberofvertices+1],...
+                                        'message',['evaporation must have md.mesh.numberofvertices+1 rows in order to force a climatology']);
+                        end
                 end % }}}
                 function disp(self) % {{{
@@ -85 +97 @@
                         fielddisplay(self,'melt','amount of ice melt in ice column [m/yr ice eq]');
                         fielddisplay(self,'refreeze','amount of ice melt refrozen in ice column [m/yr ice eq]');
+                        fielddisplay(self,'isclimatology','repeat all forcings when past last forcing time (default false)');
                         fielddisplay(self,'requested_outputs','additional outputs requested');
                 end % }}}
@@ -105 +118 @@
                         end
                         WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray');
+                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','isclimatology','format','Boolean');
 
                 end % }}}

issm/trunk-jpl/src/m/classes/SMBmeltcomponents.py ¶

@@ -16 +16 @@
                 self.runoff = float('NaN')
                 self.evaporation = float('NaN')
+                self.isclimatology = 0
                 self.requested_outputs      = []
                 #}}}
@@ -24 +25 @@
                 string="%s\n%s"%(string,fielddisplay(self,'melt','amount of ice melt in the ice column [m/yr ice eq]'))
                 string="%s\n%s"%(string,fielddisplay(self,'refreeze','amount of ice melt refrozen in the ice column [m/yr ice eq]'))
+                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,'requested_outputs','additional outputs requested'))
                 return string
@@ -85 +87 @@
 
                 md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1)
+                md = checkfield(md,'fieldname','smb.isclimatology','values',[0,1])
                 return md
         # }}}
@@ -104 +107 @@
                         outputs    =outputscopy
                 WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray')
+                WriteData(fid,prefix,'object',self,'class','smb','fieldname','isclimatology','format','Boolean')
+                if (self.isclimatology>0):
+                        md = checkfield(md,'fieldname', 'smb.accumulation', 'size',[md.mesh.numberofvertices+1],'message','accumulation must have md.mesh.numberofvertices+1 rows in order to force a climatology')
+                        md = checkfield(md,'fieldname', 'smb.melt', 'size',[md.mesh.numberofvertices+1],'message','melt must have md.mesh.numberofvertices+1 rows in order to force a climatology')
+                        md = checkfield(md,'fieldname', 'smb.refreeze', 'size',[md.mesh.numberofvertices+1],'message','refreeze must have md.mesh.numberofvertices+1 rows in order to force a climatology')
+                        md = checkfield(md,'fieldname', 'smb.evaporation', 'size',[md.mesh.numberofvertices+1],'message','evaporation must have md.mesh.numberofvertices+1 rows in order to force a climatology')
 
         # }}}

issm/trunk-jpl/src/m/classes/model.m ¶

@@ -167 +167 @@
                                 else
                                         md.settings.results_on_nodes = {'all'};
+                                end
+                        end
+                        %2019 Mar 28
+                        if isa(md.smb.isclimatology,'double')
+                                if prod(size(md.smb.isclimatology)) ~= 1
+                                        md.smb.isclimatology = 0;
                                 end
                         end