#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,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/

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

	/*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(i,iomodel,analysis_counter,analysis_type,P1Enum);
			counter++;
		}
	}

	/*Figure out smb model: */
	iomodel->FindConstant(&smb_model,"md.smb.model");
	switch(smb_model){
		case SMBforcingEnum:
			iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
			break;
		case SMBgembEnum:
			iomodel->FetchDataToInput(elements,"md.smb.Ta",SmbTaEnum);
			iomodel->FetchDataToInput(elements,"md.smb.V",SmbVEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dswrf",SmbDswrfEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dlwrf",SmbDlwrfEnum);
			iomodel->FetchDataToInput(elements,"md.smb.P",SmbPEnum);
			iomodel->FetchDataToInput(elements,"md.smb.eAir",SmbEAirEnum);
			iomodel->FetchDataToInput(elements,"md.smb.pAir",SmbPAirEnum);
			iomodel->FetchDataToInput(elements,"md.smb.zTop",SmbZTopEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dzTop",SmbDzTopEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dzMin",SmbDzMinEnum);
			iomodel->FetchDataToInput(elements,"md.smb.zY",SmbZYEnum);
			iomodel->FetchDataToInput(elements,"md.smb.zMax",SmbZMaxEnum);
			iomodel->FetchDataToInput(elements,"md.smb.zMin",SmbZMinEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Tmean",SmbTmeanEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Vmean",SmbVmeanEnum);
			iomodel->FetchDataToInput(elements,"md.smb.C",SmbCEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Tz",SmbTzEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Vz",SmbVzEnum);
			InputUpdateFromConstantx(elements,0.,SmbIsInitializedEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Dzini",SmbDziniEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Dini",SmbDiniEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Reini",SmbReiniEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Gdnini",SmbGdniniEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Gspini",SmbGspiniEnum);
			iomodel->FetchDataToInput(elements,"md.smb.ECini",SmbECiniEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Wini",SmbWiniEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Aini",SmbAiniEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Tini",SmbTiniEnum);
			iomodel->FetchDataToInput(elements,"md.smb.Sizeini",SmbSizeiniEnum);
			iomodel->FetchDataToInput(elements,"md.smb.aValue",SmbAValueEnum);
			iomodel->FetchDataToInput(elements,"md.smb.teValue",SmbTeValueEnum);
			break;
		case SMBpddEnum:
			iomodel->FindConstant(&isdelta18o,"md.smb.isdelta18o");
			iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
			iomodel->FetchDataToInput(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
			iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
			iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum);
			if(isdelta18o || ismungsm){
				iomodel->FetchDataToInput(elements,"md.smb.temperatures_lgm",SmbTemperaturesLgmEnum);
				iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
				iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
				iomodel->FetchDataToInput(elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum);
			}else{
				iomodel->FetchDataToInput(elements,"md.smb.precipitation",SmbPrecipitationEnum);
				iomodel->FetchDataToInput(elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
			}
			break;
		case SMBpddSicopolisEnum:
			iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
			iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum);
			iomodel->FindConstant(&isfirnwarming,"md.smb.isfirnwarming");
			iomodel->FetchDataToInput(elements,"md.smb.smb_corr",SmbSmbCorrEnum);
			iomodel->FetchDataToInput(elements,"md.smb.precipitation",SmbPrecipitationEnum);
			iomodel->FetchDataToInput(elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
			iomodel->FetchDataToInput(elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum);
			iomodel->FetchDataToInput(elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum);
			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(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
			iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
			iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum);
			if(isd18opd){
				iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
				iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
				if(!istemperaturescaled){
					iomodel->FetchDataToInput(elements,"md.smb.temperatures_reconstructed",SmbTemperaturesReconstructedEnum);
				}
				if(!isprecipscaled){
					iomodel->FetchDataToInput(elements,"md.smb.precipitations_reconstructed",SmbPrecipitationsReconstructedEnum);
				}
			}
			if(issetpddfac){
				iomodel->FetchDataToInput(elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
				iomodel->FetchDataToInput(elements,"md.smb.pddfac_ice",SmbPddfacIceEnum,-1.);
			}
			break;
		case SMBgradientsEnum:
			iomodel->FetchDataToInput(elements,"md.smb.href",SmbHrefEnum);
			iomodel->FetchDataToInput(elements,"md.smb.smbref",SmbSmbrefEnum);
			iomodel->FetchDataToInput(elements,"md.smb.b_pos",SmbBPosEnum);
			iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum);
			break;
		case SMBgradientselaEnum:
			iomodel->FetchDataToInput(elements,"md.smb.ela",SmbElaEnum);
			iomodel->FetchDataToInput(elements,"md.smb.b_pos",SmbBPosEnum);
			iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum);
			iomodel->FetchDataToInput(elements,"md.smb.b_max",SmbBMaxEnum);
			iomodel->FetchDataToInput(elements,"md.smb.b_min",SmbBMinEnum);
			break;
		case SMBhenningEnum:
			iomodel->FetchDataToInput(elements,"md.smb.smbref",SmbSmbrefEnum,0.);
			break;
		case SMBcomponentsEnum:
			iomodel->FetchDataToInput(elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
			iomodel->FetchDataToInput(elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
			iomodel->FetchDataToInput(elements,"md.smb.runoff",SmbRunoffEnum,0.);
			break;
		case SMBmeltcomponentsEnum:
			iomodel->FetchDataToInput(elements,"md.smb.accumulation",SmbAccumulationEnum,0.);
			iomodel->FetchDataToInput(elements,"md.smb.evaporation",SmbEvaporationEnum,0.);
			iomodel->FetchDataToInput(elements,"md.smb.melt",SmbMeltEnum,0.);
			iomodel->FetchDataToInput(elements,"md.smb.refreeze",SmbRefreezeEnum,0.);
			break;
		case SMBgradientscomponentsEnum:
			iomodel->FetchDataToInput(elements,"md.smb.accualti",SmbAccualtiEnum);
			iomodel->FetchDataToInput(elements,"md.smb.accugrad",SmbAccugradEnum);
			iomodel->FetchDataToInput(elements,"md.smb.runoffalti",SmbRunoffaltiEnum);
			iomodel->FetchDataToInput(elements,"md.smb.runoffgrad",SmbRunoffgradEnum);
			break;
		case SMBsemicEnum:
			iomodel->FetchDataToInput(elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
			iomodel->FetchDataToInput(elements,"md.smb.s0gcm",SmbS0gcmEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dailysnowfall",SmbDailysnowfallEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dailyrainfall",SmbDailyrainfallEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dailydsradiation",SmbDailydsradiationEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dailydlradiation",SmbDailydlradiationEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dailywindspeed",SmbDailywindspeedEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dailypressure",SmbDailypressureEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dailyairdensity",SmbDailyairdensityEnum);
			iomodel->FetchDataToInput(elements,"md.smb.dailyairhumidity",SmbDailyairhumidityEnum);
			iomodel->FetchDataToInput(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){/*{{{*/

	printf("start param update\n");
	int     numoutputs;
	char**  requestedoutputs = NULL;
	bool    isdelta18o,ismungsm,isd18opd,issetpddfac,interp,isfirnwarming;
	int     smb_model, smbslices;
	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_forcings");

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

	switch(smb_model){
		case SMBforcingEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
			break;
		case SMBgembEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.aIdx",SmbAIdxEnum));
			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.isclimatology",SmbIsclimatologyEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
			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,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,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,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,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,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,M));
				iomodel->DeleteData(temp,"md.smb.delta18o");
			}
			break;
		case SMBgradientsEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBgradientselaEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBhenningEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBcomponentsEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
			break;
		case SMBmeltcomponentsEnum:
			parameters->AddObject(iomodel->CopyConstantObject("md.smb.isclimatology",SmbIsclimatologyEnum));
			break;
		case SMBgradientscomponentsEnum:
				parameters->AddObject(iomodel->CopyConstantObject("md.smb.accualti",SmbAccualtiEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.smb.accugrad",SmbAccugradEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.smb.runoffalti",SmbRunoffaltiEnum));
				parameters->AddObject(iomodel->CopyConstantObject("md.smb.runoffgrad",SmbRunoffgradEnum));

			  iomodel->FetchData(&temp,&N,&M,"md.smb.accuref"); _assert_(N==2);
			  parameters->AddObject(new TransientParam(SmbAccurefEnum,&temp[0],&temp[M],interp,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,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 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");
	}

}/*}}}*/
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_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;
}/*}}}*/