%SLR class definition % % Usage: % slr=slr(); classdef slr properties (SetAccess=public) deltathickness = NaN; sealevel = NaN; spcthickness = NaN; maxiter = 0; reltol = 0; abstol = 0; love_h = 0; %provided by PREM model love_k = 0; %ideam love_l = 0; %ideam tide_love_k = 0; %ideam tide_love_h = 0; %ideam fluid_love = 0; equatorial_moi = 0; polar_moi = 0; angular_velocity = 0; rigid = 0; elastic = 0; rotation = 0; ocean_area_scaling = 0; hydro_rate = 0; %rate of steric expansion from hydrological effects. geodetic_run_frequency = 1; %how many time steps we skip before we run the geodetic part of the solver during transient geodetic = 0; %compute geodetic SLR? (in addition to steric?) degacc = 0; loop_increment = 0; horiz = 0; Ngia = NaN; Ugia = NaN; requested_outputs = {}; transitions = {}; end methods function self = slr(varargin) % {{{ switch nargin case 0 self=setdefaultparameters(self); otherwise error('constructor not supported'); end end % }}} function self = setdefaultparameters(self) % {{{ %Convergence criterion: absolute, relative and residual self.reltol=0.01; % 1 per cent self.abstol=NaN; % default %maximum of non-linear iterations. self.maxiter=5; self.loop_increment=200; %computational flags: self.geodetic=0; self.rigid=1; self.elastic=1; self.ocean_area_scaling=0; self.rotation=1; %tidal love numbers: self.tide_love_h=0.6149; %degree 2 self.tide_love_k=0.3055; % degree 2 %secular fluid love number: self.fluid_love=0.942; %moment of inertia: self.equatorial_moi=8.0077*10^37; % [kg m^2] self.polar_moi =8.0345*10^37; % [kg m^2] % mean rotational velocity of earth self.angular_velocity=7.2921*10^-5; % [s^-1] %numerical discretization accuracy self.degacc=.01; %hydro self.hydro_rate=0; %how many time steps we skip before we run SLR solver during transient self.geodetic_run_frequency=1; %output default: self.requested_outputs={'default'}; %transitions should be a cell array of vectors: self.transitions={}; %horizontal displacement? (not by default) self.horiz=0; end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ if ~ismember('SealevelriseAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslr==0), return; end md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1); md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.love_l','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.fluid_love','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.equatorial_moi','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.polar_moi','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.reltol','size',[1 1]); md = checkfield(md,'fieldname','slr.abstol','size',[1 1]); md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1); md = checkfield(md,'fieldname','slr.geodetic_run_frequency','size',[1 1],'>=',1); md = checkfield(md,'fieldname','slr.hydro_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10); md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1); md = checkfield(md,'fieldname','slr.loop_increment','NaN',1,'Inf',1,'>=',1); md = checkfield(md,'fieldname','slr.horiz','NaN',1,'Inf',1,'values',[0 1]); md = checkfield(md,'fieldname','slr.Ngia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','slr.Ugia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); %check that love numbers are provided at the same level of accuracy: if (size(self.love_h,1)~=size(self.love_k,1) | size(self.love_h,1)~=size(self.love_l,1)), error('slr error message: love numbers should be provided at the same level of accuracy'); end %cross check that whereever we have an ice load, the mask is <0 on each vertex: pos=find(self.deltathickness); maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:)); [els,vertices]=find(maskpos>0); if length(els), warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!'); end %check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not, %a coupler to a planet model is provided. if self.geodetic, if md.transient.iscoupler, %we are good; else if strcmpi(class(md.mesh),'mesh3dsurface'), %we are good else error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!'); end end end end % }}} function list=defaultoutputs(self,md) % {{{ list = {'Sealevel'}; end % }}} function disp(self) % {{{ disp(sprintf(' slr parameters:')); fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]'); fielddisplay(self,'sealevel','current sea level (prior to computation) [m]'); fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]'); fielddisplay(self,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)'); fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied'); fielddisplay(self,'maxiter','maximum number of nonlinear iterations'); fielddisplay(self,'love_h','load Love number for radial displacement'); fielddisplay(self,'love_k','load Love number for gravitational potential perturbation'); fielddisplay(self,'love_l','load Love number for horizontal displacements'); fielddisplay(self,'tide_love_k','tidal load Love number (deg 2)'); fielddisplay(self,'tide_love_h','tidal load Love number (deg 2)'); fielddisplay(self,'fluid_love','secular fluid Love number'); fielddisplay(self,'equatorial_moi','mean equatorial moment of inertia [kg m^2]'); fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]'); fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]'); fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); fielddisplay(self,'hydro_rate','rate of hydrological expansion (in mm/yr)'); fielddisplay(self,'Ngia','rate of viscous (GIA) geoid expansion (in mm/yr)'); fielddisplay(self,'Ugia','rate of viscous (GIA) bedrock uplift (in mm/yr)'); fielddisplay(self,'loop_increment','vector assembly (in the convolution) framentation'); fielddisplay(self,'geodetic','compute geodetic SLR? ( in addition to steric?) default 0'); fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)'); fielddisplay(self,'rigid','rigid earth graviational potential perturbation'); fielddisplay(self,'elastic','elastic earth graviational potential perturbation'); fielddisplay(self,'rotation','earth rotational potential perturbation'); fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'); fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); fielddisplay(self,'requested_outputs','additional outputs requested'); end % }}} function marshall(self,prefix,md,fid) % {{{ %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2); WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1); WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','abstol','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','maxiter','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','love_h','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','love_k','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','love_l','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','fluid_love','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','equatorial_moi','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','polar_moi','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','hydro_rate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3/md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','Ngia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','Ugia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray'); WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','geodetic','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.slr.requested_outputs','format','StringArray'); end % }}} function savemodeljs(self,fid,modelname) % {{{ writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness); writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel); writejs1Darray(fid,[modelname '.slr.spcthickness'],self.spcthickness); writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter); writejsdouble(fid,[modelname '.slr.reltol'],self.reltol); writejsdouble(fid,[modelname '.slr.abstol'],self.abstol); writejs1Darray(fid,[modelname '.slr.love_h'],self.love_h); writejs1Darray(fid,[modelname '.slr.love_k'],self.love_k); writejs1Darray(fid,[modelname '.slr.love_l'],self.love_l); writejsdouble(fid,[modelname '.slr.tide_love_k'],self.tide_love_k); writejsdouble(fid,[modelname '.slr.tide_love_h'],self.tide_love_h); writejsdouble(fid,[modelname '.slr.fluid_love'],self.fluid_love); writejsdouble(fid,[modelname '.slr.equatorial_moi'],self.equatorial_moi); writejsdouble(fid,[modelname '.slr.polar_moi'],self.polar_moi); writejsdouble(fid,[modelname '.slr.angular_velocity'],self.angular_velocity); writejsdouble(fid,[modelname '.slr.rigid'],self.rigid); writejsdouble(fid,[modelname '.slr.elastic'],self.elastic); writejsdouble(fid,[modelname '.slr.rotation'],self.rotation); writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling); writejsdouble(fid,[modelname '.slr.geodetic_run_frequency'],self.geodetic_run_frequency); writejs1Darray(fid,[modelname '.slr.hydro_rate'],self.hydro_rate); writejsdouble(fid,[modelname '.slr.degacc'],self.degacc); writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs); writejscellarray(fid,[modelname '.slr.transitions'],self.transitions); end % }}} function self = extrude(self,md) % {{{ self.sealevel=project3d(md,'vector',self.sealevel,'type','node'); end % }}} end end