%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
