Index: /issm/trunk-jpl/src/m/classes/SMBarma.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBarma.m (revision 27482) +++ /issm/trunk-jpl/src/m/classes/SMBarma.m (revision 27483) @@ -10,16 +10,16 @@ num_params = 0; arma_timestep = 0; - ar_order = 0; + ar_order = 0; arlag_coefs = NaN; - ma_order = 0; + ma_order = 0; malag_coefs = NaN; polynomialparams = NaN; datebreaks = NaN; - basin_id = NaN; + basin_id = NaN; lapserates = NaN; elevationbins = NaN; refelevation = NaN; steps_per_step = 1; - averaging = 0; + averaging = 0; requested_outputs = {}; end @@ -47,5 +47,5 @@ if (self.ma_order==0) self.ma_order = 1; %dummy 1 value for moving-average - self.arlag_coefs = zeros(self.num_basins,self.ma_order); %moving-average coefficients all set to 0 + self.malag_coefs = zeros(self.num_basins,self.ma_order); %moving-average coefficients all set to 0 disp(' smb.ma_order (order of moving-average model) not specified: order of moving-average model set to 0'); end @@ -76,5 +76,5 @@ md = checkfield(md,'fieldname','smb.num_params','numel',1,'NaN',1,'Inf',1,'>',0); md = checkfield(md,'fieldname','smb.num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0); - md = checkfield(md,'fieldname','smb.basin_id','Inf',1,'>=',0,'<=',md.smb.num_basins,'size',[md.mesh.numberofelements,1]); + md = checkfield(md,'fieldname','smb.basin_id','Inf',1,'>=',0,'<=',nbas,'size',[md.mesh.numberofelements,1]); if(nbas>1 && nbrk>=1 && nprm>1) md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm); @@ -89,6 +89,6 @@ md = checkfield(md,'fieldname','smb.ma_order','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','smb.arma_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %arma time step cannot be finer than ISSM timestep - md = checkfield(md,'fieldname','smb.arlag_coefs','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]); - md = checkfield(md,'fieldname','smb.malag_coefs','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ma_order]); + md = checkfield(md,'fieldname','smb.arlag_coefs','NaN',1,'Inf',1,'size',[nbas,md.smb.ar_order]); + md = checkfield(md,'fieldname','smb.malag_coefs','NaN',1,'Inf',1,'size',[nbas,md.smb.ma_order]); if(nbrk>0) @@ -100,5 +100,5 @@ end if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1) - md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); + md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,nbas],'numel',nbas); end nbas = size(md.smb.lapserates,1); @@ -109,6 +109,6 @@ 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); + md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[nbas,nbins,ntmlapse],'numel',nbas*nbins*ntmlapse); + md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[nbas,max(1,nbins-1),ntmlapse],'numel',nbas*max(1,nbins-1)*ntmlapse); if(issorted(md.smb.elevationbins,2)==0) error('md.smb.elevationbins should have rows in order of increasing elevation'); @@ -117,9 +117,9 @@ %elevationbins specified but not lapserates: this will inevitably lead to inconsistencies nbas = size(md.smb.elevationbins,1); - nbins = size(md.smb.elevationbins,2); + 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); + ntmlapse = size(md.smb.elevationbins,3); + md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[nbas,max(1,nbins-1),ntmlapse],'numel',nbas*nbins*ntmlapse); + md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[nbas,max(1,nbins-1),ntmlapse],'numel',nbas*max(1,nbins-1)*ntmlapse); end end @@ -164,6 +164,6 @@ temprefelevation = md.smb.refelevation; nbas = size(md.smb.lapserates,1); - nbins = size(md.smb.lapserates,2); - ntmlapse = size(md.smb.lapserates,3); + 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); @@ -175,6 +175,6 @@ end nbas = size(templapserates,1); - nbins = size(templapserates,2); - ntmlapse = size(templapserates,3); + nbins = size(templapserates,2); + ntmlapse = size(templapserates,3); if(any(isnan(md.smb.refelevation))) temprefelevation = zeros(1,md.smb.num_basins); @@ -205,42 +205,42 @@ polyparams2dScaled = zeros(nbas,nper*nprm); if(nprm>1) - % Case 3D % - if(nbas>1 && nper>1) - for(ii=[1:nprm]) - polyparamsScaled(:,:,ii) = polyparamsScaled(:,:,ii)*((1/yts)^(ii)); - end - % Fit in 2D array % - for(ii=[1:nprm]) - jj = 1+(ii-1)*nper; - polyparams2dScaled(:,jj:jj+nper-1) = polyparamsScaled(:,:,ii); - end - % Case 2D and higher-order params at increasing row index % - elseif(nbas==1) - for(ii=[1:nprm]) - polyparamsScaled(ii,:) = polyparamsScaled(ii,:)*((1/yts)^(ii)); - end - % Fit in row array % - for(ii=[1:nprm]) - jj = 1+(ii-1)*nper; - polyparams2dScaled(1,jj:jj+nper-1) = polyparamsScaled(ii,:); - end - % Case 2D and higher-order params at incrasing column index % - elseif(nper==1) - for(ii=[1:nprm]) - polyparamsScaled(:,ii) = polyparamsScaled(:,ii)*((1/yts)^(ii)); - end - % 2D array is already in correct format % - polyparams2dScaled = polyparamsScaled; - end - else + % Case 3D % + if(nbas>1 && nper>1) + for(ii=[1:nprm]) + polyparamsScaled(:,:,ii) = polyparamsScaled(:,:,ii)*((1/yts)^(ii)); + end + % Fit in 2D array % + for(ii=[1:nprm]) + jj = 1+(ii-1)*nper; + polyparams2dScaled(:,jj:jj+nper-1) = polyparamsScaled(:,:,ii); + end + % Case 2D and higher-order params at increasing row index % + elseif(nbas==1) + for(ii=[1:nprm]) + polyparamsScaled(ii,:) = polyparamsScaled(ii,:)*((1/yts)^(ii)); + end + % Fit in row array % + for(ii=[1:nprm]) + jj = 1+(ii-1)*nper; + polyparams2dScaled(1,jj:jj+nper-1) = polyparamsScaled(ii,:); + end + % Case 2D and higher-order params at incrasing column index % + elseif(nper==1) + for(ii=[1:nprm]) + polyparamsScaled(:,ii) = polyparamsScaled(:,ii)*((1/yts)^(ii)); + end + % 2D array is already in correct format % + polyparams2dScaled = polyparamsScaled; + end + else polyparamsScaled = polyparamsScaled*(1/yts); - % 2D array is already in correct format % - polyparams2dScaled = polyparamsScaled; - end + % 2D array is already in correct format % + polyparams2dScaled = polyparamsScaled; + end if(nper==1) %a single period (no break date) - dbreaks = zeros(nbas,1); %dummy - else - dbreaks = md.smb.datebreaks; - end + dbreaks = zeros(nbas,1); %dummy + else + dbreaks = md.smb.datebreaks; + end WriteData(fid,prefix,'name','md.smb.model','data',13,'format','Integer'); @@ -267,6 +267,6 @@ pos = find(ismember(outputs,'default')); if ~isempty(pos), - outputs(pos) = []; %remove 'default' from outputs - outputs = [outputs defaultoutputs(self,md)]; %add defaults + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(self,md)]; %add defaults end WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray'); Index: /issm/trunk-jpl/src/m/classes/SMBarma.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBarma.py (revision 27482) +++ /issm/trunk-jpl/src/m/classes/SMBarma.py (revision 27483) @@ -3,7 +3,7 @@ from checkfield import * from fielddisplay import fielddisplay +from GetAreas import * from project3d import * from WriteData import * -from GetAreas import * class SMBarma(object): @@ -21,5 +21,8 @@ self.ar_order = 0 self.arlag_coefs = np.nan + self.ma_order = 0 self.malag_coefs = np.nan + self.polynomialparams = np.nan + self.datebreaks = np.nan self.basin_id = np.nan self.lapserates = np.nan @@ -79,4 +82,8 @@ self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autoregression coefficients all set to 0 print(' smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0') + if self.ma_order == 0: + self.ma_order = 1 # Dummy 1 value for moving-average + self.malag_coefs = np.zeros((self.num_basins, self.ma_order)) # Moving-average coefficients all set to 0 + print(' smb.ma_order (order of moving-average model) not specified: order of moving-average model set to 0') if self.arma_timestep == 0: self.arma_timestep = md.timestepping.time_step # ARMA model has no prescribed time step @@ -92,39 +99,44 @@ def checkconsistency(self, md, solution, analyses): # {{{ + """ + TODO: + - Ensure that checks on shape of self.lapserates are same as those under MATLAB as matrix addressing is quite different here + """ if 'MasstransportAnalysis' in analyses: - nbas = md.smb.num_basins; - nprm = md.smb.num_params; - nbrk = md.smb.num_breaks; + nbas = md.smb.num_basins + nprm = md.smb.num_params + nbrk = md.smb.num_breaks md = checkfield(md, 'fieldname', 'smb.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) md = checkfield(md, 'fieldname', 'smb.num_params', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) md = checkfield(md, 'fieldname', 'smb.num_breaks', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) md = checkfield(md, 'fieldname', 'smb.basin_id', 'Inf', 1, '>=', 0, '<=', md.smb.num_basins, 'size', [md.mesh.numberofelements]) - if len(np.shape(self.polynomialparams)) == 1: - self.polynomialparams = np.array([[self.polynomialparams]]) - if(nbas>1 and nbrk>=1 and nprm>1): - md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm) - elif(nbas==1): - md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm) - elif(nbrk==0): - md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm) - elif(nprm==1): - md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk],'numel',nbas*(nbrk+1)*nprm) + # if len(np.shape(self.polynomialparams)) == 1: + # self.polynomialparams = np.array([[self.polynomialparams]]) + if nbas > 1 and nbrk >= 1 and nprm > 1: + md = checkfield(md, 'fieldname', 'smb.polynomialparams', 'NaN', 1, 'Inf', 1, 'size', [nbas, nbrk + 1, nprm], 'numel', nbas * (nbrk + 1) * nprm) + elif nbas == 1: + md = checkfield(md, 'fieldname', 'smb.polynomialparams', 'NaN', 1, 'Inf', 1, 'size', [nprm, nbrk + 1], 'numel', nbas * (nbrk + 1) * nprm) + elif nbrk == 0: + md = checkfield(md, 'fieldname', 'smb.polynomialparams', 'NaN', 1, 'Inf', 1, 'size', [nbas, nprm], 'numel', nbas * (nbrk + 1) * nprm) + elif nprm == 1: + md = checkfield(md, 'fieldname', 'smb.polynomialparams', 'NaN', 1, 'Inf', 1, 'size', [nbas, nbrk], 'numel', nbas * (nbrk + 1) * nprm) md = checkfield(md, 'fieldname', 'smb.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) + md = checkfield(md, 'fieldname', 'smb.ma_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) md = checkfield(md, 'fieldname', 'smb.arma_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep - md = checkfield(md, 'fieldname', 'smb.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order]) - md = checkfield(md, 'fieldname', 'smb.malag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ma_order]) - if(nbrk>0): + md = checkfield(md, 'fieldname', 'smb.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [nbas, md.smb.ar_order]) + md = checkfield(md, 'fieldname', 'smb.malag_coefs', 'NaN', 1, 'Inf', 1, 'size', [nbas, md.smb.ma_order]) + if nbrk > 0: md = checkfield(md, 'fieldname', 'smb.datebreaks', 'NaN', 1, 'Inf', 1, 'size', [nbas,nbrk]) - elif(np.size(md.smb.datebreaks)==0 or np.all(np.isnan(md.smb.datebreaks))): + elif np.size(md.smb.datebreaks) == 0 or np.all(np.isnan(md.smb.datebreaks)): pass else: raise RuntimeError('md.smb.num_breaks is 0 but md.smb.datebreaks is not empty') - if(np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1): + if np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1: if len(np.shape(self.refelevation)) == 1: self.refelevation = np.array([self.refelevation]) - md = checkfield(md, 'fieldname', 'smb.refelevation', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) - - if(np.any(np.isnan(self.lapserates) is False) or np.size(self.lapserates) > 1): + md = checkfield(md, 'fieldname', 'smb.refelevation', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [1, nbas], 'numel', nbas) + + 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: @@ -140,10 +152,10 @@ 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): + md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [nbas,max(1,nbins-1),ntmlapse], 'numel', nbas*max(1,nbins-1)*ntmlapse) + for rr in range(nbas): if(np.all(self.elevationbins[rr,0:-1]<=self.elevationbins[rr,1:])==False): raise TypeError('md.smb.elevationbins should have rows in order of increasing elevation') - 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 + 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: @@ -155,8 +167,8 @@ elif(len(np.shape(self.lapserates)) == 3): nbins = np.shape(self.lapserates)[1] - nbins = nbins-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.lapserates', 'NaN', 1, 'Inf', 1, 'size', [nbas, nbins * ntmlapse], 'numel', nbas * nbins * ntmlapse) + md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [nbas, max(1, nbins - 1) * ntmlapse], 'numel', nbas * max(1, nbins - 1) * ntmlapse) md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) @@ -168,65 +180,65 @@ def marshall(self, prefix, md, fid): # {{{ yts = md.constants.yts - nbas = md.smb.num_basins; - nprm = md.smb.num_params; - nper = md.smb.num_breaks+1; + nbas = md.smb.num_basins + nprm = md.smb.num_params + nper = md.smb.num_breaks + 1 if(np.any(np.isnan(md.smb.lapserates))): - templapserates = np.zeros((md.smb.num_basins,2,12)) + templapserates = np.zeros((nbas, 2, 12)) print(' smb.lapserates not specified: set to 0') - tempelevationbins = np.zeros((md.smb.num_basins,1,12)) #dummy elevation bins + tempelevationbins = np.zeros((nbas, 1, 12)) # Dummy elevation bins nbins = 2 ntmlapse = 12 else: - if(len(np.shape(md.smb.lapserates))==1): + if len(np.shape(md.smb.lapserates)) == 1: nbins = 1 ntmlapse = 1 - elif(len(np.shape(md.smb.lapserates))==2): + 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): + 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]) + 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 # + # Scale the parameters polyparamsScaled = np.copy(md.smb.polynomialparams) - polyparams2dScaled = np.zeros((nbas,nper*nprm)) - if(nprm>1): - # Case 3D # - if(nbas>1 and nper>1): - for ii in range(nprm): - polyparamsScaled[:,:,ii] = polyparamsScaled[:,:,ii]*(1/yts)**(ii+1) - # Fit in 2D array # - for ii in range(nprm): - polyparams2dScaled[:,ii*nper:(ii+1)*nper] = 1*polyparamsScaled[:,:,ii] - # Case 2D and higher-order params at increasing row index # - elif(nbas==1): - for ii in range(nprm): - polyparamsScaled[ii,:] = polyparamsScaled[ii,:]*(1/yts)**(ii+1) - # Fit in row array # - for ii in range(nprm): - polyparams2dScaled[0,ii*nper:(ii+1)*nper] = 1*polyparamsScaled[ii,:] - # Case 2D and higher-order params at incrasing column index # - elif(nper==1): - for ii in range(nprm): - polyparamsScaled[:,ii] = polyparamsScaled[:,ii]*(1/yts)**(ii+1) - # 2D array is already in correct format # + polyparams2dScaled = np.zeros((nbas, nper * nprm)) + if nprm > 1: + # Case 3D + if nbas > 1 and nper > 1: + for ii in range(nprm): + polyparamsScaled[:, :, ii] = polyparamsScaled[:, :, ii] * (1 / yts) ** (ii + 1) + # Fit in 2D array + for ii in range(nprm): + polyparams2dScaled[:, ii * nper:(ii + 1) * nper] = 1 * polyparamsScaled[:, :, ii] + # Case 2D and higher-order params at increasing row index + elif nbas == 1: + for ii in range(nprm): + polyparamsScaled[ii, :] = polyparamsScaled[ii, :] * (1 / yts) ** (ii + 1) + # Fit in row array + for ii in range(nprm): + polyparams2dScaled[0, ii * nper:(ii + 1) * nper] = 1 * polyparamsScaled[ii, :] + # Case 2D and higher-order params at increasing column index + elif nper == 1: + for ii in range(nprm): + polyparamsScaled[:, ii] = polyparamsScaled[:, ii] * (1 / yts) ** (ii + 1) + # 2D array is already in correct format polyparams2dScaled = np.copy(polyparamsScaled) else: - polyparamsScaled = polyparamsScaled*(1/yts) - # 2D array is already in correct format # + polyparamsScaled = polyparamsScaled * (1 / yts) + # 2D array is already in correct format polyparams2dScaled = np.copy(polyparamsScaled) - - if(nper==1): - dbreaks = np.zeros((nbas,1)) + + if nper == 1: + dbreaks = np.zeros((nbas, 1)) else: dbreaks = np.copy(md.smb.datebreaks) - if(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) + if 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((nbas)).reshape(1, nbas) areas = GetAreas(md.mesh.elements, md.mesh.x, md.mesh.y) for ii, bid in enumerate(np.unique(md.smb.basin_id)): @@ -239,9 +251,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))) + 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] + 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')