/* file:  PddSurfaceMassBlanceSicopolis.cpp
   Calculating the surface mass balance using the adapted PDD routine from SICOPOLIS.
 */

#include "./elements.h"
#include "../Numerics/numerics.h"
#include "../Exceptions/exceptions.h"
#include <cmath>

IssmDouble PddSurfaceMassBalanceSicopolis(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec,
				 IssmDouble* melt, IssmDouble* accu, IssmDouble* melt_star, IssmDouble* t_ampl, IssmDouble* p_ampl,
				 IssmDouble yts, IssmDouble s, IssmDouble desfac,
				 IssmDouble s0t, IssmDouble s0p, IssmDouble rlaps,
				 IssmDouble rho_water,IssmDouble rho_ice){

  int			imonth;				// month counter
  IssmDouble B;					// output: surface mass balance (m/a IE), melt+accumulation
  IssmDouble frac_solid, snowfall, rainfall, runoff; 
  IssmDouble saccu;				// yearly surface accumulation (m/a IE)
  IssmDouble smelt;				// yearly melt (m/a IE)
  IssmDouble smelt_star;		// yearly ...
  IssmDouble precip;				// total precipitation during 1 year
  IssmDouble sconv;				//rhow_rain/rhoi / 12 months
  IssmDouble st;					// elevation between altitude of the temp record and current altitude
  IssmDouble sp;					// elevation between altitude of the prec record and current altitude
  IssmDouble q;					// q is desert/elev. fact
  IssmDouble pdd;					// pdd factor (a * degC)
  IssmDouble tstar;				// monthly temp. after lapse rate correction (degC)
  IssmDouble precip_star;		// monthly precip after correction (m/a IE)
  IssmDouble beta1 = 2.73;		// 3 mm IE/(d*deg C),  ablation factor for snow per positive degree day.
  IssmDouble beta2 = 7.28;		// 8 mm IE/(d*deg C),  ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
  IssmDouble Pmax = 0.6;
  IssmDouble inv_twelve=1./12.; 

  sconv=(rho_water/rho_ice);		//rhow_rain/rhoi

  /* FIXME betas shoud be user input */
  beta1=beta1*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C)
  beta2=beta2*(0.001*365)*sconv; // (mm WE)/(d*deg C) --> (m IE)/(a*deg C)

  /* initalize fields */
  precip=0.0;
  tstar=0.0;
  snowfall=0.0;
  pdd=0.0;
  /* seasonal loop */
  for(imonth=0;imonth<12;imonth++){

    /********* Surface temperature correction *******/    
    st=(s-s0t)/1000.;

    // FIXME rlaps ??
	 rlaps=-6.309e-03+(-5.426e-03-(-6.309e-03))*sin((imonth+1-4)*PI/6.0)*1000.0;
    monthlytemperatures[imonth]=monthlytemperatures[imonth]-rlaps*st;//*max(st,1e-3);
    tstar=monthlytemperatures[imonth]+t_ampl[0];

    /********* Precipitation correction *************/
    /* Ref: Vizcaino et al 2010; DOI 10.1007/s00382-009-0591-y */
    if(s0p<2000.0)
      q=exp(desfac*(max(s,2000.0)-2000.0));
	 else
      q=exp(desfac*(max(s,2000.0)-s0p));

    precip_star=q*monthlyprec[imonth]*sconv*p_ampl[0]*yts; // convert precip from m/s -> m/a
    precip=precip+precip_star*inv_twelve;

    /********* compute PDD **************************/
    /* Ref: Calov & Greve 2005 Journal of Glaciology, Vol. 51, No. 172, 2005, Correspondence */
	 IssmDouble s_stat=5.0;
    IssmDouble inv_sqrt2pi =1.0/sqrt(2.0*PI);
    IssmDouble inv_s_stat  =1.0/s_stat;
    IssmDouble inv_sqrt2   =1.0/sqrt(2.0);

	 #if !defined(_HAVE_ADOLC_)
    pdd=pdd+(s_stat*inv_sqrt2pi*exp(-0.5*pow(tstar*inv_s_stat,2))
				 +0.5*tstar*erfc(-tstar*inv_s_stat*inv_sqrt2))*inv_twelve;
	 #else
	 _error_("Cannot differentiate erfc, talk to ADOLC folks (http://functions.wolfram.com/GammaBetaErf/Erfc/20/01/)");
	 #endif

	 /*Partition of precip in solid and liquid parts, Bales et al. (2009) */
	 IssmDouble temp_rain=7.2;		// Threshold monthly mean temperature for
											// precipitation = 101% rain, in deg C
	 IssmDouble temp_snow=-11.6;  // Threshold monthly mean temperature for
											// precipitation = 100% snow, in deg C

	 IssmDouble coeff1=5.4714e-01;	// Coefficients
	 IssmDouble coeff2=-9.1603e-02;	// of
	 IssmDouble coeff3=-3.314e-03;	// the
	 IssmDouble coeff4= 4.66e-04;		// fifth-order
	 IssmDouble coeff5=3.8e-05;		// polynomial
	 IssmDouble coeff6=6.0e-07;		// fit

	 if(tstar>=temp_rain)
	  frac_solid = 0.0;
	 else if(tstar<=temp_snow)
	  frac_solid = 1.0;
	 else{ 
		 frac_solid=coeff1+tstar*(coeff2
					 +tstar*(coeff3+tstar*(coeff4+tstar*(coeff5+tstar*coeff6))));
	 }

	 snowfall=snowfall+precip_star*frac_solid*inv_twelve;
  } 
  /* end of seasonal loop */ 

  rainfall=precip-snowfall;
  if(snowfall<0.0) snowfall=0.0;   // correction of
  if(rainfall<0.0) rainfall=0.0;   // negative values

  if(rainfall<=(Pmax*snowfall)){
	  if((rainfall+beta1*pdd)<=(Pmax*snowfall)) {
		  smelt_star = rainfall+beta1*pdd;
		  smelt      = 0.0;
		  runoff     = smelt;
	  }
	  else{
		  smelt_star = Pmax*snowfall;
		  smelt      = beta2*(pdd-(smelt_star-rainfall)/beta1);
		  runoff     = smelt;
	  }
  } 
  else{
	  smelt_star = Pmax*snowfall;
	  smelt      = beta2*pdd;
	  runoff     = smelt+rainfall-Pmax*snowfall;
  }

  saccu = precip;	

  /* asign output*/
  melt[0]=runoff/yts;
  accu[0]=saccu/yts;
  melt_star[0]=smelt_star/yts;
  B=(saccu-runoff)/yts;

  return B;
}