/*!\file: hydrology_core.cpp
 * \brief: core of the hydrology solution 
 */ 

#include "./analyses.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"

void hydrology_core(FemModel* femmodel){

	int i;

	/*intermediary*/
	int        step,nsteps;
	int        output_frequency,hydrology_model;
	bool       save_results;
	bool       modify_loads=true;
	bool       isefficientlayer;
	IssmDouble starttime,final_time;
	IssmDouble time,dt;

	/*first recover parameters common to all solutions*/
	femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
	femmodel->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
	femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum);
	femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);

	/*first compute slopes: */
	if (hydrology_model==HydrologyshreveEnum){
		surfaceslope_core(femmodel);
		bedslope_core(femmodel);
	}

	/*Compute number of time steps: */
	if((dt==0)|| (final_time==0)){
		dt=0;
		nsteps=1;
	}
	else nsteps=reCast<int,IssmDouble>((final_time-starttime)/dt);

	/*initialize: */
	step=0;
	time=starttime;

	/*Loop through time: */
	for(i=0;i<nsteps;i++){

		if(nsteps)if(VerboseSolution()) _printf0_("time step:" << i+1 << "/" << nsteps << "\n");
		time+=dt;
		step+=1;
		femmodel->parameters->SetParam(time,TimeEnum);
		femmodel->parameters->SetParam(step,StepEnum);

		if (hydrology_model==HydrologyshreveEnum){
			if(VerboseSolution()) _printf0_("   computing water column\n");
			femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum);
			solutionsequence_nonlinear(femmodel,modify_loads);

			/*transfer water column thickness to old water column thickness: */

			InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WatercolumnEnum,WaterColumnOldEnum);

			if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
				if(VerboseSolution()) _printf0_("   saving results \n");
				InputToResultx(femmodel,WatercolumnEnum);
				InputToResultx(femmodel,HydrologyWaterVxEnum);
				InputToResultx(femmodel,HydrologyWaterVyEnum);

				/*unload results*/
				if(VerboseSolution()) _printf0_("   saving temporary results\n");
				OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results);
			}
		}

		else if (hydrology_model==HydrologydcEnum){
			InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SedimentHeadEnum,SedimentHeadOldEnum);
			femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
			if (isefficientlayer){
				InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EplHeadEnum,EplHeadOldEnum);
			}

			if(VerboseSolution()) _printf0_("   computing water head\n");
			solutionsequence_hydro_nonlinear(femmodel);
			if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
				if(VerboseSolution()) _printf0_("   saving results \n");
				InputToResultx(femmodel,SedimentHeadEnum);
				InputToResultx(femmodel,SedimentHeadResidualEnum);
				if(isefficientlayer){
					InputToResultx(femmodel,EplHeadEnum);
					InputToResultx(femmodel,HydrologydcMaskEplactiveEnum);
				}
				/*unload results*/
				if(VerboseSolution()) _printf0_("   saving temporary results\n");
				OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results);
			}
		}

	}
}
