/*!\file GEMB module from Alex Gardner.
 * \brief: calculates SMB 
 */

#include "./SurfaceMassBalancex.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "../modules.h"
#include "../../classes/Inputs/TransientInput.h"

const double Pi = 3.141592653589793;
const double CtoK = 273.15;             // Kelvin to Celcius conversion/ice melt. point T in K
const double dts = 86400.0;              // Number of seconds in a day

/* Tolerances have to be defined for if loops involving some specific values
	like densitiy of ice or melting point temp because values are "rounded"
	(e.g rho ice = 909.99... instead of 910.0) and we don't always go into the right loop */
const double Ttol = 1e-10;
const double Dtol = 1e-11;
const double Gdntol = 1e-10;
const double Wtol = 1e-13;
const double Ptol = 1e-6;

const double CI = 2102.0;                       // heat capacity of snow/ice (J kg-1 k-1)
const double LF = 0.3345e6;             // latent heat of fusion (J kg-1)
const double LV = 2.495e6;               // latent heat of vaporization (J kg-1)
const double LS = 2.8295e6;             // latent heat of sublimation (J kg-1)
const double SB = 5.67e-8;                // Stefan-Boltzmann constant (W m-2 K-4)
const double CA = 1005.0;                    // heat capacity of air (J kg-1 K-1)
const double R = 8.314;                      // gas constant (J mol-1 K-1)

const double Delflag=-99999;

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

	int        count=0;
	IssmDouble time,dt,finaltime;
	IssmDouble timeclim=0.0;
	IssmDouble t,smb_dt;
   IssmDouble delta;
	bool       linear_interp=true;

	femmodel->parameters->FindParam(&linear_interp,TimesteppingInterpForcingEnum); /*is interpolation requested*/
	femmodel->parameters->FindParam(&time,TimeEnum);                        /*transient core time at which we run the smb core*/
   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
	femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
   femmodel->parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/

	//before starting loop, realize that the transient core runs this smb_core at time = time +deltaT.
	//go back to time - deltaT:
	time-=dt;

	IssmDouble timeinputs = time;

	/*Start loop: */
	count=1;
	for (t=time;t<=time+dt-smb_dt;t=t+smb_dt){

		for(Object* & object : femmodel->elements->objects){
			Element* element=xDynamicCast<Element*>(object);

			timeclim=time;
			if (linear_interp) timeinputs = t-time+timeclim;
			else timeinputs = t-time+timeclim+smb_dt/2;
			element->SmbGemb(timeinputs,count);
		}
		count=count+1;
	}

} /*}}}*/
void GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY){ /*{{{*/

	/* This file sets up the initial grid spacing and total grid depth.  The
	grid structure is set as constant grid length 'dzTop' for the top
	'zTop' meters of the model grid. Bellow 'zTop' the gid length increases
	linearly with depth */

	/*intermediary:*/
	IssmDouble dgpTop=0.0;
	int gpTop=0;
	int gpBottom=0;
	int i=0;
	IssmDouble gp0=0.0;
	IssmDouble z0=0.0;
	IssmDouble* dzT=NULL;
	IssmDouble* dzB=NULL;

	/*output: */
	IssmDouble* dz=NULL;

	//----------------------Calculate Grid Lengths------------------------------
	//calculate number of top grid points
	dgpTop = zTop/dzTop;

	//check to see if the top grid cell structure length (dzTop) goes evenly 
	//into specified top structure depth (zTop). Also make sure top grid cell
	//structure length (dzTop) is greater than 5 cm
	#ifndef _HAVE_AD_  //avoid the round operation check!
	if (dgpTop != round(dgpTop)){ 
		_error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop\n");
	}
	#endif
	if(dzTop < 0.05){
		_printf_("initial top grid cell length (dzTop) is < 0.05 m\n");
	}
	gpTop=reCast<int,IssmDouble>(dgpTop);

	//initialize top grid depth vector
	dzT = xNew<IssmDouble>(gpTop); 
	for (i=0;i<gpTop;i++)dzT[i]=dzTop;

	//bottom grid cell depth = x*zY^(cells from to structure)
	//figure out the number of grid points in the bottom vector (not known a priori)
	gp0 = dzTop;
	z0 = zTop;
	gpBottom = 0;
	while (zMax > z0){
		gp0= gp0 * zY;
		z0 = z0 + gp0;
		gpBottom++;
	}
	//initialize bottom vectors
	dzB = xNewZeroInit<IssmDouble>(gpBottom);
	gp0 = dzTop;
	z0 = zTop;
	for(i=0;i<gpBottom;i++){
		gp0=gp0*zY;
		dzB[i]=gp0;
	}

	//combine top and bottom dz vectors
	dz = xNew<IssmDouble>(gpTop+gpBottom);
	for(i=0;i<gpTop;i++){
		dz[i]=dzT[i];
	}
	for(i=0;i<gpBottom;i++){
		dz[gpTop+i]=dzB[i];
	}

	/*Free resouces:*/
	xDelete(dzT);
	xDelete(dzB);

	//---------NEED TO IMPLEMENT A PROPER GRID STRECHING ALGORITHM------------

	/*assign ouput pointers: */
	*pdz=dz; 
	*psize=gpTop+gpBottom;
} /*}}}*/ 
IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT){ /*{{{*/

	// calculates grain growth according to Fig. 9 of Marbouty, 1980
	// ------NO GRAIN GROWTH FOR d > 400 kg m-3 because H is set to zero------
	// ---------------this is a major limitation of the model-------------------

	// initialize
	IssmDouble F = 0.0, H=0.0, G=0.0;
	const IssmDouble E = 0.09;        //[mm d-1] model time growth constant E
	// convert T from K to degC
	T = T - CtoK;
	// convert dT from degC/m to degC/cm
	dT = dT/100;

	// temperature coefficient F
	if(T> -6.0+Ttol) F =  0.7 + ((T/-6.0) * 0.3);
	if(T<= -6.0+Ttol && T> -22.0+Ttol) F =  1.0 - ((T+6.0)/-16.0 * 0.8);
	if(T<= -22.0+Ttol && T> -40.0+Ttol) F =  0.2 - ((T+22.0)/-18.0 * 0.2);

	// density coefficient F
	if(d< 150.0-Dtol) H=1.0;

	if(d>= 150.0-Dtol && d <400.0-Dtol) H = 1.0 - ((d-150.0)/250.0);

	// temperature gradient coefficient G
	if(dT >= 0.16-Ttol && dT < 0.25-Ttol) G = ((dT - 0.16)/0.09) * 0.1;
	if(dT >= 0.25-Ttol && dT < 0.4-Ttol)  G = 0.1 + (((dT - 0.25)/0.15) * 0.57);
	if(dT >= 0.4-Ttol && dT < 0.5-Ttol)  G = 0.67 + (((dT - 0.4)/0.1) * 0.23);
	if(dT >= 0.5-Ttol && dT < 0.7-Ttol)  G = 0.9 + (((dT - 0.5)/0.2) * 0.1);
	if(dT >= 0.7-Ttol)  G = 1.0;

	// grouped coefficient Q
	return F*H*G*E;

} /*}}}*/
IssmDouble gardnerAlb(IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble dPHC, int m){ /*{{{*/
	//gardnerAlb(S1, c1, SZA, t, z1, S2, c2)
	//This is an implementation of the snow and ice broadband albedo
	//  parameterization developed by Alex Gardner.
	//Created By: Alex S. Gardner, Jet Propulsion Laboratory [alex.s.gardner@jpl.nasa.gov]
	//  Last Modified: June, 2014
	//Full Reference: Gardner, A. S., and Sharp, M. J.: A review of snow and
	//  ice albedo and the development of a new physically based broadband albedo
	//  parameterization, J. Geophys. Res., 115, F01009, 10.1029/2009jf001444,
	//  2010.

	//INPUTS
	// ONE LAYER
	//  - S1    : specific surface area of the snow or ice [cm^2 g-1]
	//  - c1    : concentration of light absorbing carbon  [ppm1]
	//  - SZA   : solar zenith angle of the incident radiation [deg]
	//  - t     : cloud optical thickness
	// TWO LAYER
	//  - z1    : depth of snow suface layer [mm w.e.]
	//  - S2    : specific surface area of bottom ice layer [cm^2 g-1]
	//  - c2    : concentration of light absorbing carbon of bottom ice
	//             layer [ppm1]
	IssmDouble c1=clabSnow;
	IssmDouble c2=clabIce;
	IssmDouble t=COT;
	IssmDouble a=0.0;

	//Single layer albedo parameterization
	//convert effective radius to specific surface area [cm2 g-1]
	IssmDouble S1 = 3.0 / (0.091 * re[0]);

	//effective solar zenith angle
	IssmDouble x = min(pow(t/(3*cos(Pi*SZA/180.0)),0.5), 1.0);
	IssmDouble u = 0.64*x + (1-x)*cos(Pi*SZA/180.0);

	// pure snow albedo
	IssmDouble as = 1.48 - pow(S1,-0.07);

	//change in pure snow albedo due to soot loading
	IssmDouble dac = max(0.04 - as, pow(-c1,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c1,0.6))*pow(S1,-0.25)));

	//Two layer albedo parameterization
	//  do two layer calculation if there is more than 1 layer
	IssmDouble z1=0.0;
	int lice=0;
	for(int l=0;(l<m & d[l]<dPHC-Dtol);l++){
		z1=z1+dz[l]*d[l]; //mm
		lice=l+1;
	}
	if (m>0 & lice<m & z1 > Dtol){
		// determine albedo values for bottom layer
		IssmDouble S2 = 3.0 / (0.091 * re[lice]);

		// pure snow albedo
		IssmDouble as2 = 1.48 - pow(S2,-0.07);

		// change in pure snow albedo due to soot loading
		IssmDouble dac2 = max(0.04 - as2, pow(-c2,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c2,0.6))*pow(S2,-0.25)));

		// determine the effective change due to finite depth and soot loading
		IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1-as) - 0.1*c1 - 0.13)));

		dac =  (as2 + dac2 - as) + A*((as + dac) - (as2 + dac2));
	}

	// change in albedo due to solar zenith angle
	IssmDouble dasz = 0.53*as*(1 - (as + dac))*pow(1 - u,1.2);

	// change in albedo due to cloud (apart from change in diffuse fraction)
	IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1 + 1.5*t,as));

	// Broadband albedo
	a = as + dac + dasz + dat;

	return a;
}   /*}}}*/
void grainGrowth(IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx,int sid){ /*{{{*/

	/*Created by: Alex S. Gardner, University of Alberta

	 *Description*: models the effective snow grain size

	 *Reference:*
	 DENDRITIC SNOW METAMORPHISM:
	 Brun, E., P. David, M. Sudul, and G. Brunot, 1992: A numerical model to
	 simulate snow-cover stratigraphy for operational avalanche forecasting.
	 Journal of Glaciology, 38, 13-22.

	 NONDENDRITIC SNOW METAMORPHISM:
	 Dry snow metamorphism:
	 Marbouty, D., 1980: An experimental study of temperature-gradient
	 metamorphism. Journal of Glaciology, 26, 303-312.

	 WET SNOW METAMORPHISM:
	 Brun, E., 1989: Investigation on wet-snow metamorphism in respect of
	 liquid-water content. Annals of Glaciology, 13, 22-26.

	 INPUTS
	 * T: grid cell temperature [K]
	 * dz: grid cell depth [m]
	 * d: grid cell density [kg m-3]
	 * W: water content [kg/m^2]
	 * re: effective grain size [mm]
	 * gdn: grain dentricity
	 * gsp: grain sphericity
	 * dt: time step of input data [s]

	 OUTPUTS
	 * re: effective grain size [mm]
	 * gdn: grain dentricity
	 * gsp: grain sphericity*/

	/*intermediary: */
	IssmDouble  dt=0.0;
	IssmDouble  Ti=0.0;
	IssmDouble* gsz=NULL;
	IssmDouble* dT=NULL;
	IssmDouble* zGPC=NULL;
	IssmDouble* lwc=NULL;
	IssmDouble  Q=0.0;

	/*output: */
	IssmDouble* re=NULL;
	IssmDouble* gdn=NULL;
	IssmDouble* gsp=NULL;

	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   grain growth module\n");

	/*Recover pointers: */
	re=*pre;
	gdn=*pgdn;
	gsp=*pgsp;

	/*only when aIdx = 1 or 2 do we run grainGrowth: */
	if(aIdx!=1 && aIdx!=2){
		/*come out as we came in: */
		return;
	}

	/*Figure out grain size from effective grain radius: */
	gsz=xNew<IssmDouble>(m); for(int i=0;i<m;i++)gsz[i]=re[i]*2.0;

	/*Convert dt from seconds to day: */
	dt=smb_dt/dts;

	/*Determine liquid-water content in percentage: */
	lwc=xNew<IssmDouble>(m); for(int i=0;i<m;i++)lwc[i]= W[i] / (d[i]*dz[i])*100.0;

	//set maximum water content by mass to 9 percent (Brun, 1980)
	for(int i=0;i<m;i++)if(lwc[i]>9.0+Wtol) lwc[i]=9.0;

	/* Calculate temperature gradiant across grid cells 
	 * Returns the average gradient across the upper and lower grid cell */

	//initialize
	dT=xNewZeroInit<IssmDouble>(m); 

	//depth of grid point center from surface
	zGPC=xNewZeroInit<IssmDouble>(m);  
	for(int i=0;i<m;i++){
		for (int j=0;j<=i;j++) zGPC[i]+=dz[j];
		zGPC[i]-=(dz[i]/2.0);
	}

	// Take forward differences on left and right edges
	if(m>2){
		dT[0] = (T[2] - T[0])/(zGPC[2]-zGPC[0]);
		dT[m-1] = (T[m-1] - T[m-3])/(zGPC[m-1]-zGPC[m-3]);
	}
	else if(m>1){
		dT[0] = (T[1] - T[0])/(zGPC[1]-zGPC[0]);
		dT[m-1] = (T[m-1] - T[m-2])/(zGPC[m-1]-zGPC[m-2]);
	}

	//Take centered differences on interior points
	for(int i=1;i<m-1;i++) dT[i] = (T[i+1]-T[i-1])/(zGPC[i+1]-zGPC[i-1]);

	// take absolute value of temperature gradient
	for(int i=0;i<m;i++)dT[i]=fabs(dT[i]);

	/*Snow metamorphism. Depends on value of snow dendricity and wetness of the snowpack: */
	for(int i=0;i<m;i++){

		// T for this layer
		Ti = T[i];

		if (gdn[i]>0.0+Gdntol){

			if(W[i]<=0.0+Wtol){
				//_printf_("Dendritic dry snow metamorphism\n");
				//index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1
				if(fabs(dT[i])<=5.0+Ttol){
					//determine coefficients
					IssmDouble A = - 2e8 * exp(-6e3 / Ti) * dt;
					IssmDouble B = 1e9 * exp(-6e3 / Ti) * dt;
					//new dentricity and sphericity for dT < 5 degC m-1
					gdn[i] += A;
					gsp[i] += B;
				}
				else{
					// new dendricity and sphericity for dT >= 5 degC m-1

					//determine coefficients
					IssmDouble C = (-2e8 * exp(-6e3 / Ti) * dt) * pow(fabs(dT[i]),.4);
					gdn[i] +=C;
					gsp[i] +=C;
				}
			}
			else{
				// wet snow metamorphism
				//_printf_("Dendritic wet snow metamorphism\n");

				//determine coefficient
				IssmDouble D = (1.0/16.0) * pow(lwc[i],3.0) * dt;

				// new dendricity for wet snow
				gdn[i] -= D;
				// new sphericity for wet snow
				gsp[i] += D;
			}
			// dendricity and sphericity can not be > 1 or < 0
         if (gdn[i]<=0.0+Gdntol)gdn[i]=0.0;
         if (gsp[i]<=0.0+Gdntol)gsp[i]=0.0;
         if (gdn[i]>=1.0-Gdntol)gdn[i]=1.0;
         if (gsp[i]>=1.0-Gdntol)gsp[i]=1.0;

         // determine new grain size (mm)
			gsz[i] = 1e-1*(gdn[i]/.99+(1-1*gdn[i]/.99)*(gsp[i]/.99*3+(1-gsp[i]/.99)*4));

		}
		else{

			//When wet-snow grains (class 6) are submitted to a
			// temperature gradient higher than 5 degC m-1, their sphericity
			// decreases according to Equations (4). When sphericity
			// reaches 0, their size increases according to the functions
			// determined by Marbouty. (Brun et al., 1992)
			if(gsp[i]>0.0+Gdntol && gsp[i]<1.0-Gdntol){

				IssmDouble F = 0.0;

				if (fabs(dT[i])>5.0+Ttol){
					F = (-2e8 * exp(-6e3 / Ti) * dt) * pow(fabs(dT[i]),.4);
				}
				else if (W[i]>0.0+Wtol){
					F = (1.0/16.0) * pow(lwc[i],3.0) * dt;
				}
				else{
					F = 1e9 * exp(-6e3 / Ti) * dt;
				}
				gsp[i] +=F;

			}
			if (gsp[i]<=0.0+Gdntol)gsp[i]=0.0;
			if (gsp[i]>=1.0-Gdntol)gsp[i]=1.0;

			/*Dry snow metamorphism (Marbouty, 1980) grouped model coefficents
			 *from Marbouty, 1980: Figure 9*/
			if(W[i]<=0.0+Wtol || (gsp[i]<=0.0+Gdntol && fabs(dT[i])>5.0+Ttol)){
				//_printf_("Nondendritic snow metamorphism\n");
				Q = Marbouty(Ti,d[i],dT[i]);

				// calculate grain growth
				gsz[i] += (Q*dt);
			}
			//Wet snow metamorphism (Brun, 1989)
			else{
				//_printf_("Nondendritic wet snow metamorphism\n");
				//wet rate of change coefficient
				IssmDouble E = (1.28e-8 + 4.22e-10 * pow(lwc[i],3.0))* (dt *dts);   // [mm^3 s^-1]

				// calculate change in grain volume and convert to grain size
				gsz[i] = 2.0 * pow(3.0/(Pi * 4.0)*((4.0/ 3.0)*Pi*pow(gsz[i]/2.0,3.0) + E),1.0/3.0);
			}

			// grains with sphericity == 1 can not have grain sizes > 2 mm (Brun, 1992)
			if (fabs(gsp[i]-1.0)<Wtol && gsz[i]>2.0-Wtol) gsz[i]=2.0;

			// grains with sphericity == 0 can not have grain sizes > 5 mm (Brun, 1992)
			if (fabs(gsp[i]-1.0)>=Wtol && gsz[i]>5.0-Wtol) gsz[i]=5.0;
		}

		//convert grain size back to effective grain radius:
		re[i]=gsz[i]/2.0;
	}

	/*Free resources:*/
	xDelete<IssmDouble>(gsz);
	xDelete<IssmDouble>(dT);
	xDelete<IssmDouble>(zGPC);
	xDelete<IssmDouble>(lwc);

	/*Assign output pointers:*/
	*pre=re;
	*pgdn=gdn;
	*pgsp=gsp;

}  /*}}}*/
void albedo(IssmDouble** pa, IssmDouble** padiff, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/

	//// Calculates Snow, firn and ice albedo as a function of:
	//   0 : direct input from aValue parameter
	//   1 : effective grain radius (Gardner & Sharp, 2009)
	//   2 : effective grain radius (Brun et al., 1992, Lefebre et al., 2003)
	//   3 : density and cloud amount (Greuell & Konzelmann, 1994)
	//   4 : exponential time decay & wetness (Bougamont & Bamber, 2005)

	//// Inputs
	// aIdx      = albedo method to use

	// Method 0
	//  aValue   = direct input value for albedo, override all changes to albedo

	// adThresh
	//  Apply below method to all areas with densities below this value, 
	//  or else apply direct input value, allowing albedo to be altered.  
	//  Default value is rho water (1023 kg m-3).

	// Methods 1 & 2
	//   re      = surface effective grain radius [mm]
	// Method 1, optional
	//  clabSnow = concentration of light absorbing carbon  [ppm1], default 0
	//  SZA      = solar zenith angle of the incident radiation [deg], default 0
	//  COT      = cloud optical thickness, default 0
	//  For TWO LAYER
	//  clabIce  = concentration of light absorbing carbon of first ice layer [ppm1], default 0

	// Method 3
	//   d       = snow surface density [kg m-3]
	//   n       = cloud amount
	//   aIce    = albedo of ice
	//   aSnow   = albedo of fresh snow

	// Method 4
	//   aIce    = albedo of ice
	//   aSnow   = albedo of fresh snow
	//   a       = grid cell albedo from prevous time step;
	//   T       = grid cell temperature [k]
	//   W       = pore water [kg]
	//   P       = precipitation [mm w.e.] or [kg m-3]
	//   EC      = surface evaporation (-) condensation (+) [kg m-2]
	//   t0wet   = time scale for wet snow (15-21.9) [d]
	//   t0dry   = warm snow timescale [15] [d]
	//   K       = time scale temperature coef. (7) [d]
	//   dt      = time step of input data [s]

	//// Output
	//   a       = grid cell albedo 

	//// Usage 
	// Method 1
	// a = albedo(1, 0.1); 

	// Method 4
	// a = albedo(4, [], [], [], 0.48, 0.85, [0.8 0.5 ... 0.48], ...
	//   [273 272.5 ... 265], [0 0.001 ... 0], 0, 0.01, 15, 15, 7, 3600)

	/*output: */
	IssmDouble* a=NULL;
	IssmDouble* adiff=NULL;

	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   albedo module\n");

	/*Recover pointers: */
	a=*pa;
	adiff=*padiff;

	//some constants:
	const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
	const IssmDouble dPHC = 830;  //Pore closeoff density
	const IssmDouble ai_max = 0.58;  //maximum ice albedo, from Lefebre,2003
	const IssmDouble ai_min = aIce;  //minimum ice albedo
	const IssmDouble as_min = 0.65;  //minimum snow albedo, from Alexander 2014

	if(aIdx==0 || (adThresh - d[0])<Dtol){
		a[0] = aValue;
	}
	else{
		if(aIdx==1){ 
			//function of effective grain radius
			// clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, int m
			a[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, SZA, COT, dPHC, m);
			adiff[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, 50, COT, dPHC, m);
		}
		else if(aIdx==2){

			// Spectral fractions  (Lefebre et al., 2003)
			// [0.3-0.8um 0.8-1.5um 1.5-2.8um]

			IssmDouble sF[3] = {0.606, 0.301, 0.093};

			// convert effective radius to grain size in meters
			IssmDouble gsz = (re[0] * 2.0) / 1000.0;

			// spectral range:
			// 0.3 - 0.8um
			IssmDouble a0 = min(0.98, 0.95 - 1.58 *pow(gsz,0.5));
			// 0.8 - 1.5um
			IssmDouble a1 = max(0., 0.95 - 15.4 *pow(gsz,0.5));
			// 1.5 - 2.8um
			IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5));

			// broadband surface albedo
			a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2;

		}
		else if(aIdx==3){

			// a as a function of density

			// calculate albedo
			a[0] = aIce + (d[0] - dIce)*(aSnow - aIce) / (dSnow - dIce) + (0.05 * (cldFrac - 0.5));
		}
		else if(aIdx==4){

			// exponential time decay & wetness

			// change in albedo with time:
			//   (d_a) = (a - a_old)/(t0)
			// where: t0 = timescale for albedo decay

			dt = dt / dts;    // convert from [s] to [d]

			// initialize variables
			IssmDouble* t0=xNew<IssmDouble>(m);
			IssmDouble* T=xNew<IssmDouble>(m);
			IssmDouble* t0warm=xNew<IssmDouble>(m);
			IssmDouble* d_a=xNew<IssmDouble>(m);

			// specify constants
			// a_wet = 0.15;        // water albedo (0.15)
			// a_new = aSnow        // new snow albedo (0.64 - 0.89)
			// a_old = aIce;        // old snow/ice albedo (0.27-0.53)
			// t0_wet = t0wet;      // time scale for wet snow (15-21.9) [d]
			// t0_dry = t0dry;      // warm snow timescale [15] [d]
			// K = 7                // time scale temperature coef. (7) [d]
			// W0 = 300;            // 200 - 600 [mm]
			const IssmDouble z_snow = 15.0;            // 16 - 32 [mm]

			// determine timescale for albedo decay
			for(int i=0;i<m;i++)if(W[i]>0.0+Wtol)t0[i]=t0wet; // wet snow timescale
			for(int i=0;i<m;i++)T[i]=TK[i] - CtoK; // change T from K to degC
			for(int i=0;i<m;i++) t0warm[i]= fabs(T[i]) * K + t0dry; //// 'warm' snow timescale
			for(int i=0;i<m;i++)if(fabs(W[i])<Wtol && T[i]>=-10.0-Ttol)t0[i]= t0warm[i];
			for(int i=0;i<m;i++)if(T[i]<-10.0-Ttol) t0[i] =  10.0 * K + t0dry; // 'cold' snow timescale

			// calculate new albedo
			for(int i=0;i<m;i++)d_a[i] = (a[i] - aIce) / t0[i] * dt;           // change in albedo
			for(int i=0;i<m;i++)a[i] -= d_a[i];                            // new albedo

			// modification of albedo due to thin layer of snow or solid
			// condensation (deposition) at the surface 

			// check if condensation occurs & if it is deposited in solid phase
			if ( EC > 0.0 + Dtol && T[0] < 0.0-Ttol) P = P + (EC/dSnow) * 1000.0;  // add cond to precip [mm]

			a[0] = aSnow - (aSnow - a[0]) * exp(-P/z_snow);

			//----------THIS NEEDS TO BE IMPLEMENTED AT A LATER DATE------------
			// modification of albedo due to thin layer of water on the surface
			// a_surf = a_wet - (a_wet - a_surf) * exp(-W_surf/W0);

			/*Free resources:*/
			xDelete<IssmDouble>(t0);
			xDelete<IssmDouble>(T);
			xDelete<IssmDouble>(t0warm);
			xDelete<IssmDouble>(d_a);

		}
		else _error_("albedo method switch should range from 0 to 4!");

		//If we do not have fresh snow
		if (aIdx<3 && aIdx>0 && (adThresh - d[0])>=Dtol){
			// In a snow layer < 10cm, account for mix of ice and snow,
			// after P. Alexander et al., 2014
			IssmDouble depthsnow=0.0;
			IssmDouble aice=0.0;
			int lice=0;
			for(int l=0;(l<m & d[l]<dPHC-Dtol);l++){
				depthsnow=depthsnow+dz[l];
				lice=l+1;
			}
			if (depthsnow<=0.1+Dtol & lice<m & d[lice]>=dPHC-Dtol){
				aice = ai_max + (as_min - ai_max)*(d[lice]-dIce)/(dPHC-dIce);
				a[0]= aice + max(a[0]-aice,0.0)*(depthsnow/0.1);
			}

			if (d[0]>=dPHC-Dtol){
				if (d[0]<dIce-Dtol){ //For continuity of albedo in firn i.e. P. Alexander et al., 2014

					//ai=ai_max + (as_min - ai_max)*(dI-dIce)/(dPHC-dIce);
					//dPHC is pore close off (830 kg m^-3)
					//dI is density of the upper firn layer

					a[0] = ai_max + (as_min - ai_max)*(d[0]-dIce)/(dPHC-dIce);

				}
				else{ //surface layer is density of ice

					//When density is > dIce (typically 910 kg m^-3, 920 is used by Alexander in MAR),
					//ai=ai_min + (ai_max - ai_min)*e^(-1*(Msw(t)/K))
					//K is a scale factor (set to 200 kg m^-2)
					//Msw(t) is the time-dependent accumulated amount of excessive surface meltwater
					//  before run-off in kg m^-2 (melt per GEMB timestep, i.e. 3 hourly)
					IssmDouble M = Msurf+W[0];
					a[0]=max(ai_min + (ai_max - ai_min)*exp(-1*(M/200)), ai_min);

				}
			}
		}
	}

	// Check for erroneous values
	if (a[0] > 1) _printf_("albedo > 1.0\n");
	else if (a[0] < 0) _printf_("albedo is negative\n");
	else if (xIsNan(a[0])) _error_("albedo == NAN\n");

	/*Assign output pointers:*/
	*pa=a;
	*padiff=adiff;

}  /*}}}*/
void thermo(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble** pT, IssmDouble* pulwrf, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble tcIdx, IssmDouble eIdx, IssmDouble teValue, IssmDouble dulwrfValue, IssmDouble teThresh, IssmDouble Ws, IssmDouble dt0, IssmDouble dzMin, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT, bool isdeltaLWup) { /*{{{*/

	/* ENGLACIAL THERMODYNAMICS*/

	/* Description: 
	   computes new temperature profile accounting for energy absorption and 
	   thermal diffusion.*/

	// INPUTS
	//  T: grid cell temperature [k]
	//  dz: grid cell depth [m]
	//  d: grid cell density [kg m-3]
	//  swf: shortwave radiation fluxes [W m-2]
	//  dlwrf: downward longwave radiation fluxes [W m-2]
	//  Ta: 2 m air temperature
	//  V:  wind velocity [m s-1]
	//  eAir: screen level vapor pressure [Pa]
	//  Ws: surface water content [kg]
	//  dt0: time step of input data [s]
	//  elev: surface elevation [m a.s.l.] 
	//  Vz: air temperature height above surface [m]
	//  Tz: wind height above surface [m]
	//  thermo_scaling: scaling factor to multiply the thermal diffusion timestep (delta t) 

	// OUTPUTS
	// T: grid cell temperature [k]
	// EC: evaporation/condensation [kg]
	// ulwrf: upward longwave radiation flux [W m-2]

	/*intermediary: */
	IssmDouble* K = NULL;
	IssmDouble* KU = NULL;
	IssmDouble* KD = NULL;
	IssmDouble* KP = NULL;
	IssmDouble* Au = NULL;
	IssmDouble* Ad = NULL;
	IssmDouble* Ap = NULL;
	IssmDouble* Nu = NULL;
	IssmDouble* Nd = NULL;
	IssmDouble* Np = NULL;
	IssmDouble* dzU = NULL;
	IssmDouble* dzD = NULL;
	IssmDouble* sw = NULL;
	IssmDouble* dT_sw = NULL;
	IssmDouble* T0 = NULL;
	IssmDouble* Tu = NULL;
	IssmDouble* Td = NULL;

	IssmDouble z0=0.0;	
	IssmDouble zT=0.0;
	IssmDouble zQ=0.0;
	IssmDouble zratio=1;
	IssmDouble dt=0.0;
	IssmDouble max_fdt=0.0;
	IssmDouble Ts=0.0;
	IssmDouble L=0.0;
	IssmDouble eS=0.0;
	IssmDouble Ri=0.0;
	IssmDouble coefM=0.0;
	IssmDouble coefH=0.0;
	IssmDouble coefHT=0.0;
	IssmDouble coefHQ=0.0;
	IssmDouble An_num=0.0;
	IssmDouble An_den_T=0.0;
	IssmDouble An_den_Q=0.0;
	IssmDouble An=0.0;
	IssmDouble C=0.0;
	IssmDouble shf=0.0;
	IssmDouble ds=0.0;
	IssmDouble dAir=0.0;
	IssmDouble TCs=0.0;
	IssmDouble lhf=0.0;
	IssmDouble EC_day=0.0;
	IssmDouble lhf_cum=0.0;
	IssmDouble shf_cum=0.0;
	IssmDouble dT_turb=0.0;
	IssmDouble turb=0.0;
	IssmDouble ulw=0.0;
	IssmDouble dT_ulw=0.0;
	IssmDouble dlw=0.0;
	IssmDouble dT_dlw=0.0;

	/*outputs:*/
	IssmDouble EC=0.0;
	IssmDouble* T=*pT;
	IssmDouble ulwrf=0.0;

	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   thermal module\n");

	ds = d[0];      // density of top grid cell

	// calculated air density [kg/m3]
	dAir = 0.029 * pAir /(R * Ta);

	// thermal capacity of top grid cell [J/k]
	TCs = d[0]*dz[0]*CI; 

	//initialize Evaporation - Condenstation 
	EC = 0.0;

	// check if all SW applied to surface or distributed throught subsurface
	// swIdx = length(swf) > 1

	// SURFACE ROUGHNESS (Bougamont, 2006)
	// wind/temperature surface roughness height [m]
	if (ds < dIce-Dtol && Ws < Wtol) z0 = 0.00012;        // 0.12 mm for dry snow
	else if (ds >= dIce-Dtol) z0 = 0.0032;             // 3.2 mm for ice
	else z0 = 0.0013;                            // 1.3 mm for wet snow

	//zT and zQ are percentage of z0 (Foken 2008)
	zratio=10;
	zT=z0/zratio;
	zQ=z0/zratio;

	// if V = 0 goes to infinity therfore if V = 0 change
	if(V<0.01-Dtol)V=0.01;

	// Bulk-transfer coefficient for turbulent fluxes
	An =  pow(0.4,2); // Bulk-transfer coefficient
	C = An*V;  // shf & lhf common coefficient
	An_den_T = (log(Tz/zT)*log(Vz/z0)); 
	An_den_Q = (log(Tz/zQ)*log(Vz/z0)); 

	// THERMAL CONDUCTIVITY (Sturm, 1997: J. Glaciology)
	// calculate new K profile [W m-1 K-1]

	// initialize conductivity
	K= xNewZeroInit<IssmDouble>(m);

	// for snow and firn (density < 910 kg m-3) (Sturm et al, 1997) or (Calonne et al., 2011)
	if (tcIdx == 2){
		for(int i=0;i<m;i++) if(d[i]<dIce-Dtol) K[i] = 0.024 - 1.23e-4 * d[i] + 2.5e-6 * (pow(d[i],2));
	}
	else{ //default (Sturm et al, 1997)
		for(int i=0;i<m;i++) if(d[i]<dIce-Dtol) K[i] = 0.138 - 1.01e-3 * d[i] + 3.233e-6 * (pow(d[i],2));
	}

	// for ice (density >= 910 kg m-3)
	for(int i=0;i<m;i++) if(d[i]>=dIce-Dtol) K[i] = 9.828 * exp(-5.7e-3*T[i]);

	// THERMAL DIFFUSION COEFFICIENTS

	// A discretization scheme which truncates the Taylor-Series expansion
	// after the 3rd term is used. See Patankar 1980, Ch. 3&4

	// discretized heat equation:

	//                 Tp = (Au*Tuo+ Ad*Tdo+ (Ap-Au-Ad)Tpo+ S) / Ap

	// where neighbor coefficients Au, Ap, & Ad are

	//                   Au = [dz_u/2KU + dz/2KP]^-1
	//                   Ad = [dz_d/2KD + dz/2KP]^-1
	//                   Ap = d*CI*dz/Dt 

	// and u & d represent grid points up and down from the center grid point 
	// point p and o identifies previous time step values. S is a source term.

	// u, d, and p conductivities
	KU = xNew<IssmDouble>(m);
	KD = xNew<IssmDouble>(m);
	KP = xNew<IssmDouble>(m);

	KU[0] = Delflag; //Thermal conductivity of air = 0.025 W/m/K
	KD[m-1] = Delflag;
	for(int i=1;i<m;i++) KU[i]= K[i-1];
	for(int i=0;i<m-1;i++) KD[i] = K[i+1];
	for(int i=0;i<m;i++) KP[i] = K[i];

	// determine u, d & p cell widths
	dzU = xNew<IssmDouble>(m);
	dzD = xNew<IssmDouble>(m);
	dzU[0]=Delflag;
	dzD[m-1]=Delflag;

	for(int i=1;i<m;i++) dzU[i]= dz[i-1];
	for(int i=0;i<m-1;i++) dzD[i] = dz[i+1];

	// determine minimum acceptable delta t (diffusion number > 1/2) [s]
	// NS: 2.16.18 divided dt by scaling factor, default set to 1/11 for stability
	dt=1e12; 
	for(int i=0;i<m;i++) dt = min(dt,CI * pow(dz[i],2) * d[i]  / (3 * K[i]) * thermo_scaling);

	// smallest possible even integer of 60 min where diffusion number > 1/2
	// must go evenly into one hour or the data frequency if it is smaller

	// all integer factors of the number of second in a day (86400 [s])
	int f[45] = {1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 30, 36, 40, 45, 48, 50, 60,
    72, 75, 80, 90, 100, 120, 144, 150, 180, 200, 225, 240, 300, 360, 400, 450, 600, 720, 900, 1200, 1800, 3600};

	// return the min integer factor that is < dt
	max_fdt=f[0];
	for(int i=0;i<45;i++){
		if (f[i]<dt-Dtol)if(f[i]>=max_fdt-Dtol)max_fdt=f[i];
	}
	dt=max_fdt;

	// determine mean (harmonic mean) of K/dz for u, d, & p
	Au = xNew<IssmDouble>(m);
	Ad = xNew<IssmDouble>(m);
	Ap = xNew<IssmDouble>(m);
	for(int i=0;i<m;i++){
		Au[i] = pow((dzU[i]/2.0/KU[i] + dz[i]/2.0/KP[i]),-1.0);
		Ad[i] = pow((dzD[i]/2.0/KD[i] + dz[i]/2.0/KP[i]),-1.0);
		Ap[i] = (d[i]*dz[i]*CI)/dt;
	}

	// create "neighbor" coefficient matrix
	Nu = xNewZeroInit<IssmDouble>(m);
	Nd = xNewZeroInit<IssmDouble>(m);
	Np = xNewZeroInit<IssmDouble>(m);
	for(int i=0;i<m;i++){
		Nu[i] = Au[i] / Ap[i];
		Nd[i] = Ad[i] / Ap[i];
		Np[i]= 1.0 - Nu[i] - Nd[i];
	}

	// specify boundary conditions: constant flux at bottom
	Nu[m-1] = 0.0;
	Np[m-1] = 1.0;

	// zero flux at surface
	Np[0] = 1.0 - Nd[0];

	// Create neighbor arrays for diffusion calculations instead of a tridiagonal matrix
	Nu[0] = 0.0;
	Nd[m-1] = 0.0;

	/* RADIATIVE FLUXES*/

	// energy supplied by shortwave radiation [J]
	sw = xNew<IssmDouble>(m);
	for(int i=0;i<m;i++) sw[i]= swf[i] * dt;

	// temperature change due to SW
	dT_sw = xNew<IssmDouble>(m);
	for(int i=0;i<m;i++) dT_sw[i]= sw[i] / (CI * d[i] * dz[i]);

	// Upward longwave radiation flux is calculated from the snow surface
	// temperature which is set equal to the average temperature of the
	// top grid cells.

	// energy supplied by downward longwave radiation to the top grid cell [J]
	dlw = dlwrf * dt;

	// temperature change due to dlw_surf
	dT_dlw = dlw / TCs;

	// PREALLOCATE ARRAYS BEFORE LOOP FOR IMPROVED PERFORMANCE
	T0 = xNewZeroInit<IssmDouble>(m+2);
	Tu=xNew<IssmDouble>(m);
	Td=xNew<IssmDouble>(m);

	/* CALCULATE ENERGY SOURCES AND DIFFUSION FOR EVERY TIME STEP [dt]*/
	for (IssmDouble i=1;i<=dt0;i+=dt){

		// PART OF ENERGY CONSERVATION CHECK
		// store initial temperature
		//T_init = T;

		// calculate temperature of snow surface (Ts)
		Ts = T[0];
		Ts = min(CtoK,Ts);    // don't allow Ts to exceed 273.15 K (0 degC)

		//TURBULENT HEAT FLUX

		// Monin-Obukhov Stability Correction
		// Reference:
		// Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
		// Journal of Climatology, 2, 65-84.

		// calculate the Bulk Richardson Number (Ri)
		Ri = pow(100000/pAir,0.286)*(2.0*9.81*(Ta - Ts)) / (Tz*(Ta + Ts)* pow(V/(Vz),2.0));

		IssmDouble PhiM;
		IssmDouble PhiH;

		// calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
		if (false){
			// do not allow Ri to exceed 0.16
			Ri = min(Ri, 0.16); //Ohmura, 1982

			// calculate momentum 'coefM' stability factor
			if (Ri > 0.0+Ttol){
				// if stable
				coefM = 1.0/(1.0-5.2*Ri);
			}
			else {
				coefM =pow (1.0-18*Ri,-0.25);
			}

			// calculate heat/wind 'coef_H' stability factor
			if (Ri <= -0.03+Ttol) coefH = coefM/1.3;
			else coefH = coefM;

         coefHT = coefH*An_den_T;
         coefHQ = coefH*An_den_Q;

		}
		else if(false){

			// do not allow Ri to exceed 0.19
			Ri = min(Ri, 0.19); //Ohmura, 1982

			// calculate momentum 'coefM' stability factor
			if (Ri > 0.0+Ttol){
				// if stable
				//coefM = pow(1.0-5.0*Ri,2.0); //Fitzpatrick et al., 2017, from Brock et al., 2010
				coefM=1+5.3*min((Ri/(1-5*Ri)),0.5);
				coefH=1+8*min((Ri/(1-5*Ri)),0.5);
			}
			else {
				//coefM =pow(1.0-16.0*max(Ri,-1.0),0.75); //Fitzpatrick et al., 2017, from Brock et al., 2010
				coefM=pow(1-19*max(Ri/1.5,-2.0),-0.25);
				coefH=0.95*pow(1-11.6*max(Ri/1.5,-2.0),-0.5);
			}

			coefHT = coefH*An_den_T;
			coefHQ = coefH*An_den_Q;

		}
      else if (false){
         // Greuell and Konzelman, 1994
         // calculate momentum 'coefM' stability factor

         if (Ri > 0.0+Ttol){
            // if stable
            coefM=1+15*Ri*pow(1+Ri,1/2);
            coefH=1;
         }
         else {
            coefM=pow(1-15*Ri/(1+75*pow(0.4/log(Tz/zT),2)*pow(Tz/zT*fabs(Ri),1/2)),-1);
            coefH=1;
         }

         coefHT = coefH*An_den_T;
         coefHQ = coefH*An_den_Q;

      }
      else {
         IssmDouble a1=1.0;
         IssmDouble b1=2/3;
         IssmDouble c1=5.0;
         IssmDouble d1=0.35;
         IssmDouble PhiMz;
         IssmDouble PhiHz;
         IssmDouble PhiMz0=0.0;
         IssmDouble PhiHzT=0.0;
         IssmDouble PhiHzQ=0.0;
         IssmDouble zL;
         IssmDouble zLT;
         IssmDouble zLM;

         if (Ri > 0.0+Ttol){
            // if stable

            if(Ri < 0.2-Ttol){
               zL = Ri/(1-5*Ri);
            }
            else{
               zL=Ri;
            }
            //zL = min(zL, 0.5); //Sjoblom, 2014
            zLM=max(zL/Vz*z0,1e-3);
            zLT=max(zL/Tz*zT,1e-3);

            //Ding et al. 2020, from Beljaars and Holtslag (1991)
            PhiMz=-1*(a1*zL + b1*(zL-c1/d1)*exp(-1*d1*zL) + b1*c1/d1);
            PhiHz=-1*(pow(1+2*a1*zL/3,1.5) + b1*(zL-c1/d1)*exp(-1*d1*zL) + b1*c1/d1 - 1);
            PhiMz0=-1*(a1*zLM + b1*(zLM-c1/d1)*exp(-1*d1*zLM) + b1*c1/d1);
            PhiHzT=-1*(pow(1+2*a1*zLT/3,1.5) + b1*(zLT-c1/d1)*exp(-1*d1*zLT) + b1*c1/d1 - 1);

            PhiHzQ=PhiHzT;
         }
         else {
            IssmDouble xm;
            IssmDouble xh;
            IssmDouble xmT;
            IssmDouble xmM;

            zL = Ri/1.5; //max(Ri, -0.5+Ttol)/1.5; //Hogstrom (1996)
            //zL = max(zL, -2.0); //Sjoblom, 2014
            zLM=min(zL/Vz*z0,-1e-3);
            zLT=min(zL/Tz*zT,-1e-3);

            if (true){ //Sjoblom, 2014
               xm=pow(1.0-19*zL,-0.25);
               PhiMz=2*log((1+xm)/2.0) + log((1+pow(xm,2))/2.0) - 2*atan(xm) + Pi/2;

               xh=0.95*pow(1.0-11.6*zL,-0.5);
               PhiHz=2*log((1+pow(xh,2))/2.0);
            }
            else{ //Ding et al., 2020
               xm=pow(1.0-16*zL,0.25);
               xmM=pow(1.0-16*zLM,0.25);
               xmT=pow(1.0-16*zLT,0.25);
               PhiMz=2*log((1+xm)/2.0) + log((1+pow(xm,2))/2.0) - 2*atan(xm) + Pi/2;
               PhiMz0=2*log((1+xmM)/2.0) + log((1+pow(xmM,2))/2.0) - 2*atan(xmM) + Pi/2;

               PhiHz=2*log((1+pow(xm,2))/2.0);
               PhiHzT=2*log((1+pow(xmT,2))/2.0);

               PhiHzQ=PhiHzT;
            }
         }

         PhiM=PhiMz;
         PhiH=PhiHz;
         coefM = log(Vz/z0) - PhiMz + PhiMz0; //Ding et al., 2019
         coefHT = log(Tz/zT) - PhiHz + PhiHzT; //Sjoblom, 2014, after Foken 2008
         coefHQ = log(Tz/zQ) - PhiHz + PhiHzQ; //Sjoblom, 2014, after Foken 2008

      }

      //// Sensible Heat
      // calculate the sensible heat flux [W m-2](Patterson, 1998)
      shf = dAir * C * CA * (Ta - Ts) * pow(100000/pAir,0.286);

      // adjust using Monin-Obukhov stability theory
      shf = shf/(coefM*coefHT);

      //// Latent Heat
      // determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water
      if (Ts >= CtoK-Ttol){
         L = LV; //for liquid water at 273.15 k to vapor

         //for liquid surface (assume liquid on surface when Ts == 0 deg C)
         // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A
         //eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
         // Murray 1967, https://cran.r-project.org/web/packages/humidity/vignettes/humidity-measures.html
         eS = 610.78 * exp(17.2693882 * (Ts - CtoK - 0.01) / (Ts - 35.86));
      }
      else{
         L = LS; // latent heat of sublimation

         // for an ice surface Murphy and Koop, 2005 [Equation 7]
         //eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
         // for an ice surface Ding et al., 2019 after Bolton, 1980
         eS = 610.78 * exp(21.8745584 * (Ts - CtoK - 0.01) / (Ts - 7.66));
      }

      // Latent heat flux [W m-2]
      lhf = C * L * (eAir - eS) / (461.9*(Ta+Ts)/2);

      // adjust using Monin-Obukhov stability theory (if lhf '+' then there is energy and mass gained at the surface,
      // if '-' then there is mass and energy loss at the surface.
      lhf = lhf/(coefM*coefHQ);

		//mass loss (-)/acreation(+) due to evaporation/condensation [kg]
		EC_day = lhf * dts / L;

		// temperature change due turbulent fluxes
		turb = (shf + lhf)* dt;
		dT_turb = turb  / TCs;

		// upward longwave contribution
		IssmDouble deltaULW=0.0;
		IssmDouble emissivity=1.0;
		//If user wants to set a upward long wave bias
		if(isdeltaLWup) deltaULW = dulwrfValue;
		//If user wants to directly set emissivity, or grain radius is larger than the
		// threshold, or eIdx is 2 and we have wet snow or ice, use prescribed emissivity
		if(eIdx==0 || (teThresh - re[0])<Gdntol || (eIdx==2 && z0>0.001)) emissivity = teValue;
		ulw = - (SB * pow(Ts,4.0)* emissivity + deltaULW) * dt; 
		ulwrf = ulwrf - ulw/dt0;

		dT_ulw = ulw / TCs;

		// new grid point temperature

		//SW penetrates surface
		if(!isconstrainsurfaceT){
			for(int j=0;j<m;j++) T[j] = T[j] + dT_sw[j];
			T[0] = T[0] + dT_dlw + dT_ulw + dT_turb;
		}

		// temperature diffusion
		for(int j=0;j<m;j++) T0[1+j]=T[j];
		T0[0]=Ta;
		T0[m+1]=T[m-1];
		for(int j=0;j<m;j++) Tu[j] = T0[j];
		for(int j=0;j<m;j++) Td[j] = T0[2+j];
		for(int j=0;j<m;j++) T[j] = (Np[j] * T[j]) + (Nu[j] * Tu[j]) + (Nd[j] * Td[j]);

		// calculate cumulative evaporation (+)/condensation(-)
		if(!isconstrainsurfaceT) EC = EC + (EC_day/dts)*dt;
		lhf_cum=lhf_cum+lhf*dt/dt0;
		shf_cum=shf_cum+shf*dt/dt0;

		/* CHECK FOR ENERGY (E) CONSERVATION [UNITS: J]
		//energy flux across lower boundary (energy supplied by underling ice)
		base_flux = Ad(-1)*(T_init()-T_init(-1)) * dt;

		E_used = sum((T - T_init) * (d*dz*CI));
		E_sup = ((sum(swf)  * dt) + dlw + ulw + turb + base_flux);

		E_diff = E_used - E_sup;

		if fabs(E_diff) > 1E-6 || isnan(E_diff)
		disp(T(1))
		_error_("energy not conserved in thermodynamics equations");
		*/
	}

	/*Free resources:*/
	xDelete<IssmDouble>(K);
	xDelete<IssmDouble>(KU);
	xDelete<IssmDouble>(KD);
	xDelete<IssmDouble>(KP);
	xDelete<IssmDouble>(Au);
	xDelete<IssmDouble>(Ad);
	xDelete<IssmDouble>(Ap);
	xDelete<IssmDouble>(Nu);
	xDelete<IssmDouble>(Nd);
	xDelete<IssmDouble>(Np);
	xDelete<IssmDouble>(dzU);
	xDelete<IssmDouble>(dzD);
	xDelete<IssmDouble>(sw);
	xDelete<IssmDouble>(dT_sw);
	xDelete<IssmDouble>(T0);
	xDelete<IssmDouble>(Tu);
	xDelete<IssmDouble>(Td);

	/*Assign output pointers:*/
	*pEC=EC;
	*plhf=lhf_cum;
	*pshf=shf_cum;
	*pT=T;
	*pulwrf=ulwrf;

}  /*}}}*/
void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble dswdiff, IssmDouble as, IssmDouble asdiff, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid){ /*{{{*/

	// DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE

	// swIdx = 0 : all absorbed SW energy is assigned to the top grid cell

	// swIdx = 1 : absorbed SW is distributed with depth as a function of:
	//   default   : snow density (taken from Bassford, 2002)
	//   if aIdx=2 : grain size in 3 spectral bands (Brun et al., 1992)

	// Inputs
	//   swIdx   = shortwave allowed to penetrate surface (0 = No, 1 = Yes)
	//   aIdx    = method for calculating albedo (1-4)
	//   dsw     = downward shortwave radiative flux [w m-2]
	//   dswdiff  = downward shortwave diffuse radiative flux [w m-2]
	//   as      = surface albedo
	//   asdiff  = surface albedo for diffuse radiation
	//   d       = grid cell density [kg m-3]
	//   dz      = grid cell depth [m]
	//   re      = grid cell effective grain radius [mm]

	// Outputs
	//   swf     = absorbed shortwave radiation [W m-2]
	//

	/*outputs: */
	IssmDouble* swf=NULL;

	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   shortwave module\n");

	/*Initialize and allocate: */
	swf=xNewZeroInit<IssmDouble>(m);

	// SHORTWAVE FUNCTION
	if (swIdx == 0 | (dIce - d[0])<Dtol) {// all sw radation is absorbed in by the top grid cell

		// calculate surface shortwave radiation fluxes [W m-2]
		if (aIdx == 1){
			swf[0] = (1.0 - as) * max(0.0,(dsw - dswdiff)) +  (1.0 - asdiff) * dswdiff;
		}
		else{
			swf[0] = (1.0 - as) * dsw;
		}
	}
	else{ // sw radation is absorbed at depth within the glacier

		if (aIdx == 2){    // function of effective radius (3 spectral bands)

			IssmDouble * gsz=NULL;
			IssmDouble * B1_cum=NULL;
			IssmDouble * B2_cum=NULL;
			IssmDouble* h =NULL;
			IssmDouble* B1 =NULL;
			IssmDouble* B2 =NULL;
			IssmDouble* exp1 = NULL;
			IssmDouble* exp2 = NULL;
			IssmDouble*  Qs1 = NULL;
			IssmDouble*  Qs2 = NULL;

			// convert effective radius [mm] to grain size [m]
			gsz=xNew<IssmDouble>(m);
			for(int i=0;i<m;i++) gsz[i]= (re[i] * 2.0) / 1000.0;

			// Spectral fractions [0.3-0.8um 0.8-1.5um 1.5-2.8um]
			// (Lefebre et al., 2003)
			IssmDouble sF[3] = {0.606, 0.301, 0.093};

			// initialize variables
			B1_cum=xNew<IssmDouble>(m+1);
			B2_cum=xNew<IssmDouble>(m+1);
			for(int i=0;i<m+1;i++){
				B1_cum[i]=1.0;
				B2_cum[i]=1.0;
			}

			// spectral albedos:
			// 0.3 - 0.8um
			IssmDouble a0 = min(0.98, 0.95 - 1.58 *pow(gsz[0],0.5));
			// 0.8 - 1.5um
			IssmDouble a1 = max(0.0, 0.95 - 15.4 *pow(gsz[0],0.5));
			// 1.5 - 2.8um
			IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5));

			// separate net shortwave radiative flux into spectral ranges
			IssmDouble swfS[3];
			swfS[0] = (sF[0] * dsw) * (1.0 - a0);
			swfS[1] = (sF[1] * dsw) * (1.0 - a1);
			swfS[2] = (sF[2] * dsw) * (1.0 - a2);

			// absorption coefficient for spectral range:
			h =xNew<IssmDouble>(m);
			B1 =xNew<IssmDouble>(m);
			B2 =xNew<IssmDouble>(m);
			for(int i=0;i<m;i++) h[i]= d[i] /(pow(gsz[i],0.5));
			for(int i=0;i<m;i++) B1[i] = 0.0192 * h[i];                 // 0.3 - 0.8um
			for(int i=0;i<m;i++) B2[i]= 0.1098 * h[i];                 // 0.8 - 1.5um
			// B3 = +inf                     // 1.5 - 2.8um

			// cumulative extinction factors
			exp1 = xNew<IssmDouble>(m); 
			exp2 = xNew<IssmDouble>(m); 
			for(int i=0;i<m;i++) exp1[i]=exp(-B1[i]*dz[i]);
			for(int i=0;i<m;i++) exp2[i]=exp(-B2[i]*dz[i]);

			for(int i=0;i<m;i++){
				IssmDouble cum1=exp1[0];
				IssmDouble cum2=exp2[0];
				for(int j=1;j<=i;j++){
					cum1 = cum1*exp1[j];
					cum2 = cum2*exp2[j];
				}
				B1_cum[i+1]=cum1;
				B2_cum[i+1]=cum2;
			}

			// flux across grid cell boundaries
			Qs1 = xNew<IssmDouble>(m+1);
			Qs2 = xNew<IssmDouble>(m+1);
			for(int i=0;i<m+1;i++){
				Qs1[i] = swfS[0] * B1_cum[i];
				Qs2[i] = swfS[1] * B2_cum[i];
			}

			// net energy flux to each grid cell
			for(int i=0;i<m;i++) swf[i]= (Qs1[i]-Qs1[i+1]) + (Qs2[i]-Qs2[i+1]);

			// add flux absorbed at surface
			swf[0] = swf[0]+ swfS[2];

			/*Free resources: */
			xDelete<IssmDouble>(gsz);
			xDelete<IssmDouble>(B1_cum);
			xDelete<IssmDouble>(B2_cum);
			xDelete<IssmDouble>(h);
			xDelete<IssmDouble>(B1);
			xDelete<IssmDouble>(B2);
			xDelete<IssmDouble>(exp1);
			xDelete<IssmDouble>(exp2);
			xDelete<IssmDouble>(Qs1);
			xDelete<IssmDouble>(Qs2);

		}
		else{  //function of grid cell density

			/*intermediary: */
			IssmDouble* B_cum = NULL;
			IssmDouble* exp_B = NULL;
			IssmDouble* Qs    = NULL;
			IssmDouble* B    = NULL;

			// fraction of sw radiation absorbed in top grid cell (wavelength > 0.8um)
			IssmDouble SWs = 0.36;

			// SWs and SWss coefficients need to be better constranted. Greuell
			// and Konzelmann 1994 used SWs = 0.36 and SWss = 0.64 as this the
			// the % of SW radiation with wavelengths > and < 800 nm
			// respectively.  This, however, may not account for the fact that
			// the albedo of wavelengths > 800 nm has a much lower albedo.

			// calculate surface shortwave radiation fluxes [W m-2]
			IssmDouble swf_s = SWs * (1.0 - as) * dsw;

			// calculate surface shortwave radiation fluxes [W m-2]
			IssmDouble swf_ss = (1.0-SWs) * (1.0 - as) * dsw;

			// SW allowed to penetrate into snowpack
			IssmDouble Bs = 10.0;    // snow SW extinction coefficient [m-1] (Bassford,2006)
			IssmDouble Bi = 1.3;   // ice SW extinction coefficient [m-1] (Bassford,2006)

			// calculate extinction coefficient B [m-1] vector
			B=xNew<IssmDouble>(m);
			for(int i=0;i<m;i++) B[i] = Bs + (300.0 - d[i]) * ((Bs - Bi)/(dIce - 300.0));

			// cumulative extinction factor
			B_cum = xNew<IssmDouble>(m+1);
			exp_B = xNew<IssmDouble>(m);
			for(int i=0;i<m;i++)exp_B[i]=exp(-B[i]*dz[i]);

			B_cum[0]=1.0;
			for(int i=0;i<m;i++){
				IssmDouble cum_B=exp_B[0];
				for(int j=1;j<=i;j++) cum_B=cum_B*exp_B[j];
				B_cum[i+1]=  cum_B;
			}

			// flux across grid cell boundaries
			Qs=xNew<IssmDouble>(m+1);
			for(int i=0;i<m+1;i++) Qs[i] = swf_ss * B_cum[i];

			// net energy flux to each grid cell
			for(int i=0;i<m;i++) swf[i] = (Qs[i]-Qs[i+1]);

			// add flux absorbed at surface
			swf[0] += swf_s;

			/*Free resources:*/
			xDelete<IssmDouble>(B_cum);
			xDelete<IssmDouble>(exp_B);
			xDelete<IssmDouble>(Qs);
			xDelete<IssmDouble>(B);
		}
	}
	/*Assign output pointers: */
	*pswf=swf;

} /*}}}*/ 
void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* pRa, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid){ /*{{{*/

	// Adds precipitation and deposition to the model grid

	// Author: Alex Gardner, University of Alberta
	// Date last modified: JAN, 2008

	/* Description:
	   adjusts the properties of the top grid cell to account for accumulation
	   T_air & T = Air and top grid cell temperatures [K]
		Tmean =  average surface temperature [K]
		Vmean =  average wind velocity [m s-1]
		V =  wind velocity [m s-1]
		C =  average accumulation rate [kg m-2 yr-1]
	   dz = topgrid cell length [m]
	   d = density of top grid gell [kg m-3]
		P = precipitation [mm w.e.] or [kg m-3]
		Ra = rainfall [mm w.e.] or [kg m-3]
		re = effective grain radius [mm]
		gdn = grain dentricity
		gsp = grain sphericity*/

	// MAIN FUNCTION
	// specify constants
	IssmDouble dSnow = 150;    // density of snow [kg m-3]
	IssmDouble reNew = 0.05;    // new snow grain size [mm]
	IssmDouble gdnNew = 1.0;     // new snow dendricity 
	IssmDouble gspNew = 0.5;   // new snow sphericity 

	/*intermediary: */
	IssmDouble* mInit=NULL;
	bool        top=true;
	IssmDouble  mass=0;
	IssmDouble  massinit=0;
	IssmDouble  mass_diff=0;

	/*output: */
	IssmDouble* T=NULL;
	IssmDouble* dz=NULL;
	IssmDouble* d=NULL;
	IssmDouble* W=NULL;
	IssmDouble* a=NULL;
	IssmDouble* adiff=NULL;
	IssmDouble* re=NULL;
	IssmDouble* gdn=NULL;
	IssmDouble* gsp=NULL;
	IssmDouble  ra=0.0;
	int         m=0;

	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   accumulation module\n");

	/*Recover pointers: */
	T=*pT;
	dz=*pdz;
	d=*pd;
	W=*pW;
	a=*pa;
	adiff=*padiff;
	re=*pre;
	gdn=*pgdn;
	gsp=*pgsp;
	m=*pm;
	ra=*pRa;

	//Density of fresh snow [kg m-3]
	switch (dsnowIdx){
		case 0: // Default value defined above
			break;

		case 1: // Density of Antarctica snow
			dSnow = 350;
			//dSnow = 360; //FirnMICE Lundin et al., 2017
			break;

		case 2: // Density of Greenland snow, Fausto et al., 2018
			dSnow = 315;
			//From Vionnet et al., 2012 (Crocus)
			gdnNew = min(max(1.29 - 0.17*V,0.20),1.0);
			gspNew = min(max(0.08*V + 0.38,0.5),0.9);
			reNew=max(1e-1*(gdnNew/.99+(1-1*gdnNew/.99)*(gspNew/.99*3+(1-gspNew/.99)*4))/2,Gdntol);
			break;

		case 3: //Surface snow accumulation density from Kaspers et al., 2004, Antarctica
			//dSnow = alpha1 + beta1*T + delta1*C + epsilon1*W
			//     7.36x10-2  1.06x10-3  6.69x10-2  4.77x10-3 
			dSnow=(7.36e-2 + 1.06e-3*min(Tmean,CtoK-Ttol) + 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000;
			break;

		case 4: // Kuipers Munneke and others (2015), Greenland
			dSnow = 481.0 + 4.834*(Tmean-CtoK);
			break;
	}

	// determine initial mass
	mInit=xNew<IssmDouble>(m);
	for(int i=0;i<m;i++) mInit[i]= d[i] * dz[i];
	massinit=0.0; 
	for(int i=0;i<m;i++)massinit+=mInit[i];

	if (P > 0.0+Ptol){

		if (T_air <= CtoK+Ttol){ // if snow

			IssmDouble  z_snow = P/dSnow;               // depth of snow
			IssmDouble dfall = gdnNew;
			IssmDouble sfall = gspNew;
			IssmDouble refall = reNew;

			// if snow depth is greater than specified min dz, new cell created
			if (z_snow > dzMin+Dtol){

				newcell(&T,T_air,top,m); //new cell T
				newcell(&dz,z_snow,top,m); //new cell dz
				newcell(&d,dSnow,top,m); //new cell d
				newcell(&W,0.0,top,m); //new cell W
				newcell(&a,aSnow,top,m); //new cell a
				newcell(&adiff,aSnow,top,m); //new cell a
				newcell(&re,refall,top,m); //new cell grain size
				newcell(&gdn,dfall,top,m); //new cell grain dendricity
				newcell(&gsp,sfall,top,m); //new cell grain sphericity
				m=m+1;
			}
			else { // if snow depth is less than specified minimum dz snow

				mass = mInit[0] + P;         // grid cell adjust mass

				dz[0] = dz[0] + P/dSnow;    // adjust grid cell depth      
				d[0] = mass / dz[0];    // adjust grid cell density

				// adjust variables as a linearly weighted function of mass
				// adjust temperature (assume P is same temp as air)
				T[0] = (T_air * P + T[0] * mInit[0])/mass;

				// adjust a, re, gdn & gsp
				if(aIdx>0)a[0] = (aSnow * P + a[0] * mInit[0])/mass;
				gdn[0] = dfall; 
				gsp[0] = sfall; 
				re[0] = max(1e-1*(gdn[0]/.99+(1-1*gdn[0]/.99)*(gsp[0]/.99*3+(1-gsp[0]/.99)*4))/2,Gdntol);
			}
		}
		else{ // if rain    

			/*rain is added by increasing the mass and temperature of the ice
			  of the top grid cell.  Temperatures are set artifically high to
			  account for the latent heat of fusion.  This is the same as
			  directly adding liquid water to the the snow pack surface but
			  makes the numerics easier.*/

			// grid cell adjust mass
			mass = mInit[0] + P;

			// adjust temperature
			// liquid: must account for latent heat of fusion
			T[0] = (P *(T_air + LF/CI) + T[0] * mInit[0]) / mass;

			// adjust grid cell density
			d[0] = mass / dz[0];

			// if d > the density of ice, d = dIce
			if (d[0] > dIce-Dtol){
				d[0] = dIce;           // adjust d
				dz[0] = mass / d[0];    // dz is adjusted to conserve mass
			}

			ra=P;
		}

		// check for conservation of mass
		mass=0; for(int i=0;i<m;i++)mass+=d[i]*dz[i]; 

		mass_diff = mass - massinit - P;

		#ifndef _HAVE_AD_  //avoid round operation. only check in forward mode.
		mass_diff = round(mass_diff * 100.0)/100.0;
		if (mass_diff > 0) _error_("mass not conserved in accumulation function");
		#endif

	}
	/*Free resources:*/
	if(mInit)xDelete<IssmDouble>(mInit);

	/*Assign output pointers:*/
	*pT=T;
	*pdz=dz;
	*pd=d;
	*pW=W;
	*pa=a;
	*padiff=adiff;
	*pre=re;
	*pgdn=gdn;
	*pgsp=gsp;
	*pm=m;
	*pRa=ra;
} /*}}}*/
void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pF, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble Ra, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid){ /*{{{*/

	//// MELT ROUTINE

	// Description:
	// computes the quantity of meltwater due to snow temperature in excess of
	// 0 deg C, determines pore water content and adjusts grid spacing

	/*intermediary:*/
	IssmDouble* m=NULL;
	IssmDouble* maxF=NULL;
	IssmDouble* dW=NULL;
	IssmDouble* exsW=NULL;
	IssmDouble* exsT=NULL;
	IssmDouble* surpT=NULL;
	IssmDouble* surpE=NULL;
	IssmDouble* flxDn=NULL;
	IssmDouble  ER=0.0;
	IssmDouble* EI=NULL;
	IssmDouble* EW=NULL;
	IssmDouble* M=NULL;
	int*        D=NULL;

	IssmDouble sumM=0.0;
	IssmDouble sumER=0.0;
	IssmDouble addE=0.0;
	IssmDouble mSum0=0.0;
	IssmDouble sumE0=0.0;
	IssmDouble mSum1=0.0;
	IssmDouble sumE1=0.0;
	IssmDouble dE=0.0;
	IssmDouble dm=0.0;
	IssmDouble X=0.0;
	IssmDouble Wi=0.0;

	int        D_size = 0;

	/*outputs:*/
	IssmDouble  Msurf = 0.0;
	IssmDouble  mAdd = 0.0;
	IssmDouble  surplusE = 0.0;
	IssmDouble  surplusT = 0.0;
	IssmDouble  dz_add = 0.0;
	IssmDouble  Rsum = 0.0;
	IssmDouble  Fsum = 0.0;
	IssmDouble* T=*pT;
	IssmDouble* d=*pd;
	IssmDouble* dz=*pdz;
	IssmDouble* W=*pW;
	IssmDouble* a=*pa;
	IssmDouble* adiff=*padiff;
	IssmDouble* re=*pre;
	IssmDouble* gdn=*pgdn;
	IssmDouble* gsp=*pgsp;
	int         n=*pn;
	IssmDouble* R=NULL;

	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   melt module\n");

	//// INITIALIZATION

	/*Allocations: */
	M=xNewZeroInit<IssmDouble>(n); 
	maxF=xNewZeroInit<IssmDouble>(n); 
	dW=xNewZeroInit<IssmDouble>(n); 

	// store initial mass [kg] and energy [J]
	m=xNew<IssmDouble>(n); for(int i=0;i<n;i++) m[i] = dz[i]* d[i];                    // grid cell mass [kg]
	EI=xNew<IssmDouble>(n); for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI;               // initial enegy of snow/ice
	EW=xNew<IssmDouble>(n); for(int i=0;i<n;i++)EW[i]= W[i] * (LF + CtoK * CI);     // initial enegy of water

	mSum0 = cellsum(W,n) + cellsum(m,n);        // total mass [kg]
	sumE0 = cellsum(EI,n) + cellsum(EW,n);      // total energy [J]

	// initialize melt and runoff scalars
	Rsum = 0.0;       // runoff [kg]
	Fsum = 0.0;       // refreeze [kg]
	sumM = 0.0;       // total melt [kg]
	mAdd = 0.0;       // mass added/removed to/from base of model [kg]
	addE = 0.0;       // energy added/removed to/from base of model [J]
	dz_add=0.0;      // thickness of the layer added/removed to/from base of model [m]

	// calculate temperature excess above 0 deg C
	exsT=xNewZeroInit<IssmDouble>(n);
	for(int i=0;i<n;i++) exsT[i]= max(0.0, T[i] - CtoK);        // [K] to [degC]

	// new grid point center temperature, T [K]
	// for(int i=0;i<n;i++) T[i]-=exsT[i];
	for(int i=0;i<n;i++) T[i]=min(T[i],CtoK);

	// specify irreducible water content saturation [fraction]
	const IssmDouble Swi = 0.07;                     // assumed constant after Colbeck, 1974
	const IssmDouble dPHC = 830;                     //Pore closeoff density

	//// REFREEZE PORE WATER
	// check if any pore water
	if (cellsum(W,n) >0.0+Wtol){
		if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("      pore water refreeze\n");
		// calculate maximum freeze amount, maxF [kg]
		for(int i=0;i<n;i++) maxF[i] = max(0.0, -((T[i] - CtoK) * m[i] * CI) / LF);

		// freeze pore water and change snow/ice properties
		for(int i=0;i<n;i++) dW[i] = min(maxF[i], W[i]);    // freeze mass [kg]   
		for(int i=0;i<n;i++) W[i] -= dW[i];                                            // pore water mass [kg]
		for(int i=0;i<n;i++) m[i] += dW[i];                                            // new mass [kg]
		for(int i=0;i<n;i++) d[i] = m[i] / dz[i];                                    // density [kg m-3]   
		for(int i=0;i<n;i++) if(m[i]>Wtol) T[i] = T[i] + (dW[i]*(LF+(CtoK - T[i])*CI)/(m[i]*CI));      // temperature [K]

		// if pore water froze in ice then adjust d and dz thickness
		for(int i=0;i<n;i++)if(d[i]> dIce-Dtol)d[i]=dIce;
		for(int i=0;i<n;i++) dz[i]= m[i]/d[i];

	}

	// squeeze water from snow pack
	exsW=xNew<IssmDouble>(n); 
	for(int i=0;i<n;i++){
		Wi= (dIce - d[i]) * Swi * (m[i] / d[i]);        // irreducible water content [kg]
		exsW[i] = max(0.0, W[i] - Wi);                  // water "squeezed" from snow [kg]
	}

	//// MELT, PERCOLATION AND REFREEZE

	// initialize refreeze, runoff, flxDn and dW vectors [kg]
	IssmDouble* F = xNewZeroInit<IssmDouble>(n);

	// Add previous refreeze to F and reset dW
	for(int i=0;i<n;i++){
		F[i]=F[i]+dW[i];
		dW[i] = 0.0;
	}

	// run melt algorithm if there is melt water or excess pore water
	if ((cellsum(exsT,n) > 0.0+Ttol) || (cellsum(exsW,n) > 0.0+Wtol)){
		// _printf_(""MELT OCCURS");
		// check to see if thermal energy exceeds energy to melt entire cell
		// if so redistribute temperature to lower cells (temperature surplus)
		// (maximum T of snow before entire grid cell melts is a constant
		// LF/CI = 159.1342)
		surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = max(0.0, exsT[i]- LF/CI);

		if (cellsum(surpT,n) > 0.0+Ttol ){
			// _printf_("T Surplus");
			// calculate surplus energy
			surpE=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpE[i] = surpT[i] * CI * m[i];

			int i = 0;
			while (cellsum(surpE,n) > 0.0+Ttol && i<n){

				if (i<n-1){
					// use surplus energy to increase the temperature of lower cell
					T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];

					exsT[i+1] = max(0.0, T[i+1] - CtoK) + exsT[i+1];
					T[i+1] = min(CtoK, T[i+1]);

					surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
					surpE[i+1] = surpT[i+1] * CI * m[i+1];
				}
				else{
					surplusT=max(0.0, exsT[i] - LF/CI);
					surplusE=surpE[i];
					if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0){
						_printf0_(" WARNING: surplus energy at the base of GEMB column\n");
					}
				}

				// adjust current cell properties (again 159.1342 is the max T)
				exsT[i] = LF/CI;
				surpE[i] = 0.0;
				i = i + 1;
			}
		}

		// convert temperature excess to melt [kg]
		IssmDouble Mmax=0.0;
		for(int i=0;i<n;i++){
			Mmax=exsT[i] * d[i] * dz[i] * CI / LF;
			M[i] = min(Mmax, m[i]);  // melt
		}
		Msurf = M[0];
		sumM = max(0.0,cellsum(M,n)-Ra);  // total melt [kg] minus the liquid rain that had been added

		// calculate maximum refreeze amount, maxF [kg]
		for(int i=0;i<n;i++)maxF[i] = max(0.0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);

		// initialize refreeze, runoff, flxDn and dW vectors [kg]
		IssmDouble* R = xNewZeroInit<IssmDouble>(n);

		flxDn=xNewZeroInit<IssmDouble>(n+1);

		// determine the deepest grid cell where melt/pore water is generated
		X = 0;
		for(int i=n-1;i>=0;i--){
			if(M[i]> 0.0+Wtol || exsW[i]> 0.0+Wtol){
				X=i;
				break;
			}
		}

		IssmDouble depthice=0.0;
		int Xi=0;
		//// meltwater percolation
		for(int i=0;i<n;i++){
			// calculate total melt water entering cell
			IssmDouble inM = M[i]+ flxDn[i];

			depthice=0.0;
			if (d[i] >= dPHC-Dtol){
				for(int l=i;(l<n && d[l]>=dPHC-Dtol);l++) depthice=depthice+dz[l];
			}

			// break loop if there is no meltwater and if depth is > mw_depth
			if (fabs(inM) < Wtol && i > X){
				break;
			}
			// if reaches impermeable ice layer all liquid water runs off (R)
			else if (d[i] >= dIce-Dtol || (d[i] >= dPHC-Dtol && depthice>0.1)){  // dPHC = pore hole close off [kg m-3]
				// _printf_("ICE LAYER");
				// no water freezes in this cell
				// no water percolates to lower cell
				// cell ice temperature & density do not change

				m[i] = m[i] - M[i];                     // mass after melt
				Wi = (dIce-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
				dW[i] = max(min(inM, Wi - W[i]),-1*W[i]);            // change in pore water
				R[i] = max(0.0, inM - dW[i]);             // runoff
			}
			// check if no energy to refreeze meltwater
			else if (fabs(maxF[i]) < Dtol){
				// _printf_("REFREEZE == 0");
				// no water freezes in this cell
				// cell ice temperature & density do not change

				m[i] = m[i] - M[i];                     // mass after melt
				Wi = (dIce-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
				dW[i] = max(min(inM, Wi - W[i]),-1*W[i]);              // change in pore water
				flxDn[i+1] = max(0.0, inM - dW[i]);         // meltwater out
				R[i] = 0.0;
			}
			// some or all meltwater refreezes
			else{
				// change in density and temperature
				// _printf_("MELT REFREEZE");
				//-----------------------melt water-----------------------------
				m[i] = m[i] - M[i];
				IssmDouble dz_0 = m[i]/d[i];          
				IssmDouble dMax = (dIce - d[i])*dz_0;              // d max = dIce
				IssmDouble F1 = min(min(inM,dMax),maxF[i]);         // maximum refreeze               
				m[i] = m[i] + F1;                       // mass after refreeze
				d[i] = m[i]/dz_0;

				//-----------------------pore water-----------------------------
				Wi = (dIce-d[i])* Swi * dz_0;            // irreducible water 
				dW[i] = max(min(inM - F1, Wi-W[i]),-1*W[i]);         // change in pore water
				IssmDouble F2 = 0.0;                                 

				if (dW[i] < 0.0-Wtol){                         // excess pore water
					dMax = (dIce - d[i])*dz_0;          // maximum refreeze                                             
					IssmDouble maxF2 = min(dMax, maxF[i]-F1);      // maximum refreeze
					F2 = min(-1*dW[i], maxF2);            // pore water refreeze
					m[i] = m[i] + F2;                   // mass after refreeze
					d[i] = m[i]/dz_0;
				}

				F[i] = F[i] + F1 + F2;

				flxDn[i+1] = max(0.0,inM - F1 - dW[i]);     // meltwater out        
				if (m[i]>Wtol){
					T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
				}

				// check if an ice layer forms 
				if (fabs(d[i] - dIce) < Dtol){
					// _printf_("ICE LAYER FORMS");
					// excess water runs off
					R[i] = flxDn[i+1];
					// no water percolates to lower cell
					flxDn[i+1] = 0.0;
				}
			}
			Xi=Xi+1;
		}

		//// GRID CELL SPACING AND MODEL DEPTH
		for(int i=0;i<n;i++)if (W[i] < 0.0-Wtol) _error_("negative pore water generated in melt equations");

		// delete all cells with zero mass
		// adjust pore water
		for(int i=0;i<n;i++)W[i] += dW[i];

		//calculate Rsum:
		Rsum=cellsum(R,n) + flxDn[Xi];

		// delete all cells with zero mass
		D_size=0; for(int i=0;i<n;i++)if(m[i]> (0.0+Wtol))D_size++; 
		D=xNew<int>(D_size);
		D_size=0; for(int i=0;i<n;i++)if(m[i]> (0.0+Wtol)){ D[D_size] = i; D_size++;}

		celldelete(&m,n,D,D_size);
		celldelete(&W,n,D,D_size);
		celldelete(&d,n,D,D_size);
		celldelete(&dz,n,D,D_size);
		celldelete(&T,n,D,D_size);
		celldelete(&a,n,D,D_size);
		celldelete(&adiff,n,D,D_size);
		celldelete(&re,n,D,D_size);
		celldelete(&gdn,n,D,D_size);
		celldelete(&gsp,n,D,D_size);
		celldelete(&EI,n,D,D_size);
		celldelete(&EW,n,D,D_size);
		n=D_size;
		xDelete<int>(D);

		// calculate new grid lengths
		for(int i=0;i<n;i++)dz[i] = m[i] / d[i];

		/*Free resources:*/
		xDelete<IssmDouble>(R);
	}

	//calculate Fsum:
	Fsum=cellsum(F,n);

	/*Free resources:*/
	xDelete<IssmDouble>(F);

	//Manage the layering to match the user defined requirements
	managelayers(&mAdd, &dz_add, &addE, &m, &EI, &EW, &T, &d, &dz, &W, &a, &adiff, &re, &gdn, &gsp, &n, dzMin, zMax, zMin, zTop, zY);

	//// CHECK FOR MASS AND ENERGY CONSERVATION

	// calculate final mass [kg] and energy [J]
	sumER = Rsum * (LF + CtoK * CI);
	for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI;
	for(int i=0;i<n;i++)EW[i] = W[i] * (LF + CtoK * CI);

	mSum1 = cellsum(W,n) + cellsum(m,n) + Rsum;
	sumE1 = cellsum(EI,n) + cellsum(EW,n);

	/*checks: */
	for(int i=0;i<n;i++) if (W[i]<0.0-Wtol) _error_("negative pore water generated in melt equations\n");

	/*only in forward mode! avoid round in AD mode as it is not differentiable: */
	#ifndef _HAVE_AD_
	dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
	dE = round(sumE0 - sumE1 - sumER +  addE - surplusE);
	if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
			<< "dm: " << dm << " dE: " << dE << "\n");
	#endif

	/*Free resources:*/
	if(m)xDelete<IssmDouble>(m);
	if(EI)xDelete<IssmDouble>(EI);
	if(EW)xDelete<IssmDouble>(EW);
	if(maxF)xDelete<IssmDouble>(maxF);
	if(dW)xDelete<IssmDouble>(dW);
	if(exsW)xDelete<IssmDouble>(exsW);
	if(exsT)xDelete<IssmDouble>(exsT);
	if(surpT)xDelete<IssmDouble>(surpT);
	if(surpE)xDelete<IssmDouble>(surpE);
	if(flxDn)xDelete<IssmDouble>(flxDn);
	if(D)xDelete<int>(D);
	if(M)xDelete<IssmDouble>(M);

	/*Assign output pointers:*/
	*pMs=Msurf;
	*pM=sumM;
	*pR=Rsum;
	*pF=Fsum;
	*pmAdd=mAdd;
	*pdz_add=dz_add;

	*pT=T;
	*pd=d;
	*pdz=dz;
	*pW=W;
	*pa=a;
	*padiff=adiff;
	*pre=re;
	*pgdn=gdn;
	*pgsp=gsp;
	*pn=n;

} /*}}}*/ 
void managelayers(IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble* paddE, IssmDouble** pm, IssmDouble** pEI, IssmDouble** pEW, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY){ /*{{{*/

	/*intermediary:*/
	IssmDouble* Zcum=NULL;
	IssmDouble* dzMin2=NULL;
	IssmDouble* dzMax2=NULL;
	int*        D=NULL;

	IssmDouble zY2=zY;
	IssmDouble X=0.0;
	int X1=0;
	int X2=0;
	int D_size = 0;

	IssmDouble Ztot=0.0;
	IssmDouble T_bot=0.0;
	IssmDouble m_bot=0.0;
	IssmDouble dz_bot=0.0;
	IssmDouble d_bot=0.0;
	IssmDouble W_bot=0.0;
	IssmDouble a_bot=0.0;
	IssmDouble adiff_bot=0.0;
	IssmDouble re_bot=0.0;
	IssmDouble gdn_bot=0.0;
	IssmDouble gsp_bot=0.0;
	IssmDouble EI_bot=0.0;
	IssmDouble EW_bot=0.0;
	bool       top=false;

	/*outputs:*/
	IssmDouble  mAdd = 0.0;
	IssmDouble  addE = 0.0;
	IssmDouble  dz_add = 0.0;
	IssmDouble* T=*pT;
	IssmDouble* d=*pd;
	IssmDouble* dz=*pdz;
	IssmDouble* W=*pW;
	IssmDouble* a=*pa;
	IssmDouble* adiff=*padiff;
	IssmDouble* re=*pre;
	IssmDouble* gdn=*pgdn;
	IssmDouble* gsp=*pgsp;
	IssmDouble* m=*pm;
	IssmDouble* EI=*pEI;
	IssmDouble* EW=*pEW;
	int         n=*pn;

	//Merging of cells as they are burried under snow.
	Zcum=xNew<IssmDouble>(n);
	dzMin2=xNew<IssmDouble>(n);
	dzMax2=xNew<IssmDouble>(n);

	X=0;
	Zcum[0]=dz[0]; // Compute a cumulative depth vector
	for (int i=0;i<n;i++){
		if (i==0){
			dzMin2[i]=dzMin;
		}
		else{
			Zcum[i]=Zcum[i-1]+dz[i];
			if (Zcum[i]<=zTop+Dtol){
				dzMin2[i]=dzMin;
				X=i;
			}
			else{
				dzMin2[i]=zY2*dzMin2[i-1];
			}
		}
	}

	// Check to see if any cells are too small and need to be merged
	for (int i=0; i<n; i++){
		if ( (i<=X && dz[i]<dzMin-Dtol) || (i>X && dz[i]<dzMin2[i]-Dtol) ) {

			if (i==n-1){
				X2=i;
				//find closest cell to merge with
				for(int j=n-2;j>=0;j--){
					if(m[j]!=Delflag){
						X1=j;
						break;
					}
				}
			}
			else{
				X1=i+1;
				X2=i;
			}

			// adjust variables as a linearly weighted function of mass
			IssmDouble m_new = m[X2] + m[X1];
			T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
			a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new;
			adiff[X1] = (adiff[X2]*m[X2] + adiff[X1]*m[X1]) / m_new;
         //use grain properties from lower cell
			re[X1] = re[X2]; 
			gdn[X1] = gdn[X2]; 
			gsp[X1] = gsp[X2]; 

			// merge with underlying grid cell and delete old cell
			dz[X1] = dz[X2] + dz[X1];                 // combine cell depths
			d[X1] = m_new / dz[X1];                   // combine top densities
			W[X1] = W[X1] + W[X2];                     // combine liquid water
			m[X1] = m_new;                             // combine top masses

			// set cell to -99999 for deletion
			m[X2] = Delflag;
		}
	}

	// delete combined cells
	D_size=0; for(int i=0;i<n;i++)if(m[i]> Delflag+Wtol)D_size++; 
	D=xNew<int>(D_size); 
	D_size=0; for(int i=0;i<n;i++)if(m[i]> Delflag+Wtol){ D[D_size] = i; D_size++;}

	celldelete(&m,n,D,D_size);
	celldelete(&W,n,D,D_size);
	celldelete(&dz,n,D,D_size);
	celldelete(&d,n,D,D_size);
	celldelete(&T,n,D,D_size);
	celldelete(&a,n,D,D_size);
	celldelete(&adiff,n,D,D_size);
	celldelete(&re,n,D,D_size);
	celldelete(&gdn,n,D,D_size);
	celldelete(&gsp,n,D,D_size);
	celldelete(&EI,n,D,D_size);
	celldelete(&EW,n,D,D_size);
	n=D_size;
	xDelete<int>(D);

	// check if any of the cell depths are too large
	X=0;
	Zcum[0]=dz[0]; // Compute a cumulative depth vector
	for (int i=0;i<n;i++){
		if (i==0){
			dzMax2[i]=dzMin*2;
		}
		else{
			Zcum[i]=Zcum[i-1]+dz[i];
			if (Zcum[i]<=zTop+Dtol){
				dzMax2[i]=dzMin*2;
				X=i;
			}
			else{
				dzMax2[i]=max(zY2*dzMin2[i-1],dzMin*2);
			}
		}
	}

	for (int j=n-1;j>=0;j--){
		if ((j<X && dz[j] > dzMax2[j]+Dtol) || (dz[j] > dzMax2[j]*zY2+Dtol)){
			// _printf_("dz > dzMin * 2");
			// split in two
			cellsplit(&dz, n, j,.5);
			cellsplit(&W, n, j,.5);
			cellsplit(&m, n, j,.5);
			cellsplit(&T, n, j,1.0);
			cellsplit(&d, n, j,1.0);
			cellsplit(&a, n, j,1.0);
			cellsplit(&adiff, n, j,1.0);
			cellsplit(&EI, n, j,.5);
			cellsplit(&EW, n, j,.5);
			cellsplit(&re, n, j,1.0);
			cellsplit(&gdn, n, j,1.0);
			cellsplit(&gsp, n, j,1.0);
			n++;
		}
	}

	//// CORRECT FOR TOTAL MODEL DEPTH
	// WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT
	// INTERPRETATION
	// calculate total model depth
	Ztot = cellsum(dz,n);

	if (Ztot < zMin-Dtol){
		// printf("Total depth < zMin %f \n", Ztot);
		// mass and energy to be added
		mAdd = m[n-1]+W[n-1];
		addE = T[n-1]*m[n-1]*CI + W[n-1]*(LF+CtoK*CI);

		// add a grid cell of the same size and temperature to the bottom
		dz_bot=dz[n-1];
		T_bot=T[n-1];
		W_bot=W[n-1];
		m_bot=m[n-1];
		d_bot=d[n-1];
		a_bot=a[n-1];
		adiff_bot=adiff[n-1];
		re_bot=re[n-1];
		gdn_bot=gdn[n-1];
		gsp_bot=gsp[n-1];
		EI_bot=EI[n-1];
		EW_bot=EW[n-1];

		dz_add=dz_bot;

		newcell(&dz,dz_bot,top,n);
		newcell(&T,T_bot,top,n);
		newcell(&W,W_bot,top,n);
		newcell(&m,m_bot,top,n);
		newcell(&d,d_bot,top,n);
		newcell(&a,a_bot,top,n);
		newcell(&adiff,adiff_bot,top,n);
		newcell(&re,re_bot,top,n);
		newcell(&gdn,gdn_bot,top,n);
		newcell(&gsp,gsp_bot,top,n);
		newcell(&EI,EI_bot,top,n);
		newcell(&EW,EW_bot,top,n);
		n=n+1;
	}
	else if (Ztot > zMax+Dtol){
		// printf("Total depth > zMax %f \n", Ztot);
		// mass and energy loss
		mAdd = -(m[n-1]+W[n-1]);
		addE = -(T[n-1]*m[n-1]*CI) - (W[n-1]*(LF+CtoK*CI));
		dz_add=-(dz[n-1]);

		// remove a grid cell from the bottom
		D_size=n-1;
		D=xNew<int>(D_size);

		for(int i=0;i<n-1;i++) D[i]=i;
		celldelete(&dz,n,D,D_size);
		celldelete(&T,n,D,D_size);
		celldelete(&W,n,D,D_size);
		celldelete(&m,n,D,D_size);
		celldelete(&d,n,D,D_size);
		celldelete(&a,n,D,D_size);
		celldelete(&adiff,n,D,D_size);
		celldelete(&re,n,D,D_size);
		celldelete(&gdn,n,D,D_size);
		celldelete(&gsp,n,D,D_size);
		celldelete(&EI,n,D,D_size);
		celldelete(&EW,n,D,D_size);
		n=D_size;
		xDelete<int>(D);
	}

	/*Free resources:*/
 	xDelete<IssmDouble>(Zcum);
	xDelete<IssmDouble>(dzMin2);
	xDelete<IssmDouble>(dzMax2);
	if(D)xDelete<int>(D);

	/*Assign output pointers:*/
	*pT=T;
	*pd=d;
	*pdz=dz;
	*pW=W;
	*pa=a;
	*padiff=adiff;
	*pre=re;
	*pgdn=gdn;
	*pgsp=gsp;
	*pn=n;
	*pm=m;
	*pEI=EI;
	*pEW=EW;

	*pmAdd=mAdd;
	*paddE=addE;
	*pdz_add=dz_add;

} /*}}}*/ 
void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, int aIdx, int swIdx, IssmDouble adThresh, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid){ /*{{{*/

	//// THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEAMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPACTION IS COMPENSATED FOR BY TRACES OF SNOW???]

	//// FUNCTION INFO

	// Author: Alex Gardner, University of Alberta
	// Date last modified: FEB, 2008 

	// Description: 
	//   computes the densification of snow/firn using the emperical model of
	//   Herron and Langway (1980) or the semi-emperical model of Anthern et al.
	//   (2010)

	// Inputs:
	//   denIdx = densification model to use:
	//       1 = emperical model of Herron and Langway (1980)
	//       2 = semi-imerical model of Anthern et al. (2010)
	//       3 = physical model from Appendix B of Anthern et al. (2010)
	//   d   = initial snow/firn density [kg m-3]
	//   T   = temperature [K]
	//   dz  = grid cell size [m]
	//   C   = average accumulation rate [kg m-2 yr-1]
	//   dt  = time lapsed [s]
	//   re  = effective grain radius [mm];
	//   Ta  = mean annual temperature                                          

	// Reference: 
	// Herron and Langway (1980), Anthern et al. (2010)

	//// FOR TESTING
	// denIdx = 2;
	// d = 800;
	// T = 270;
	// dz = 0.005;
	// C = 200;
	// dt = 60*60;
	// re = 0.7;
	// Tmean = 273.15-18;

	//// MAIN FUNCTION
	// specify constants
	dt      = dt / dts;  // convert from [s] to [d]
	// R     = 8.314        // gas constant [mol-1 K-1]
	// Ec    = 60           // activation energy for self-diffusion of water
	//                      // molecules through the ice tattice [kJ mol-1]
	// Eg    = 42.4         // activation energy for grain growth [kJ mol-1]

	/*intermediary: */
	IssmDouble c0=0.0;
	IssmDouble c1=0.0;
   IssmDouble c2=0.0;
	IssmDouble H=0.0;
	IssmDouble M0=0.0;
	IssmDouble M1=0.0;
   IssmDouble M2=0.0;
	IssmDouble c0arth=0.0;
	IssmDouble c1arth=0.0;

	/*output: */
	IssmDouble* dz=NULL;
	IssmDouble* d=NULL;

	if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   densification module\n");

	/*Recover pointers: */
	dz=*pdz;
	d=*pd;

	// initial mass
	IssmDouble* mass_init = xNew<IssmDouble>(m);for(int i=0;i<m;i++) mass_init[i]=d[i] * dz[i];

	/*allocations and initialization of overburden pressure and factor H: */
	IssmDouble* cumdz = xNew<IssmDouble>(m-1);
	cumdz[0]=dz[0];
	for(int i=1;i<m-1;i++)cumdz[i]=cumdz[i-1]+dz[i];

	IssmDouble* obp = xNew<IssmDouble>(m);
	obp[0]=0.0;
	for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];

	// calculate new snow/firn density for:
	//   snow with densities <= 550 [kg m-3]
	//   snow with densities > 550 [kg m-3]

	for(int i=0;i<m;i++){
		switch (denIdx){
			case 1: // Herron and Langway (1980)
				c0 = (11.0 * exp(-10160.0 / (T[i] * R))) * C/1000.0;
				c1 = (575.0 * exp(-21400.0 / (T[i]* R))) * pow(C/1000.0,.5);
				break;

			case 2: // Arthern et al. (2010) [semi-emperical]
				// common variable
				// NOTE: Ec=60000, Eg=42400 (i.e. should be in J not kJ)
				H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
				c0 = 0.07 * H;
				c1 = 0.03 * H;
				break;

			case 3: // Arthern et al. (2010) [physical model eqn. B1]

				// common variable
				H = exp((-60000.0/(T[i] * R))) * obp[i] / pow(re[i]/1000.0,2.0);
				c0 = 9.2e-9 * H;
				c1 = 3.7e-9 * H;
				break;

			case 4: // Li and Zwally (2004)
				c0 = (C/dIce) * (139.21 - 0.542*Tmean)*8.36*pow(CtoK - T[i],-2.061);
				c1 = c0;
				break;

			case 5: // Helsen et al. (2008)
				// common variable
				c0 = (C/dIce) * (76.138 - 0.28965*Tmean)*8.36*pow(CtoK - T[i],-2.061);
				c1 = c0;
				break;

			case 6: // Ligtenberg and others (2011) [semi-emperical], Antarctica 
				// common variable
				// From literature: H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
				H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
				c0arth = 0.07 * H;
				c1arth = 0.03 * H;
				//ERA-5 old
				//M0 = max(2.3128 - (0.2480 * log(C)),0.25);
				//M1 = max(2.7950 - (0.3318 * log(C)),0.25);
				// ERA5 new aIdx=1, swIdx=0
				if (aIdx==1 && swIdx==0){
					if (fabs(adThresh - 820) < Dtol){
                  // ERA5 v4
                  M0 = max(1.5131 - (0.1317 * log(C)),0.25);
                  M1 = max(1.8819 - (0.2158 * log(C)),0.25);
					}
					else{
						// ERA5 new aIdx=1, swIdx=0
						//M0 = max(1.8785 - (0.1811 * log(C)),0.25);
						//M1 = max(2.0005 - (0.2346 * log(C)),0.25);
						// ERA5 new aIdx=1, swIdx=0, bare ice
						M0 = max(1.8422 - (0.1688 * log(C)),0.25);
						M1 = max(2.4979 - (0.3225 * log(C)),0.25);
					}
				}
				// ERA5 new aIdx=2, swIdx=1
				else if (aIdx<3 && swIdx>0){
					M0 = max(2.2191 - (0.2301 * log(C)),0.25);
					M1 = max(2.2917 - (0.2710 * log(C)),0.25);
				}
				// ERA5 new aIdx=2, swIdx=0
				//else if (aIdx==2){
				//}
				//From Ligtenberg
				//H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
				//M0 = max(1.435 - (0.151 * log(C)),0.25);
				//M1 = max(2.366 - (0.293 * log(C)),0.25);
            //RACMO
            M0 = max(1.6383 - (0.1691 * log(C)),0.25);
            M1 = max(1.9991 - (0.2414 * log(C)),0.25);
				c0 = M0*c0arth;
				c1 = M1*c1arth;
				c2 = M2*c1arth;
				break;

			case 7: // Kuipers Munneke and others (2015) [semi-emperical], Greenland
				// common variable
				// From literature: H = exp((-60000.0/(T[i] * R)) + (42400.0/(T[i] * R))) * (C * 9.81);
				H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
				c0arth = 0.07 * H;
				c1arth = 0.03 * H;
				// ERA5 old
				//M0 = max(1.8554 - (0.1316 * log(C)),0.25);
				//M1 = max(2.8901 - (0.3014 * log(C)),0.25);
				// ERA5 new aIdx=1, swIdx=0
				if (aIdx==1 && swIdx==0){
					if (fabs(adThresh - 820) < Dtol){
						// ERA5 v4
						M0 = max(1.3566 - (0.1350 * log(C)),0.25);
						M1 = max(1.8705 - (0.2290 * log(C)),0.25);
					}
					else{
						// ERA5 new aIdx=1, swIdx=0
						//M0 = max(1.4574 - (0.1123 * log(C)),0.25);
						//M1 = max(2.0238 - (0.2070 * log(C)),0.25);
						// ERA5 new aIdx=1, swIdx=0, bare ice
						M0 = max(1.4318 - (0.1055 * log(C)),0.25);
						M1 = max(2.0453 - (0.2137 * log(C)),0.25);
					}
				}
				// ERA5 new aIdx=2, swIdx=1
				else if (aIdx<3 && swIdx>0){
					M0 = max(1.7834 - (0.1409 * log(C)),0.25);
					M1 = max(1.9260 - (0.1527 * log(C)),0.25);
				}
				// ERA5 new aIdx=2, swIdx=0
				//else if (aIdx==2){
				//}
				// From Kuipers Munneke
				//M0 = max(1.042 - (0.0916 * log(C)),0.25);
				//M1 = max(1.734 - (0.2039 * log(C)),0.25);
				// RACMO
            M0 = max(1.2691 - (0.1184 * log(C)),0.25);
            M1 = max(1.9983 - (0.2511 * log(C)),0.25);
            c0 = M0*c0arth;
            c1 = M1*c1arth;
            c2 = M2*c1arth;
				break;
		}

      // new snow density
      if(d[i] <= 550.0+Dtol) d[i] = d[i] + (c0 * (dIce - d[i]) / 365.0 * dt);
      else if(d[i] <= 830.0+Dtol | fabs(c2)<Dtol) d[i] = d[i] + (c1 * (dIce - d[i]) / 365.0 * dt);
      else d[i] = d[i] + (c2 * (dIce - d[i]) / 365.0 * dt);

		// do not allow densities to exceed the density of ice
		if(d[i] > dIce-Ptol) d[i]=dIce;

		// calculate new grid cell length
		dz[i] = mass_init[i] / d[i];
	}

	/*Free resources:*/
	xDelete<IssmDouble>(mass_init);
	xDelete<IssmDouble>(cumdz);
	xDelete<IssmDouble>(obp);

	/*Assign output pointers:*/
	*pdz=dz;
	*pd=d;

} /*}}}*/