/*!\file: CreateParameters.cpp
 * \brief general driver for creating parameters dataset
 */ 

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "../../toolkits/toolkits.h"
#include "../../classes/classes.h"
#include "../../shared/shared.h"
#include "../MeshPartitionx/MeshPartitionx.h"
#include "../ParseToolkitsOptionsx/ParseToolkitsOptionsx.h"
#include "./ModelProcessorx.h"

void CreateParameters(Parameters* parameters,IoModel* iomodel,char* rootpath,FILE* toolkitsoptionsfid,const int solution_type){

	int         i,j,m,k;
	int         numoutputs,basalforcing_model,timestepping_type;
	char**      requestedoutputs = NULL;
	char**      outputonnodes = NULL;
	char*       fieldname = NULL;
	IssmDouble  time;

	/*parameters for mass flux:*/
	int          mass_flux_num_profiles     = 0;
	bool         qmu_mass_flux_present      = false;
	bool         autodiff_mass_flux_present = false;
	bool         mass_flux_present          = false;
	bool         interp;
	IssmDouble **array                      = NULL;
	int         *mdims_array                = NULL;
	int         *ndims_array                = NULL;
	IssmDouble  *temp_matrix                = NULL;
	int          temp_m,temp_n;
	IssmDouble  *matrix                     = NULL;
	int          count;

	IssmDouble *temp = NULL;
	IssmDouble *transparam = NULL;
	IssmDouble  yts;
	int         N,M;

	/*Copy some constants from iomodel */
	parameters->AddObject(iomodel->CopyConstantObject("md.mesh.domain_type",DomainTypeEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.mesh.domain_dimension",DomainDimensionEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.settings.output_frequency",SettingsOutputFrequencyEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.settings.sb_coupling_frequency",SettingsSbCouplingFrequencyEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.settings.recording_frequency",SettingsRecordingFrequencyEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.constants.yts",ConstantsYtsEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.debug.profiling",DebugProfilingEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.mesh.average_vertex_connectivity",MeshAverageVertexConnectivityEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.settings.waitonlock",SettingsWaitonlockEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.mesh.numberofvertices",MeshNumberofverticesEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.settings.io_gather",SettingsIoGatherEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.settings.solver_residue_threshold",SettingsSolverResidueThresholdEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.isautodiff",AutodiffIsautodiffEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.qmu.isdakota",QmuIsdakotaEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.inversion.iscontrol",InversionIscontrolEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum));
	parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.parameterization",FrontalForcingsParamEnum));
	parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1));  

	  {/*This is specific to ice...*/
		parameters->AddObject(iomodel->CopyConstantObject("md.mesh.elementtype",MeshElementtypeEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.steadystate.reltol",SteadystateReltolEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.steadystate.maxiter",SteadystateMaxiterEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.migration",GroundinglineMigrationEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.friction_interpolation",GroundinglineFrictionInterpolationEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.groundingline.melt_interpolation",GroundinglineMeltInterpolationEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.isstressbalance",TransientIsstressbalanceEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.ismasstransport",TransientIsmasstransportEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.issmb",TransientIssmbEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.isthermal",TransientIsthermalEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.isgroundingline",TransientIsgroundinglineEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.isgia",TransientIsgiaEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.isesa",TransientIsesaEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.isdamageevolution",TransientIsdamageevolutionEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.ishydrology",TransientIshydrologyEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.ismovingfront",TransientIsmovingfrontEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.isslr",TransientIsslrEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.iscoupler",TransientIscouplerEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.isoceancoupling",TransientIsoceancouplingEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum));

		/*For stress balance only*/
		parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isFS",FlowequationIsFSEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.rift_penalty_threshold",StressbalanceRiftPenaltyThresholdEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.maxiter",StressbalanceMaxiterEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.restol",StressbalanceRestolEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.reltol",StressbalanceReltolEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.abstol",StressbalanceAbstolEnum));
		if(iomodel->domaintype==Domain3DEnum)
		 parameters->AddObject(iomodel->CopyConstantObject("md.mesh.numberoflayers",MeshNumberoflayersEnum));
	  }

	/*amr properties*/	
	int amrtype,amr_frequency;
	iomodel->FindConstant(&amr_frequency,"md.transient.amr_frequency");
	if(solution_type==TransientSolutionEnum && amr_frequency){
		/*Load common amr parameters*/
		parameters->AddObject(iomodel->CopyConstantObject("md.amr.type",AmrTypeEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_distance",AmrGroundingLineDistanceEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_distance",AmrIceFrontDistanceEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_threshold",AmrThicknessErrorThresholdEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_groupthreshold",AmrThicknessErrorGroupThresholdEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_maximum",AmrThicknessErrorMaximumEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_threshold",AmrDeviatoricErrorThresholdEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_groupthreshold",AmrDeviatoricErrorGroupThresholdEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_maximum",AmrDeviatoricErrorMaximumEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.amr.restart",AmrRestartEnum));
		/*Load specific amr parameters*/
		iomodel->FindConstant(&amrtype,"md.amr.type");
		switch(amrtype){
			#ifdef _HAVE_NEOPZ_
			case AmrNeopzEnum:
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.level_max",AmrLevelMaxEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.gradation",AmrGradationEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.lag",AmrLagEnum));
				break;
			#endif

			#ifdef _HAVE_BAMG_
			case AmrBamgEnum:
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.hmin",AmrHminEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.hmax",AmrHmaxEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.err",AmrErrEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.keepmetric",AmrKeepMetricEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.gradation",AmrGradationEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.groundingline_resolution",AmrGroundingLineResolutionEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.icefront_resolution",AmrIceFrontResolutionEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.thicknesserror_resolution",AmrThicknessErrorResolutionEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.amr.deviatoricerror_resolution",AmrDeviatoricErrorResolutionEnum));
				/*Convert fieldname to enum and put it in params*/
				iomodel->FindConstant(&fieldname,"md.amr.fieldname");
				parameters->AddObject(new IntParam(AmrFieldEnum,StringToEnumx(fieldname)));
				xDelete<char>(fieldname);
				break;
			#endif

			default:
				_error_("Adaptive mesh refinement "<<EnumToStringx(amrtype)<<" not implemented yet");
		}
	}

	/*Basal forcing parameters*/
	parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.model",BasalforcingsEnum));
	iomodel->FindConstant(&basalforcing_model,"md.basalforcings.model");
	switch(basalforcing_model){
		case FloatingMeltRateEnum:
			/*Nothing to add to parameters*/
			break;
		case LinearFloatingMeltRateEnum:
			iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
			iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_melting_rate"); 
			if(N==1){
				_assert_(M==1);
				parameters->AddObject(new DoubleParam(BasalforcingsDeepwaterMeltingRateEnum,transparam[0]));
			}
			else{
				_assert_(N==2);
				parameters->AddObject(new TransientParam(BasalforcingsDeepwaterMeltingRateEnum,&transparam[0],&transparam[M],interp,M));
			}
			xDelete<IssmDouble>(transparam);
			iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_melting_rate"); 
			if(N==1){
				_assert_(M==1);
				parameters->AddObject(new DoubleParam(BasalforcingsUpperwaterMeltingRateEnum,transparam[0]));
			}
			else{
				_assert_(N==2);
				parameters->AddObject(new TransientParam(BasalforcingsUpperwaterMeltingRateEnum,&transparam[0],&transparam[M],interp,M));
			}
			xDelete<IssmDouble>(transparam);
			iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_elevation"); 
			if(N==1){
				_assert_(M==1);
				parameters->AddObject(new DoubleParam(BasalforcingsDeepwaterElevationEnum,transparam[0]));
			}
			else{
				_assert_(N==2);
				parameters->AddObject(new TransientParam(BasalforcingsDeepwaterElevationEnum,&transparam[0],&transparam[M],interp,M)); 
			}
			xDelete<IssmDouble>(transparam);
			iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_elevation"); 
			if(N==1){
				_assert_(M==1);
				parameters->AddObject(new DoubleParam(BasalforcingsUpperwaterElevationEnum,transparam[0]));
			}
			else{
				_assert_(N==2);
				parameters->AddObject(new TransientParam(BasalforcingsUpperwaterElevationEnum,&transparam[0],&transparam[M],interp,M));
			}
			xDelete<IssmDouble>(transparam);
			break;
		case SpatialLinearFloatingMeltRateEnum:
			/*Nothing to add to parameters:*/
			break;
		case MismipFloatingMeltRateEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.threshold_thickness",BasalforcingsThresholdThicknessEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.upperdepth_melt",BasalforcingsUpperdepthMeltEnum));
			break;
		case MantlePlumeGeothermalFluxEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.mantleconductivity",BasalforcingsMantleconductivityEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.nusselt",BasalforcingsNusseltEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.dtbg",BasalforcingsDtbgEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.plumeradius",BasalforcingsPlumeradiusEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.topplumedepth",BasalforcingsTopplumedepthEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.bottomplumedepth",BasalforcingsBottomplumedepthEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.plumex",BasalforcingsPlumexEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.plumey",BasalforcingsPlumeyEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.crustthickness",BasalforcingsCrustthicknessEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.uppercrustthickness",BasalforcingsUppercrustthicknessEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.uppercrustheat",BasalforcingsUppercrustheatEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.lowercrustheat",BasalforcingsLowercrustheatEnum));
			break;
		case BasalforcingsPicoEnum:
			iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsPicoNumBasinsEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.maxboxcount",BasalforcingsPicoMaxboxcountEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_T",BasalforcingsPicoGammaTEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.isplume",BasalforcingsPicoIsplumeEnum));
			iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_temperature");
			_assert_(M>=1 && N>=1); 
			parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceantemperatureEnum,transparam,&transparam[N*(M-1)],interp,N,M));
			xDelete<IssmDouble>(transparam);
			iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_salinity");
			_assert_(M>=1 && N>=1); 
			parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceansalinityEnum,transparam,&transparam[N*(M-1)],interp,N,M));
			xDelete<IssmDouble>(transparam);
			break;
		case BasalforcingsIsmip6Enum:
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsIsmip6NumBasinsEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_0",BasalforcingsIsmip6Gamma0Enum));
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmip6IsLocalEnum));	
			iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.delta_t");
			parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6DeltaTEnum,transparam,N));
			xDelete<IssmDouble>(transparam);
			iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.tf_depths");
			parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6TfDepthsEnum,transparam,N));
			xDelete<IssmDouble>(transparam);
			break;
		case BeckmannGoosseFloatingMeltRateEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum));
			break;
		default:
			_error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");
	}

	/*some parameters that did not come with the iomodel: */
	parameters->AddObject(new IntParam(SolutionTypeEnum,solution_type));

	/*Time stepping*/
	parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.type",TimesteppingTypeEnum));
	iomodel->FindConstant(&timestepping_type,"md.timestepping.type");
	switch(timestepping_type){
		case FixedTimesteppingEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.start_time",TimesteppingStartTimeEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.final_time",TimesteppingFinalTimeEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.time_step",TimesteppingTimeStepEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.interp_forcings",TimesteppingInterpForcingsEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.coupling_time",TimesteppingCouplingTimeEnum));
			break;
		case AdaptiveTimesteppingEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.start_time",TimesteppingStartTimeEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.final_time",TimesteppingFinalTimeEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.time_step_min",TimesteppingTimeStepMinEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.time_step_max",TimesteppingTimeStepMaxEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.cfl_coefficient",TimesteppingCflCoefficientEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.interp_forcings",TimesteppingInterpForcingsEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.timestepping.coupling_time",TimesteppingCouplingTimeEnum));
			break;
		default:
			_error_("Time stepping \""<<EnumToStringx(timestepping_type)<<"\" not supported yet");
	}
	iomodel->FindConstant(&time,"md.timestepping.start_time");
	parameters->AddObject(new DoubleParam(TimeEnum,time));  
	parameters->AddObject(new IntParam(StepEnum,0));  

	/*By default, save all results*/
	parameters->AddObject(new BoolParam(SaveResultsEnum,true));

	/*Should we output results on nodes?*/
	iomodel->FindConstant(&outputonnodes,&numoutputs,"md.settings.results_on_nodes");
	parameters->AddObject(new IntParam(SettingsNumResultsOnNodesEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new StringArrayParam(SettingsResultsOnNodesEnum,outputonnodes,numoutputs));
	iomodel->DeleteData(&outputonnodes,numoutputs,"md.settings.results_on_nodes");

	/*Requested outputs */
	iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.transient.requested_outputs");
	parameters->AddObject(new IntParam(TransientNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new StringArrayParam(TransientRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(&requestedoutputs,numoutputs,"md.transient.requested_outputs");

	iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.steadystate.requested_outputs");
	parameters->AddObject(new IntParam(SteadystateNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new StringArrayParam(SteadystateRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(&requestedoutputs,numoutputs,"md.steadystate.requested_outputs");

	int materialstype;
	iomodel->FindConstant(&materialstype,"md.materials.type");

	switch(materialstype){
		case MaticeEnum:
		case MatdamageiceEnum:
		case MatenhancediceEnum:
		case MatestarEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_ice",MaterialsRhoIceEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_water",MaterialsRhoSeawaterEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_freshwater",MaterialsRhoFreshwaterEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.mu_water",MaterialsMuWaterEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.heatcapacity",MaterialsHeatcapacityEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermalconductivity",MaterialsThermalconductivityEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.temperateiceconductivity",MaterialsTemperateiceconductivityEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.effectiveconductivity_averaging",MaterialsEffectiveconductivityAveragingEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.latentheat",MaterialsLatentheatEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.beta",MaterialsBetaEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.meltingpoint",MaterialsMeltingpointEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.constants.referencetemperature",ConstantsReferencetemperatureEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.mixed_layer_capacity",MaterialsMixedLayerCapacityEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermal_exchange_velocity",MaterialsThermalExchangeVelocityEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.constants.g",ConstantsGEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.rheology_law",MaterialsRheologyLawEnum));

			/*gia: */
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.lithosphere_shear_modulus",MaterialsLithosphereShearModulusEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.lithosphere_density",MaterialsLithosphereDensityEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.mantle_shear_modulus",MaterialsMantleShearModulusEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.mantle_density",MaterialsMantleDensityEnum));

			/*slr:*/
			parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum));
			break;
		case MaterialsEnum:{
			int nnat,dummy;
			int* nature=NULL;
			iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature");
			for(int i=0;i<nnat;i++){
				switch(IoCodeToEnumMaterials(nature[i])){
					case MatlithoEnum:
						break;
					case MaticeEnum:
					case MatdamageiceEnum:
					case MatenhancediceEnum:
					case MatestarEnum:
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_ice",MaterialsRhoIceEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_water",MaterialsRhoSeawaterEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_freshwater",MaterialsRhoFreshwaterEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.mu_water",MaterialsMuWaterEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.heatcapacity",MaterialsHeatcapacityEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermalconductivity",MaterialsThermalconductivityEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.temperateiceconductivity",MaterialsTemperateiceconductivityEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.effectiveconductivity_averaging",MaterialsEffectiveconductivityAveragingEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.latentheat",MaterialsLatentheatEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.beta",MaterialsBetaEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.meltingpoint",MaterialsMeltingpointEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.constants.referencetemperature",ConstantsReferencetemperatureEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.mixed_layer_capacity",MaterialsMixedLayerCapacityEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermal_exchange_velocity",MaterialsThermalExchangeVelocityEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.constants.g",ConstantsGEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.rheology_law",MaterialsRheologyLawEnum));

						/*gia: */
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.lithosphere_shear_modulus",MaterialsLithosphereShearModulusEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.lithosphere_density",MaterialsLithosphereDensityEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.mantle_shear_modulus",MaterialsMantleShearModulusEnum));
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.mantle_density",MaterialsMantleDensityEnum));

						/*slr:*/
						parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum));
						break;
				}
			}
			break;

	}
		default:
			_error_("Material "<< EnumToStringx(materialstype) <<" not supported yet");
	}

	int smb_model;
	iomodel->FindConstant(&smb_model,"md.smb.model");
	switch(smb_model){
		case SMBforcingEnum:
		case SMBgradientsEnum:
		case SMBgradientselaEnum:
		case SMBhenningEnum:
		case SMBcomponentsEnum:
		case SMBmeltcomponentsEnum:
		case SMBgradientscomponentsEnum:
			/*Nothing to add*/
			break;
		case SMBgembEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.aSnow",SmbASnowEnum));
			break;
		case SMBpddEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlapslgm",SmbRlapslgmEnum));
			break;
		case SMBpddSicopolisEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
			break;
		case SMBd18opddEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlapslgm",SmbRlapslgmEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.dpermil",SmbDpermilEnum));
			break;
		case SMBsemicEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfac",SmbDesfacEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.rdl",SmbRdlEnum));
			break;
		default:
			_error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
	}

	int hydrology_model;
	iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
	if(hydrology_model==HydrologydcEnum){
		/*FIXME: this cshould go to Analysis!!!*/
		parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_compressibility",HydrologydcSedimentCompressibilityEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_porosity",HydrologydcSedimentPorosityEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_thickness",HydrologydcSedimentThicknessEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.water_compressibility",HydrologydcWaterCompressibilityEnum));
		parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.isefficientlayer",HydrologydcIsefficientlayerEnum));

		bool isefficientlayer;
		iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
		if(isefficientlayer){
			parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_compressibility",HydrologydcEplCompressibilityEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_porosity",HydrologydcEplPorosityEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_initial_thickness",HydrologydcEplInitialThicknessEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_colapse_thickness",HydrologydcEplColapseThicknessEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_max_thickness",HydrologydcEplMaxThicknessEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_conductivity",HydrologydcEplConductivityEnum));
		}
	}
	else if(hydrology_model==HydrologyshreveEnum){
		/*Nothing to add*/
	}
	else if(hydrology_model==HydrologyshaktiEnum){
		/*Nothing to add*/
	}
	else if(hydrology_model==HydrologypismEnum){
		/*Nothing to add*/
	}
	else if(hydrology_model==HydrologyGlaDSEnum){
		/*Nothing to add*/
	}
	else{
		_error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");
	}

	if(materialstype==MatdamageiceEnum){
		iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.damage.requested_outputs");
		parameters->AddObject(new IntParam(DamageEvolutionNumRequestedOutputsEnum,numoutputs));
		if(numoutputs)parameters->AddObject(new StringArrayParam(DamageEvolutionRequestedOutputsEnum,requestedoutputs,numoutputs));
		iomodel->DeleteData(&requestedoutputs,numoutputs,"md.damage.requested_outputs");
	}

	/*Deal with mass flux segments: {{{*/
	iomodel->FetchData(&qmu_mass_flux_present,"md.qmu.mass_flux_segments_present");
	iomodel->FetchData(&autodiff_mass_flux_present,"md.autodiff.mass_flux_segments_present");

	if(qmu_mass_flux_present || autodiff_mass_flux_present)mass_flux_present=true;
	else mass_flux_present=false;
	parameters->AddObject(new BoolParam(MassFluxSegmentsPresentEnum,mass_flux_present));

	if(mass_flux_present){

		/*Fetch the mass flux segments necessary to compute the mass fluxes.  Build a DoubleMatArrayParam object out of them: */ 
		iomodel->FetchData(&array,&mdims_array,&ndims_array,&mass_flux_num_profiles,"md.qmu.mass_flux_segments");
		if(mass_flux_num_profiles==0)_error_("mass_flux_num_profiles is 0, when MassFlux computations were requested!");

		/*Go through segments, and extract those that belong to this cpu: */
		for(i=0;i<mass_flux_num_profiles;i++){
			temp_matrix=array[i];
			temp_m=mdims_array[i];
			temp_n=ndims_array[i];
			_assert_(temp_n==5);

			m=0;
			for(j=0;j<temp_m;j++){
				if (  iomodel->my_elements[reCast<int>(*(temp_matrix+5*j+4))-1] )m++;
			}
			if(m){
				matrix=xNewZeroInit<IssmDouble>(5*m);
				count=0;
				for(j=0;j<temp_m;j++){
					if (iomodel->my_elements[reCast<int>(*(temp_matrix+5*j+4))-1]){
						for(k=0;k<5;k++)*(matrix+5*count+k)=*(temp_matrix+5*j+k);
						count++;
					}
				}
			}
			else{
				matrix=NULL;
			}

			/*Assign: */
			array[i]=matrix;
			mdims_array[i]=m;
			ndims_array[i]=5;

			/*Free temporary matrix: */
			xDelete<IssmDouble>(temp_matrix);
		}

		/*Ok, we have an array of segments, different on every cpu. Create a DoubleMatArrayParam object with it: */
		parameters->AddObject(new DoubleMatArrayParam(MassFluxSegmentsEnum,array,mass_flux_num_profiles,mdims_array,ndims_array));

		/*Free data: */
		for(i=0;i<mass_flux_num_profiles;i++){
			IssmDouble* matrix=array[i];
			xDelete<IssmDouble>(matrix);
		}
		xDelete<int>(mdims_array); 
		xDelete<int>(ndims_array);
		xDelete<IssmDouble*>(array);
	}
	/*}}}*/

	/*Before returning, create parameters in case we are running Qmu or control types runs: */
	CreateParametersControl(parameters,iomodel,solution_type);

	#ifdef _HAVE_DAKOTA_
	CreateParametersDakota(parameters,iomodel,rootpath);
	#endif

	/*Now, deal with toolkits options, which need to be put into the parameters dataset: */
	ParseToolkitsOptionsx(parameters,toolkitsoptionsfid);

	#ifdef _HAVE_AD_
	if(VerboseMProcessor()) _printf0_("   starting autodiff parameters \n");
	CreateParametersAutodiff(parameters,iomodel);
	if(VerboseMProcessor()) _printf0_("   ending autodiff parameters \n");
	#endif
}
