Index: /issm/trunk-jpl/src/m/classes/spatiallinearbasalforcings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/spatiallinearbasalforcings.m (revision 22945) +++ /issm/trunk-jpl/src/m/classes/spatiallinearbasalforcings.m (revision 22945) @@ -0,0 +1,110 @@ +%SPATIAL LINEAR BASAL FORCINGS class definition +% +% Usage: +% spatiallinearbasalforcings=spatiallinearbasalforcings(); + +classdef spatiallinearbasalforcings + properties (SetAccess=public) + groundedice_melting_rate = NaN; + deepwater_melting_rate = NaN; + deepwater_elevation = NaN; + upperwater_elevation = NaN; + geothermalflux = NaN; + end + methods + function self = spatiallinearbasalforcings(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + case 1 + lb=varargin{1}; + if strcmpi(class(lb),'linearbasalforcings'), + nvertices=length(lb.groundedice_melting_rate); + self.groundedice_melting_rate=lb.groundedice_melting_rate; + self.geothermalflux=lb.geothermalflux; + self.deepwater_elevation=lb.deepwater_elevation*ones(nvertices,1); + self.deepwater_melting_rate=lb.deepwater_melting_rate*ones(nvertices,1); + self.upperwater_elevation=lb.upperwater_elevation*ones(nvertices,1); + else + self=structtoobj(spatiallinearbasalforcings(),varargin{1}); + end + 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_melting_rate=project3d(md,'vector',self.deepwater_melting_rate,'type','node','layer',1); + self.deepwater_elevation=project3d(md,'vector',self.deepwater_elevation,'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) % {{{ + + %default values for melting parameterization: can't, we don't have the vertices information! + %self.deepwater_melting_rate = 50; + %self.deepwater_elevation = -800; + %self.upperwater_elevation = -400; + + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ismember('MasstransportAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.ismasstransport==0), + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','basalforcings.deepwater_melting_rate','NaN',1,'Inf',1,'timeseries',1,'>=',0); + md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','NaN',1,'Inf',1,'timeseries',1,'<',0); + end + if ismember('BalancethicknessAnalysis',analyses), + error('not implemented yet!'); + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','basalforcings.deepwater_melting_rate','>=',0,'numel',1); + md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','<','basalforcings.upperwater_elevation','numel',1); + md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','<=',0,'numel',1); + end + if ismember('ThermalAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isthermal==0), + error('not implemented yet!'); + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','basalforcings.deepwater_melting_rate','>=',0,'numel',1); + md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','<','basalforcings.upperwater_elevation','numel',1); + md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','<',0,'numel',1); + md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'timeseries',1,'>=',0); + end + end % }}} + function disp(self) % {{{ + disp(sprintf(' basal forcings parameters:')); + + fielddisplay(self,'groundedice_melting_rate','basal melting rate (positive if melting) [m/yr]'); + fielddisplay(self,'deepwater_melting_rate','basal melting rate (positive if melting applied for floating ice whith base < deepwater_elevation) [m/yr]'); + fielddisplay(self,'deepwater_elevation','elevation of ocean deepwater [m]'); + fielddisplay(self,'upperwater_elevation','elevation of ocean upperwater [m]'); + fielddisplay(self,'geothermalflux','geothermal heat flux [W/m^2]'); + + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + yts=md.constants.yts; + + floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); + pos=find(md.geometry.base<=md.basalforcings.deepwater_elevation); + floatingice_melting_rate(pos)=md.basalforcings.deepwater_melting_rate(pos); + pos=find(md.geometry.base>md.basalforcings.deepwater_elevation & md.geometry.base