Changeset 27465

Timestamp:

12/21/22 05:46:25 (2 years ago)

Author:

vverjans

Message:

NEW allowing monthly varying lapse rates in SMBarma

Location:

issm/trunk-jpl

Files:

• src/c/classes/Elements/Element.cpp (modified) (2 diffs)
• src/m/classes/SMBarma.m (modified) (5 diffs)
• src/m/classes/SMBarma.py (modified) (7 diffs)
• test/Archives/Archive257.arch (modified) ( previous)
• test/NightlyRun/test257.m (modified) (3 diffs)
• test/NightlyRun/test257.py (modified) (3 diffs)

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

@@ -2431 +2431 @@
    const int numvertices = this->GetNumberOfVertices();
    bool isadjustsmb = false;
-        int basinid,bb1,bb2;
-   IssmDouble ela,refelevation_b;
+        int basinid,bb1,bb2,mindex;
+        IssmDouble ela,refelevation_b,time,dt,fracyear,yts;
+   IssmDouble monthsteps[12]  = {0.,1./12,2./12,3./12,4./12,5./12,6./12,7./12,8./12,9./12,10./12,11./12};
    IssmDouble* surfacelist  = xNew<IssmDouble>(numvertices);
    IssmDouble* smblist      = xNew<IssmDouble>(numvertices);
-   /* numelevbins values of lapse rates */
+   /* numelevbins values of lapse rates at current month */
         IssmDouble* lapserates_b = xNew<IssmDouble>(numelevbins);
-   /* (numelevbins-1) limits between elevation bins (be cautious with indexing) */
+   /* (numelevbins-1) limits between elevation bins at current month (be cautious with indexing) */
         IssmDouble* elevbins_b   = xNew<IssmDouble>(numelevbins-1);
+
+        /*Find month of current time step*/
+        this->parameters->FindParam(&yts,ConstantsYtsEnum);
+   this->parameters->FindParam(&time,TimeEnum);
+   this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); _assert_(dt<yts);
+   fracyear     = time/yts-floor(time/yts);
+   for(int i=1;i<12;i++){
+                if(fracyear>=monthsteps[i-1]) mindex = i-1;
+        }
+   if(fracyear>=monthsteps[11]) mindex = 11;
 
    /*Retrieve SMB values non-adjusted for SMB lapse rate*/
@@ -2447 +2458 @@
    this->GetInputValue(&basinid,SmbBasinsIdEnum);
    refelevation_b = refelevation[basinid];
-        for(int ii=0;ii<(numelevbins-1);ii++) elevbins_b[ii] = elevbins[basinid*(numelevbins-1)+ii];
+        /*Retrieve bins and laps rates for this basin at this month*/
+        for(int ii=0;ii<(numelevbins-1);ii++) elevbins_b[ii] = elevbins[basinid*(numelevbins-1)*12+mindex*(numelevbins-1)+ii];
         for(int ii=0;ii<numelevbins;ii++){
-                lapserates_b[ii] = lapserates[basinid*numelevbins+ii];
+                lapserates_b[ii] = lapserates[basinid*numelevbins*12+mindex*numelevbins+ii];
                 if(lapserates_b[ii]!=0) isadjustsmb=true;
         }
+
         /*Adjust SMB if any lapse rate value is non-zero*/
         if(isadjustsmb){

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

@@ -102 +102 @@
                                         md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins);
                                 end
-                                [nbas,nbins] = size(md.smb.lapserates);
-                                if (any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins]))==0) || numel(md.smb.lapserates)>1)
-                                        md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins],'numel',md.smb.num_basins*nbins);
-                                        md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins-1],'numel',md.smb.num_basins*(nbins-1));
+                                nbas     = size(md.smb.lapserates,1);
+                                nbins    = size(md.smb.lapserates,2);
+                                ntmlapse = size(md.smb.lapserates,3);
+                                if(ntmlapse>1 && ntmlapse~=12)
+                                        error('3rd dimension of md.smb.lapserates must be of size 1 or 12 (for monthly lapse rates)');
+                                end
+                                if (any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins*ntmlapse]))==0) || numel(md.smb.lapserates)>1)
+                                        md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins,ntmlapse],'numel',md.smb.num_basins*nbins*ntmlapse);
+                                        md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,max(1,nbins-1),ntmlapse],'numel',md.smb.num_basins*max(1,nbins-1)*ntmlapse);
                                         if(issorted(md.smb.elevationbins,2)==0)
                                                 error('md.smb.elevationbins should have rows in order of increasing elevation');
                                         end
-                                elseif (isnan(md.smb.elevationbins(1,1))==0 || numel(md.smb.elevationbins)>1)
+                                elseif (isnan(md.smb.elevationbins(1,1,1))==0 || numel(md.smb.elevationbins)>1)
                                         %elevationbins specified but not lapserates: this will inevitably lead to inconsistencies
-                                        [nbas,nbins] = size(md.smb.elevationbins);
-                                        nbins        = nbins+1;
-                                        md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins],'numel',md.smb.num_basins*nbins);
-                                        md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins-1],'numel',md.smb.num_basins*(nbins-1));
+                                        nbas     = size(md.smb.elevationbins,1);
+               nbins    = size(md.smb.elevationbins,2);
+                                        nbins    = nbins+1;
+               ntmlapse = size(md.smb.elevationbins,3);
+                                        md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,max(1,nbins-1),ntmlapse],'numel',md.smb.num_basins*nbins*ntmlapse);
+                                        md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,max(1,nbins-1),ntmlapse],'numel',md.smb.num_basins*max(1,nbins-1)*ntmlapse);
                                 end
                         end
@@ -135 +142 @@
                         fielddisplay(self,'arlag_coefs','basin-specific vectors of AR lag coefficients [unitless]');
                         fielddisplay(self,'malag_coefs','basin-specific vectors of MA lag coefficients [unitless]');
-                        fielddisplay(self,'lapserates','basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)');
-                        fielddisplay(self,'elevationbins','basin-specific separations between elevation bins, 1 row per basin, 1 column per limit between bins [m] (default: no basin separation)');
+                        fielddisplay(self,'lapserates','basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin, dimension 3 can be of size 12 to prescribe monthly varying values [m ice eq yr^-1 m^-1] (default: no lapse rate)');
+                        fielddisplay(self,'elevationbins','basin-specific separations between elevation bins, 1 row per basin, 1 column per limit between bins, dimension 3 can be of size 12 to prescribe monthly varying values [m] (default: no basin separation)');
                         fielddisplay(self,'refelevation','basin-specific reference elevations at which SMB is calculated, and from which SMB is downscaled using lapserates (default: basin mean elevation) [m]');
                         fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
@@ -156 +163 @@
                         tempelevationbins = md.smb.elevationbins;
                         temprefelevation  = md.smb.refelevation;
-                        [nbas,nbins]      = size(md.smb.lapserates);
-                        if(any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins]))))
-                                templapserates = zeros(md.smb.num_basins,2);
+                        nbas     = size(md.smb.lapserates,1);
+         nbins    = size(md.smb.lapserates,2);
+         ntmlapse = size(md.smb.lapserates,3);
+                        if(any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins*ntmlapse]))))
+                                templapserates = zeros(md.smb.num_basins,2,12);
                                 disp('      smb.lapserates not specified: set to 0');
-                           tempelevationbins = zeros(md.smb.num_basins,1); %dummy elevation bins
+                           tempelevationbins = zeros(md.smb.num_basins,1,12); %dummy elevation bins
+                        elseif(ntmlapse==1)
+                                templapserates    = repmat(templapserates,1,1,12); %same values each month
+                                tempelevationbins = repmat(tempelevationbins,1,1,12); %same values each month
                         end
                         if(any(isnan(md.smb.refelevation)))
@@ -173 +185 @@
                                         temprefelevation(ii) = sum(areas(indices).*elemsh)/sum(areas(indices));
                                 end
-                                if(any(reshape(md.smb.lapserates,[1,nbas*nbins])~=0))
+                                if(any(reshape(md.smb.lapserates,[1,nbas*nbins*12])~=0))
                                         disp('      smb.refelevation not specified: Reference elevations set to mean surface elevation of basins');
                                 end
                         end
-                        [nbas,nbins] = size(templapserates);
+                        nbas     = size(templapserates,1);
+         nbins    = size(templapserates,2);
+                        temp2dlapserates    = zeros(nbas,nbins*12);
+                        temp2delevationbins = zeros(nbas,max(1,nbins-1)*12);
+                        for(ii=[1:12])
+                                jj = 1+(ii-1)*nbins;
+                                temp2dlapserates(:,jj:jj+nbins-1)    = templapserates(:,:,ii);
+                                kk = 1+(ii-1)*(nbins-1);
+                                temp2delevationbins(:,kk:kk+nbins-2) = tempelevationbins(:,:,ii);
+                        end
 
                         % Scale the parameters %
@@ -234 +255 @@
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','malag_coefs','format','DoubleMat','name','md.smb.malag_coefs','yts',yts);
                         WriteData(fid,prefix,'data',dbreaks,'name','md.smb.datebreaks','format','DoubleMat','scale',yts);
-                        WriteData(fid,prefix,'data',templapserates,'format','DoubleMat','name','md.smb.lapserates','scale',1./yts,'yts',yts);
-                        WriteData(fid,prefix,'data',tempelevationbins,'format','DoubleMat','name','md.smb.elevationbins');
+                        WriteData(fid,prefix,'data',temp2dlapserates,'format','DoubleMat','name','md.smb.lapserates','scale',1./yts,'yts',yts);
+                        WriteData(fid,prefix,'data',temp2delevationbins,'format','DoubleMat','name','md.smb.elevationbins');
                         WriteData(fid,prefix,'data',temprefelevation,'format','DoubleMat','name','md.smb.refelevation');
                         WriteData(fid,prefix,'data',nbins,'format','Integer','name','md.smb.num_bins');

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

@@ -49 +49 @@
         s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of AR lag coefficients [unitless]'))
         s += '{}\n'.format(fielddisplay(self, 'malag_coefs', 'basin-specific vectors of MA lag coefficients [unitless]'))
-        s += '{}\n'.format(fielddisplay(self, 'lapserates', 'basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)'))
-        s += '{}\n'.format(fielddisplay(self, 'elevationbins', 'basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate)'))
+        s += '{}\n'.format(fielddisplay(self, 'lapserates', 'basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin, dimension 3 can be of size 12 to prescribe monthly varying values [m ice eq yr^-1 m^-1] (default: no lapse rate)'))
+        s += '{}\n'.format(fielddisplay(self, 'elevationbins', 'basin-specific separations between elevation bins, 1 row per basin, 1 column per limit between bins, dimension 3 can be of size 12 to prescribe monthly varying values [m] (default: no basin separation)'))
         s += '{}\n'.format(fielddisplay(self, 'refelevation', 'basin-specific reference elevations at which SMB is calculated, and from which SMB is downscaled using lapserates (default: basin mean elevation) [m]'))
         s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'))
@@ -127 +127 @@
 
             if(np.any(np.isnan(self.lapserates) is False) or np.size(self.lapserates) > 1):
+                nbas = md.smb.num_basins
                 if len(np.shape(self.lapserates)) == 1:
-                    self.lapserates = np.array([self.lapserates])
                     nbins = 1
-                else:
+                    self.lapserates = np.reshape(self.lapserates,[nbas,nbins,1])
+                elif(len(np.shape(self.lapserates)) == 2):
                     nbins = np.shape(self.lapserates)[1]
-                if len(np.shape(self.elevationbins)) == 1:
-                    self.elevationbins = np.array([self.elevationbins])
-                md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins], 'numel', md.smb.num_basins*nbins)
-                md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins-1], 'numel', md.smb.num_basins*(nbins-1))
+                    self.lapserates = np.reshape(self.lapserates,[nbas,nbins,1])
+                elif(len(np.shape(self.lapserates)) == 3):
+                    nbins = np.shape(self.lapserates)[1]
+                ntmlapse = np.shape(self.lapserates)[2]
+                if len(np.shape(self.elevationbins)) < 3:
+                    self.elevationbins = np.reshape(self.elevationbins,[nbas,max(1,nbins-1),ntmlapse])
+                md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [nbas,nbins,ntmlapse], 'numel', md.smb.num_basins*nbins*ntmlapse)
+                md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [nbas,max(1,nbins-1),ntmlapse], 'numel', md.smb.num_basins*max(1,nbins-1)*ntmlapse)
                 for rr in range(md.smb.num_basins):
                     if(np.all(self.elevationbins[rr,0:-1]<=self.elevationbins[rr,1:])==False):
@@ -141 +146 @@
             elif(np.any(np.isnan(self.elevationbins) is False) or np.size(self.elevationbins) > 1):
                 #elevationbins specified but not lapserates: this will inevitably lead to inconsistencies
+                nbas = md.smb.num_basins
                 if len(np.shape(self.elevationbins)) == 1:
-                    self.elevationbins = np.array([self.elevationbins])
                     nbins = 1
-                else:
-                    nbins = np.shape(self.elevationbins)[1]+1
-                md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins], 'numel', md.smb.num_basins*nbins)
-                md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins-1], 'numel', md.smb.num_basins*(nbins-1))
+                    self.elevationbins = np.reshape(self.elevationbins,[nbas,nbins,1])
+                elif(len(np.shape(self.lapserates)) == 2):
+                    nbins = np.shape(self.elevationbins)[1]
+                    self.elevationbins = np.reshape(self.elevationbins,[nbas,nbins,1])
+                elif(len(np.shape(self.lapserates)) == 3):
+                    nbins = np.shape(self.lapserates)[1]
+                nbins = nbins-1
+                ntmlapse = np.shape(self.lapserates)[2]
+                md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins*ntmlapse], 'numel', md.smb.num_basins*nbins*ntmlapse)
+                md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins,max(1,nbins-1)*ntmlapse], 'numel', md.smb.num_basins*max(1,nbins-1)*ntmlapse)
 
         md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1])
@@ -160 +171 @@
         nprm = md.smb.num_params;
         nper = md.smb.num_breaks+1;
-        templapserates    = np.copy(md.smb.lapserates)
-        tempelevationbins = np.copy(md.smb.elevationbins)
+        if(len(np.shape(md.smb.lapserates))==1):
+            nbins    = 1
+            ntmlapse = 1
+        elif(len(np.shape(md.smb.lapserates))==2):
+            nbins    = np.shape(md.smb.lapserates)[1]
+            ntmlapse = 1
+        elif(len(np.shape(md.smb.lapserates))==3):
+            nbins    = np.shape(md.smb.lapserates)[1]
+            ntmlapse = np.shape(md.smb.lapserates)[2]
+        templapserates    = np.reshape(md.smb.lapserates,[nbas,nbins,ntmlapse])
+        tempelevationbins = np.reshape(md.smb.elevationbins,[nbas,max(1,nbins-1),ntmlapse])
         temprefelevation  = np.copy(md.smb.refelevation)
         # Scale the parameters #
@@ -198 +218 @@
 
         if(np.any(np.isnan(md.smb.lapserates))):
-            templapserates = np.zeros((md.smb.num_basins,2))
+            templapserates = np.zeros((md.smb.num_basins,2,12))
             print('      smb.lapserates not specified: set to 0')
-            tempelevationbins = np.zeros((md.smb.num_basins,1)) #dummy elevation bins
+            tempelevationbins = np.zeros((md.smb.num_basins,1,12)) #dummy elevation bins
+        elif(ntmlapse==1):
+            templapserates    = np.repeat(templapserates,12,axis=2) 
+            tempelevationbins = np.repeat(tempelevationbins,12,axis=2) 
         if(np.any(np.isnan(md.smb.refelevation))):
             temprefelevation = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins)
@@ -213 +236 @@
                 print('      smb.refelevation not specified: Reference elevations set to mean surface elevation of basins')
         nbins = np.shape(templapserates)[1]
+        temp2dlapserates    = np.zeros((nbas,nbins*12))
+        temp2delevationbins = np.zeros((nbas,max(12,(nbins-1)*12)))
+        for ii in range(12):
+            temp2dlapserates[:,ii*nbins:(ii+1)*nbins] = templapserates[:,:,ii]
+            temp2delevationbins[:,ii*(nbins-1):(ii+1)*(nbins-1)] = tempelevationbins[:,:,ii]
 
         WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 13, 'format', 'Integer')
@@ -226 +254 @@
         WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'malag_coefs', 'format', 'DoubleMat', 'name', 'md.smb.malag_coefs', 'yts', yts)
         WriteData(fid, prefix, 'data', dbreaks, 'name', 'md.smb.datebreaks', 'format', 'DoubleMat','scale',yts)
-        WriteData(fid, prefix, 'data', templapserates, 'name', 'md.smb.lapserates', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)
-        WriteData(fid, prefix, 'data', tempelevationbins, 'name', 'md.smb.elevationbins', 'format', 'DoubleMat')
+        WriteData(fid, prefix, 'data', temp2dlapserates, 'name', 'md.smb.lapserates', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)
+        WriteData(fid, prefix, 'data', temp2delevationbins, 'name', 'md.smb.elevationbins', 'format', 'DoubleMat')
         WriteData(fid, prefix, 'data', temprefelevation, 'name', 'md.smb.refelevation', 'format', 'DoubleMat')
         WriteData(fid, prefix, 'data', nbins, 'name', 'md.smb.num_bins', 'format', 'Integer')

issm/trunk-jpl/test/NightlyRun/test257.m

@@ -41 +41 @@
 
 md.timestepping.start_time = 0;
-md.timestepping.time_step  = 1;
-md.timestepping.final_time = 7;
+md.timestepping.time_step  = 1/12;
+md.timestepping.final_time = 2;
 md.smb                     = SMBarma();
 md.smb.num_basins          = 3; %number of basins
@@ -55 +55 @@
 md.smb.arlag_coefs         = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]];
 md.smb.malag_coefs         = [1.0;0;0.2];
-md.smb.lapserates          = [0.01,0.0;0.01,-0.01;0.0,-0.01];
-md.smb.elevationbins       = [100;150;100];
+lm0                        = [1e-4*[1,-0.1,-1];1e-6*[1,-0.1,-1];1e-5*[1,-0.1,-1]];
+lm1                        = [1e-4*[2,-0.2,-2];1e-6*[2,-0.2,-2];1e-5*[2,-0.2,-2]];
+lm2                        = [1e-4*[3,-0.3,-3];1e-6*[3,-0.3,-3];1e-5*[3,-0.3,-3]];
+lm3                        = [1e-4*[4,-0.4,-4];1e-6*[4,-0.4,-4];1e-5*[4,-0.4,-4]];
+lm4                        = [1e-4*[5,-0.5,-5];1e-6*[5,-0.5,-5];1e-5*[5,-0.5,-5]];
+lm5                        = [1e-4*[6,-0.6,-6];1e-6*[6,-0.6,-6];1e-5*[6,-0.6,-6]];
+lm6                        = [1e-4*[7,-0.7,-7];1e-6*[7,-0.7,-7];1e-5*[7,-0.7,-7]];
+lm7                        = [1e-4*[8,-0.8,-8];1e-6*[8,-0.8,-8];1e-5*[8,-0.8,-8]];
+lm8                        = [1e-4*[9,-0.9,-9];1e-6*[9,-0.9,-9];1e-5*[9,-0.9,-9]];
+lm9                        = [1e-4*[10,-1,-10];1e-6*[10,-1.0,-10];1e-5*[10,-1.0,-10]];
+lm10                       = [1e-4*[11,-1.1,-11];1e-6*[11,-1.1,-11];1e-5*[11,-1.1,-11]];
+lm11                       = [1e-4*[12,-1.2,-12];1e-6*[12,-1.2,-12];1e-5*[12,-1.2,-12]];
+md.smb.lapserates          = cat(3,lm0,lm1,lm2,lm3,lm4,lm5,lm6,lm7,lm8,lm9,lm10,lm11);
+md.smb.elevationbins       = repmat([100,300;200,400;250,450],1,1,12);
 
 %Stochastic forcing
@@ -79 +91 @@
         (md.results.TransientSolution(1).IceVolume),...
         (md.results.TransientSolution(1).SmbMassBalance),...
-        (md.results.TransientSolution(2).Vx),...
-        (md.results.TransientSolution(2).Vy),...
-        (md.results.TransientSolution(2).Vel),...
-        (md.results.TransientSolution(2).Thickness),...
-        (md.results.TransientSolution(2).IceVolume),...
-        (md.results.TransientSolution(2).SmbMassBalance),...
-        (md.results.TransientSolution(7).Vx),...
-        (md.results.TransientSolution(7).Vy),...
-        (md.results.TransientSolution(7).Vel),...
-        (md.results.TransientSolution(7).Thickness),...
-        (md.results.TransientSolution(7).IceVolume),...
-        (md.results.TransientSolution(7).SmbMassBalance),...
+        (md.results.TransientSolution(12).Vx),...
+        (md.results.TransientSolution(12).Vy),...
+        (md.results.TransientSolution(12).Vel),...
+        (md.results.TransientSolution(12).Thickness),...
+        (md.results.TransientSolution(12).IceVolume),...
+        (md.results.TransientSolution(12).SmbMassBalance),...
+        (md.results.TransientSolution(24).Vx),...
+        (md.results.TransientSolution(24).Vy),...
+        (md.results.TransientSolution(24).Vel),...
+        (md.results.TransientSolution(24).Thickness),...
+        (md.results.TransientSolution(24).IceVolume),...
+        (md.results.TransientSolution(24).SmbMassBalance),...
         };

issm/trunk-jpl/test/NightlyRun/test257.py

@@ -46 +46 @@
 
 md.timestepping.start_time = 0
-md.timestepping.time_step = 1
-md.timestepping.final_time = 8
+md.timestepping.time_step = 1/12
+md.timestepping.final_time = 2
 md.smb = SMBarma()
 md.smb.num_basins = 3  # number of basins
@@ -60 +60 @@
 md.smb.arlag_coefs = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1], [0.4, -0.4, 0.1, -0.1]])
 md.smb.malag_coefs = np.array([[1.0],[0],[0.2]])
-md.smb.lapserates        = np.array([[0.01,0.0],[0.01,-0.01],[0.0,-0.01]])
-md.smb.elevationbins  = np.array([100,150,100]).reshape(md.smb.num_basins,1)
+
+lm0                   = np.array([1e-4*np.array([1,-0.1,-1]),1e-6*np.array([1,-0.1,-1]),1e-5*np.array([1,-0.1,-1])])
+lm1                   = np.array([1e-4*np.array([2,-0.2,-2]),1e-6*np.array([2,-0.2,-2]),1e-5*np.array([2,-0.2,-2])])
+lm2                   = np.array([1e-4*np.array([3,-0.3,-3]),1e-6*np.array([3,-0.3,-3]),1e-5*np.array([3,-0.3,-3])])
+lm3                   = np.array([1e-4*np.array([4,-0.4,-4]),1e-6*np.array([4,-0.4,-4]),1e-5*np.array([4,-0.4,-4])])
+lm4                   = np.array([1e-4*np.array([5,-0.5,-5]),1e-6*np.array([5,-0.5,-5]),1e-5*np.array([5,-0.5,-5])])
+lm5                   = np.array([1e-4*np.array([6,-0.6,-6]),1e-6*np.array([6,-0.6,-6]),1e-5*np.array([6,-0.6,-6])])
+lm6                   = np.array([1e-4*np.array([7,-0.7,-7]),1e-6*np.array([7,-0.7,-7]),1e-5*np.array([7,-0.7,-7])])
+lm7                   = np.array([1e-4*np.array([8,-0.8,-8]),1e-6*np.array([8,-0.8,-8]),1e-5*np.array([8,-0.8,-8])])
+lm8                   = np.array([1e-4*np.array([9,-0.9,-9]),1e-6*np.array([9,-0.9,-9]),1e-5*np.array([9,-0.9,-9])])
+lm9                   = np.array([1e-4*np.array([10,-1.0,-10]),1e-6*np.array([10,-1.0,-10]),1e-5*np.array([10,-1.0,-10])])
+lm10                  = np.array([1e-4*np.array([11,-1.1,-11]),1e-6*np.array([11,-1.1,-11]),1e-5*np.array([11,-1.1,-11])])
+lm11                  = np.array([1e-4*np.array([12,-1.2,-12]),1e-6*np.array([12,-1.2,-12]),1e-5*np.array([12,-1.2,-12])])
+md.smb.lapserates     = np.stack((lm0,lm1,lm2,lm3,lm4,lm5,lm6,lm7,lm8,lm9,lm10,lm11),axis=2)
+ebins                 = np.array([[100,300],[200,400],[250,450]])
+md.smb.elevationbins  = np.stack([ebins for ii in range(12)],axis=2)
+
 
 # Stochastic forcing
@@ -90 +105 @@
     md.results.TransientSolution[0].IceVolume,
     md.results.TransientSolution[0].SmbMassBalance,
-    md.results.TransientSolution[1].Vx,
-    md.results.TransientSolution[1].Vy,
-    md.results.TransientSolution[1].Vel,
-    md.results.TransientSolution[1].Thickness,
-    md.results.TransientSolution[1].IceVolume,
-    md.results.TransientSolution[1].SmbMassBalance,
-    md.results.TransientSolution[6].Vx,
-    md.results.TransientSolution[6].Vy,
-    md.results.TransientSolution[6].Vel,
-    md.results.TransientSolution[6].Thickness,
-    md.results.TransientSolution[6].IceVolume,
-    md.results.TransientSolution[6].SmbMassBalance
+    md.results.TransientSolution[11].Vx,
+    md.results.TransientSolution[11].Vy,
+    md.results.TransientSolution[11].Vel,
+    md.results.TransientSolution[11].Thickness,
+    md.results.TransientSolution[11].IceVolume,
+    md.results.TransientSolution[11].SmbMassBalance,
+    md.results.TransientSolution[23].Vx,
+    md.results.TransientSolution[23].Vy,
+    md.results.TransientSolution[23].Vel,
+    md.results.TransientSolution[23].Thickness,
+    md.results.TransientSolution[23].IceVolume,
+    md.results.TransientSolution[23].SmbMassBalance
 ]