Index: /issm/trunk-jpl/src/c/shared/Random/random.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Random/random.cpp (revision 27275) +++ /issm/trunk-jpl/src/c/shared/Random/random.cpp (revision 27276) @@ -44,4 +44,5 @@ for(int i=0;i',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.const','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); %scheme fails if passed as column vector - md = checkfield(md,'fieldname','smb.trend','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); %scheme fails if passed as column vector - md = checkfield(md,'fieldname','smb.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); - md = checkfield(md,'fieldname','smb.ar_initialtime','numel',1,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','smb.ar_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]); - - 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); - 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)); - 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) - %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)); - end - end - md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]); - md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]); - md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); - end % }}} - function disp(self) % {{{ - disp(sprintf(' surface forcings parameters:')); - fielddisplay(self,'num_basins','number of different basins [unitless]'); - fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); - fielddisplay(self,'const','basin-specific constant values [m ice eq./yr]'); - fielddisplay(self,'trend','basin-specific trend values [m ice eq. yr^(-2)]'); - fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); - fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]'); - fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); - fielddisplay(self,'arlag_coefs','basin-specific vectors of 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,'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'); - fielddisplay(self, 'averaging', 'averaging methods from short to long steps'); - disp(sprintf('%51s 0: Arithmetic (default)',' ')); - disp(sprintf('%51s 1: Geometric',' ')); - disp(sprintf('%51s 2: Harmonic',' ')); - fielddisplay(self,'requested_outputs','additional outputs requested'); - - end % }}} - function marshall(self,prefix,md,fid) % {{{ - - yts=md.constants.yts; - - templapserates = md.smb.lapserates; - 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); - disp(' smb.lapserates not specified: set to 0'); - tempelevationbins = zeros(md.smb.num_basins,1); %dummy elevation bins - end - if(any(isnan(md.smb.refelevation))) - temprefelevation = zeros(1,md.smb.num_basins); - areas = GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y); - for ii=1:md.smb.num_basins - indices = find(md.smb.basin_id==ii); - elemsh = zeros(numel(indices),1); - for jj=1:numel(indices) - elemsh(jj) = mean(md.geometry.surface(md.mesh.elements(indices(jj),:))); - end - temprefelevation(ii) = sum(areas(indices).*elemsh)/sum(areas(indices)); - end - if(any(reshape(md.smb.lapserates,[1,nbas*nbins])~=0)) - disp(' smb.refelevation not specified: Reference elevations set to mean surface elevation of basins'); - end - end - [nbas,nbins] = size(templapserates); - - WriteData(fid,prefix,'name','md.smb.model','data',13,'format','Integer'); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','num_basins','format','Integer'); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_order','format','Integer'); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_initialtime','format','Double','scale',yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','ar_timestep','format','Double','scale',yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','basin_id','data',self.basin_id-1,'name','md.smb.basin_id','format','IntMat','mattype',2); %0-indexed - WriteData(fid,prefix,'object',self,'class','smb','fieldname','const','format','DoubleMat','name','md.smb.const','scale',1./yts,'yts',yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','trend','format','DoubleMat','name','md.smb.trend','scale',1./(yts^2),'yts',yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','arlag_coefs','format','DoubleMat','name','md.smb.arlag_coefs','yts',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',temprefelevation,'format','DoubleMat','name','md.smb.refelevation'); - WriteData(fid,prefix,'data',nbins,'format','Integer','name','md.smb.num_bins'); - WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer'); - - %process requested outputs - outputs = self.requested_outputs; - pos = find(ismember(outputs,'default')); - if ~isempty(pos), - 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'); - - end % }}} - end -end Index: sm/trunk-jpl/src/m/classes/SMBautoregression.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBautoregression.py (revision 27275) +++ (revision ) @@ -1,188 +1,0 @@ -import numpy as np - -from checkfield import * -from fielddisplay import fielddisplay -from project3d import * -from WriteData import * -from GetAreas import * - -class SMBautoregression(object): - """SMBAUTOREGRESSION class definition - - Usage: - SMBautoregression = SMBautoregression() - """ - - def __init__(self, *args): # {{{ - self.num_basins = 0 - self.const = np.nan - self.trend = np.nan - self.ar_order = 0 - self.ar_initialtime = 0 - self.ar_timestep = 0 - self.arlag_coefs = np.nan - self.basin_id = np.nan - self.lapserates = np.nan - self.elevationbins = np.nan - self.refelevation = np.nan - self.steps_per_step = 1 - self.averaging = 0 - self.requested_outputs = [] - - nargs = len(args) - if nargs == 0: - self.setdefaultparameters() - else: - raise Exception('constructor not supported') - # }}} - - def __repr__(self): # {{{ - s = ' surface forcings parameters:\n' - s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'const', 'basin-specific constant values [m ice eq./yr]')) - s += '{}\n'.format(fielddisplay(self, 'trend', 'basin-specific trend values [m ice eq. yr^(-2)]')) - s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]')) - s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) - s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of 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, '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')) - s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) - s += '\t\t{}\n'.format('0: Arithmetic (default)') - s += '\t\t{}\n'.format('1: Geometric') - s += '\t\t{}\n'.format('2: Harmonic') - s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested')) - return s - # }}} - - def setdefaultparameters(self): #{{{ - self.ar_order = 0.0 # Autoregression model of order 0 - # }}} - - def extrude(self, md): # {{{ - return self # Nothing for now - # }}} - - def defaultoutputs(self, md): # {{{ - return [] - # }}} - - def initialize(self, md): # {{{ - if np.all(np.isnan(self.trend)): - self.trend = np.zeros((1, self.num_basins)) # No trend in SMB - print(' smb.trend (trend) not specified: value set to 0') - if self.ar_order == 0: - self.ar_order = 1 # Dummy 1 value for autoregression - self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autorgression coefficients all set to 0 - print(' smb.ar_order (order of autoregressive model) not specified: order of autoregressive model set to 0') - if self.ar_initialtime == 0: - self.ar_initialtime = md.timestepping.start_time # Autoregression model has no prescribed initial time - print(' smb.ar_initialtime (initial time in the autoregressive model parameterization) not specified: set to md.timestepping.start_time') - if self.ar_timestep == 0: - self.ar_timestep = md.timestepping.time_step # Autoregression model has no prescribed time step - print(' smb.ar_timestep (timestep of autoregressive model) not specified: set to md.timestepping.time_step') - if np.all(np.isnan(self.arlag_coefs)): - self.arlag_coefs = np.zeros((self.num_basins, self.ar_order)) # Autoregression model of order 0 - print(' smb.arlag_coefs (lag coefficients) not specified: order of autoregressive model set to 0') - return self - # }}} - - def checkconsistency(self, md, solution, analyses): # {{{ - if 'MasstransportAnalysis' in analyses: - md = checkfield(md, 'fieldname', 'smb.num_basins', '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.const)) == 1: - self.const = np.array([self.const]) - self.trend = np.array([self.trend]) - md = checkfield(md, 'fieldname', 'smb.const', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) # Scheme fails if passed as column vector - md = checkfield(md, 'fieldname', 'smb.trend', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) # Scheme fails if passed as column vector; NOTE: As opposed to MATLAB implementation, pass list - md = checkfield(md, 'fieldname', 'smb.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) - md = checkfield(md, 'fieldname', 'smb.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) - md = checkfield(md, 'fieldname', 'smb.ar_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]) - - 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): - if len(np.shape(self.lapserates)) == 1: - self.lapserates = np.array([self.lapserates]) - nbins = 1 - else: - 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)) - for rr in range(md.smb.num_basins): - 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 - 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)) - - md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) - md = checkfield(md, 'fieldname', 'smb.averaging', 'numel', [1], 'values', [0, 1, 2]) - md = checkfield(md, 'fieldname', 'smb.requested_outputs', 'stringrow', 1) - return md - # }}} - - def marshall(self, prefix, md, fid): # {{{ - yts = md.constants.yts - - templapserates = np.copy(md.smb.lapserates) - tempelevationbins = np.copy(md.smb.elevationbins) - temprefelevation = np.copy(md.smb.refelevation) - if(np.any(np.isnan(md.smb.lapserates))): - templapserates = np.zeros((md.smb.num_basins,2)) - print(' smb.lapserates not specified: set to 0') - tempelevationbins = np.zeros((md.smb.num_basins,1)) #dummy elevation bins - if(np.any(np.isnan(md.smb.refelevation))): - temprefelevation = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins) - areas = GetAreas(md.mesh.elements, md.mesh.x, md.mesh.y) - for ii, bid in enumerate(np.unique(md.smb.basin_id)): - indices = np.where(md.smb.basin_id == bid)[0] - elemsh = np.zeros((len(indices))) - for jj in range(len(indices)): - elemsh[jj] = np.mean(md.geometry.surface[md.mesh.elements[indices[jj], :] - 1]) - temprefelevation[0, ii] = np.sum(areas[indices] * elemsh) / np.sum(areas[indices]) - if(np.any(templapserates != 0)): - print(' smb.refelevation not specified: Reference elevations set to mean surface elevation of basins') - nbins = np.shape(templapserates)[1] - - WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 13, 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'num_basins', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ar_order', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ar_initialtime', 'format', 'Double', 'scale', yts) - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'ar_timestep', 'format', 'Double', 'scale', yts) - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.smb.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'const', 'format', 'DoubleMat', 'name', 'md.smb.const', 'scale', 1 / yts, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'trend', 'format', 'DoubleMat', 'name', 'md.smb.trend', 'scale', 1 / (yts ** 2), 'yts', yts) - WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.smb.arlag_coefs', 'yts', 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', temprefelevation, 'name', 'md.smb.refelevation', 'format', 'DoubleMat') - WriteData(fid, prefix, 'data', nbins, 'name', 'md.smb.num_bins', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer') - - # Process requested outputs - outputs = self.requested_outputs - indices = [i for i, x in enumerate(outputs) if x == 'default'] - if len(indices) > 0: - outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:] - outputs = outputscopy - WriteData(fid, prefix, 'data', outputs, 'name', 'md.smb.requested_outputs', 'format', 'StringArray') - - # }}} Index: sm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.m (revision 27275) +++ (revision ) @@ -1,117 +1,0 @@ -%AUTOREGRESSION LINEAR BASAL FORCINGS class definition -% -% Usage: -% autoregressionlinearbasalforcings=spatiallinearbasalforcings(); - -classdef autoregressionlinearbasalforcings - - %VV: not verified yet (27Jan2022) - - properties (SetAccess=public) - num_basins = 0; - const = NaN; - trend = NaN; - ar_order = 0; - ar_initialtime = 0; - ar_timestep = 0; - arlag_coefs = NaN; - basin_id = NaN; - groundedice_melting_rate = NaN; - deepwater_elevation = NaN; - upperwater_melting_rate = NaN; - upperwater_elevation = NaN; - geothermalflux = NaN; - end - methods - function self = autoregressionlinearbasalforcings(varargin) % {{{ - switch nargin - case 0 - self=setdefaultparameters(self); - otherwise - error('constructor not supported'); - end - end % }}} - function self = extrude(self,md) % {{{ - self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1); - self.deepwater_elevation=project3d(md,'vector',self.deepwater_elevation,'type','node','layer',1); - self.upperwater_melting_rate=project3d(md,'vector',self.upperwater_melting_rate,'type','node','layer',1); - self.upperwater_elevation=project3d(md,'vector',self.upperwater_elevation,'type','node','layer',1); - self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux - end % }}} - function self = initialize(self,md) % {{{ - - if isnan(self.groundedice_melting_rate), - self.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - disp(' no basalforcings.groundedice_melting_rate specified: values set as zero'); - end - - end % }}} - function self = setdefaultparameters(self) % {{{ - %Nothing for now - end % }}} - function md = checkconsistency(self,md,solution,analyses) % {{{ - - if ismember('MasstransportAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.ismasstransport==0), - md = checkfield(md,'fieldname','basalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); - md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); - md = checkfield(md,'fieldname','basalforcings.upperwater_melting_rate','NaN',1,'Inf',1,'>=',0,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); - md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','NaN',1,'Inf',1,'<=',0,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); - md = checkfield(md,'fieldname','basalforcings.basin_id','Inf',1,'>=',0,'<=',md.basalforcings.num_basins,'size',[md.mesh.numberofelements,1]); - md = checkfield(md,'fieldname','basalforcings.const','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); - md = checkfield(md,'fieldname','basalforcings.trend','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); - md = checkfield(md,'fieldname','basalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); - md = checkfield(md,'fieldname','basalforcings.ar_initialtime','numel',1,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','basalforcings.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep - md = checkfield(md,'fieldname','basalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins,md.basalforcings.ar_order]); - end - if ismember('BalancethicknessAnalysis',analyses), - error('not implemented yet!'); - end - if ismember('ThermalAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isthermal==0), - error('not implemented yet!'); - end - if numel(md.basalforcings.geothermalflux)>1 - md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'timeseries',1,'>=',0); - end - end % }}} - function disp(self) % {{{ - disp(sprintf(' autoregression linear basal forcings parameters:')); - disp(sprintf(' autoregressive model is applied for deepwater_melting_rate')); - - fielddisplay(self,'num_basins','number of different basins [unitless]'); - fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); - fielddisplay(self,'const','basin-specific constant values [m/yr]'); - fielddisplay(self,'trend','basin-specific trend values [m yr^(-2)]'); - fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); - fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]'); - fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); - fielddisplay(self,'arlag_coefs','basin-specific vectors of lag coefficients [unitless]'); - fielddisplay(self,'deepwater_elevation','basin-specific elevation of ocean deepwater [m]'); - fielddisplay(self,'upperwater_melting_rate','basin-specific basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]'); - fielddisplay(self,'upperwater_elevation','basin-specific elevation of ocean upperwater [m]'); - fielddisplay(self,'groundedice_melting_rate','node-specific basal melting rate (positive if melting) [m/yr]'); - fielddisplay(self,'geothermalflux','node-specific geothermal heat flux [W/m^2]'); - - end % }}} - function marshall(self,prefix,md,fid) % {{{ - - yts=md.constants.yts; - - WriteData(fid,prefix,'name','md.basalforcings.model','data',9,'format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','name','md.basalforcings.geothermalflux','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'fieldname','num_basins','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','ar_order','format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','ar_initialtime','format','Double','scale',yts); - WriteData(fid,prefix,'object',self,'fieldname','ar_timestep','format','Double','scale',yts); - WriteData(fid,prefix,'object',self,'fieldname','basin_id','data',self.basin_id-1,'name','md.basalforcings.basin_id','format','IntMat','mattype',2); %0-indexed - WriteData(fid,prefix,'object',self,'fieldname','const','format','DoubleMat','name','md.basalforcings.const','scale',1./yts,'yts',yts); - WriteData(fid,prefix,'object',self,'fieldname','trend','format','DoubleMat','name','md.basalforcings.trend','scale',1./(yts^2),'yts',yts); - WriteData(fid,prefix,'object',self,'fieldname','arlag_coefs','format','DoubleMat','name','md.basalforcings.arlag_coefs','yts',yts); - WriteData(fid,prefix,'object',self,'fieldname','deepwater_elevation','format','DoubleMat','name','md.basalforcings.deepwater_elevation'); - WriteData(fid,prefix,'object',self,'fieldname','upperwater_melting_rate','format','DoubleMat','name','md.basalforcings.upperwater_melting_rate','scale',1./yts); - WriteData(fid,prefix,'object',self,'fieldname','upperwater_elevation','format','DoubleMat','name','md.basalforcings.upperwater_elevation'); - end % }}} - end -end Index: sm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.py (revision 27275) +++ (revision ) @@ -1,122 +1,0 @@ -import numpy as np - -from checkfield import * -from fielddisplay import fielddisplay -from project3d import * -from WriteData import * - -class autoregressionlinearbasalforcings(object): - """autoregressionlinearbasalforcings class definition - - Usage: - autoregressionlinearbasalforcings = autoregressionlinearbasalforcings() - """ - - def __init__(self, *args): # {{{ - self.num_basins = 0 - self.const = np.nan - self.trend = np.nan - self.ar_order = 0 - self.ar_initialtime = 0 - self.ar_timestep = 0 - self.arlag_coefs = np.nan - self.basin_id = np.nan - self.groundedice_melting_rate = np.nan - self.deepwater_elevation = np.nan - self.upperwater_melting_rate = np.nan - self.upperwater_elevation = np.nan - self.geothermalflux = np.nan - - nargs = len(args) - if nargs == 0: - self.setdefaultparameters() - else: - raise Exception('constructor not supported') - # }}} - - def __repr__(self): # {{{ - s = ' surface forcings parameters:\n' - s += ' autoregressive model is applied for deepwater_melting_rate\n' - s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'const', 'basin-specific constant values [m ice eq./yr]')) - s += '{}\n'.format(fielddisplay(self, 'trend', 'basin-specific trend values [m ice eq. yr^(-2)]')) - s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]')) - s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) - s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of lag coefficients [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'deepwater_elevation', 'basin-specific elevation of ocean deepwater [m]')) - s += '{}\n'.format(fielddisplay(self, 'upperwater_melting_rate', 'basin-specic basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]')) - s += '{}\n'.format(fielddisplay(self, 'upperwater_elevation', 'basin-specific elevation of ocean upperwater [m]')) - s += '{}\n'.format(fielddisplay(self, 'groundedice_melting_rate','node-specific basal melting rate (positive if melting) [m/yr]')) - s += '{}\n'.format(fielddisplay(self, 'geothermalflux','node-specific geothermal heat flux [W/m^2]')) - return s - # }}} - - def setdefaultparameters(self): #{{{ - #Nothing for now - return self - # }}} - - def extrude(self, md): # {{{ - return self # Nothing for now - # }}} - - def initialize(self, md): # {{{ - if np.all(np.isnan(self.groundedice_melting_rate)): - self.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices)) - print(" no basalforcings.groundedice_melting_rate specified: values set as zero") - return self - # }}} - - def checkconsistency(self, md, solution, analyses): # {{{ - if 'MasstransportAnalysis' in analyses: - md = checkfield(md, 'fieldname', 'basalforcings.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) - md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1) - - if len(np.shape(self.deepwater_elevation)) == 1: - self.deepwater_elevation = np.array([self.deepwater_elevation]) - self.upperwater_elevation = np.array([self.upperwater_elevation]) - self.upperwater_melting_rate = np.array([self.upperwater_melting_rate]) - md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) - md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', 'NaN', 1, 'Inf', 1, '<=', 0, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) - md = checkfield(md, 'fieldname', 'basalforcings.upperwater_melting_rate', 'NaN', 1, 'Inf', 1,'>=', 0, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) - md = checkfield(md, 'fieldname', 'basalforcings.basin_id', 'Inf', 1, '>=', 0, '<=', md.basalforcings.num_basins, 'size', [md.mesh.numberofelements]) - - if len(np.shape(self.const)) == 1: - self.const = np.array([self.const]) - self.trend = np.array([self.trend]) - - md = checkfield(md, 'fieldname', 'basalforcings.const', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) # Scheme fails if passed as column vector - md = checkfield(md, 'fieldname', 'basalforcings.trend', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) # Scheme fails if passed as column vector; NOTE: As opposed to MATLAB implementation, pass list - md = checkfield(md, 'fieldname', 'basalforcings.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) - md = checkfield(md, 'fieldname', 'basalforcings.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) - md = checkfield(md, 'fieldname', 'basalforcings.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep - md = checkfield(md, 'fieldname', 'basalforcings.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.basalforcings.num_basins, md.basalforcings.ar_order]) - if 'BalancethicknessAnalysis' in analyses: - raise Exception('not implemented yet!') - if 'ThermalAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isthermal: - raise Exception('not implemented yet!') - - return md - # }}} - - def marshall(self, prefix, md, fid): # {{{ - yts = md.constants.yts - - WriteData(fid, prefix, 'name', 'md.basalforcings.model', 'data', 9, 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'groundedice_melting_rate', 'name', 'md.basalforcings.groundedice_melting_rate', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'geothermalflux', 'name', 'md.basalforcings.geothermalflux', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'num_basins', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_order', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_initialtime', 'format', 'Double', 'scale', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_timestep', 'format', 'Double', 'scale', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.basalforcings.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed - WriteData(fid, prefix, 'object', self, 'fieldname', 'const', 'format', 'DoubleMat', 'name', 'md.basalforcings.const', 'scale', 1 / yts, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'trend', 'format', 'DoubleMat', 'name', 'md.basalforcings.trend', 'scale', 1 / (yts ** 2), 'yts', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.basalforcings.arlag_coefs', 'yts', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'deepwater_elevation', 'format', 'DoubleMat', 'name', 'md.basalforcings.deepwater_elevation') - WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.upperwater_melting_rate', 'scale', 1 / yts, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_elevation', 'format', 'DoubleMat', 'name', 'md.basalforcings.upperwater_elevation') - - # }}} Index: sm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.m (revision 27275) +++ (revision ) @@ -1,89 +1,0 @@ -%FRONTAL FORCINGS rignot autoregression class definition -% -% Usage: -% frontalforcingsrignotautoregression=frontalforcingsrignotautoregression(); - -classdef frontalforcingsrignotautoregression - properties (SetAccess=public) - num_basins = 0; - const = NaN; - trend = NaN; - ar_order = 0; - ar_initialtime = 0; - ar_timestep = 0; - arlag_coefs = NaN; - basin_id = NaN; - subglacial_discharge = NaN; - end - methods - function self = frontalforcingsrignot(varargin) % {{{ - switch nargin - case 0 - self=setdefaultparameters(self); - case 1 - inputstruct=varargin{1}; - list1 = properties('frontalforcingsrignotautoregression'); - list2 = fieldnames(inputstruct); - for i=1:length(list1) - fieldname = list1{i}; - if ismember(fieldname,list2), - self.(fieldname) = inputstruct.(fieldname); - end - end - otherwise - error('constructor not supported'); - end - end % }}} - function self = extrude(self,md) % {{{ - % nothing for now - end % }}} - function self = setdefaultparameters(self) % {{{ - - self.basin_id = NaN; - self.num_basins = 0; - self.subglacial_discharge = NaN; - self.ar_order = 0.0; %autoregression model of order 0 - - end % }}} - function md = checkconsistency(self,md,solution,analyses) % {{{ - %Early return - if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end - - md = checkfield(md,'fieldname','frontalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); - md = checkfield(md,'fieldname','frontalforcings.basin_id','Inf',1,'>=',0,'<=',md.frontalforcings.num_basins,'size',[md.mesh.numberofelements 1]); - md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'NaN',1,'Inf',1,'timeseries',1); - md = checkfield(md,'fieldname','frontalforcings.const','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins); - md = checkfield(md,'fieldname','frontalforcings.trend','NaN',1,'Inf',1,'size',[1,md.frontalforcings.num_basins],'numel',md.frontalforcings.num_basins); - md = checkfield(md,'fieldname','frontalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); - md = checkfield(md,'fieldname','frontalforcings.ar_initialtime','numel',1,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','frontalforcings.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep - md = checkfield(md,'fieldname','frontalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]); - - end % }}} - function disp(self) % {{{ - disp(sprintf(' Frontalforcings parameters:')); - fielddisplay(self,'num_basins','number of different basins [unitless]'); - fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); - fielddisplay(self,'subglacial_discharge','sum of subglacial discharge for each basin [m/d]'); - fielddisplay(self,'const','basin-specific constant term [∘C]'); - fielddisplay(self,'trend','basin-specific trend values [∘C yr^(-1)]'); - fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); - fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]'); - fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); - fielddisplay(self,'arlag_coefs','basin-specific vectors of lag coefficients [unitless]'); - end % }}} - function marshall(self,prefix,md,fid) % {{{ - yts=md.constants.yts; - WriteData(fid,prefix,'name','md.frontalforcings.parameterization','data',3,'format','Integer'); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_basins','format','Integer'); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','subglacial_discharge','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_order','format','Integer'); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_initialtime','format','Double','scale',yts); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','ar_timestep','format','Double','scale',yts); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','basin_id','data',self.basin_id-1,'name','md.frontalforcings.basin_id','format','IntMat','mattype',2); %0-indexed - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','const','format','DoubleMat','name','md.frontalforcings.const'); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','trend','format','DoubleMat','name','md.frontalforcings.trend','scale',1./yts,'yts',yts); - WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','arlag_coefs','format','DoubleMat','name','md.frontalforcings.arlag_coefs','yts',yts); - end % }}} - end -end Index: sm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotautoregression.py (revision 27275) +++ (revision ) @@ -1,91 +1,0 @@ -# -*- coding: utf-8 -*- -import numpy as np -from checkfield import checkfield -from fielddisplay import fielddisplay -from MatlabFuncs import * -from WriteData import WriteData - - -class frontalforcingsrignotautoregression(object): - """FRONTALFORCINGSRIGNOTAUTOREGRESSION class definition - - Usage: - frontalforcingsrignotautoregression = frontalforcingsrignotautoregression() - """ - - def __init__(self, *args): # {{{ - self.num_basins = 0 - self.const = np.nan - self.trend = np.nan - self.ar_order = 0 - self.ar_initialtime = 0 - self.ar_timestep = 0 - self.arlag_coefs = np.nan - self.basin_id = np.nan - self.subglacial_discharge = np.nan - - if len(args) == 0: - self.setdefaultparameters() - else: - error('constructor not supported') - - def __repr__(self): # {{{ - s = ' Frontalforcings parameters:\n' - s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'subglacial_discharge', 'sum of subglacial discharge for each basin [m/d]')) - s += '{}\n'.format(fielddisplay(self, 'const', 'basin-specific constant values [°C]')) - s += '{}\n'.format(fielddisplay(self, 'trend', 'basin-specific trend values [°C yr^(-1)]')) - s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]')) - s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]')) - s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) - s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of lag coefficients [unitless]')) - return s - #}}} - - def setdefaultparameters(self): # {{{ - self.basin_id = np.nan - self.num_basins = 0 - self.subglacial_discharge = np.nan - self.ar_order = 0.0 # Autoregression model of order 0 - return self - #}}} - - def checkconsistency(self, md, solution, analyses): # {{{ - # Early return - if not (solution == 'TransientSolution') or not md.transient.ismovingfront: - return md - - md = checkfield(md, 'fieldname', 'frontalforcings.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) - md = checkfield(md, 'fieldname', 'frontalforcings.basin_id', 'Inf', 1, '>=', 0, '<=', md.frontalforcings.num_basins, 'size', [md.mesh.numberofelements]) - md = checkfield(md, 'fieldname', 'frontalforcings.subglacial_discharge', '>=', 0, 'NaN', 1, 'Inf', 1, 'timeseries', 1) - if len(np.shape(self.const)) == 1: - self.const = np.array([self.const]) - self.trend = np.array([self.trend]) - md = checkfield(md, 'fieldname', 'frontalforcings.const', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins) - md = checkfield(md, 'fieldname', 'frontalforcings.trend', 'NaN', 1, 'Inf', 1, 'size', [1, md.frontalforcings.num_basins], 'numel', md.frontalforcings.num_basins) - md = checkfield(md, 'fieldname', 'frontalforcings.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) - md = checkfield(md, 'fieldname', 'frontalforcings.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) - md = checkfield(md, 'fieldname', 'frontalforcings.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep - md = checkfield(md, 'fieldname', 'frontalforcings.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ar_order]) - return md - # }}} - - def extrude(self, md): # {{{ - # Nothing for now - return self - # }}} - - def marshall(self, prefix, md, fid): # {{{ - yts = md.constants.yts - WriteData(fid, prefix, 'name', 'md.frontalforcings.parameterization', 'data', 3, 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'num_basins', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'subglacial_discharge', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ar_order', 'format', 'Integer') - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ar_initialtime', 'format', 'Double', 'scale', yts) - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'ar_timestep', 'format', 'Double', 'scale', yts) - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'basin_id', 'data', self.basin_id - 1, 'name', 'md.frontalforcings.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'const', 'format', 'DoubleMat', 'name', 'md.frontalforcings.const') - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'trend', 'format', 'DoubleMat', 'name', 'md.frontalforcings.trend', 'scale', 1 / yts, 'yts', yts) - WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.arlag_coefs', 'yts', yts) - # }}} Index: /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m (revision 27275) +++ /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m (revision 27276) @@ -45,4 +45,5 @@ ['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'],... ['./' basename '-' ncdate '.nc'],... + '/media/vincent/TOSH4TB/GeorgiaTech/DataSearch/BedMachine/BedMachineGreenland-2021-04-20.nc',... }; @@ -89,7 +90,9 @@ if strcmp(string,'mask') | strcmp(string,'source'), %Need nearest neighbor to avoid interpolation between 0 and 2 - output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest'); + %output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest'); + output = InterpFromGridToMesh(xdata,flipud(ydata),flipud(data),double(X),double(Y),NaN); %VV else - output = InterpFromGrid(xdata,ydata,data,double(X),double(Y)); + %output = InterpFromGrid(xdata,ydata,data,double(X),double(Y)); + output = InterpFromGridToMesh(xdata,flipud(ydata),flipud(data),double(X),double(Y),NaN); %VV end Index: /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpSeaRISE.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpSeaRISE.m (revision 27275) +++ /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpSeaRISE.m (revision 27276) @@ -58,4 +58,8 @@ searisenc='/totten_1/ModelData/SeaRISE/Antarctica_5km_dev1.0.nc'; end + case {'vincentxps159500'} + if hemisphere==1 + searisenc = '/media/vincent/TOSH4TB/GeorgiaTech/ISSM/fromSeroussi/GrIS04Feb2022/seariseData/Greenland_5km_dev1.2.nc'; + end otherwise error('hostname not supported yet'); @@ -80,7 +84,9 @@ if verbose, disp([' -- SeaRISE: interpolating ' string]); end if strcmpi(string,'LandMask'); - dataout = InterpFromGrid(xdata,ydata,data,xproj,yproj,'nearest'); + %dataout = InterpFromGrid(xdata,ydata,data,xproj,yproj,'nearest'); + dataout = InterpFromGridToMesh(xdata,ydata,data,xproj,yproj,NaN); %VV else - dataout = InterpFromGrid(xdata,ydata,data,xproj,yproj); + %dataout = InterpFromGrid(xdata,ydata,data,xproj,yproj); + dataout = InterpFromGridToMesh(xdata,ydata,data,xproj,yproj,NaN); %VV end dataout = reshape(dataout,size(X,1),size(X,2));