/*!\file FloatingiceMeltingRatex
 * \brief: calculates Floating ice melting rate
 */

#include "./FloatingiceMeltingRatex.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "./../../classes/Inputs2/DatasetInput2.h"

void FloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/

	/*Intermediaties*/
	int  basalforcing_model;

	/*First, get BMB model from parameters*/
	femmodel->parameters->FindParam(&basalforcing_model,BasalforcingsEnum);

	/*branch to correct module*/
	switch(basalforcing_model){
		case FloatingMeltRateEnum:
		case MantlePlumeGeothermalFluxEnum:
			/*Nothing to be done*/
			break;
		case LinearFloatingMeltRateEnum:
			if(VerboseSolution())_printf0_("	  call Linear Floating melting rate module\n");
			LinearFloatingiceMeltingRatex(femmodel);
			break;
		case MismipFloatingMeltRateEnum:
			if(VerboseSolution())_printf0_("	  call Mismip Floating melting rate module\n");
			MismipFloatingiceMeltingRatex(femmodel);
			break;
		case SpatialLinearFloatingMeltRateEnum:
			if(VerboseSolution())_printf0_("	  call Spatial Linear Floating melting rate module\n");
			SpatialLinearFloatingiceMeltingRatex(femmodel);
			break;
		case BasalforcingsPicoEnum:
			if(VerboseSolution())_printf0_("   call Pico Floating melting rate module\n");
			FloatingiceMeltingRatePicox(femmodel);
			break;
		case BasalforcingsIsmip6Enum:
			if(VerboseSolution())_printf0_("   call ISMIP 6 Floating melting rate module\n");
			FloatingiceMeltingRateIsmip6x(femmodel);
			break;
		case BeckmannGoosseFloatingMeltRateEnum:
			if(VerboseSolution())_printf0_("        call BeckmannGoosse Floating melting rate module\n");
			BeckmannGoosseFloatingiceMeltingRatex(femmodel);
			break;
		default:
			_error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");
	}

}/*}}}*/

void LinearFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/

	for(int i=0;i<femmodel->elements->Size();i++){
		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
		element->LinearFloatingiceMeltingRate();
	}

}/*}}}*/
void SpatialLinearFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/

	for(int i=0;i<femmodel->elements->Size();i++){
		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
		element->SpatialLinearFloatingiceMeltingRate();
	}

}/*}}}*/
void MismipFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/

	for(int i=0;i<femmodel->elements->Size();i++){
		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
		element->MismipFloatingiceMeltingRate();
	}
}
/*}}}*/
void FloatingiceMeltingRateIsmip6x(FemModel* femmodel){/*{{{*/

	int         num_basins, basinid, num_depths, domaintype;
	IssmDouble  area, tf, base, time;
	bool        islocal;
	IssmDouble* tf_depths = NULL;

	femmodel->parameters->FindParam(&num_basins,BasalforcingsIsmip6NumBasinsEnum);
	femmodel->parameters->FindParam(&tf_depths,&num_depths,BasalforcingsIsmip6TfDepthsEnum); _assert_(tf_depths); 
	femmodel->parameters->FindParam(&islocal,BasalforcingsIsmip6IsLocalEnum);

	/*Binary search works for vectors that are sorted in increasing order only, make depths positive*/
	for(int i=0;i<num_depths;i++) tf_depths[i] = -tf_depths[i];

	IssmDouble* tf_weighted_avg     = xNewZeroInit<IssmDouble>(num_basins);
	IssmDouble* tf_weighted_avg_cpu = xNewZeroInit<IssmDouble>(num_basins);
	IssmDouble* areas_summed        = xNewZeroInit<IssmDouble>(num_basins);
	IssmDouble* areas_summed_cpu    = xNewZeroInit<IssmDouble>(num_basins);

	/*Get TF at each ice shelf point - linearly intepolate in depth and time*/
	for(int i=0;i<femmodel->elements->Size();i++){
		Element* element     = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
		int      numvertices = element->GetNumberOfVertices();

		/*Set melt to 0 if non floating*/
		if(!element->IsIceInElement() || !element->IsFloating() || !element->IsOnBase()){
			IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
			element->AddInput2(BasalforcingsFloatingiceMeltingRateEnum,values,P1DGEnum);
			element->AddInput2(BasalforcingsIsmip6TfShelfEnum,values,P1DGEnum);
			xDelete<IssmDouble>(values);
			continue;
		}

		/*Get TF on all vertices*/
		IssmDouble*    tf_test        = xNew<IssmDouble>(numvertices);
		IssmDouble*    depth_vertices = xNew<IssmDouble>(numvertices);
		DatasetInput2* tf_input = element->GetDatasetInput2(BasalforcingsIsmip6TfEnum); _assert_(tf_input);

		element->GetInputListOnVertices(&depth_vertices[0],BaseEnum);

		Gauss* gauss=element->NewGauss();
		for(int iv=0;iv<numvertices;iv++){
			gauss->GaussVertex(iv);

			/*Find out where the ice shelf base is within tf_depths*/
			IssmDouble depth = -depth_vertices[iv]; /*NOTE: make sure we are dealing with depth>0*/
			int offset;
			int found=binary_search(&offset,depth,tf_depths,num_depths); 
			if(!found) _error_("depth not found");

			if (offset==-1){
				/*get values for the first depth: */
				_assert_(depth<=tf_depths[0]);
				tf_input->GetInputValue(&tf_test[iv],gauss,0);
			}
			else if(offset==num_depths-1){
				/*get values for the last time: */
				_assert_(depth>=tf_depths[num_depths-1]);
				tf_input->GetInputValue(&tf_test[iv],gauss,num_depths-1);
			}
			else {
				/*get values between two times [offset:offset+1], Interpolate linearly*/
				_assert_(depth>=tf_depths[offset] && depth<tf_depths[offset+1]);
				IssmDouble deltaz=tf_depths[offset+1]-tf_depths[offset];
				IssmDouble alpha2=(depth-tf_depths[offset])/deltaz;
				IssmDouble alpha1=(1.-alpha2);
				IssmDouble tf1,tf2;
				tf_input->GetInputValue(&tf1,gauss,offset);
				tf_input->GetInputValue(&tf2,gauss,offset+1);
				tf_test[iv] = alpha1*tf1 + alpha2*tf2;
			}
		}

		element->AddInput2(BasalforcingsIsmip6TfShelfEnum,tf_test,P1DGEnum);
		xDelete<IssmDouble>(tf_test);
		xDelete<IssmDouble>(depth_vertices);
		delete gauss;
	}

	if(!islocal) {
		/*Compute sums of tf*area and shelf-area per cpu*/
		for(int i=0;i<femmodel->elements->Size();i++){
			Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
			if(!element->IsOnBase() || !element->IsIceInElement() || !element->IsFloating()) continue;
			/*Spawn basal element if on base to compute element area*/
			Element* basalelement = element->SpawnBasalElement();
			Input2* tf_input=basalelement->GetInput2(BasalforcingsIsmip6TfShelfEnum); _assert_(tf_input);
			basalelement->GetInputValue(&basinid,BasalforcingsIsmip6BasinIdEnum);
			Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
			tf_input->GetInputValue(&tf,gauss);
			delete gauss;
			area=basalelement->GetHorizontalSurfaceArea();
			tf_weighted_avg[basinid]+=tf*area;
			areas_summed[basinid]   +=area;	
			/*Delete spawned element if we are in 3D*/
			basalelement->FindParam(&domaintype,DomainTypeEnum);
			if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
		}

		/*Syncronize across cpus*/
		ISSM_MPI_Allreduce(tf_weighted_avg,tf_weighted_avg_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
		ISSM_MPI_Allreduce(areas_summed,areas_summed_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());

		/*Compute weighted means and save*/
		for(int k=0;k<num_basins;k++){tf_weighted_avg_cpu[k] = tf_weighted_avg_cpu[k]/areas_summed_cpu[k];}
		femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6AverageTfEnum,tf_weighted_avg_cpu,num_basins));
	}

   /*Compute meltrates*/
	for(int i=0;i<femmodel->elements->Size();i++){
		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
		element->Ismip6FloatingiceMeltingRate();
	}

	/*Cleanup and return */
	xDelete<IssmDouble>(tf_weighted_avg);
	xDelete<IssmDouble>(tf_weighted_avg_cpu);
	xDelete<IssmDouble>(areas_summed);
	xDelete<IssmDouble>(areas_summed_cpu);
	xDelete<IssmDouble>(tf_depths);
}
/*}}}*/
void BeckmannGoosseFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/

	for(int i=0;i<femmodel->elements->Size();i++){
		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
		element->BeckmannGoosseFloatingiceMeltingRate();
	}
}
/*}}}*/
