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

#include "./cores.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);

		/*Using the Shreve based Model*/
		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,WatercolumnEnum,WaterColumnOldEnum);

			if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps || i==0)){
				if(VerboseSolution()) _printf0_("   saving results \n");
				int outputs[3] = {WatercolumnEnum,HydrologyWaterVxEnum,HydrologyWaterVyEnum};
				femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);

				/*unload results*/
				if(VerboseSolution()) _printf0_("   saving temporary results\n");
				OutputResultsx(femmodel);
			}
		}

		/*Using the double continuum model*/
		else if (hydrology_model==HydrologydcEnum){
			InputDuplicatex(femmodel,SedimentHeadEnum,SedimentHeadOldEnum);
			femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
			if (isefficientlayer){
				InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum);
				InputDuplicatex(femmodel,HydrologydcEplThicknessEnum,HydrologydcEplThicknessOldEnum);
			}
			
			/*Proceed now to heads computations*/
			if(VerboseSolution()) _printf0_("   computing water head\n");
			solutionsequence_hydro_nonlinear(femmodel);
			if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps || i==0)){
				if(VerboseSolution()) _printf0_("   saving results \n");
				if(isefficientlayer){
					int outputs[8] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveNodeEnum,HydrologydcMaskEplactiveEltEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum};
					femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],8);
				}
				else{
					int outputs[2] = {SedimentHeadEnum,SedimentHeadResidualEnum};
					femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
				}
				/*unload results*/
				if(VerboseSolution()) _printf0_("   saving temporary results\n");
				OutputResultsx(femmodel);
			}
		}

	}
}
