%SMBsemic Class definition
%
%   Usage:
%      SMBsemic=SMBsemic();

classdef SMBsemic
	properties (SetAccess=public)
		dailysnowfall		= NaN;
		dailyrainfall		= NaN;
		dailydsradiation	= NaN;
		dailydlradiation	= NaN;
		dailywindspeed		= NaN;
		dailypressure		= NaN;
		dailyairdensity		= NaN;
		dailyairhumidity	= NaN;
		dailytemperature	= NaN;
		desfac				= 0;
		rlaps				= 0;
		rdl					= 0;
		s0gcm				= NaN;
		steps_per_step		= 1;
		averaging			= 0;
		requested_outputs	= {};
	end
	methods
		function self = SMBsemic(varargin) % {{{
			switch nargin
				case 0
					self=setdefaultparameters(self);
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function self = extrude(self,md) % {{{
			self.dailysnowfall=project3d(md,'vector',self.dailysnowfall,'type','node');
			self.dailyrainfall=project3d(md,'vector',self.dailyrainfall,'type','node');
			self.dailydsradiation=project3d(md,'vector',self.dailydsradiation,'type','node');
			self.dailydlradiation=project3d(md,'vector',self.dailydlradiation,'type','node');
			self.dailywindspeed=project3d(md,'vector',self.dailywindspeed,'type','node');
			self.dailypressure=project3d(md,'vector',self.dailypressure,'type','node');
			self.dailyairdensity=project3d(md,'vector',self.dailyairdensity,'type','node');
			self.dailyairhumidity=project3d(md,'vector',self.dailyairhumidity,'type','node');
			self.dailytemperature=project3d(md,'vector',self.dailytemperature,'type','node');
			self.s0gcm=project3d(md,'vector',self.s0gcm,'type','node');

		end % }}}
		function list = defaultoutputs(self,md) % {{{
			list = {'SmbMassBalance'};
		end % }}}
		function self = initialize(self,md) % {{{

			if isnan(self.s0gcm),
				self.s0gcm=zeros(md.mesh.numberofvertices,1);
				disp('      no SMBsemic.s0gcm specified: values set as zero');
			end

		end % }}}
		function self = setdefaultparameters(self) % {{{

			self.desfac		= -log(2.0)/1000;
			self.rlaps		= 7.4;
			self.rdl			= 0.29;
			self.requested_outputs={'default'};

		end % }}}zo
		function md = checkconsistency(self,md,solution,analyses) % {{{

			if ismember('MasstransportAnalysis',analyses),
				md = checkfield(md,'fieldname','smb.desfac','<=',1,'numel',1);
				md = checkfield(md,'fieldname','smb.s0gcm','>=',0,'NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
				md = checkfield(md,'fieldname','smb.rlaps','>=',0,'numel',1);
				md = checkfield(md,'fieldname','smb.rdl','>=',0,'numel',1);
				md = checkfield(md,'fieldname','smb.dailysnowfall','timeseries',1,'NaN',1,'Inf',1);
				md = checkfield(md,'fieldname','smb.dailyrainfall','timeseries',1,'NaN',1,'Inf',1);
				md = checkfield(md,'fieldname','smb.dailydsradiation','timeseries',1,'NaN',1,'Inf',1);
				md = checkfield(md,'fieldname','smb.dailydlradiation','timeseries',1,'NaN',1,'Inf',1);
				md = checkfield(md,'fieldname','smb.dailywindspeed','timeseries',1,'NaN',1,'Inf',1);
				md = checkfield(md,'fieldname','smb.dailypressure','timeseries',1,'NaN',1,'Inf',1);
				md = checkfield(md,'fieldname','smb.dailyairdensity','timeseries',1,'NaN',1,'Inf',1);
				md = checkfield(md,'fieldname','smb.dailyairhumidity','timeseries',1,'NaN',1,'Inf',1);
				md = checkfield(md,'fieldname','smb.dailytemperature','timeseries',1,'NaN',1,'Inf',1);
			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:'));

			disp(sprintf('   Interface for coupling GCM data to the energy balance model SEMIC (Krapp et al (2017) https://doi.org/10.5194/tc-11-1519-2017).'));
			disp(sprintf('   The implemented coupling uses daily mean GCM input to calculate yearly mean smb, accumulation, ablation, and surface temperature.'));
			disp(sprintf('   smb and temperatures are updated every year'));
			disp(sprintf('\n   SEMIC parameters:'));
			fielddisplay(self,'dailysnowfall','daily surface dailysnowfall [m/s]');
			fielddisplay(self,'dailyrainfall','daily surface dailyrainfall [m/s]');
			fielddisplay(self,'dailydsradiation','daily downwelling shortwave radiation [W/m2]');
			fielddisplay(self,'dailydlradiation','daily downwelling longwave radiation [W/m2]');
			fielddisplay(self,'dailywindspeed','daily surface wind speed [m/s]');
			fielddisplay(self,'dailypressure','daily surface pressure [Pa]');
			fielddisplay(self,'dailyairdensity','daily air density [kg/m3]');
			fielddisplay(self,'dailyairhumidity','daily air specific humidity [kg/kg]');
			fielddisplay(self,'dailytemperature','daily surface air temperature [K]');
			fielddisplay(self,'rlaps','present day lapse rate (default is 7.4 [degree/km]; Erokhina et al. 2017)');
			fielddisplay(self,'desfac','desertification elevation factor (default is -log(2.0)/1000 [1/km]; Vizcaino et al. 2010)');
			fielddisplay(self,'rdl','longwave downward radiation decrease (default is 0.29 [W/m^2/km]; Marty et al. 2002)');
			fielddisplay(self,'s0gcm','GCM reference elevation; (default is 0) [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;

			WriteData(fid,prefix,'name','md.smb.model','data',12,'format','Integer');

			WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfac','format','Double');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0gcm','format','DoubleMat','mattype',1);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','rlaps','format','Double');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','rdl','format','Double');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailysnowfall','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyrainfall','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailydsradiation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailydlradiation','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailywindspeed','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailypressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyairdensity','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyairhumidity','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailytemperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
			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