Index: /issm/trunk-jpl/src/m/classes/snowpack.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/snowpack.m	(revision 15914)
+++ /issm/trunk-jpl/src/m/classes/snowpack.m	(revision 15914)
@@ -0,0 +1,367 @@
+%SNOWPACK class definition
+%
+%   Usage:
+%      snowpack=snowpack();
+
+classdef snowpack
+	properties (SetAccess=public) 
+
+		%first, the configuration fields, by category: 
+		%snowpack:  %{{{
+		snowpack_meas_tss = 0;
+		snowpack_enforce_measured_snow_heights = 0;
+		snowpack_sw_mode = 0;
+		snowpack_incoming_longwave = 0;
+		snowpack_height_of_wind_value = 0;
+		snowpack_height_of_meteo_values = 0;
+		snowpack_neutral = 0;
+		snowpack_roughness_length = 0;
+		snowpack_number_slopes = 0;
+		snowpack_snow_redistribution = 0;
+		snowpack_calculation_step_length = 0;
+		snowpack_change_bc = 0;
+		snowpack_thresh_change_bc = 0;
+		snowpack_snp_soil = 0;
+		snowpack_soil_flux = 0;
+		snowpack_geo_heat = 0;
+		snowpack_canopy = 0;
+		%}}}
+		%snowpackadvanced:  %{{{
+		snowpackadvanced_variant = ''; % use 320 kg m-3 for fixed density
+		snowpackadvanced_hn_density = '';
+		%}}}
+		%general:  %{{{
+		general_pluginpath = '';
+		general_buff_chunk_size = 0;
+		general_buff_before = 0;
+		%}}}
+		%input {{{
+		input_coordsys = '';
+		input_coordparam = '';
+		input_time_zone = 0;
+		input_meteo = '';
+		input_meteopath = '';
+		input_station1 = '';
+		input_snowfile1 = '';
+		%}}}
+		%output {{{
+		output_coordsys = '';
+		output_coordparam = '';
+		output_time_zone = 0;
+		output_meteopath = '';
+		output_experiment = '';
+		output_ts_write = 0;
+		output_ts_start = 0;
+		output_ts_days_between = 0;
+		output_profile = '';
+		output_prof_write = 0;
+		output_prof_start = 0;
+		output_prof_days_between = 0;
+		%}}}
+		%interpolations1d %{{{
+		interpolations1d_window_size = 0; %that is 5 d and 2 h; 1 d = 86400
+		interpolations1d_hnw_resample = '';
+		interpolations1d_hs_resample = '';
+		interpolations1d_tsg_resample = '';
+		interpolations1d_rho_hn_resample = '';
+		interpolations1d_vw_resample = '';
+		interpolations1d_vw_args = '';
+		%}}}
+		%filters {{{
+		filters_ta_filter1 = '';
+		filters_ta_arg1 = NaN;
+		filters_rh_filter1 = '';
+		filters_rh_arg1 = NaN;
+		filters_rh_filter2 = '';
+		filters_rh_arg2 = NaN;
+		filters_iswr_filter1 = '';
+		filters_iswr_arg1 = NaN;
+		filters_iswr_filter2 = '';
+		filters_iswr_arg2 = NaN;
+		filters_rswr_filter1 = '';
+		filters_rswr_arg1 = NaN;
+		filters_rswr_filter2 = '';
+		filters_rswr_arg2 = NaN;
+
+		%for ta between 190 and 280 k;
+		filters_ilwr_filter1 = '';
+		filters_ilwr_arg1 = NaN;
+		filters_ilwr_filter2 = '';
+		filters_ilwr_arg2 = NaN;
+		filters_tss_filter1 = '';
+		filters_tss_arg1 = NaN;
+		filters_tsg_filter1 = '';
+		filters_tsg_arg1 = NaN;
+		filters_vw_filter1 = '';
+		filters_vw_arg1 = NaN;
+		filters_vw_filter2 = '';
+		filters_vw_arg2 = NaN;
+		%}}}
+
+	end
+	methods
+		function obj = snowpack(varargin) % {{{
+			switch nargin
+				case 0
+					obj=setdefaultparameters(obj);
+				case 1
+					inputstruct=varargin{1};
+					list1 = properties('snowpack');
+					list2 = fieldnames(inputstruct);
+					for i=1:length(list1)
+						fieldname = list1{i};
+						if ismember(fieldname,list2),
+							obj.(fieldname) = inputstruct.(fieldname);
+						end
+					end
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function obj = setdefaultparameters(obj) % {{{
+	
+		%snowpack:  %{{{
+		obj.snowpack_meas_tss = 1;
+		obj.snowpack_enforce_measured_snow_heights = 0;
+		obj.snowpack_sw_mode = 0;
+		obj.snowpack_incoming_longwave = 1;
+		obj.snowpack_height_of_wind_value = 12.;
+		obj.snowpack_height_of_meteo_values = 12.;
+		obj.snowpack_neutral = 0;
+		obj.snowpack_roughness_length = 0.002;
+		obj.snowpack_number_slopes = 1;
+		obj.snowpack_snow_redistribution = 1;
+		obj.snowpack_calculation_step_length = 15.0;
+		obj.snowpack_change_bc = 0;
+		obj.snowpack_thresh_change_bc = -1.0;
+		obj.snowpack_snp_soil = 0;
+		obj.snowpack_soil_flux = 0;
+		obj.snowpack_geo_heat = 0.06;
+		obj.snowpack_canopy = 0;
+		%}}}
+		%snowpackadvanced:  %{{{
+		obj.snowpackadvanced_variant = 'ANTARCTICA'; % use 320 kg m-3 for fixed density
+		obj.snowpackadvanced_hn_density = 'EVENT';
+		%}}}
+		%general:  %{{{
+		obj.general_pluginpath = '/usr/local/lib/meteoio/plugins/';
+		obj.general_buff_chunk_size = 90;
+		obj.general_buff_before = 1.5;
+		%}}}
+		%input {{{
+		obj.input_coordsys = 'ch1903';
+		obj.input_coordparam = 'null';
+		obj.input_time_zone = 8;
+		obj.input_meteo = 'smet';
+		obj.input_meteopath = './input';
+		obj.input_station1 = 'domec.smet';
+		obj.input_snowfile1 = 'domec.sno';
+		%}}}
+		%output {{{
+		obj.output_coordsys = 'ch1903';
+		obj.output_coordparam = 'null';
+		obj.output_time_zone = 8;
+		obj.output_meteopath = './output';
+		obj.output_experiment = 'smet';
+		obj.output_ts_write = 1;
+		obj.output_ts_start = 0.0;
+		obj.output_ts_days_between = 0.04166667;
+		obj.output_profile = 'ascii';
+		obj.output_prof_write = 1;
+		obj.output_prof_start = 0.0;
+		obj.output_prof_days_between = 0.04166667;
+		%}}}
+		%interpolations1d %{{{
+		obj.interpolations1d_window_size = 439200 %that is 5 d and 2 h; 1 d = 86400
+		obj.interpolations1d_hnw_resample = 'none';
+		obj.interpolations1d_hs_resample = 'linear';
+		obj.interpolations1d_tsg_resample = 'linear';
+		obj.interpolations1d_rho_hn_resample = 'none';
+		obj.interpolations1d_vw_resample = 'nearest_neighbour';
+		obj.interpolations1d_vw_args = 'extrapolate';
+		%}}}
+		%filters {{{
+		obj.filters_ta_filter1 = 'min_max';
+		obj.filters_ta_arg1 = [190 280];
+		obj.filters_rh_filter1 = 'min_max';
+		obj.filters_rh_arg1 = [0.01 1.2];
+		obj.filters_rh_filter2 = 'min_max';
+		obj.filters_rh_arg2 = {'soft' 0.01 1.0};
+		obj.filters_iswr_filter1 = 'min_max';
+		obj.filters_iswr_arg1 = [-10 1500];
+		obj.filters_iswr_filter2 = 'min_max';
+		obj.filters_iswr_arg2 = {'soft' 0 1500};
+		obj.filters_rswr_filter1 = 'min_max';
+		obj.filters_rswr_arg1 = [-10 1500];
+		obj.filters_rswr_filter2 = 'min_max';
+		obj.filters_rswr_arg2 = {'soft' 0 1500};
+
+		%for ta between 190 and 280 k;
+		obj.filters_ilwr_filter1 = 'min_max';
+		obj.filters_ilwr_arg1 = [30 355];
+		obj.filters_ilwr_filter2 = 'min_max';
+		obj.filters_ilwr_arg2 = {'soft' 35 350};
+		obj.filters_tss_filter1 = 'min_max';
+		obj.filters_tss_arg1 = [180 275];
+		obj.filters_tsg_filter1 = 'min_max';
+		obj.filters_tsg_arg1 = [200 275];
+		obj.filters_vw_filter1 = 'min_max';
+		obj.filters_vw_arg1 = [-2 70];
+		obj.filters_vw_filter2 = 'min_max';
+		obj.filters_vw_arg2 = {'soft' 0 50};
+		%}}}
+
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+			%snowpack:  %{{{
+			md=checkfield(md,'snowpack.snowpack_meas_tss','values',[0 1]);
+			md=checkfield(md,'snowpack.snowpack_enforce_measured_snow_heights','values',[0 1]);
+			md=checkfield(md,'snowpack.snowpack_sw_mode','values',[0 1 2]);
+			md=checkfield(md,'snowpack.snowpack_incoming_longwave','values',[0 1]);
+			md=checkfield(md,'snowpack.snowpack_height_of_wind_value,'>=',0);
+			md=checkfield(md,'snowpack.snowpack_height_of_meteo_values','>=',0);
+			md=checkfield(md,'snowpack.snowpack_neutral','values',[-1 0 1]);
+			md=checkfield(md,'snowpack.snowpack_roughness_length','>=',0);
+			md=checkfield(md,'snowpack.snowpack_number_slopes','values',[1 3 5 9]);
+			md=checkfield(md,'snowpack.snowpack_snow_redistribution','values',[0 1]);
+			md=checkfield(md,'snowpack.snowpack_calculation_step_length','>',0);
+			md=checkfield(md,'snowpack.snowpack_change_bc','values',[0 1]);
+			md=checkfield(md,'snowpack.snowpack_thresh_change_bc','<=',0);
+			md=checkfield(md,'snowpack.snowpack_snp_soil','values',[0 1]);
+			md=checkfield(md,'snowpack.snowpack_soil_flux,'values',[0 1]);
+			md=checkfield(md,'snowpack.snowpack_geo_heat,'>=',0);
+			md=checkfield(md,'snowpack.snowpack_canopy,'values',[0 1]);
+			%}}}
+			%snowpackadvanced:  %{{{
+			md=checkfield(md,'snowpack.snowpackadvanced_variant = 'ANTARCTICA'; % use 320 kg m-3 for fixed density
+			md=checkfield(md,'snowpack.snowpackadvanced_hn_density = 'EVENT';
+			%}}}
+			%general:  %{{{
+			md=checkfield(md,'snowpack.general_pluginpath = '/usr/local/lib/meteoio/plugins/';
+			md=checkfield(md,'snowpack.general_buff_chunk_size = 90;
+			md=checkfield(md,'snowpack.general_buff_before = 1.5;
+			%}}}
+			%input {{{
+			md=checkfield(md,'snowpack.input_coordsys = 'ch1903';
+			md=checkfield(md,'snowpack.input_coordparam = 'null';
+			md=checkfield(md,'snowpack.input_time_zone = 8;
+			md=checkfield(md,'snowpack.input_meteo = 'smet';
+			md=checkfield(md,'snowpack.input_meteopath = './input';
+			md=checkfield(md,'snowpack.input_station1 = 'domec.smet';
+			md=checkfield(md,'snowpack.input_snowfile1 = 'domec.sno';
+			%}}}
+			%output {{{
+			md=checkfield(md,'snowpack.output_coordsys = 'ch1903';
+			md=checkfield(md,'snowpack.output_coordparam = 'null';
+			md=checkfield(md,'snowpack.output_time_zone = 8;
+			md=checkfield(md,'snowpack.output_meteopath = './output';
+			md=checkfield(md,'snowpack.output_experiment = 'smet';
+			md=checkfield(md,'snowpack.output_ts_write = 1;
+			md=checkfield(md,'snowpack.output_ts_start = 0.0;
+			md=checkfield(md,'snowpack.output_ts_days_between = 0.04166667;
+			md=checkfield(md,'snowpack.output_profile = 'ascii';
+			md=checkfield(md,'snowpack.output_prof_write = 1;
+			md=checkfield(md,'snowpack.output_prof_start = 0.0;
+			md=checkfield(md,'snowpack.output_prof_days_between = 0.04166667;
+			%}}}
+			%interpolations1d %{{{
+			md=checkfield(md,'snowpack.interpolations1d_window_size = 439200 %that is 5 d and 2 h; 1 d = 86400
+			md=checkfield(md,'snowpack.interpolations1d_hnw_resample = 'none';
+			md=checkfield(md,'snowpack.interpolations1d_hs_resample = 'linear';
+			md=checkfield(md,'snowpack.interpolations1d_tsg_resample = 'linear';
+			md=checkfield(md,'snowpack.interpolations1d_rho_hn_resample = 'none';
+			md=checkfield(md,'snowpack.interpolations1d_vw_resample = 'nearest_neighbour';
+			md=checkfield(md,'snowpack.interpolations1d_vw_args = 'extrapolate';
+			%}}}
+			%filters {{{
+			md=checkfield(md,'snowpack.filters_ta_filter1 = 'min_max';
+			md=checkfield(md,'snowpack.filters_ta_arg1 = [190 280];
+			md=checkfield(md,'snowpack.filters_rh_filter1 = 'min_max';
+			md=checkfield(md,'snowpack.filters_rh_arg1 = [0.01 1.2];
+			md=checkfield(md,'snowpack.filters_rh_filter2 = 'min_max';
+			md=checkfield(md,'snowpack.filters_rh_arg2 = {'soft' 0.01 1.0};
+			md=checkfield(md,'snowpack.filters_iswr_filter1 = 'min_max';
+			md=checkfield(md,'snowpack.filters_iswr_arg1 = [-10 1500];
+			md=checkfield(md,'snowpack.filters_iswr_filter2 = 'min_max';
+			md=checkfield(md,'snowpack.filters_iswr_arg2 = {'soft' 0 1500};
+			md=checkfield(md,'snowpack.filters_rswr_filter1 = 'min_max';
+			md=checkfield(md,'snowpack.filters_rswr_arg1 = [-10 1500];
+			md=checkfield(md,'snowpack.filters_rswr_filter2 = 'min_max';
+			md=checkfield(md,'snowpack.filters_rswr_arg2 = {'soft' 0 1500};
+
+			%for ta between 190 and 280 k;
+			md=checkfield(md,'snowpack.filters_ilwr_filter1 = 'min_max';
+			md=checkfield(md,'snowpack.filters_ilwr_arg1 = [30 355];
+			md=checkfield(md,'snowpack.filters_ilwr_filter2 = 'min_max';
+			md=checkfield(md,'snowpack.filters_ilwr_arg2 = {'soft' 35 350};
+			md=checkfield(md,'snowpack.filters_tss_filter1 = 'min_max';
+			md=checkfield(md,'snowpack.filters_tss_arg1 = [180 275];
+			md=checkfield(md,'snowpack.filters_tsg_filter1 = 'min_max';
+			md=checkfield(md,'snowpack.filters_tsg_arg1 = [200 275];
+			md=checkfield(md,'snowpack.filters_vw_filter1 = 'min_max';
+			md=checkfield(md,'snowpack.filters_vw_arg1 = [-2 70];
+			md=checkfield(md,'snowpack.filters_vw_filter2 = 'min_max';
+			md=checkfield(md,'snowpack.filters_vw_arg2 = {'soft' 0 50};
+			%}}}
+
+		end % }}}
+		function disp(obj) % {{{
+
+			disp(sprintf('   StressBalance solution parameters:'));
+
+			disp(sprintf('\n      %s','Convergence criteria:'));
+			fielddisplay(obj,'restol','mechanical equilibrium residual convergence criterion');
+			fielddisplay(obj,'reltol','velocity relative convergence criterion, NaN: not applied');
+			fielddisplay(obj,'abstol','velocity absolute convergence criterion, NaN: not applied');
+			fielddisplay(obj,'isnewton','0: Picard''s fixed point, 1: Newton''s method, 2: hybrid');
+			fielddisplay(obj,'maxiter','maximum number of nonlinear iterations');
+			fielddisplay(obj,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)');
+
+			disp(sprintf('\n      %s','boundary conditions:'));
+			fielddisplay(obj,'spcvx','x-axis velocity constraint (NaN means no constraint) [m/yr]');
+			fielddisplay(obj,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]');
+			fielddisplay(obj,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]');
+
+			disp(sprintf('\n      %s','Rift options:'));
+			fielddisplay(obj,'rift_penalty_threshold','threshold for instability of mechanical constraints');
+			fielddisplay(obj,'rift_penalty_lock','number of iterations before rift penalties are locked');
+
+			disp(sprintf('\n      %s','Penalty options:'));
+			fielddisplay(obj,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset');
+			fielddisplay(obj,'vertex_pairing','pairs of vertices that are penalized');
+
+			disp(sprintf('\n      %s','Other:'));
+			fielddisplay(obj,'shelf_dampening','use dampening for floating ice ? Only for FS model');
+			fielddisplay(obj,'FSreconditioning','multiplier for incompressibility equation. Only for FS model');
+			fielddisplay(obj,'referential','local referential');
+			fielddisplay(obj,'loadingforce','loading force applied on each point [N/m^3]');
+			fielddisplay(obj,'requested_outputs','additional outputs requested');
+
+		end % }}}
+		function marshall(obj,md,fid) % {{{
+
+			yts=365.0*24.0*3600.0;
+
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','spcvx','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','spcvy','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','spcvz','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','restol','format','Double');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','reltol','format','Double');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','abstol','format','Double');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','isnewton','format','Integer');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','FSreconditioning','format','Double');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','viscosity_overshoot','format','Double');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','maxiter','format','Integer');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','shelf_dampening','format','Integer');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','vertex_pairing','format','DoubleMat','mattype',3);
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','penalty_factor','format','Double');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','rift_penalty_lock','format','Integer');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','rift_penalty_threshold','format','Integer');
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','referential','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'class','snowpack','fieldname','requested_outputs','format','DoubleMat','mattype',3);
+			WriteData(fid,'data',obj.loadingforce(:,1),'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum);
+			WriteData(fid,'data',obj.loadingforce(:,2),'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum);
+			WriteData(fid,'data',obj.loadingforce(:,3),'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum);
+		end % }}}
+	end
+end