#include "./SmbAnalysis.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"

// FIX
#include "./shared/io/Print/Print.h"

/*Model processing*/
void SmbAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
	/*No constraints*/
}/*}}}*/
void SmbAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
	/*No loads*/
}/*}}}*/
void SmbAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
	::CreateNodes(nodes,iomodel,SmbAnalysisEnum,P1Enum);
}/*}}}*/
int  SmbAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
	return 1;
}/*}}}*/
void SmbAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/

	int    smb_model;
	bool   isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled,isfirnwarming,isstochastic;

	/*Update elements: */
	int counter=0;
	for(int i=0;i<iomodel->numberofelements;i++){
		if(iomodel->my_elements[i]){
			Element* element=(Element*)elements->GetObjectByOffset(counter);
			element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum);
			counter++;
		}
	}

	/*Figure out smb model: */
	iomodel->FindConstant(&smb_model,"md.smb.model");
	iomodel->FindConstant(&isstochastic,"md.stochasticforcing.isstochasticforcing");
	switch(smb_model){
		case SMBforcingEnum:
			iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
			if(isstochastic){
				iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
				iomodel->FetchDataToInput(inputs,elements,"md.smb.mass_balance",BaselineSmbMassBalanceEnum,0.);
			}
			break;
		case SMBgembEnum:
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Ta",SmbTaEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.V",SmbVEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dswrf",SmbDswrfEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dswdiffrf",SmbDswdiffrfEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dlwrf",SmbDlwrfEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.P",SmbPEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.eAir",SmbEAirEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.pAir",SmbPAirEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.zTop",SmbZTopEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dzTop",SmbDzTopEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dzMin",SmbDzMinEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.zY",SmbZYEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.zMax",SmbZMaxEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.zMin",SmbZMinEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Tmean",SmbTmeanEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Vmean",SmbVmeanEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.C",SmbCEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Tz",SmbTzEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Vz",SmbVzEnum);
			InputUpdateFromConstantx(inputs,elements,false,SmbIsInitializedEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Dzini",SmbDziniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Dini",SmbDiniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Reini",SmbReiniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Gdnini",SmbGdniniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Gspini",SmbGspiniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.ECini",SmbECiniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Wini",SmbWiniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Aini",SmbAiniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Adiffini",SmbAdiffiniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Tini",SmbTiniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.Sizeini",SmbSizeiniEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.aValue",SmbAValueEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dulwrfValue",SmbDulwrfValueEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.teValue",SmbTeValueEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.szaValue",SmbSzaValueEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.cotValue",SmbCotValueEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.ccsnowValue",SmbCcsnowValueEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.cciceValue",SmbCciceValueEnum);
			break;
		case SMBpddEnum:
			iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
			iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
			iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
			iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
			if(isdelta18o || ismungsm){
				iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_lgm",SmbTemperaturesLgmEnum);
				iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
				iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
				iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum);
			}else{
				iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitation",SmbPrecipitationEnum);
				iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
			}
			if(issetpddfac){
				iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
				iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
			}
			break;
		case SMBpddSicopolisEnum:
			iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
			iomodel->FindConstant(&isfirnwarming,"md.smb.isfirnwarming");
			iomodel->FetchDataToInput(inputs,elements,"md.smb.smb_corr",SmbSmbCorrEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum);
			iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
			iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitation",SmbPrecipitationEnum);
			break;
		case SMBd18opddEnum:
			iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
			iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
			iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
			iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
			iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
			iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.s0p",SmbS0pEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
			if(isd18opd){
				iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
				iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
				if(!istemperaturescaled){
					/*Fetch array*/
					IssmDouble* doublearray = NULL;
					int         M,N;
					iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
					if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
					if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");

					/*Build indices*/
					int* ids = xNew<int>(N); for(int i=0;i<N;i++) ids[i] = i;

					for(Object* & object : elements->objects){
						Element* element=xDynamicCast<Element*>(object);
						element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs,iomodel,SmbTemperaturesReconstructedEnum);
					}
					xDelete<int>(ids);
					iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
				}
				if(!isprecipscaled){
					/*Fetch array*/
					IssmDouble* doublearray = NULL;
					int         M,N;
					iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
					if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
					if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");

					/*Build indices*/
					int* ids = xNew<int>(N); for(int i=0;i<N;i++) ids[i] = i;

					for(Object* & object : elements->objects){
						Element* element=xDynamicCast<Element*>(object);
						element->DatasetInputCreate(doublearray,M-1,N,ids,N,inputs,iomodel,SmbPrecipitationsReconstructedEnum);
					}
					xDelete<int>(ids);
					iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
				}
			}
			if(issetpddfac){
				iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
				iomodel->FetchDataToInput(inputs,elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
			}
			break;
		case SMBgradientsEnum:
			iomodel->FetchDataToInput(inputs,elements,"md.smb.href",SmbHrefEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.b_pos",SmbBPosEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.b_neg",SmbBNegEnum);
			break;
		case SMBgradientselaEnum:
			iomodel->FetchDataToInput(inputs,elements,"md.smb.ela",SmbElaEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.b_pos",SmbBPosEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.b_neg",SmbBNegEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.b_max",SmbBMaxEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.b_min",SmbBMinEnum);
			break;
		case SMBautoregressionEnum:
         iomodel->FetchDataToInput(inputs,elements,"md.smb.basin_id",SmbBasinsIdEnum);
         break;
		case SMBhenningEnum:
			iomodel->FetchDataToInput(inputs,elements,"md.smb.smbref",SmbSmbrefEnum,0.);
			break;
		case SMBcomponentsEnum:
			iomodel->FetchDataToInput(inputs,elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.runoff",SmbRunoffEnum,0.);
			break;
		case SMBmeltcomponentsEnum:
			iomodel->FetchDataToInput(inputs,elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.melt",SmbMeltEnum,0.);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.refreeze",SmbRefreezeEnum,0.);
			break;
		case SMBgradientscomponentsEnum:
			/* Nothing to add to input */
			break;
		case SMBsemicEnum:
			iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.s0gcm",SmbS0gcmEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dailysnowfall",SmbDailysnowfallEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyrainfall",SmbDailyrainfallEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dailydsradiation",SmbDailydsradiationEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dailydlradiation",SmbDailydlradiationEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dailywindspeed",SmbDailywindspeedEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dailypressure",SmbDailypressureEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairdensity",SmbDailyairdensityEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairhumidity",SmbDailyairhumidityEnum);
			iomodel->FetchDataToInput(inputs,elements,"md.smb.dailytemperature",SmbDailytemperatureEnum);
			break;
		default:
			_error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
	}

}/*}}}*/
void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/

	int     numoutputs;
	char**  requestedoutputs = NULL;
	bool    isdelta18o,ismungsm,isd18opd,issetpddfac,interp,cycle,isfirnwarming;
	int     smb_model, smbslices, averaging;
	IssmDouble *temp = NULL;
	int         N,M;

	parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));

	iomodel->FindConstant(&smb_model,"md.smb.model");
	iomodel->FindConstant(&interp,"md.timestepping.interp_forcing");
	iomodel->FindConstant(&cycle,"md.timestepping.cycle_forcing");

	iomodel->FindConstant(&smbslices,"md.smb.steps_per_step");
	parameters->AddObject(new IntParam(SmbStepsPerStepEnum,smbslices));

	parameters->AddObject(iomodel->CopyConstantObject("md.smb.averaging",SmbAveragingEnum));

	switch(smb_model){
		case SMBforcingEnum:
			break;
		case SMBgembEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.eIdx",SmbEIdxEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.swIdx",SmbSwIdxEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.denIdx",SmbDenIdxEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.dsnowIdx",SmbDsnowIdxEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.cldFrac",SmbCldFracEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0wet",SmbT0wetEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.t0dry",SmbT0dryEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.K",SmbKEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.aSnow",SmbASnowEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIce",SmbAIceEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.dt",SmbDtEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isgraingrowth",SmbIsgraingrowthEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isalbedo",SmbIsalbedoEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isshortwave",SmbIsshortwaveEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isthermal",SmbIsthermalEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isaccumulation",SmbIsaccumulationEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismelt",SmbIsmeltEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdensification",SmbIsdensificationEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isturbulentflux",SmbIsturbulentfluxEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdeltaLWup",SmbIsdeltaLWupEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.teThresh",SmbTeThreshEnum));
			break;
		case SMBpddEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdelta18o",SmbIsdelta18oEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
			iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
			iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");

			if(ismungsm){
			  iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
			  parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,cycle,M));
			  iomodel->DeleteData(temp,"md.smb.Pfac");

			  iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
			  parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,cycle,M));
			  iomodel->DeleteData(temp,"md.smb.Tdiff");

			  iomodel->FetchData(&temp,&N,&M,"md.smb.sealev"); _assert_(N==2);
			  parameters->AddObject(new TransientParam(SmbSealevEnum,&temp[0],&temp[M],interp,cycle,M));
			  iomodel->DeleteData(temp,"md.smb.sealev");
			}
			if(isdelta18o){
				iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
				parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,cycle,M));
				iomodel->DeleteData(temp,"md.smb.delta18o");

				iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o_surface"); _assert_(N==2);
				parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,cycle,M));
				iomodel->DeleteData(temp,"md.smb.delta18o_surface");
			}

			break;
		case SMBpddSicopolisEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isfirnwarming",SmbIsfirnwarmingEnum));
			break;
		case SMBd18opddEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isd18opd",SmbIsd18opdEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.issetpddfac",SmbIssetpddfacEnum));
			iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
			iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
			iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
			if(isd18opd){
				parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.smb.istemperaturescaled",SmbIstemperaturescaledEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipscaled",SmbIsprecipscaledEnum));
				iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
				parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,cycle,M));
				iomodel->DeleteData(temp,"md.smb.delta18o");

				IssmDouble yts;
				bool istemperaturescaled,isprecipscaled;
				iomodel->FindConstant(&yts,"md.constants.yts");
				iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
				iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
				if(!istemperaturescaled){
					/*Fetch array*/
					IssmDouble* doublearray = NULL;
					int         M,N;
					iomodel->FetchData(&doublearray,&M,&N,"md.smb.temperatures_reconstructed");
					if(M!=iomodel->numberofvertices+1) _error_("md.smb.temperatures_reconstructed should have nbv+1 rows");
					if(N%12!=0) _error_("md.smb.temperatures_reconstructed should have a multiple of 12 columns (since it is monthly)");
					int numyears = N/12; _assert_(numyears*12==N);

					/*Check times*/
					#ifdef _ISSM_DEBUG_
					for(int i=0;i<numyears;i++){
						for(int j=1;j<12;j++){
							//_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
							_assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
						}
					}
					#endif

					/*Build time*/
					IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
					parameters->AddObject(new DoubleVecParam(SmbTemperaturesReconstructedYearsEnum,times,numyears));
					xDelete<IssmDouble>(times);
					iomodel->DeleteData(doublearray,"md.smb.temperatures_reconstructed");
				}
				if(!isprecipscaled){
					/*Fetch array*/
					IssmDouble* doublearray = NULL;
					int         M,N;
					iomodel->FetchData(&doublearray,&M,&N,"md.smb.precipitations_reconstructed");
					if(M!=iomodel->numberofvertices+1) _error_("md.smb.precipitations_reconstructed should have nbv+1 rows");
					if(N%12!=0) _error_("md.smb.precipitations_reconstructed should have a multiple of 12 columns (since it is monthly)");
					int numyears = N/12; _assert_(numyears*12==N);

					/*Check times*/
					#ifdef _ISSM_DEBUG_
					for(int i=0;i<numyears;i++){
						for(int j=1;j<12;j++){
							//_assert_(floor(doublearray[(M-1)*N+i*12+j]/yts)==floor(doublearray[(M-1)*N+i*12]/yts));
							_assert_(doublearray[(M-1)*N+i*12+j]>doublearray[(M-1)*N+i*12+j-1]);
						}
					}
					#endif

					/*Build time*/
					IssmDouble* times = xNew<IssmDouble>(numyears); for(int i=0;i<numyears;i++) times[i] = doublearray[(M-1)*N+i*12];
					parameters->AddObject(new DoubleVecParam(SmbPrecipitationsReconstructedYearsEnum,times,numyears));
					xDelete<IssmDouble>(times);
					iomodel->DeleteData(doublearray,"md.smb.precipitations_reconstructed");
				}
			}

			break;
		case SMBgradientsEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBgradientselaEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBautoregressionEnum:
         /*Nothing to add to parameters*/
         break;
		case SMBhenningEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBcomponentsEnum:
			break;
		case SMBmeltcomponentsEnum:
			break;
		case SMBgradientscomponentsEnum:
				parameters->AddObject(iomodel->CopyConstantObject("md.smb.accualti",SmbAccualtiEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.smb.runoffalti",SmbRunoffaltiEnum));

				iomodel->FetchData(&temp,&N,&M,"md.smb.accugrad"); _assert_(N==2);
				parameters->AddObject(new TransientParam(SmbAccugradEnum,&temp[0],&temp[M],interp,cycle,M));
				iomodel->DeleteData(temp,"md.smb.accugrad");
				iomodel->FetchData(&temp,&N,&M,"md.smb.runoffgrad"); _assert_(N==2);
				parameters->AddObject(new TransientParam(SmbRunoffgradEnum,&temp[0],&temp[M],interp,cycle,M));
				iomodel->DeleteData(temp,"md.smb.runoffgrad");

				iomodel->FetchData(&temp,&N,&M,"md.smb.accuref"); _assert_(N==2);
				parameters->AddObject(new TransientParam(SmbAccurefEnum,&temp[0],&temp[M],interp,cycle,M));
				iomodel->DeleteData(temp,"md.smb.accuref");
				iomodel->FetchData(&temp,&N,&M,"md.smb.runoffref"); _assert_(N==2);
				parameters->AddObject(new TransientParam(SmbRunoffrefEnum,&temp[0],&temp[M],interp,cycle,M));
				iomodel->DeleteData(temp,"md.smb.runoffref");
			break;
		case SMBsemicEnum:
			/*Nothing to add to parameters*/
			break;
		default:
			_error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
	}

	iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.smb.requested_outputs");
	parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(&requestedoutputs,numoutputs,"md.smb.requested_outputs");

}/*}}}*/

/*Finite Element Analysis*/
void           SmbAnalysis::Core(FemModel* femmodel){/*{{{*/

	int    smb_model;

	/*Figure out smb model: */
	femmodel->parameters->FindParam(&smb_model,SmbEnum);

	/*branch to correct module*/
	switch(smb_model){
		case SMBforcingEnum:
			SmbForcingx(femmodel);
			break;
		case SMBgembEnum:
			Gembx(femmodel);
			break;
		case SMBpddEnum:
			bool isdelta18o,ismungsm;
			femmodel->parameters->FindParam(&isdelta18o,SmbIsdelta18oEnum);
			femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
			if(isdelta18o){
				if(VerboseSolution()) _printf0_("   call Delta18oParameterization module\n");
				Delta18oParameterizationx(femmodel);
			}
			if(ismungsm){
				if(VerboseSolution()) _printf0_("   call MungsmtpParameterization module\n");
				MungsmtpParameterizationx(femmodel);
			}
			if(VerboseSolution()) _printf0_("   call positive degree day module\n");
			PositiveDegreeDayx(femmodel);
			break;
		case SMBpddSicopolisEnum:
			if(VerboseSolution()) _printf0_("   call SICOPOLIS positive degree day module\n");
			PositiveDegreeDaySicopolisx(femmodel);
			break;
		case SMBd18opddEnum:
			bool isd18opd;
			femmodel->parameters->FindParam(&isd18opd,SmbIsd18opdEnum);
			if(isd18opd){
				if(VerboseSolution()) _printf0_("   call Delta18opdParameterization module\n");
				Delta18opdParameterizationx(femmodel);
				if(VerboseSolution()) _printf0_("   call positive degree day module\n");
				PositiveDegreeDayx(femmodel);
			}
			break;
		case SMBgradientsEnum:
			if(VerboseSolution())_printf0_("	call smb gradients module\n");
			SmbGradientsx(femmodel);
			break;
		case SMBgradientselaEnum:
			if(VerboseSolution())_printf0_("	call smb gradients ela module\n");
			SmbGradientsElax(femmodel);
			break;
		case SMBautoregressionEnum:
         if(VerboseSolution())_printf0_("    call smb autoregression module\n");
         Smbautoregressionx(femmodel);
         break;
		case SMBhenningEnum:
			if(VerboseSolution())_printf0_("  call smb Henning module\n");
			SmbHenningx(femmodel);
			break;
		case SMBcomponentsEnum:
			if(VerboseSolution())_printf0_("  call smb Components module\n");
			SmbComponentsx(femmodel);
			break;
		case SMBmeltcomponentsEnum:
			if(VerboseSolution())_printf0_("  call smb Melt Components module\n");
			SmbMeltComponentsx(femmodel);
			break;
		case SMBgcmEnum:
			/*Nothing to be done*/
			break;
		case SMBgradientscomponentsEnum:
			if(VerboseSolution())_printf0_("	call smb gradients components module\n");
			SmbGradientsComponentsx(femmodel);
			break;
		case SMBsemicEnum:
			#ifdef _HAVE_SEMIC_
			if(VerboseSolution())_printf0_("  call smb SEMIC module\n");
			SmbSemicx(femmodel);
			#else
			_error_("SEMIC not installed");
			#endif //_HAVE_SEMIC_
			break;
		default:
			_error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
	}

}/*}}}*/
void           SmbAnalysis::PreCore(FemModel* femmodel){/*{{{*/
	_error_("not implemented");
}/*}}}*/
ElementVector* SmbAnalysis::CreateDVector(Element* element){/*{{{*/
	_error_("not implemented");
}/*}}}*/
ElementMatrix* SmbAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
_error_("Not implemented");
}/*}}}*/
ElementMatrix* SmbAnalysis::CreateKMatrix(Element* element){/*{{{*/
	_error_("not implemented yet");
}/*}}}*/
ElementVector* SmbAnalysis::CreatePVector(Element* element){/*{{{*/
_error_("not implemented yet");
}/*}}}*/
void           SmbAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
	   _error_("not implemented yet");
}/*}}}*/
void           SmbAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
	_error_("Not implemented yet");
}/*}}}*/
void           SmbAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
	_error_("not implemented yet");
}/*}}}*/
void           SmbAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
	/*Default, do nothing*/
	return;
}/*}}}*/
