/*!\file FrontalForcingsx
 * \brief: compute ice frontal melting rate
 */

#include "./FrontalForcingsx.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "../../shared/Random/random.h"

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

	/*Recover melt_parameterization*/
	int melt_parameterization;
	femmodel->parameters->FindParam(&melt_parameterization,FrontalForcingsParamEnum);

	/*Calculate melting rate*/
	switch(melt_parameterization){
		case FrontalForcingsDefaultEnum:
			break;
		case FrontalForcingsRignotAutoregressionEnum:
			Thermalforcingautoregressionx(femmodel);
			/*Do not break here, call IcefrontAreax(),RignotMeltParameterizationx()*/
		case FrontalForcingsRignotEnum:
			femmodel->IcefrontAreax();
			femmodel->RignotMeltParameterizationx();
			break;
		default:
			_error_("Frontal forcings "<<EnumToStringx(melt_parameterization)<<" not supported yet");
	}
}/*}}}*/
void ThermalforcingautoregressionInitx(FemModel* femmodel){/*{{{*/

   /*Initialization step of Thermalforcingautoregressionx*/
   int M,N,Nphi,arorder,numbasins,my_rank;
   IssmDouble starttime,tstep_ar,tinit_ar;
   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
   femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
   IssmDouble* beta0    = NULL;
   IssmDouble* beta1    = NULL;
   IssmDouble* phi      = NULL;
   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
   femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum);
   femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
   femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum);    _assert_(M==numbasins);
   femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
   femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);

   /*AR model spin-up with 0 noise to initialize ThermalforcingValuesAutoregressionEnum (688 = log(0.001)/log(0.99): decaying time of inluence of phi[0]=0.99 to 0.001 of beta_0*/
	int nspin = 688;
   for(Object* &object:femmodel->elements->objects){
      Element* element      = xDynamicCast<Element*>(object); //generate element object
      element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,FrontalForcingsRignotAutoregressionEnum);
   }

   /*Cleanup*/
   xDelete<IssmDouble>(beta0);
   xDelete<IssmDouble>(beta1);
   xDelete<IssmDouble>(phi);
}/*}}}*/
void Thermalforcingautoregressionx(FemModel* femmodel){/*{{{*/

   /*Get time parameters*/
   IssmDouble time,dt,starttime,tstep_ar;
   femmodel->parameters->FindParam(&time,TimeEnum);
   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
   femmodel->parameters->FindParam(&tstep_ar,FrontalForcingsAutoregressionTimestepEnum);

   /*Initialize module at first time step*/
   if(time<=starttime+dt){ThermalforcingautoregressionInitx(femmodel);}
   /*Determine if this is a time step for the AR model*/
   bool isstepforar = false;

   #ifndef _HAVE_AD_
   if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt) isstepforar = true;
   #else
   _error_("not implemented yet");
   #endif

   /*Load parameters*/
	bool isstochastic;
   bool istfstochastic = false;
	int M,N,Nphi,arorder,numbasins,my_rank;
   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
   femmodel->parameters->FindParam(&arorder,FrontalForcingsAutoregressiveOrderEnum);
   IssmDouble tinit_ar;
   IssmDouble* beta0      = NULL;
   IssmDouble* beta1      = NULL;
   IssmDouble* phi        = NULL;

	femmodel->parameters->FindParam(&tinit_ar,FrontalForcingsAutoregressionInitialTimeEnum);
   femmodel->parameters->FindParam(&beta0,&M,FrontalForcingsBeta0Enum);    _assert_(M==numbasins);
   femmodel->parameters->FindParam(&beta1,&M,FrontalForcingsBeta1Enum);    _assert_(M==numbasins);
   femmodel->parameters->FindParam(&phi,&M,&Nphi,FrontalForcingsPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);

	femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
	if(isstochastic){
		int  numstochasticfields;
		int* stochasticfields;
		femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
		femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
		for(int i=0;i<numstochasticfields;i++){
			if(stochasticfields[i]==FrontalForcingsRignotAutoregressionEnum) istfstochastic = true;
		}
		xDelete<int>(stochasticfields);
	}
   /*Time elapsed with respect to AR model initial time*/
   IssmDouble telapsed_ar = time-tinit_ar;

   /*Loop over each element to compute Thermal Forcing at vertices*/
   for(Object* &object:femmodel->elements->objects){
      Element* element = xDynamicCast<Element*>(object);
      element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,istfstochastic,FrontalForcingsRignotAutoregressionEnum);
   }

   /*Cleanup*/
   xDelete<IssmDouble>(beta0);
   xDelete<IssmDouble>(beta1);
   xDelete<IssmDouble>(phi);
}/*}}}*/
