%SMBgemb Class definition. 
%   This is the class that hosts all the inputs for the Alberta Glacier Surface Mass Balance Model 
%   Alex Gardner, University of Alberta.
%   
%   Usage:
%      SMBgemb=SMBgemb();

classdef SMBgemb
	properties (SetAccess=public)  
	% {{{
		%each one of these properties is a transient forcing to the GEMB model, loaded from meteorological data derived 
		%from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number 
		%of time steps. )

		%solution choices
		isgraingrowth;
		isalbedo;
		isshortwave;
		isthermal;
		isaccumulation;
		ismelt;
		isdensification;
		isturbulentflux;    

		%inputs: 
		Ta    = NaN; %2 m air temperature, in Kelvin
		V     = NaN; %wind speed (m/s-1)
		dswrf = NaN; %downward shortwave radiation flux [W/m^2]
		dlwrf = NaN; %downward longwave radiation flux [W/m^2]
		P     = NaN; %precipitation [mm w.e. / m^2]
		eAir  = NaN; %screen level vapor pressure [Pa]
		pAir  = NaN; %surface pressure [Pa]
		
		Tmean = NaN; %mean annual temperature [K]
		C     = NaN; %mean annual snow accumulation [kg m-2 yr-1]
		Tz    = NaN; %height above ground at which temperature (T) was sampled [m]
		Vz    = NaN; %height above ground at which wind (V) eas sampled [m]
        
		% Initialization of snow properties
		Dzini = NaN; %cell depth (m)
		Dini = NaN; %snow density (kg m-3)
		Reini = NaN; %effective grain size (mm)
		Gdnini = NaN; %grain dendricity (0-1)
		Gspini = NaN; %grain sphericity (0-1)
		ECini = NaN; %evaporation/condensation (kg m-2)
		Wini = NaN; %Water content (kg m-2)
		Aini = NaN; %albedo (0-1)
		Tini = NaN; %snow temperature (K)
		Sizeini = NaN; %Number of layers

		%settings: 
		aIdx   = NaN; %method for calculating albedo and subsurface absorption (default is 1)
		              % 1: effective grain radius [Gardner & Sharp, 2009]
					  % 2: effective grain radius [Brun et al., 2009]
					  % 3: density and cloud amount [Greuell & Konzelmann, 1994]
					  % 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
		swIdx  = NaN; % apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1)

		denIdx = NaN; %densification model to use (default is 2):
					% 1 = emperical model of Herron and Langway (1980)
					% 2 = semi-emerical model of Anthern et al. (2010)
					% 3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010)
					% 4 = DO NOT USE: emperical model of Li and Zwally (2004)
					% 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)

		zTop  = NaN; % depth over which grid length is constant at the top of the snopack (default 10) [m]
		dzTop = NaN; % initial top vertical grid spacing (default .05) [m] 
		dzMin = NaN; % initial min vertical allowable grid spacing (default dzMin/2) [m] 
		zY    = NaN; % strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]
		zMax = NaN; %initial max model depth (default is min(thickness,500)) [m]
		zMin = NaN; %initial min model depth (default is min(thickness,30)) [m]
		outputFreq = NaN; %output frequency in days (default is monthly, 30)

		%specific albedo parameters: 
		%Method 1 and 2: 
		aSnow = NaN; % new snow albedo (0.64 - 0.89)
		aIce  = NaN; % range 0.27-0.58 for old snow
		%Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo
		cldFrac = NaN; % average cloud amount
		%Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005)
		t0wet = NaN; % time scale for wet snow (15-21.9) 
		t0dry = NaN; % warm snow timescale (30) 
		K     = NaN; % time scale temperature coef. (7) 

		%densities:
		InitDensityScaling= NaN; %initial scaling factor multiplying the density of ice, which describes the density of the snowpack.
		
		requested_outputs      = {};

		%Several fields are missing from the standard GEMB model, which are capture intrinsically by ISSM. 
		%dateN: that's the last row of the above fields. 
		%dt:    included in dateN. Not an input.  
		%elev:  this is taken from the ISSM surface itself.

	end % }}}
	methods
		function self = SMBgemb(varargin) % {{{
			switch nargin
				case 2
					mesh=varargin{1}; 
					geometry=varargin{2}; 
					self=setdefaultparameters(self,mesh,geometry);
				otherwise
					error('constructor not supported: need geometry and mesh to set defaults');
			end
		end % }}}
		function self = extrude(self,md) % {{{

			self.Ta=project3d(md,'vector',self.Ta,'type','node');
			self.V=project3d(md,'vector',self.V,'type','node');
			self.dswrf=project3d(md,'vector',self.dswrf,'type','node');
			self.dswrf=project3d(md,'vector',self.dswrf,'type','node');
			self.P=project3d(md,'vector',self.P,'type','node');
			self.eAir=project3d(md,'vector',self.eAir,'type','node');
			self.pAir=project3d(md,'vector',self.pAir,'type','node');

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

		self.isgraingrowth=1;
		self.isalbedo=1;
		self.isshortwave=1;
		self.isthermal=1;
		self.isaccumulation=1;
		self.ismelt=1;
		self.isdensification=1;
		self.isturbulentflux=1;
	
		self.aIdx = 1;
		self.swIdx = 1;
		self.denIdx = 2;
		self.zTop=10*ones(mesh.numberofelements,1);
		self.dzTop = .05* ones (mesh.numberofelements,1);
		self.dzMin = self.dzTop/2;
		self.InitDensityScaling = 1.0;
		
		self.zMax=250*ones(mesh.numberofelements,1);
		self.zMin=130*ones(mesh.numberofelements,1);
		self.zY = 1.10*ones(mesh.numberofelements,1);
		self.outputFreq = 30;
		
		%additional albedo parameters
		self.aSnow = 0.85;
		self.aIce = 0.48;
		self.cldFrac = 0.1; 
		self.t0wet = 15;
		self.t0dry = 30;
		self.K = 7;
        
		self.Dzini=0.05*ones(mesh.numberofelements,2);
		self.Dini=910.0*ones(mesh.numberofelements,2); 
		self.Reini=2.5*ones(mesh.numberofelements,2);
		self.Gdnini=0.0*ones(mesh.numberofelements,2);
		self.Gspini=0.0*ones(mesh.numberofelements,2);
		self.ECini=0.0*ones(mesh.numberofelements,1);
		self.Wini=0.0*ones(mesh.numberofelements,2);
		self.Aini=self.aSnow*ones(mesh.numberofelements,2); 
		self.Tini=273.15*ones(mesh.numberofelements,2); 
% 		/!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh). 
% 		If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp 
		self.Sizeini=2*ones(mesh.numberofelements,1);

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


			md = checkfield(md,'fieldname','smb.isgraingrowth','values',[0 1]);
			md = checkfield(md,'fieldname','smb.isalbedo','values',[0 1]);
			md = checkfield(md,'fieldname','smb.isshortwave','values',[0 1]);
			md = checkfield(md,'fieldname','smb.isthermal','values',[0 1]);
			md = checkfield(md,'fieldname','smb.isaccumulation','values',[0 1]);
			md = checkfield(md,'fieldname','smb.ismelt','values',[0 1]);
			md = checkfield(md,'fieldname','smb.isdensification','values',[0 1]);
			md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]);

			md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value
			md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<',45); %max 500 km/h
			md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400);
			md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'>=',0);
			md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',100);
			md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1);

			md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value
			md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0); 
			md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); 
			md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); 

			md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'Inf',1,'values',[1,2,3,4]);
			md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'Inf',1,'values',[0,1]);
			md = checkfield(md,'fieldname','smb.denIdx','NaN',1,'Inf',1,'values',[1,2,3,4,5]);

			md = checkfield(md,'fieldname','smb.zTop','NaN',1,'Inf',1,'>=',0);
			md = checkfield(md,'fieldname','smb.dzTop','NaN',1,'Inf',1,'>',0);
			md = checkfield(md,'fieldname','smb.dzMin','NaN',1,'Inf',1,'>',0);
			md = checkfield(md,'fieldname','smb.zY','NaN',1,'Inf',1,'>=',1);
			md = checkfield(md,'fieldname','smb.outputFreq','NaN',1,'Inf',1,'>',0,'<',10*365); %10 years max 
			md = checkfield(md,'fieldname','smb.InitDensityScaling','NaN',1,'Inf',1,'>=',0,'<=',1);

			switch self.aIdx,
				case {1 2}
					md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'>=',.64,'<=',.89);
					md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'>=',.27,'<=',.58);
				case 3
					md = checkfield(md,'fieldname','smb.cldFrac','NaN',1,'Inf',1,'>=',0,'<=',1);
				case 4
					md = checkfield(md,'fieldname','smb.t0wet','NaN',1,'Inf',1,'>=',15,'<=',21.9);
					md = checkfield(md,'fieldname','smb.t0dry','NaN',1,'Inf',1,'>=',30,'<=',30);
					md = checkfield(md,'fieldname','smb.K','NaN',1,'Inf',1,'>=',7,'<=',7);
			end

			%check zTop is < local thickness:
			he=sum(md.geometry.thickness(md.mesh.elements),2)/size(md.mesh.elements,2);
			if any(he<self.zTop),
				error('SMBgemb consistency check error: zTop should be smaller than local ice thickness');
			end
			md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);

		end % }}}
		function disp(self) % {{{
			
			disp(sprintf('   surface forcings for SMB GEMB model :'));
			
			fielddisplay(self,'isgraingrowth','run grain growth module (default true)');
			fielddisplay(self,'isalbedo','run albedo module (default true)');
			fielddisplay(self,'isshortwave','run short wave module (default true)');
			fielddisplay(self,'isthermal','run thermal module (default true)');
			fielddisplay(self,'isaccumulation','run accumulation module (default true)');
			fielddisplay(self,'ismelt','run melting  module (default true)');
			fielddisplay(self,'isdensification','run densification module (default true)');
			fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)');
			fielddisplay(self,'Ta','2 m air temperature, in Kelvin');
			fielddisplay(self,'V','wind speed (m/s-1)');
			fielddisplay(self,'dlwrf','downward shortwave radiation flux [W/m^2]');
			fielddisplay(self,'dswrf','downward longwave radiation flux [W/m^2]');
			fielddisplay(self,'P','precipitation [mm w.e. / m^2]');
			fielddisplay(self,'eAir','screen level vapor pressure [Pa]');
			fielddisplay(self,'pAir','surface pressure [Pa]');
			fielddisplay(self,'Tmean','mean annual temperature [K]');
			fielddisplay(self,'C','mean annual snow accumulation [kg m-2 yr-1]');
			fielddisplay(self,'Tz','height above ground at which temperature (T) was sampled [m]');
			fielddisplay(self,'Vz','height above ground at which wind (V) eas sampled [m]');
			fielddisplay(self,'zTop','depth over which grid length is constant at the top of the snopack (default 10) [m]');
			fielddisplay(self,'dzTop','initial top vertical grid spacing (default .05) [m] ');
			fielddisplay(self,'dzMin','initial min vertical allowable grid spacing (default dzMin/2) [m] ');
			fielddisplay(self,'zMax','initial max model depth (default is min(thickness,500)) [m]');
			fielddisplay(self,'zMin','initial min model depth (default is min(thickness,30)) [m]');
			fielddisplay(self,'zY','strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]');
			fielddisplay(self,'InitDensityScaling',{'initial scaling factor multiplying the density of ice','which describes the density of the snowpack.'});
			fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)');
			fielddisplay(self,'aIdx',{'method for calculating albedo and subsurface absorption (default is 1)',...
									'1: effective grain radius [Gardner & Sharp, 2009]',...
									'2: effective grain radius [Brun et al., 2009]',...
									'3: density and cloud amount [Greuell & Konzelmann, 1994]',...
									'4: exponential time decay & wetness [Bougamont & Bamber, 2005]'});
                                
			%snow properties init
			fielddisplay(self,'Dzini','Initial cel depth when restart [m]');
			fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]');
			fielddisplay(self,'Reini','Initial grain size when restart [mm]');
			fielddisplay(self,'Gdnini','Initial grain dendricity when restart [-]');
			fielddisplay(self,'Gspini','Initial grain sphericity when restart [-]');
			fielddisplay(self,'ECini','Initial evaporation/condensation when restart [kg m-2]');
			fielddisplay(self,'Wini','Initial snow water content when restart [kg m-2]');
			fielddisplay(self,'Aini','Initial albedo when restart [-]');
			fielddisplay(self,'Tini','Initial snow temperature when restart [K]');
			fielddisplay(self,'Sizeini','Initial number of layers when restart [K]');
            
            
			%additional albedo parameters: 
			switch self.aIdx
			case {1 2}
				fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)');
				fielddisplay(self,'aIce','albedo of ice (0.27-0.58)');
			case 3
				fielddisplay(self,'cldFrac','average cloud amount');
			case 4
				fielddisplay(self,'t0wet','time scale for wet snow (15-21.9) [d]');
				fielddisplay(self,'t0dry','warm snow timescale (30) [d]');
				fielddisplay(self,'K','time scale temperature coef. (7) [d]');
			end

			fielddisplay(self,'swIdx','apply all SW to top grid cell (0) or allow SW to penetrate surface (1) [default 1]');
			fielddisplay(self,'denIdx',{'densification model to use (default is 2):',...
									'1 = emperical model of Herron and Langway (1980)',...
									'2 = semi-emerical model of Anthern et al. (2010)',...
									'3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010)',...
									'4 = DO NOT USE: emperical model of Li and Zwally (2004)',...
									'5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)'});
			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',8,'format','Integer');
			
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','isgraingrowth','format','Boolean');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','isalbedo','format','Boolean');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','isshortwave','format','Boolean');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','isthermal','format','Boolean');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','isaccumulation','format','Boolean');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','ismelt','format','Boolean');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','isdensification','format','Boolean');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','isturbulentflux','format','Boolean');
            
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','V','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dswrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dlwrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','eAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);         

			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tmean','format','DoubleMat','mattype',2);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','C','format','DoubleMat','mattype',2);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tz','format','DoubleMat','mattype',2);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Vz','format','DoubleMat','mattype',2);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','zTop','format','DoubleMat','mattype',2);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dzTop','format','DoubleMat','mattype',2);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','dzMin','format','DoubleMat','mattype',2);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','zY','format','DoubleMat','mattype',2);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','zMax','format','DoubleMat','mattype',2);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','zMin','format','DoubleMat','mattype',2);
		
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','aIdx','format','Integer');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','swIdx','format','Integer');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','denIdx','format','Integer');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','InitDensityScaling','format','Double');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','outputFreq','format','Double');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','aSnow','format','Double');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','aIce','format','Double');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','cldFrac','format','Double');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0wet','format','Double');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double');
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double');
            
			%snow properties init
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzini','format','DoubleMat','mattype',3);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dini','format','DoubleMat','mattype',3);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Reini','format','DoubleMat','mattype',3);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gdnini','format','DoubleMat','mattype',3);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gspini','format','DoubleMat','mattype',3);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','ECini','format','DoubleMat','mattype',2);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Wini','format','DoubleMat','mattype',3);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Aini','format','DoubleMat','mattype',3);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tini','format','DoubleMat','mattype',3);
			WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizeini','format','IntMat','mattype',2);

			%figure out dt from forcings: 
			time=self.Ta(end,:); %assume all forcings are on the same time step
			dtime=diff(time,1);
			dt=min(dtime);
            
			WriteData(fid,prefix,'data',dt,'name','md.smb.dt','format','Double','scale',yts);
            
			% Check if smb_dt goes evenly into transient core time step
			if (mod(md.timestepping.time_step,dt) >= 1e-10)
                error('smb_dt/dt = %f. The number of SMB time steps in one transient core time step has to be an an integer',md.timestepping.time_step/dt);
			end
			
			%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
