Changeset 21917

Timestamp:

08/08/17 17:34:18 (8 years ago)

Author:

schlegel

Message:

CHG: add capability for more than one climatology to be passed in

Location:

issm/trunk-jpl/src

Files:

• c/classes/Elements/Element.cpp (modified) (3 diffs)
• c/classes/Inputs/TransientInput.cpp (modified) (1 diff)
• c/classes/Inputs/TransientInput.h (modified) (1 diff)
• m/classes/SMBd18opdd.m (modified) (1 diff)
• m/classes/SMBd18opdd.py (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -539 +539 @@
         IssmDouble Delta18oTime;
         IssmDouble dpermil,f;
-        IssmDouble time,yts,time_yr,month;
+        IssmDouble time,yts,time_yr,month,time_clim;
         this->parameters->FindParam(&time,TimeEnum);
         this->parameters->FindParam(&yts,ConstantsYtsEnum);
         this->parameters->FindParam(&f,SmbFEnum);
         time_yr=floor(time/yts)*yts;
+        time_clim=time_yr;
 
         /*Get some pdd parameters*/
@@ -551 +552 @@
         Input*     input=this->inputs->GetInput(SmbTemperaturesPresentdayEnum);    _assert_(input);
         Input*     input2=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input2);
+        int offset;
+
+        offset=dynamic_cast<TransientInput*>(input)->GetTimeInputOffset(time_yr);
+        if (fmod(offset,12)!=1){
+                time_clim=floor(dynamic_cast<TransientInput*>(input)->GetTimeByOffset(offset-fmod(offset,12)+1)/yts)*yts;
+        }
+
         /*loop over vertices: */
         Gauss* gauss=this->NewGauss();
@@ -556 +564 @@
                 for(int iv=0;iv<numvertices;iv++) {
                         gauss->GaussVertex(iv);
-                        input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts);
-                        input2->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month/12.*yts);
+                        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;

issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp

@@ -246 +246 @@
 
         delete input;
+}
+/*}}}*/
+int  TransientInput::GetTimeInputOffset(IssmDouble time){/*{{{*/
+
+        int        found;
+        int        offset;
+
+        /*go through the timesteps, and figure out which interval we
+         *     *fall within. Then interpolate the values on this interval: */
+        found=binary_search(&offset,time,this->timesteps,this->numtimesteps);
+        if(!found) _error_("Input not found (is TransientInput sorted ?)");
+
+        return offset;
+}
+/*}}}*/
+IssmDouble  TransientInput::GetTimeByOffset(int offset){/*{{{*/
+
+        if (offset < 0) offset=0;
+        _assert_(this->timesteps[offset]);
+
+        return this->timesteps[offset];
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h

@@ -70 +70 @@
                 void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime);
                 int  GetInputInterpolationType(){_error_("not implemented yet!");};
+                IssmDouble  GetTimeByOffset(int offset);
                 Input* GetTimeInput(IssmDouble time);
+                int  GetTimeInputOffset(IssmDouble time);
                 void GetTimeValues(IssmDouble* values,IssmDouble time){_error_("not implemented yet");};
                 void GetVectorFromInputs(Vector<IssmDouble>* vector,int* doflist);

issm/trunk-jpl/src/m/classes/SMBd18opdd.m

@@ -76 +76 @@
                                 md = checkfield(md,'fieldname','smb.rlapslgm','>=',0,'numel',1);
                                 if(self.isd18opd==1) 
-                                        md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1);
-                                        md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1,'Inf',1,'timeseries',1);
+                                        lent=size(self.temperatures_presentday,2);
+                                        lenp=size(self.precipitations_presentday,2);
+                                        multt=ceil(lent/12)*12;
+                                        multp=ceil(lenp/12)*12;
+                                        md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1 multt],'NaN',1,'Inf',1,'timeseries',1);
+                                        md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1 multp],'NaN',1,'Inf',1,'timeseries',1);
                                         md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1);
                                         md = checkfield(md,'fieldname','smb.dpermil','>=',0,'numel',1);

issm/trunk-jpl/src/m/classes/SMBd18opdd.py

@@ -97 +97 @@
 
                         if self.isd18opd:
-                                md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1)
-                                md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1,'Inf',1,'timeseries',1)
+                                lent=float(np.size(self.temperatures_presentday,1))
+                                lenp=float(np.size(self.precipitations_presentday,1))
+                                multt=np.ceil(lent/12.)*12.
+                                multp=np.ceil(lenp/12.)*12.
+                                md = checkfield(md,'fieldname','smb.temperatures_presentday','size',[md.mesh.numberofvertices+1,multt],'NaN',1,'Inf',1,'timeseries',1)
+                                md = checkfield(md,'fieldname','smb.precipitations_presentday','size',[md.mesh.numberofvertices+1,multp],'NaN',1,'Inf',1,'timeseries',1)
                                 md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,np.nan],'singletimeseries',1)
                                 md = checkfield(md,'fieldname','smb.dpermil','>=',0,'numel',[1])