Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27464) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27465) @@ -2431,12 +2431,23 @@ 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(numvertices); IssmDouble* smblist = xNew(numvertices); - /* numelevbins values of lapse rates */ + /* numelevbins values of lapse rates at current month */ IssmDouble* lapserates_b = xNew(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(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=monthsteps[i-1]) mindex = i-1; + } + if(fracyear>=monthsteps[11]) mindex = 11; /*Retrieve SMB values non-adjusted for SMB lapse rate*/ @@ -2447,9 +2458,11 @@ 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=',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,6 +142,6 @@ 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,9 +163,14 @@ 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,9 +185,18 @@ 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,6 +255,6 @@ 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'); Index: /issm/trunk-jpl/src/m/classes/SMBarma.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBarma.py (revision 27464) +++ /issm/trunk-jpl/src/m/classes/SMBarma.py (revision 27465) @@ -49,6 +49,6 @@ 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,13 +127,18 @@ 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,11 +146,17 @@ 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,6 +171,15 @@ 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,7 +218,10 @@ 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,4 +236,9 @@ 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,6 +254,6 @@ 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') Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 27464) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 27465) @@ -41,6 +41,6 @@ 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,6 +55,18 @@ 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,15 +91,15 @@ (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),... }; Index: /issm/trunk-jpl/test/NightlyRun/test257.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.py (revision 27464) +++ /issm/trunk-jpl/test/NightlyRun/test257.py (revision 27465) @@ -46,6 +46,6 @@ 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,6 +60,21 @@ 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,15 +105,15 @@ 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 ]