#include "./SmbAnalysis.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.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){/*{{{*/
	::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;
	
	/*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->Constant(&smb_model,SmbEnum);

			
	switch(smb_model){
		case SMBforcingEnum:
			iomodel->FetchDataToInput(elements,SmbMassBalanceEnum,0.);
			break;
		case SMBgembEnum:
			iomodel->FetchDataToInput(elements,SmbTaEnum);
			iomodel->FetchDataToInput(elements,SmbVEnum);
			iomodel->FetchDataToInput(elements,SmbDswrfEnum);
			iomodel->FetchDataToInput(elements,SmbDlwrfEnum);
			iomodel->FetchDataToInput(elements,SmbPEnum);
			iomodel->FetchDataToInput(elements,SmbEAirEnum);
			iomodel->FetchDataToInput(elements,SmbPAirEnum);
			iomodel->FetchDataToInput(elements,SmbZTopEnum);
			iomodel->FetchDataToInput(elements,SmbDzTopEnum);
			iomodel->FetchDataToInput(elements,SmbDzMinEnum);
			iomodel->FetchDataToInput(elements,SmbZYEnum);
			iomodel->FetchDataToInput(elements,SmbZMaxEnum);
			iomodel->FetchDataToInput(elements,SmbZMinEnum);
			iomodel->FetchDataToInput(elements,SmbTmeanEnum);
			iomodel->FetchDataToInput(elements,SmbCEnum);
			iomodel->FetchDataToInput(elements,SmbTzEnum);
			iomodel->FetchDataToInput(elements,SmbVzEnum);
			break;
		case SMBpddEnum:
			iomodel->Constant(&isdelta18o,SmbIsdelta18oEnum);
			iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
			iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
			iomodel->FetchDataToInput(elements,SmbS0pEnum);
			iomodel->FetchDataToInput(elements,SmbS0tEnum);
			if(isdelta18o || ismungsm){
				iomodel->FetchDataToInput(elements,SmbTemperaturesLgmEnum);
				iomodel->FetchDataToInput(elements,SmbTemperaturesPresentdayEnum);
				iomodel->FetchDataToInput(elements,SmbPrecipitationsPresentdayEnum);
				iomodel->FetchDataToInput(elements,SmbPrecipitationsLgmEnum);
			}
			else{
				iomodel->FetchDataToInput(elements,SmbPrecipitationEnum);
				iomodel->FetchDataToInput(elements,SmbMonthlytemperaturesEnum);
			}
			break;
		case SMBd18opddEnum:
			iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
			iomodel->Constant(&isd18opd,SmbIsd18opdEnum);
			iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
			iomodel->FetchDataToInput(elements,SmbS0pEnum);
			iomodel->FetchDataToInput(elements,SmbS0tEnum);
			if (isd18opd){
			        iomodel->FetchDataToInput(elements,SmbTemperaturesPresentdayEnum);
			        iomodel->FetchDataToInput(elements,SmbPrecipitationsPresentdayEnum);
			}

			break;
		case SMBgradientsEnum:
			iomodel->FetchDataToInput(elements,SmbHrefEnum);
			iomodel->FetchDataToInput(elements,SmbSmbrefEnum);
			iomodel->FetchDataToInput(elements,SmbBPosEnum);
			iomodel->FetchDataToInput(elements,SmbBNegEnum);
			break;
		case SMBhenningEnum:
			iomodel->FetchDataToInput(elements,SmbSmbrefEnum,0.);
			break;
		case SMBcomponentsEnum:
			iomodel->FetchDataToInput(elements,SmbAccumulationEnum,0.);
			iomodel->FetchDataToInput(elements,SmbEvaporationEnum,0.);
			iomodel->FetchDataToInput(elements,SmbRunoffEnum,0.);
			break;
		case SMBmeltcomponentsEnum:
			iomodel->FetchDataToInput(elements,SmbAccumulationEnum,0.);
			iomodel->FetchDataToInput(elements,SmbEvaporationEnum,0.);
			iomodel->FetchDataToInput(elements,SmbMeltEnum,0.);
			iomodel->FetchDataToInput(elements,SmbRefreezeEnum,0.);
			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,interp;
	int     smb_model;
	IssmDouble *temp = NULL;
	int         N,M;
	
	parameters->AddObject(iomodel->CopyConstantObject(SmbEnum));
	
	iomodel->Constant(&smb_model,SmbEnum);
	iomodel->Constant(&interp,TimesteppingInterpForcingsEnum);
	
	switch(smb_model){
		case SMBforcingEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBgembEnum:
			parameters->AddObject(iomodel->CopyConstantObject(SmbAIdxEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbSwIdxEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbDenIdxEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbOutputFreqEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbCldFracEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbT0wetEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbT0dryEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbKEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbASnowEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbAIceEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbDtEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsgraingrowthEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsalbedoEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsshortwaveEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsthermalEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsaccumulationEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsmeltEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsdensificationEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsturbulentfluxEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbInitDensityScalingEnum));
			break;
		case SMBpddEnum:
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsdelta18oEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsmungsmEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbDesfacEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbRlapsEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbRlapslgmEnum));
			iomodel->Constant(&isdelta18o,SmbIsdelta18oEnum);
			iomodel->Constant(&ismungsm,SmbIsmungsmEnum);

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

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

				iomodel->FetchData(&temp,&N,&M,SmbDelta18oSurfaceEnum); _assert_(N==2);
				parameters->AddObject(new TransientParam(SmbDelta18oSurfaceEnum,&temp[0],&temp[M],interp,M));
				iomodel->DeleteData(temp,SmbDelta18oSurfaceEnum);
			}
			break;
		case SMBd18opddEnum:
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsmungsmEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbIsd18opdEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbDesfacEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbRlapsEnum));
			parameters->AddObject(iomodel->CopyConstantObject(SmbRlapslgmEnum));
			iomodel->Constant(&ismungsm,SmbIsmungsmEnum);
			iomodel->Constant(&isd18opd,SmbIsd18opdEnum);
			if(isd18opd){
				iomodel->FetchData(&temp,&N,&M,SmbDelta18oEnum); _assert_(N==2);
				parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
				iomodel->DeleteData(temp,SmbDelta18oEnum);
				
				parameters->AddObject(iomodel->CopyConstantObject(SmbDpermilEnum));
			}
			break;
		case SMBgradientsEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBhenningEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBcomponentsEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBmeltcomponentsEnum:
			/*Nothing to add to parameters*/
			break;
		default:
			_error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
	}

	iomodel->FetchData(&requestedoutputs,&numoutputs,SmbRequestedOutputsEnum);
	parameters->AddObject(new IntParam(SmbNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new StringArrayParam(SmbRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(&requestedoutputs,numoutputs,SmbRequestedOutputsEnum);

}/*}}}*/

/*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:
			/*Nothing to be done*/
			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 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 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;
		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;
}/*}}}*/