![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 27483
- Timestamp:
- 12/28/22 10:51:15 (2 years ago)
- Location:
- issm/trunk-jpl/src/m/classes
- Files:
-
- 2 edited
-
SMBarma.m (modified) (11 diffs)
-
SMBarma.py (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/classes/SMBarma.m
r27467 r27483 10 10 num_params = 0; 11 11 arma_timestep = 0; 12 ar_order 12 ar_order = 0; 13 13 arlag_coefs = NaN; 14 ma_order 14 ma_order = 0; 15 15 malag_coefs = NaN; 16 16 polynomialparams = NaN; 17 17 datebreaks = NaN; 18 basin_id 18 basin_id = NaN; 19 19 lapserates = NaN; 20 20 elevationbins = NaN; 21 21 refelevation = NaN; 22 22 steps_per_step = 1; 23 averaging 23 averaging = 0; 24 24 requested_outputs = {}; 25 25 end … … 47 47 if (self.ma_order==0) 48 48 self.ma_order = 1; %dummy 1 value for moving-average 49 self. arlag_coefs = zeros(self.num_basins,self.ma_order); %moving-average coefficients all set to 049 self.malag_coefs = zeros(self.num_basins,self.ma_order); %moving-average coefficients all set to 0 50 50 disp(' smb.ma_order (order of moving-average model) not specified: order of moving-average model set to 0'); 51 51 end … … 76 76 md = checkfield(md,'fieldname','smb.num_params','numel',1,'NaN',1,'Inf',1,'>',0); 77 77 md = checkfield(md,'fieldname','smb.num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0); 78 md = checkfield(md,'fieldname','smb.basin_id','Inf',1,'>=',0,'<=', md.smb.num_basins,'size',[md.mesh.numberofelements,1]);78 md = checkfield(md,'fieldname','smb.basin_id','Inf',1,'>=',0,'<=',nbas,'size',[md.mesh.numberofelements,1]); 79 79 if(nbas>1 && nbrk>=1 && nprm>1) 80 80 md = checkfield(md,'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm); … … 89 89 md = checkfield(md,'fieldname','smb.ma_order','numel',1,'NaN',1,'Inf',1,'>=',0); 90 90 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 91 md = checkfield(md,'fieldname','smb.arlag_coefs','NaN',1,'Inf',1,'size',[ md.smb.num_basins,md.smb.ar_order]);92 md = checkfield(md,'fieldname','smb.malag_coefs','NaN',1,'Inf',1,'size',[ md.smb.num_basins,md.smb.ma_order]);91 md = checkfield(md,'fieldname','smb.arlag_coefs','NaN',1,'Inf',1,'size',[nbas,md.smb.ar_order]); 92 md = checkfield(md,'fieldname','smb.malag_coefs','NaN',1,'Inf',1,'size',[nbas,md.smb.ma_order]); 93 93 94 94 if(nbrk>0) … … 100 100 end 101 101 if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1) 102 md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1, md.smb.num_basins],'numel',md.smb.num_basins);102 md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,nbas],'numel',nbas); 103 103 end 104 104 nbas = size(md.smb.lapserates,1); … … 109 109 end 110 110 if (any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins*ntmlapse]))==0) || numel(md.smb.lapserates)>1) 111 md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[ md.smb.num_basins,nbins,ntmlapse],'numel',md.smb.num_basins*nbins*ntmlapse);112 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);111 md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[nbas,nbins,ntmlapse],'numel',nbas*nbins*ntmlapse); 112 md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[nbas,max(1,nbins-1),ntmlapse],'numel',nbas*max(1,nbins-1)*ntmlapse); 113 113 if(issorted(md.smb.elevationbins,2)==0) 114 114 error('md.smb.elevationbins should have rows in order of increasing elevation'); … … 117 117 %elevationbins specified but not lapserates: this will inevitably lead to inconsistencies 118 118 nbas = size(md.smb.elevationbins,1); 119 119 nbins = size(md.smb.elevationbins,2); 120 120 nbins = nbins+1; 121 122 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);123 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);121 ntmlapse = size(md.smb.elevationbins,3); 122 md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[nbas,max(1,nbins-1),ntmlapse],'numel',nbas*nbins*ntmlapse); 123 md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[nbas,max(1,nbins-1),ntmlapse],'numel',nbas*max(1,nbins-1)*ntmlapse); 124 124 end 125 125 end … … 164 164 temprefelevation = md.smb.refelevation; 165 165 nbas = size(md.smb.lapserates,1); 166 167 166 nbins = size(md.smb.lapserates,2); 167 ntmlapse = size(md.smb.lapserates,3); 168 168 if(any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins*ntmlapse])))) 169 169 templapserates = zeros(md.smb.num_basins,2,12); … … 175 175 end 176 176 nbas = size(templapserates,1); 177 178 177 nbins = size(templapserates,2); 178 ntmlapse = size(templapserates,3); 179 179 if(any(isnan(md.smb.refelevation))) 180 180 temprefelevation = zeros(1,md.smb.num_basins); … … 205 205 polyparams2dScaled = zeros(nbas,nper*nprm); 206 206 if(nprm>1) 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 207 % Case 3D % 208 if(nbas>1 && nper>1) 209 for(ii=[1:nprm]) 210 polyparamsScaled(:,:,ii) = polyparamsScaled(:,:,ii)*((1/yts)^(ii)); 211 end 212 % Fit in 2D array % 213 for(ii=[1:nprm]) 214 jj = 1+(ii-1)*nper; 215 polyparams2dScaled(:,jj:jj+nper-1) = polyparamsScaled(:,:,ii); 216 end 217 % Case 2D and higher-order params at increasing row index % 218 elseif(nbas==1) 219 for(ii=[1:nprm]) 220 polyparamsScaled(ii,:) = polyparamsScaled(ii,:)*((1/yts)^(ii)); 221 end 222 % Fit in row array % 223 for(ii=[1:nprm]) 224 jj = 1+(ii-1)*nper; 225 polyparams2dScaled(1,jj:jj+nper-1) = polyparamsScaled(ii,:); 226 end 227 % Case 2D and higher-order params at incrasing column index % 228 elseif(nper==1) 229 for(ii=[1:nprm]) 230 polyparamsScaled(:,ii) = polyparamsScaled(:,ii)*((1/yts)^(ii)); 231 end 232 % 2D array is already in correct format % 233 polyparams2dScaled = polyparamsScaled; 234 end 235 else 236 236 polyparamsScaled = polyparamsScaled*(1/yts); 237 238 239 237 % 2D array is already in correct format % 238 polyparams2dScaled = polyparamsScaled; 239 end 240 240 if(nper==1) %a single period (no break date) 241 242 243 244 241 dbreaks = zeros(nbas,1); %dummy 242 else 243 dbreaks = md.smb.datebreaks; 244 end 245 245 246 246 WriteData(fid,prefix,'name','md.smb.model','data',13,'format','Integer'); … … 267 267 pos = find(ismember(outputs,'default')); 268 268 if ~isempty(pos), 269 outputs(pos) = []; 270 outputs = [outputs defaultoutputs(self,md)]; 269 outputs(pos) = []; %remove 'default' from outputs 270 outputs = [outputs defaultoutputs(self,md)]; %add defaults 271 271 end 272 272 WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray'); -
issm/trunk-jpl/src/m/classes/SMBarma.py
r27467 r27483 3 3 from checkfield import * 4 4 from fielddisplay import fielddisplay 5 from GetAreas import * 5 6 from project3d import * 6 7 from WriteData import * 7 from GetAreas import *8 8 9 9 class SMBarma(object): … … 21 21 self.ar_order = 0 22 22 self.arlag_coefs = np.nan 23 self.ma_order = 0 23 24 self.malag_coefs = np.nan 25 self.polynomialparams = np.nan 26 self.datebreaks = np.nan 24 27 self.basin_id = np.nan 25 28 self.lapserates = np.nan … … 79 82 self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autoregression coefficients all set to 0 80 83 print(' smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0') 84 if self.ma_order == 0: 85 self.ma_order = 1 # Dummy 1 value for moving-average 86 self.malag_coefs = np.zeros((self.num_basins, self.ma_order)) # Moving-average coefficients all set to 0 87 print(' smb.ma_order (order of moving-average model) not specified: order of moving-average model set to 0') 81 88 if self.arma_timestep == 0: 82 89 self.arma_timestep = md.timestepping.time_step # ARMA model has no prescribed time step … … 92 99 93 100 def checkconsistency(self, md, solution, analyses): # {{{ 101 """ 102 TODO: 103 - Ensure that checks on shape of self.lapserates are same as those under MATLAB as matrix addressing is quite different here 104 """ 94 105 if 'MasstransportAnalysis' in analyses: 95 nbas = md.smb.num_basins ;96 nprm = md.smb.num_params ;97 nbrk = md.smb.num_breaks ;106 nbas = md.smb.num_basins 107 nprm = md.smb.num_params 108 nbrk = md.smb.num_breaks 98 109 md = checkfield(md, 'fieldname', 'smb.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) 99 110 md = checkfield(md, 'fieldname', 'smb.num_params', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) 100 111 md = checkfield(md, 'fieldname', 'smb.num_breaks', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) 101 112 md = checkfield(md, 'fieldname', 'smb.basin_id', 'Inf', 1, '>=', 0, '<=', md.smb.num_basins, 'size', [md.mesh.numberofelements]) 102 if len(np.shape(self.polynomialparams)) == 1:103 self.polynomialparams = np.array([[self.polynomialparams]])104 if (nbas>1 and nbrk>=1 and nprm>1):105 md = checkfield(md, 'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm)106 elif (nbas==1):107 md = checkfield(md, 'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm)108 elif (nbrk==0):109 md = checkfield(md, 'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm)110 elif (nprm==1):111 md = checkfield(md, 'fieldname','smb.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk],'numel',nbas*(nbrk+1)*nprm)113 # if len(np.shape(self.polynomialparams)) == 1: 114 # self.polynomialparams = np.array([[self.polynomialparams]]) 115 if nbas > 1 and nbrk >= 1 and nprm > 1: 116 md = checkfield(md, 'fieldname', 'smb.polynomialparams', 'NaN', 1, 'Inf', 1, 'size', [nbas, nbrk + 1, nprm], 'numel', nbas * (nbrk + 1) * nprm) 117 elif nbas == 1: 118 md = checkfield(md, 'fieldname', 'smb.polynomialparams', 'NaN', 1, 'Inf', 1, 'size', [nprm, nbrk + 1], 'numel', nbas * (nbrk + 1) * nprm) 119 elif nbrk == 0: 120 md = checkfield(md, 'fieldname', 'smb.polynomialparams', 'NaN', 1, 'Inf', 1, 'size', [nbas, nprm], 'numel', nbas * (nbrk + 1) * nprm) 121 elif nprm == 1: 122 md = checkfield(md, 'fieldname', 'smb.polynomialparams', 'NaN', 1, 'Inf', 1, 'size', [nbas, nbrk], 'numel', nbas * (nbrk + 1) * nprm) 112 123 md = checkfield(md, 'fieldname', 'smb.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) 124 md = checkfield(md, 'fieldname', 'smb.ma_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) 113 125 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 114 md = checkfield(md, 'fieldname', 'smb.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [ md.smb.num_basins, md.smb.ar_order])115 md = checkfield(md, 'fieldname', 'smb.malag_coefs', 'NaN', 1, 'Inf', 1, 'size', [ md.smb.num_basins, md.smb.ma_order])116 if (nbrk>0):126 md = checkfield(md, 'fieldname', 'smb.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [nbas, md.smb.ar_order]) 127 md = checkfield(md, 'fieldname', 'smb.malag_coefs', 'NaN', 1, 'Inf', 1, 'size', [nbas, md.smb.ma_order]) 128 if nbrk > 0: 117 129 md = checkfield(md, 'fieldname', 'smb.datebreaks', 'NaN', 1, 'Inf', 1, 'size', [nbas,nbrk]) 118 elif (np.size(md.smb.datebreaks)==0 or np.all(np.isnan(md.smb.datebreaks))):130 elif np.size(md.smb.datebreaks) == 0 or np.all(np.isnan(md.smb.datebreaks)): 119 131 pass 120 132 else: 121 133 raise RuntimeError('md.smb.num_breaks is 0 but md.smb.datebreaks is not empty') 122 134 123 if (np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1):135 if np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1: 124 136 if len(np.shape(self.refelevation)) == 1: 125 137 self.refelevation = np.array([self.refelevation]) 126 md = checkfield(md, 'fieldname', 'smb.refelevation', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins)127 128 if (np.any(np.isnan(self.lapserates) is False) or np.size(self.lapserates) > 1):138 md = checkfield(md, 'fieldname', 'smb.refelevation', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [1, nbas], 'numel', nbas) 139 140 if (np.any(np.isnan(self.lapserates) is False) or np.size(self.lapserates) > 1): 129 141 nbas = md.smb.num_basins 130 142 if len(np.shape(self.lapserates)) == 1: … … 140 152 self.elevationbins = np.reshape(self.elevationbins,[nbas,max(1,nbins-1),ntmlapse]) 141 153 md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [nbas,nbins,ntmlapse], 'numel', md.smb.num_basins*nbins*ntmlapse) 142 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)143 for rr in range( md.smb.num_basins):154 md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [nbas,max(1,nbins-1),ntmlapse], 'numel', nbas*max(1,nbins-1)*ntmlapse) 155 for rr in range(nbas): 144 156 if(np.all(self.elevationbins[rr,0:-1]<=self.elevationbins[rr,1:])==False): 145 157 raise TypeError('md.smb.elevationbins should have rows in order of increasing elevation') 146 elif (np.any(np.isnan(self.elevationbins) is False) or np.size(self.elevationbins) > 1):147 # elevationbins specified but not lapserates: this will inevitably lead to inconsistencies158 elif (np.any(np.isnan(self.elevationbins) is False) or np.size(self.elevationbins) > 1): 159 # Elevationbins specified but not lapserates: this will inevitably lead to inconsistencies 148 160 nbas = md.smb.num_basins 149 161 if len(np.shape(self.elevationbins)) == 1: … … 155 167 elif(len(np.shape(self.lapserates)) == 3): 156 168 nbins = np.shape(self.lapserates)[1] 157 nbins = nbins -1169 nbins = nbins - 1 158 170 ntmlapse = np.shape(self.lapserates)[2] 159 md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [ md.smb.num_basins, nbins*ntmlapse], 'numel', md.smb.num_basins*nbins*ntmlapse)160 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)171 md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [nbas, nbins * ntmlapse], 'numel', nbas * nbins * ntmlapse) 172 md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [nbas, max(1, nbins - 1) * ntmlapse], 'numel', nbas * max(1, nbins - 1) * ntmlapse) 161 173 162 174 md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) … … 168 180 def marshall(self, prefix, md, fid): # {{{ 169 181 yts = md.constants.yts 170 nbas = md.smb.num_basins ;171 nprm = md.smb.num_params ;172 nper = md.smb.num_breaks +1;182 nbas = md.smb.num_basins 183 nprm = md.smb.num_params 184 nper = md.smb.num_breaks + 1 173 185 if(np.any(np.isnan(md.smb.lapserates))): 174 templapserates = np.zeros(( md.smb.num_basins,2,12))186 templapserates = np.zeros((nbas, 2, 12)) 175 187 print(' smb.lapserates not specified: set to 0') 176 tempelevationbins = np.zeros(( md.smb.num_basins,1,12)) #dummy elevation bins188 tempelevationbins = np.zeros((nbas, 1, 12)) # Dummy elevation bins 177 189 nbins = 2 178 190 ntmlapse = 12 179 191 else: 180 if (len(np.shape(md.smb.lapserates))==1):192 if len(np.shape(md.smb.lapserates)) == 1: 181 193 nbins = 1 182 194 ntmlapse = 1 183 elif (len(np.shape(md.smb.lapserates))==2):195 elif len(np.shape(md.smb.lapserates)) == 2: 184 196 nbins = np.shape(md.smb.lapserates)[1] 185 197 ntmlapse = 1 186 elif (len(np.shape(md.smb.lapserates))==3):198 elif len(np.shape(md.smb.lapserates)) == 3: 187 199 nbins = np.shape(md.smb.lapserates)[1] 188 200 ntmlapse = np.shape(md.smb.lapserates)[2] 189 templapserates = np.reshape(md.smb.lapserates,[nbas, nbins,ntmlapse])190 tempelevationbins = np.reshape(md.smb.elevationbins, [nbas,max(1,nbins-1),ntmlapse])201 templapserates = np.reshape(md.smb.lapserates,[nbas, nbins, ntmlapse]) 202 tempelevationbins = np.reshape(md.smb.elevationbins, [nbas, max(1, nbins - 1), ntmlapse]) 191 203 temprefelevation = np.copy(md.smb.refelevation) 192 # Scale the parameters #204 # Scale the parameters 193 205 polyparamsScaled = np.copy(md.smb.polynomialparams) 194 polyparams2dScaled = np.zeros((nbas, nper*nprm))195 if (nprm>1):196 # Case 3D #197 if (nbas>1 and nper>1):198 for ii in range(nprm): 199 polyparamsScaled[:, :,ii] = polyparamsScaled[:,:,ii]*(1/yts)**(ii+1)200 # Fit in 2D array #201 for ii in range(nprm): 202 polyparams2dScaled[:, ii*nper:(ii+1)*nper] = 1*polyparamsScaled[:,:,ii]203 # Case 2D and higher-order params at increasing row index #204 elif (nbas==1):205 for ii in range(nprm): 206 polyparamsScaled[ii, :] = polyparamsScaled[ii,:]*(1/yts)**(ii+1)207 # Fit in row array #208 for ii in range(nprm): 209 polyparams2dScaled[0, ii*nper:(ii+1)*nper] = 1*polyparamsScaled[ii,:]210 # Case 2D and higher-order params at incr asing column index #211 elif (nper==1):212 for ii in range(nprm): 213 polyparamsScaled[:, ii] = polyparamsScaled[:,ii]*(1/yts)**(ii+1)214 # 2D array is already in correct format #206 polyparams2dScaled = np.zeros((nbas, nper * nprm)) 207 if nprm > 1: 208 # Case 3D 209 if nbas > 1 and nper > 1: 210 for ii in range(nprm): 211 polyparamsScaled[:, :, ii] = polyparamsScaled[:, :, ii] * (1 / yts) ** (ii + 1) 212 # Fit in 2D array 213 for ii in range(nprm): 214 polyparams2dScaled[:, ii * nper:(ii + 1) * nper] = 1 * polyparamsScaled[:, :, ii] 215 # Case 2D and higher-order params at increasing row index 216 elif nbas == 1: 217 for ii in range(nprm): 218 polyparamsScaled[ii, :] = polyparamsScaled[ii, :] * (1 / yts) ** (ii + 1) 219 # Fit in row array 220 for ii in range(nprm): 221 polyparams2dScaled[0, ii * nper:(ii + 1) * nper] = 1 * polyparamsScaled[ii, :] 222 # Case 2D and higher-order params at increasing column index 223 elif nper == 1: 224 for ii in range(nprm): 225 polyparamsScaled[:, ii] = polyparamsScaled[:, ii] * (1 / yts) ** (ii + 1) 226 # 2D array is already in correct format 215 227 polyparams2dScaled = np.copy(polyparamsScaled) 216 228 else: 217 polyparamsScaled = polyparamsScaled *(1/yts)218 # 2D array is already in correct format #229 polyparamsScaled = polyparamsScaled * (1 / yts) 230 # 2D array is already in correct format 219 231 polyparams2dScaled = np.copy(polyparamsScaled) 220 221 if (nper==1):222 dbreaks = np.zeros((nbas, 1))232 233 if nper == 1: 234 dbreaks = np.zeros((nbas, 1)) 223 235 else: 224 236 dbreaks = np.copy(md.smb.datebreaks) 225 237 226 if (ntmlapse==1):227 templapserates = np.repeat(templapserates, 12,axis=2)228 tempelevationbins = np.repeat(tempelevationbins, 12,axis=2)229 if (np.any(np.isnan(md.smb.refelevation))):230 temprefelevation = np.zeros(( md.smb.num_basins)).reshape(1,md.smb.num_basins)238 if ntmlapse == 1: 239 templapserates = np.repeat(templapserates, 12, axis = 2) 240 tempelevationbins = np.repeat(tempelevationbins, 12, axis = 2) 241 if np.any(np.isnan(md.smb.refelevation)): 242 temprefelevation = np.zeros((nbas)).reshape(1, nbas) 231 243 areas = GetAreas(md.mesh.elements, md.mesh.x, md.mesh.y) 232 244 for ii, bid in enumerate(np.unique(md.smb.basin_id)): … … 239 251 print(' smb.refelevation not specified: Reference elevations set to mean surface elevation of basins') 240 252 nbins = np.shape(templapserates)[1] 241 temp2dlapserates = np.zeros((nbas, nbins*12))242 temp2delevationbins = np.zeros((nbas, max(12,(nbins-1)*12)))253 temp2dlapserates = np.zeros((nbas, nbins * 12)) 254 temp2delevationbins = np.zeros((nbas, max(12, (nbins - 1) * 12))) 243 255 for ii in range(12): 244 temp2dlapserates[:, ii*nbins:(ii+1)*nbins] = templapserates[:,:,ii]245 temp2delevationbins[:, ii*(nbins-1):(ii+1)*(nbins-1)] = tempelevationbins[:,:,ii]256 temp2dlapserates[:, ii * nbins:(ii + 1) * nbins] = templapserates[:, :, ii] 257 temp2delevationbins[:, ii * (nbins - 1):(ii + 1) * (nbins - 1)] = tempelevationbins[:, :, ii] 246 258 247 259 WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 13, 'format', 'Integer')
Note:
See TracChangeset
for help on using the changeset viewer.
