Changeset 22249

Timestamp:

11/09/17 16:58:40 (7 years ago)

Author:

schlegel

Message:

CHG: update Gemb code to use tolerances on all absolute comparisons, make stricter archive tolerances, Smb output is now in m/yr ice so that it can be run with issm mass transport on

Location:

issm/trunk-jpl

Files:

• src/c/classes/Elements/Element.cpp (modified) (9 diffs)
• src/c/modules/SurfaceMassBalancex/Gembx.cpp (modified) (104 diffs)
• src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (2 diffs)
• src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (modified) (1 diff)
• test/Archives/Archive243.arch (modified) ( previous)
• test/Archives/Archive244.arch (modified) ( previous)
• test/NightlyRun/test243.m (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -2521 +2521 @@
 
         /*Intermediary variables: {{{*/
-        IssmDouble isinitialized;
-        IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin;
-        IssmDouble Tmean; 
-        IssmDouble C; 
-        IssmDouble Tz,Vz; 
+        IssmDouble isinitialized=0.0;
+        IssmDouble zTop=0.0;
+        IssmDouble dzTop=0.0;
+        IssmDouble zMax=0.0;
+        IssmDouble zMin=0.0;
+        IssmDouble zY=0.0;
+        IssmDouble dzMin=0.0;
+        IssmDouble Tmean=0.0;
+        IssmDouble C=0.0;
+        IssmDouble Tz,Vz=0.0;
         IssmDouble rho_ice, rho_water,aSnow,aIce;
         IssmDouble time,dt;
         IssmDouble t,smb_dt;
         IssmDouble yts;
-        IssmDouble Ta,V,dlw,dsw,P,eAir,pAir;
+        IssmDouble Ta=0.0;
+        IssmDouble V=0.0;
+        IssmDouble dlw=0.0;
+        IssmDouble dsw=0.0;
+        IssmDouble P=0.0;
+        IssmDouble eAir=0.0;
+        IssmDouble pAir=0.0;
         int        aIdx=0;
         int        denIdx=0;
         int        swIdx=0;
         IssmDouble cldFrac,t0wet, t0dry, K;
-        IssmDouble ulw;
-        IssmDouble netSW;
-        IssmDouble netLW;
-        IssmDouble lhf, shf, dayEC;
-        IssmDouble initMass;
-        IssmDouble sumR, sumM, sumEC, sumP, sumW,sumMassAdd;
-        IssmDouble sumdz_add;
-        IssmDouble sumMass,dMass;
+        IssmDouble ulw=0.0;
+        IssmDouble netSW=0.0;
+        IssmDouble netLW=0.0;
+        IssmDouble lhf=0.0;
+        IssmDouble shf=0.0;
+        IssmDouble dayEC=0.0;
+        IssmDouble initMass=0.0;
+        IssmDouble sumR=0.0;
+        IssmDouble sumM=0.0;
+        IssmDouble sumEC=0.0;
+        IssmDouble sumP=0.0;
+        IssmDouble sumW=0.0;
+        IssmDouble sumMassAdd=0.0;
+        IssmDouble sumdz_add=0.0;
+        IssmDouble sumMass=0.0;
+        IssmDouble dMass=0.0;
         bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux;
-        IssmDouble init_scaling;
+        IssmDouble init_scaling=0.0;
 
         /*}}}*/
@@ -2553 +2572 @@
         IssmDouble* gdn = NULL;
         IssmDouble* gsp = NULL;
-        IssmDouble  EC = 0;
+        IssmDouble  EC = 0.0;
         IssmDouble* W = NULL;
         IssmDouble* a = NULL;
         IssmDouble* swf=NULL;
         IssmDouble* T = NULL;
-        IssmDouble  T_bottom;
-        IssmDouble  M;
-        IssmDouble  R; 
-        IssmDouble  mAdd;
-        IssmDouble  dz_add;
-    
+        IssmDouble  T_bottom = 0.0;
+        IssmDouble  M = 0.0;
+        IssmDouble  R = 0.0;
+        IssmDouble  mAdd = 0.0;
+        IssmDouble  dz_add = 0.0;
+
         IssmDouble* dzini=NULL;
         IssmDouble* dini = NULL;
@@ -2572 +2591 @@
         IssmDouble* aini = NULL;
         IssmDouble* Tini = NULL;
-    
-        int         m;
+
+        int         m=0;
         int         count=0;
         /*}}}*/
@@ -2606 +2625 @@
         parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum);
         parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum);
-    
+
         /*}}}*/
         /*Retrieve inputs: {{{*/
@@ -2645 +2664 @@
         /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/
         if(isinitialized==0.0){
-        if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n");
-        //if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" <<
-        //dz[i] << "\n");
-        
-        DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDziniEnum)); _assert_(dz_input);
-        DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDiniEnum));_assert_(d_input);
-        DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReiniEnum));_assert_(re_input);
-        DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdniniEnum));_assert_(gdn_input);
-        DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspiniEnum));_assert_(gsp_input);
-        DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECiniEnum));_assert_(EC_input);
-        DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWiniEnum));_assert_(W_input);
-        DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAiniEnum));_assert_(a_input);
-        DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTiniEnum));_assert_(T_input);
-
-        dz_input->GetValues(&dzini,&m);
-        d_input->GetValues(&dini,&m);
-        re_input->GetValues(&reini,&m);
-        gdn_input->GetValues(&gdnini,&m);
-        gsp_input->GetValues(&gspini,&m);
-        EC_input->GetInputValue(&EC);
-        W_input->GetValues(&Wini,&m);
-        a_input->GetValues(&aini,&m);
-        T_input->GetValues(&Tini,&m);
-        
-        /*Retrive the correct value of m (without the zeroes at the end)*/
-        Input* Size_input=this->GetInput(SmbSizeiniEnum); _assert_(Size_input);
-        Size_input->GetInputValue(&m);
-        
-        if(m==2){ //Snow properties are initialized with default values. Vertical grid has to be initialized too
-//            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w DEFAULT values\n");
-            
-            /*initialize profile variables:*/
-            GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);
-            
-            d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=dini[0]; //ice density [kg m-3]
-            re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=reini[0];         //set grain size to old snow [mm]
-            gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=gdnini[0];         //set grain dentricity to old snow
-            gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=gspini[0];         //set grain sphericity to old snow
-            W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0];             //set water content to zero [kg m-2]
-            a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0];         //set albedo equal to fresh snow [fraction]
-            T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
-/*            /!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties) because don't know Tmean yet when set default values.
-            Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/
-            
-            //fixed lower temperature bounday condition - T is fixed
-            T_bottom=T[m-1];
-        }
-        else{ //Retrieve snow properties from previous run. Need to provide values for all layers
-//            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n");
-            
-            dz = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzini[i];
-            d = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=dini[i];
-            re = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=reini[i];
-            gdn = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnini[i];
-            gsp = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gspini[i];
-            W = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wini[i];
-            a = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=aini[i];
-            T = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Tini[i];
-
-            //fixed lower temperature bounday condition - T is fixed
-            T_bottom=T[m-1];
-        }
-        
-        /*Flag the initialization:*/
-        this->AddInput(new DoubleInput(SmbIsInitializedEnum,1.0));
-    }
-    else{
-        /*Recover inputs: */
-        DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum)); _assert_(dz_input);
-        DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));_assert_(d_input);
-        DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));_assert_(re_input);
-        DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));_assert_(gdn_input);
-        DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));_assert_(gsp_input);
-        DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));_assert_(EC_input);
-        DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));_assert_(W_input);
-        DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));_assert_(a_input);
-        DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));_assert_(T_input);
-
-        /*Recover arrays: */
-        dz_input->GetValues(&dz,&m);
-        d_input->GetValues(&d,&m);
-        re_input->GetValues(&re,&m);
-        gdn_input->GetValues(&gdn,&m);
-        gsp_input->GetValues(&gsp,&m);
-        EC_input->GetInputValue(&EC);
-        W_input->GetValues(&W,&m);
-        a_input->GetValues(&a,&m);
-        T_input->GetValues(&T,&m);
-
-        //fixed lower temperature bounday condition - T is fixed
-        T_bottom=T[m-1];
-
-    } /*}}}*/
+                if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n");
+                //if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" <<
+                //dz[i] << "\n");
+
+                DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDziniEnum)); _assert_(dz_input);
+                DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDiniEnum));_assert_(d_input);
+                DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReiniEnum));_assert_(re_input);
+                DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdniniEnum));_assert_(gdn_input);
+                DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspiniEnum));_assert_(gsp_input);
+                DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECiniEnum));_assert_(EC_input);
+                DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWiniEnum));_assert_(W_input);
+                DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAiniEnum));_assert_(a_input);
+                DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTiniEnum));_assert_(T_input);
+
+                dz_input->GetValues(&dzini,&m);
+                d_input->GetValues(&dini,&m);
+                re_input->GetValues(&reini,&m);
+                gdn_input->GetValues(&gdnini,&m);
+                gsp_input->GetValues(&gspini,&m);
+                EC_input->GetInputValue(&EC);
+                W_input->GetValues(&Wini,&m);
+                a_input->GetValues(&aini,&m);
+                T_input->GetValues(&Tini,&m);
+
+                /*Retrive the correct value of m (without the zeroes at the end)*/
+                Input* Size_input=this->GetInput(SmbSizeiniEnum); _assert_(Size_input);
+                Size_input->GetInputValue(&m);
+
+                if(m==2){ //Snow properties are initialized with default values. Vertical grid has to be initialized too
+                        //            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w DEFAULT values\n");
+
+                        /*initialize profile variables:*/
+                        GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);
+
+                        d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=dini[0]; //ice density [kg m-3]
+                        re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=reini[0];         //set grain size to old snow [mm]
+                        gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=gdnini[0];         //set grain dentricity to old snow
+                        gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=gspini[0];         //set grain sphericity to old snow
+                        W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0];             //set water content to zero [kg m-2]
+                        a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0];         //set albedo equal to fresh snow [fraction]
+                        T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
+                        /*            /!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties) because don't know Tmean yet when set default values.
+                                                          Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/
+
+                        //fixed lower temperature bounday condition - T is fixed
+                        T_bottom=T[m-1];
+                }
+                else{ //Retrieve snow properties from previous run. Need to provide values for all layers
+                        //            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n");
+
+                        dz = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzini[i];
+                        d = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=dini[i];
+                        re = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=reini[i];
+                        gdn = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnini[i];
+                        gsp = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gspini[i];
+                        W = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wini[i];
+                        a = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=aini[i];
+                        T = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Tini[i];
+
+                        //fixed lower temperature bounday condition - T is fixed
+                        T_bottom=T[m-1];
+                }
+
+                /*Flag the initialization:*/
+                this->AddInput(new DoubleInput(SmbIsInitializedEnum,1.0));
+        }
+        else{
+                /*Recover inputs: */
+                DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum)); _assert_(dz_input);
+                DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));_assert_(d_input);
+                DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));_assert_(re_input);
+                DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));_assert_(gdn_input);
+                DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));_assert_(gsp_input);
+                DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));_assert_(EC_input);
+                DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));_assert_(W_input);
+                DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));_assert_(a_input);
+                DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));_assert_(T_input);
+
+                /*Recover arrays: */
+                dz_input->GetValues(&dz,&m);
+                d_input->GetValues(&d,&m);
+                re_input->GetValues(&re,&m);
+                gdn_input->GetValues(&gdn,&m);
+                gsp_input->GetValues(&gsp,&m);
+                EC_input->GetInputValue(&EC);
+                W_input->GetValues(&W,&m);
+                a_input->GetValues(&a,&m);
+                T_input->GetValues(&T,&m);
+
+                //fixed lower temperature bounday condition - T is fixed
+                T_bottom=T[m-1];
+
+        } /*}}}*/
 
         // determine initial mass [kg]
         initMass=0; for(int i=0;i<m;i++) initMass += dz[i]*d[i] + W[i];
-    
-        // initialize cumulative variables
+
+        // initialize cumulative variables
         sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
         sumdz_add=0;
@@ -2771 +2790 @@
 
                 /*Snow, firn and ice albedo:*/
-                if(isalbedo)albedo(a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,m,this->Sid());
-                
-                                        
+                if(isalbedo)albedo(a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
+
+
                 /*Distribution of absorbed short wave radation with depth:*/
-                if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,m,this->Sid());
-                
+                if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
+
                 /*Calculate net shortwave [W m-2]*/
                 netSW = cellsum(swf,m);
 
                 /*Thermal profile computation:*/
-                if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,this->Sid());
+                if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,rho_ice,this->Sid());
 
                 /*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell. 
                  * need to fix this in case all or more of cell evaporates */
                 dz[0] = dz[0] + EC / d[0];
-                
+
                 /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
-                if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,this->Sid());
+                if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,rho_ice,this->Sid());
 
                 /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
                  * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
-                if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,this->Sid());
+                if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid());
 
                 /*Allow non-melt densification and determine compaction [m]*/
                 if(isdensification)densification(d,dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid());
-                
+
                 /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every 
                  * sub-time step in thermo equations*/
@@ -2803 +2822 @@
                 /*Calculate net longwave [W m-2]*/
                 netLW = dlw - ulw;
-                
+
                 /*Calculate turbulent heat fluxes [W m-2]*/
-                if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,this->Sid());
-                
+                if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());
+
                 /*Verbose some resuls in debug mode: {{{*/
                 if(VerboseSmb() && 0){ 
                         _printf_("smb log: count[" << count << "] m[" << m << "] " 
-                                << setprecision(16)   << "T[" << cellsum(T,m)  << "] " 
-                                                          << "d[" << cellsum(d,m)  << "] "
-                                                          << "dz[" << cellsum(dz,m)  << "] "
-                                                          << "a[" << cellsum(a,m)  << "] "
-                                                          << "W[" << cellsum(W,m)  << "] "
-                                                          << "re[" << cellsum(re,m)  << "] "
-                                                          << "gdn[" << cellsum(gdn,m)  << "] "
-                                                          << "gsp[" << cellsum(gsp,m)  << "] "
-                                                          << "swf[" << netSW << "] "
-                                                                          << "\n");
+                                                << setprecision(16)   << "T[" << cellsum(T,m)  << "] " 
+                                                << "d[" << cellsum(d,m)  << "] "
+                                                << "dz[" << cellsum(dz,m)  << "] "
+                                                << "a[" << cellsum(a,m)  << "] "
+                                                << "W[" << cellsum(W,m)  << "] "
+                                                << "re[" << cellsum(re,m)  << "] "
+                                                << "gdn[" << cellsum(gdn,m)  << "] "
+                                                << "gsp[" << cellsum(gsp,m)  << "] "
+                                                << "swf[" << netSW << "] "
+                                                << "\n");
                 } /*}}}*/
-                
+
                 /*Sum component mass changes [kg m-2]*/
                 sumMassAdd = mAdd + sumMassAdd;
@@ -2841 +2860 @@
                 if (dMass != 0.0) _printf_("total system mass not conserved in MB function \n");
                 #endif
-                
+
                 /*Check bottom grid cell T is unchanged:*/
                 if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
-                
+
                 /*Free ressources: */
                 xDelete<IssmDouble>(swf);
@@ -2859 +2878 @@
         this->AddInput(new DoubleArrayInput(SmbGspEnum,gsp,m));
         this->AddInput(new DoubleArrayInput(SmbTEnum,T,m));
-        this->AddInput(new DoubleInput(SmbECEnum,sumEC/yts));
+        this->AddInput(new DoubleInput(SmbECEnum,sumEC/dt/rho_ice));
         this->AddInput(new DoubleArrayInput(SmbWEnum,W,m));
         this->AddInput(new DoubleArrayInput(SmbAEnum,a,m));
-        this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/yts));
-        this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/yts));
-        this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/yts));
+        this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/dt/rho_ice));
+        this->AddInput(new DoubleInput(SmbMeltEnum,sumM/dt/rho_ice));
+        this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/dt/rho_ice));
+        this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/dt/rho_ice));
         this->AddInput(new DoubleInput(SmbDz_addEnum,sumdz_add/yts));
-        this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/yts));
+        this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/dt));
 
         /*Free allocations:{{{*/

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

@@ -8 +8 @@
 
 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 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 (mol-1 K-1)
 
 void Gembx(FemModel* femmodel){  /*{{{*/
@@ -25 +43 @@
 
         /*intermediary:*/
-        IssmDouble dgpTop;
-        int gpTop, gpBottom;
-        int i;
-        IssmDouble gp0,z0;
+        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;
@@ -101 +121 @@
 
         // initialize
-        IssmDouble F = 0, H=0, G=0;
+        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 ºC
-        T = T - 273.15;
+        T = T - CtoK;
 
         // temperature coefficient F
         if(T>-6.0) F =  0.7 + ((T/-6.0) * 0.3);
-        if(T<=-6.0 && T>-22.0) F =  1 - ((T+6.0)/-16.0 * 0.8);
+        if(T<=-6.0 && T>-22.0) F =  1.0 - ((T+6.0)/-16.0 * 0.8);
         if(T<=-22.0 && T>-40.0) F =  0.2 - ((T+22.0)/-18.0 * 0.2);
 
@@ -114 +134 @@
         if(d<150.0) H=1.0;
 
-        if(d>=150.0 & d<400.0) H = 1 - ((d-150.0)/250.0);
+        if(d>=150.0 & d<400.0) H = 1.0 - ((d-150.0)/250.0);
 
         // temperature gradient coefficient G
@@ -121 +141 @@
         if(dT >= 0.4  && dT < 0.5)  G = 0.67 + (((dT - 0.4)/0.1) * 0.23);
         if(dT >= 0.5  && dT < 0.7)  G = 0.9 + (((dT - 0.5)/0.2) * 0.1);
-        if(dT >= .7              )  G = 1.0;
+        if(dT >= 0.7              )  G = 1.0;
 
         // grouped coefficient Q
@@ -164 +184 @@
 
         /*intermediary: */
-        IssmDouble  dt;
+        IssmDouble  dt=0.0;
         IssmDouble* gsz=NULL;
         IssmDouble* dT=NULL;
         IssmDouble* zGPC=NULL;
         IssmDouble* lwc=NULL;
-        IssmDouble  Q=0;
+        IssmDouble  Q=0.0;
 
         if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   grain growth module\n");
@@ -183 +203 @@
 
         /*Convert dt from seconds to day: */
-        dt=smb_dt/86400.0;
+        dt=smb_dt/dts;
 
         /*Determine liquid-water content in percentage: */
@@ -189 +209 @@
 
         //set maximum water content by mass to 9 percent (Brun, 1980)
-        for(int i=0;i<m;i++)if(lwc[i]>9.0) lwc[i]=9.0;
+        for(int i=0;i<m;i++)if(lwc[i]>9.0-Wtol) lwc[i]=9.0;
 
 
@@ -202 +222 @@
         for(int i=0;i<m;i++){
                 for (int j=0;j<=i;j++) zGPC[i]+=dz[j];
-                zGPC[i]-=dz[i]/2;
+                zGPC[i]-=(dz[i]/2.0);
         }
 
@@ -217 +237 @@
         /*Snow metamorphism. Depends on value of snow dendricity and wetness of the snowpack: */
         for(int i=0;i<m;i++){
-                if (gdn[i]>0){
+                if (gdn[i]>0.0+Gdntol){
 
                         //_printf_("Dendritic dry snow metamorphism\n");
 
                         //index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1
-                        if(dT[i]<5.0){
+                        if(dT[i]<5.0-Ttol){
                                 //determine coefficients
                                 IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt;
@@ -240 +260 @@
 
                         // wet snow metamorphism
-                        if(W[i]>0.0){
+                        if(W[i]>0.0+Wtol){
 
                                 //_printf_("D}ritic wet snow metamorphism\n");
@@ -253 +273 @@
 
                         // dendricity and sphericity can not be > 1 or < 0
-                        if (gdn[i]<0.0)gdn[i]=0.0;
-                        if (gsp[i]<0.0)gsp[i]=0.0;
-                        if (gdn[i]>1.0)gdn[i]=1.0;
-                        if (gsp[i]>1.0)gsp[i]=1.0;
+                        if (gdn[i]<0.0+Wtol)gdn[i]=0.0;
+                        if (gsp[i]<0.0+Wtol)gsp[i]=0.0;
+                        if (gdn[i]>1.0-Wtol)gdn[i]=1.0;
+                        if (gsp[i]>1.0-Wtol)gsp[i]=1.0;
 
                         // determine new grain size (mm)
-                        gsz[i] = 0.1 + (1-gdn[i])*0.25 + (0.5-gsp[i])*0.1;
+                        gsz[i] = 0.1 + (1.0-gdn[i])*0.25 + (0.5-gsp[i])*0.1;
 
                 }
@@ -271 +291 @@
 
                         // calculate grain growth
-                        gsz[i] += Q* dt;
+                        gsz[i] += (Q*dt);
 
                         //Wet snow metamorphism (Brun, 1989)
-                        if(W[i]>0.0){
+                        if(W[i]>0.0+Wtol){
                                 //_printf_("Nond}ritic wet snow metamorphism\n");
                                 //wet rate of change coefficient
-                                IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *86400.0);   // [mm^3 s^-1]
+                                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
@@ -285 +305 @@
 
                         // grains with sphericity == 1 can not have grain sizes > 2 mm (Brun, 1992)
-                        if (gsp[i]==1.0 && gsz[i]>2.0) gsz[i]=2.0;
+                        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 (gsp[i]!=1.0 && gsz[i]>5.0) gsz[i]=5.0;
+                        if (fabs(gsp[i]-1.0)>Wtol && gsz[i]>5.0-Wtol) gsz[i]=5.0;
                 }
 
@@ -302 +322 @@
 
 }  /*}}}*/
-void albedo(IssmDouble* a, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, int m,int sid) { /*{{{*/
+void albedo(IssmDouble* a, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
 
         //// Calculates Snow, firn and ice albedo as a function of:
@@ -350 +370 @@
         //some constants:
         const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
-        const IssmDouble dIce = 910;    // density of ice [kg m-3]
 
         if(aIdx==1){ 
@@ -356 +375 @@
         
         //convert effective radius to specific surface area [cm2 g-1]
-        IssmDouble S = 3.0 / (.091 * re[0]);
+        IssmDouble S = 3.0 / (0.091 * re[0]);
         
         //determine broadband albedo
-        a[0]= 1.48 - pow(S,-.07);
+        a[0]= 1.48 - pow(S,-0.07);
         }
         else if(aIdx==2){
@@ -369 +388 @@
         
         // convert effective radius to grain size in meters
-        IssmDouble gsz = (re[0] * 2) / 1000.0;
+        IssmDouble gsz = (re[0] * 2.0) / 1000.0;
         
         // spectral range:
@@ -398 +417 @@
         // where: t0 = timescale for albedo decay
         
-        dt = dt / 86400;    // convert from [s] to [d]
+        dt = dt / dts;    // convert from [s] to [d]
         
         // initialize variables
@@ -414 +433 @@
         // K = 7                // time scale temperature coef. (7) [d]
         // W0 = 300;            // 200 - 600 [mm]
-        const IssmDouble z_snow = 15;            // 16 - 32 [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)t0[i]=t0wet; // wet snow timescale
-        for(int i=0;i<m;i++)T[i]=TK[i] - 273.15; // change T from K to degC
+        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(W[i]==0.0 && T[i]>=-10)t0[i]= t0warm[i];
-        for(int i=0;i<m;i++)if(T[i]<-10) t0[i] =  10 * K + t0dry; // 'cold' 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
@@ -431 +450 @@
         
         // check if condensation occurs & if it is deposited in solid phase
-        if ( EC > 0 && T[0] < 0) P = P + (EC/dSnow) * 1000;  // add cond to precip [mm]
+        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);
@@ -453 +472 @@
         else if (xIsNan(a[0])) _error_("albedo == NAN\n");
 }  /*}}}*/
-void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz,int sid) { /*{{{*/
+void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid) { /*{{{*/
 
         /* ENGLACIAL THERMODYNAMICS*/
@@ -495 +514 @@
         IssmDouble* sw = NULL;
         IssmDouble* dT_sw = NULL;
-        IssmDouble* lw = NULL;
         IssmDouble* T0 = NULL;
         IssmDouble* Tu = NULL;
         IssmDouble* Td = NULL;
 
-        IssmDouble  z0; 
-        IssmDouble  dt;
-        IssmDouble max_fdt=0;
-        IssmDouble  Ts=0;
-        IssmDouble  L;
-        IssmDouble  eS;
-        IssmDouble  Ri=0;
-        IssmDouble  coefM;
-        IssmDouble  coefH;
-        IssmDouble An;
-        IssmDouble C;
-        IssmDouble shf;
-        IssmDouble SB;
-        IssmDouble CI; 
-        IssmDouble ds;
-        IssmDouble dAir;
-        IssmDouble TCs;
-        IssmDouble lhf;
-        IssmDouble EC_day;
-        IssmDouble dT_turb;
-        IssmDouble turb;
-        IssmDouble ulw;
-        IssmDouble dT_ulw;
-        IssmDouble dlw;
-        IssmDouble dT_dlw;
+        IssmDouble  z0=0.0;     
+        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 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 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;
+        IssmDouble EC=0.0;
 
         if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   thermal module\n");
 
-        // INITIALIZE
-        CI = 2102;          // heat capacity of snow/ice (J kg-1 k-1)
-        // CA = 1005;                  // heat capacity of air (J kg-1 k-1)
-        // LF = 0.3345E6;              // latent heat of fusion(J kg-1)
-        // LV = 2.495E6;               // latent heat of vaporization(J kg-1)
-        // dIce = 910;                 // density of ice [kg m-3]
-        // dSnow = 300;                // density of snow [kg m-3]
-        SB = 5.67E-8;       // Stefan-Boltzmann constant [W m-2 K-4]
-
         ds = d[0];      // density of top grid cell
 
         // calculated air density [kg/m3]
-        dAir = 0.029 * pAir /(8.314 * Ta);
+        dAir = 0.029 * pAir /(R * Ta);
 
         // thermal capacity of top grid cell [J/k]
@@ -549 +556 @@
 
         //initialize Evaporation - Condenstation 
-        EC = 0;
+        EC = 0.0;
         
         // check if all SW applied to surface or distributed throught subsurface
@@ -556 +563 @@
         // SURFACE ROUGHNESS (Bougamont, 2006)
         // wind/temperature surface roughness height [m]
-        if (ds < 910 && Ws == 0) z0 = 0.00012;       // 0.12 mm for dry snow
-        else if (ds >= 910) z0 = 0.0032;             // 3.2 mm for ice
+        if (ds < dIce-Dtol && fabs(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
 
         // if V = 0 goes to infinity therfore if V = 0 change
-        if(V<.01)V=.01;
+        if(V<0.01+Dtol)V=0.01;
         
         // Bulk-transfer coefficient for turbulent fluxes
@@ -574 +581 @@
 
         // for snow and firn (density < 910 kg m-3) (Sturn et al, 1997)
-        for(int i=0;i<m;i++) if(d[i]<910) K[i] = 0.138 - 1.01E-3 * d[i] + 3.233E-6 * (pow(d[i],2));
+        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]>=910) K[i] = 9.828 * exp(-5.7E-3*T[i]);
+        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
@@ -632 +639 @@
         max_fdt=f[0];
         for(int i=0;i<45;i++){
-                if (f[i]<dt)if(f[i]>=max_fdt)max_fdt=f[i];
+                if (f[i]<dt-Dtol)if(f[i]>=max_fdt-Dtol)max_fdt=f[i];
         }
         dt=max_fdt;
@@ -641 +648 @@
         Ap = xNew<IssmDouble>(m);
         for(int i=0;i<m;i++){
-                Au[i] = pow((dzU[i]/2/KP[i] + dz[i]/2/KU[i]),-1);
-                Ad[i] = pow((dzD[i]/2/KP[i] + dz[i]/2/KD[i]),-1);
+                Au[i] = pow((dzU[i]/2.0/KP[i] + dz[i]/2.0/KU[i]),-1.0);
+                Ad[i] = pow((dzD[i]/2.0/KP[i] + dz[i]/2.0/KD[i]),-1.0);
                 Ap[i] = (d[i]*dz[i]*CI)/dt;
         }
@@ -653 +660 @@
                 Nu[i] = Au[i] / Ap[i];
                 Nd[i] = Ad[i] / Ap[i];
-                Np[i]= 1 - Nu[i] - Nd[i];
+                Np[i]= 1.0 - Nu[i] - Nd[i];
         }
         
         // specify boundary conditions: constant flux at bottom
-        Nu[m-1] = 0;
-        Np[m-1] = 1;
+        Nu[m-1] = 0.0;
+        Np[m-1] = 1.0;
         
         // zero flux at surface
-        Np[0] = 1 - Nd[0];
+        Np[0] = 1.0 - Nd[0];
         
         // Create neighbor arrays for diffusion calculations instead of a tridiagonal matrix
-        Nu[0] = 0;
-        Nd[m-1] = 0;
+        Nu[0] = 0.0;
+        Nd[m-1] = 0.0;
         
         /* RADIATIVE FLUXES*/
@@ -704 +711 @@
                 // calculated.  The estimated enegy balance & melt are significanly
                 // less when Ts is taken as the mean of the x top grid cells.
-                Ts = (T[0] + T[1])/2;
-                Ts = fmin(273.15,Ts);    // don't allow Ts to exceed 273.15 K (0°C)
+                Ts = (T[0] + T[1])/2.0;
+                Ts = fmin(CtoK,Ts);    // don't allow Ts to exceed 273.15 K (0 degC)
                 
                 //TURBULENT HEAT FLUX
     
-                // Monin–Obukhov Stability Correction
+                // Monin-Obukhov Stability Correction
                 // Reference:
                 // Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
@@ -715 +722 @@
 
                 // calculate the Bulk Richardson Number (Ri)
-                Ri = (2*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
+                Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
                 
-                // calculate Monin-–Obukhov stability factors 'coefM' and 'coefH'
+                // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
     
                 // do not allow Ri to exceed 0.19
@@ -723 +730 @@
 
                 // calculate momentum 'coefM' stability factor
-                if (Ri > 0){
+                if (Ri > 0.0+Ttol){
                         // if stable
-                        coefM = 1/(1-5.2*Ri);
+                        coefM = 1.0/(1.0-5.2*Ri);
                 }
                 else {
-                        coefM =pow (1-18*Ri,-0.25);
+                        coefM =pow (1.0-18*Ri,-0.25);
                 }
                 
                 // calculate heat/wind 'coef_H' stability factor
-                if (Ri < -0.03) coefH = 1.3 * coefM;
+                if (Ri < -0.03-Ttol) coefH = 1.3 * coefM;
                 else coefH = coefM;
                 
                 //// Sensible Heat
                 // calculate the sensible heat flux [W m-2](Patterson, 1998)
-                shf = C * 1005 * (Ta - Ts);
+                shf = C * CA * (Ta - Ts);
 
                 // adjust using Monin–Obukhov stability theory
@@ -744 +751 @@
                 //// Latent Heat
                 // determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water
-                if (Ts >= 273.15){
-                        L = 2.495E6;
+                if (Ts >= CtoK-Ttol){
+                        L = LV;
 
                         // for an ice surface Murphy and Koop, 2005 [Equation 7]
@@ -751 +758 @@
                 }
                 else{
-                        L = 2.8295E6; // latent heat of sublimation for liquid surface (assume liquid on surface when Ts == 0 deg C)
+                        L = LS; // latent heat of sublimation 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 - 273.15) / (240.97 + Ts - 273.15));
+                        eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
                 }
                 
@@ -759 +766 @@
                 lhf = C * L * (eAir - eS) * 0.622 / pAir;
                 
-                // adjust using Monin–Obukhov stability theory (if lhf '+' then there is energy and mass gained at the surface, 
+                // 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 * coefH);
 
                 //mass loss (-)/acreation(+) due to evaporation/condensation [kg]
-                EC_day = lhf * 86400 / L;
+                EC_day = lhf * dts / L;
 
                 // temperature change due turbulent fluxes
@@ -787 +794 @@
                 
                 // calculate cumulative evaporation (+)/condensation(-)
-                EC = EC + (EC_day/86400)*dt;
+                EC = EC + (EC_day/dts)*dt;
     
                 /* CHECK FOR ENERGY (E) CONSERVATION [UNITS: J]
@@ -819 +826 @@
         xDelete<IssmDouble>(sw);
         xDelete<IssmDouble>(dT_sw);
-        xDelete<IssmDouble>(lw);
         xDelete<IssmDouble>(T0);
         xDelete<IssmDouble>(Tu);
@@ -829 +835 @@
 
 }  /*}}}*/
-void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid){ /*{{{*/
+void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid){ /*{{{*/
 
         // DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE
@@ -865 +871 @@
         
                 // calculate surface shortwave radiation fluxes [W m-2]
-                swf[0] = (1 - as) * dsw;
+                swf[0] = (1.0 - as) * dsw;
         }
         else{ // sw radation is absorbed at depth within the glacier
@@ -884 +890 @@
                         // convert effective radius [mm] to grain size [m]
                         gsz=xNew<IssmDouble>(m);
-                        for(int i=0;i<m;i++) gsz[i]= (re[i] * 2) / 1000;
+                        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]
@@ -894 +900 @@
                         B2_cum=xNew<IssmDouble>(m+1);
                         for(int i=0;i<m+1;i++){
-                                B1_cum[i]=1;
-                                B2_cum[i]=1;
+                                B1_cum[i]=1.0;
+                                B2_cum[i]=1.0;
                         }
 
@@ -901 +907 @@
                         // spectral albedos:
                         // 0.3 - 0.8um
-                        IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz[0],0.5));
+                        IssmDouble a0 = fmin(0.98, 1.0 - 1.58 *pow(gsz[0],0.5));
                         // 0.8 - 1.5um
-                        IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz[0],0.5));
+                        IssmDouble a1 = fmax(0.0, 0.95 - 15.4 *pow(gsz[0],0.5));
                         // 1.5 - 2.8um
                         IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5));
@@ -909 +915 @@
                         // separate net shortwave radiative flux into spectral ranges
                         IssmDouble swfS[3];
-                        swfS[0] = (sF[0] * dsw) * (1 - a0);
-                        swfS[1] = (sF[1] * dsw) * (1 - a1);
-                        swfS[2] = (sF[2] * dsw) * (1 - a2);
+                        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 coeficient for spectral range:
@@ -918 +924 @@
                         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] = .0192 * h[i];                 // 0.3 - 0.8um
-                        for(int i=0;i<m;i++) B2[i]= .1098 * h[i];                 // 0.8 - 1.5um
+                        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
 
@@ -986 +992 @@
 
                         // calculate surface shortwave radiation fluxes [W m-2]
-                        IssmDouble swf_s = SWs * (1 - as) * dsw;
+                        IssmDouble swf_s = SWs * (1.0 - as) * dsw;
 
                         // calculate surface shortwave radiation fluxes [W m-2]
-                        IssmDouble swf_ss = (1-SWs) * (1 - as) * dsw;
+                        IssmDouble swf_ss = (1.0-SWs) * (1.0 - as) * dsw;
 
                         // SW allowed to penetrate into snowpack
-                        IssmDouble Bs = 10;    // snow SW extinction coefficient [m-1] (Bassford,2006)
+                        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 - d[i]) * ((Bs - Bi)/(910 - 300));
+                        for(int i=0;i<m;i++) B[i] = Bs + (300.0 - d[i]) * ((Bs - Bi)/(dIce - 300.0));
 
                         // cumulative extinction factor
@@ -1004 +1010 @@
                         for(int i=0;i<m;i++)exp_B[i]=exp(-B[i]*dz[i]);
 
-                        B_cum[0]=1;
+                        B_cum[0]=1.0;
                         for(int i=0;i<m;i++){
                                 IssmDouble cum_B=exp_B[0];
@@ -1032 +1038 @@
 
 } /*}}}*/ 
-void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, int sid){ /*{{{*/
+void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid){ /*{{{*/
 
         // Adds precipitation and deposition to the model grid
@@ -1051 +1057 @@
         // MAIN FUNCTION
         // specify constants
-        const IssmDouble dIce = 910;     // density of ice [kg m-3]
         const IssmDouble dSnow = 150;    // density of snow [kg m-3]
         const IssmDouble reNew = 0.1;    // new snow grain size [mm]
-        const IssmDouble gdnNew = 1;     // new snow dendricity 
+        const IssmDouble gdnNew = 1.0;     // new snow dendricity 
         const IssmDouble gspNew = 0.5;   // new snow sphericity 
 
@@ -1060 +1065 @@
         IssmDouble* mInit=NULL;
         bool        top=true;
-        IssmDouble  mass, massinit, mass_diff;
+        IssmDouble  mass=0;
+        IssmDouble  massinit=0;
+        IssmDouble  mass_diff=0;
 
         /*output: */
@@ -1071 +1078 @@
         IssmDouble* gdn=NULL;
         IssmDouble* gsp=NULL;
-        int         m;
+        int         m=0;
 
         if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   accumulation module\n");
@@ -1089 +1096 @@
         mInit=xNew<IssmDouble>(m);
         for(int i=0;i<m;i++) mInit[i]= d[i] * dz[i];
-        massinit=0; for(int i=0;i<m;i++)massinit+=mInit[i];
-
-        if (P > 0){
+        massinit=0.0; 
+        for(int i=0;i<m;i++)massinit+=mInit[i];
+
+        if (P > 0.0+Dtol){
                         
 
-                if (T_air <= 273.15){ // if snow
+                if (T_air <= CtoK+Ttol){ // if snow
 
                         IssmDouble  z_snow = P/dSnow;               // depth of snow
 
                         // if snow depth is greater than specified min dz, new cell created
-                        if (z_snow > dzMin){
+                        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,top,m); //new cell W
+                                newcell(&W,0.0,top,m); //new cell W
                                 newcell(&a,aSnow,top,m); //new cell a
                                 newcell(&re,reNew,top,m); //new cell grain size
@@ -1137 +1145 @@
                           makes the numerics easier.*/
 
-                        IssmDouble LF = 0.3345E6;  // latent heat of fusion(J kg-1)
-                        IssmDouble CI = 2102;      // specific heat capacity of snow/ice (J kg-1 k-1)
-
                         // grid cell adjust mass
                         mass = mInit[0] + P;
@@ -1151 +1156 @@
 
                         // if d > the density of ice, d = dIce
-                        if (d[0] > dIce){
+                        if (d[0] > dIce+Dtol){
                                 d[0] = dIce;           // adjust d
                                 dz[0] = mass / d[0];    // dz is adjusted to conserve mass
@@ -1163 +1168 @@
                 
                 #ifndef _HAVE_ADOLC_  //avoid round operation. only check in forward mode.
-                mass_diff = round(mass_diff * 100)/100;
+                mass_diff = round(mass_diff * 100.0)/100.0;
                 if (mass_diff > 0) _error_("mass not conserved in accumulation function");
                 #endif
@@ -1182 +1187 @@
         *pm=m;
 } /*}}}*/
-void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, int sid){ /*{{{*/
+void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid){ /*{{{*/
 
         //// MELT ROUTINE
@@ -1200 +1205 @@
         IssmDouble* F=NULL;
         IssmDouble* flxDn=NULL;
-        IssmDouble  ER=0;
+        IssmDouble  ER=0.0;
         IssmDouble* EI=NULL;
         IssmDouble* EW=NULL;
@@ -1206 +1211 @@
         int*       D=NULL;
         
-        IssmDouble sumM;
-        IssmDouble sumER;
-        IssmDouble addE;
-        IssmDouble mSum0;
-        IssmDouble sumE0;
-        IssmDouble mSum1;
-        IssmDouble sumE1;
-        IssmDouble dE;
-        IssmDouble dm;
-        IssmDouble X;
-        IssmDouble Wi;
+        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;
     
-        IssmDouble Ztot;
-        IssmDouble T_bot;
-        IssmDouble m_bot;
-        IssmDouble dz_bot;
-        IssmDouble d_bot;
-        IssmDouble W_bot;
-        IssmDouble a_bot;
-        IssmDouble re_bot;
-        IssmDouble gdn_bot;
-        IssmDouble gsp_bot;
-        IssmDouble EI_bot;
-        IssmDouble EW_bot;
+        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 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;
     
@@ -1235 +1240 @@
         IssmDouble* dzMin2=NULL;
         IssmDouble zY2=1.025;
-        bool lastCellFlag;
+        bool lastCellFlag = false;
         int X1=0;
         int X2=0;
     
-        int        D_size;
+        int        D_size = 0;
 
         /*outputs:*/
-        IssmDouble  mAdd;
-        IssmDouble dz_add;
-        IssmDouble  Rsum;
+        IssmDouble  mAdd = 0.0;
+        IssmDouble dz_add = 0.0;
+        IssmDouble  Rsum = 0.0;
         IssmDouble* T=*pT;
         IssmDouble* d=*pd;
@@ -1254 +1259 @@
         IssmDouble* gsp=*pgsp;
         int         n=*pn;
-        IssmDouble* R=0;
+        IssmDouble* R=NULL;
         
         if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   melt module\n");
@@ -1265 +1270 @@
         dW=xNew<IssmDouble>(n); 
 
-        // specify constants
-        const IssmDouble CtoK = 273.15;  // clecius to Kelvin conversion
-        const IssmDouble CI = 2102;      // specific heat capacity of snow/ice (J kg-1 k-1)
-        const IssmDouble LF = 0.3345E6;  // latent heat of fusion(J kg-1)
-        const IssmDouble dPHC = 830;     // pore hole close off density[kg m-3]
-        const IssmDouble dIce = 910;     // density of ice [kg m-3]
-
         // 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]
@@ -1281 +1279 @@
 
         // initialize melt and runoff scalars
-        Rsum = 0;       // runoff [kg]
-        sumM = 0;       // total melt [kg]
-        mAdd = 0;       // mass added/removed to/from base of model [kg]
-        addE = 0;       // energy added/removed to/from base of model [J]
-        dz_add=0;      // thickness of the layer added/removed to/from base of model [m]
+        Rsum = 0.0;       // runoff [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]= fmax(0, T[i] - CtoK);        // [K] to [°C]
+        for(int i=0;i<n;i++) exsT[i]= fmax(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]-=exsT[i];
+        for(int i=0;i<n;i++) T[i]=fmin(T[i],CtoK);
+        
         // specify irreducible water content saturation [fraction]
         const IssmDouble Swi = 0.07;                     // assumed constant after Colbeck, 1974
@@ -1299 +1298 @@
         //// REFREEZE PORE WATER
         // check if any pore water
-        if (cellsum(W,n) > 0){
+        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] = fmax(0, -((T[i] - CtoK) * m[i] * CI) / LF);
+                for(int i=0;i<n;i++) maxF[i] = fmax(0.0, -((T[i] - CtoK) * m[i] * CI) / LF);
 
                 // freeze pore water and change snow/ice properties
@@ -1319 +1318 @@
         exsW=xNew<IssmDouble>(n); 
         for(int i=0;i<n;i++){
-                Wi= (910 - d[i]) * Swi * (m[i] / d[i]);        // irreducible water content [kg]
-                exsW[i] = fmax(0, W[i] - Wi);                  // water "squeezed" from snow [kg]
+                Wi= (dIce - d[i]) * Swi * (m[i] / d[i]);        // irreducible water content [kg]
+                exsW[i] = fmax(0.0, W[i] - Wi);                  // water "squeezed" from snow [kg]
         }
 
@@ -1326 +1325 @@
 
         // run melt algorithm if there is melt water or excess pore water
-        if ((cellsum(exsT,n) > 0) || (cellsum(exsW,n) > 0)){
+        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
@@ -1332 +1331 @@
                 // (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] = fmax(0, exsT [i]- 159.1342);
-        
-                if (cellsum(surpT,n) > 0 ){
+                surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0.0, exsT[i]- 159.1342);
+        
+                if (cellsum(surpT,n) > 0.0 + Ttol ){
                         // _printf_("T Surplus");
                         // calculate surplus energy
@@ -1340 +1339 @@
             
                         int i = 0;
-                        while (cellsum(surpE,n) > 0){
+                        while (cellsum(surpE,n) > 0.0 + Ttol){
                                 // 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] = fmax(0, T[i+1] - CtoK) + exsT[i+1];
+                                exsT[i+1] = fmax(0.0, T[i+1] - CtoK) + exsT[i+1];
                                 T[i+1] = fmin(CtoK, T[i+1]);
                 
-                                surpT[i+1] = fmax(0, exsT[i+1] - 159.1342);
+                                surpT[i+1] = fmax(0.0, exsT[i+1] - 159.1342);
                                 surpE[i+1] = surpT[i+1] * CI * m[i+1];
                 
                                 // adjust current cell properties (again 159.1342 is the max T)
                                 exsT[i] = 159.1342;
-                                surpE[i] = 0;
+                                surpE[i] = 0.0;
                                 i = i + 1;
                         }
@@ -1362 +1361 @@
 
                 // calculate maximum refreeze amount, maxF [kg]
-                for(int i=0;i<n;i++)maxF[i] = fmax(0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);
+                for(int i=0;i<n;i++)maxF[i] = fmax(0.0, -((T[i] - CtoK) * d[i] * dz[i] * CI)/ LF);
 
                 // initialize refreeze, runoff, flxDn and dW vectors [kg]
@@ -1368 +1367 @@
                 IssmDouble* R=xNewZeroInit<IssmDouble>(n);
 
-                for(int i=0;i<n;i++)dW[i] = 0;
+                for(int i=0;i<n;i++)dW[i] = 0.0;
                 flxDn=xNewZeroInit<IssmDouble>(n+1); for(int i=0;i<n;i++)flxDn[i+1]=F[i];
 
@@ -1374 +1373 @@
                 X = 0;
                 for(int i=n-1;i>=0;i--){
-                        if(M[i]>0 || reCast<int,IssmDouble>(exsW[i])){
+                        if(M[i]>0.0+Wtol || exsW[i]>0.0+Wtol){
                                 X=i;
                                 break;
@@ -1386 +1385 @@
 
                         // break loop if there is no meltwater and if depth is > mw_depth
-                        if (inM == 0 && i > X){
+                        if (inM == 0.0 && i > X){
                                 break;
                         }
 
                         // if reaches impermeable ice layer all liquid water runs off (R)
-                        else if (d[i] >= dIce){  // dPHC = pore hole close off [kg m-3]
+                        else if (d[i] >= dIce-Dtol){  // dPHC = pore hole close off [kg m-3]
                                 // _printf_("ICE LAYER");
                                 // no water freezes in this cell
@@ -1398 +1397 @@
 
                                 m[i] = m[i] - M[i];                     // mass after melt
-                                Wi = (910-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
+                                Wi = (dIce-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
                                 dW[i] = fmin(inM, Wi - W[i]);            // change in pore water
-                                R[i] = fmax(0, inM - dW[i]);             // runoff
+                                R[i] = fmax(0.0, inM - dW[i]);             // runoff
                         }
                         // check if no energy to refreeze meltwater
-                        else if (maxF[i] == 0){
+                        else if (fabs(maxF[i]) < Dtol){
                                 // _printf_("REFREEZE == 0");
                                 // no water freezes in this cell
@@ -1409 +1408 @@
 
                                 m[i] = m[i] - M[i];                     // mass after melt
-                                Wi = (910-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
+                                Wi = (dIce-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
                                 dW[i] = fmin(inM, Wi-W[i]);              // change in pore water
-                                flxDn[i+1] = fmax(0, inM-dW[i]);         // meltwater out
-                                F[i] = 0;                               // no freeze 
+                                flxDn[i+1] = fmax(0.0, inM-dW[i]);         // meltwater out
+                                F[i] = 0.0;                               // no freeze 
                         }
                         // some or all meltwater refreezes
@@ -1426 +1425 @@
 
                                 //-----------------------pore water-----------------------------
-                                Wi = (910-d[i])* Swi * dz_0;            // irreducible water 
+                                Wi = (dIce-d[i])* Swi * dz_0;            // irreducible water 
                                 dW[i] = fmin(inM - F1, Wi-W[i]);         // change in pore water
-                                if (-dW[i]>W[i] ){
-                                        dW[i]= W[i];
+                                if (-1*dW[i]>W[i]-Wtol ){
+                                        dW[i]= -1*W[i];
                                 }
-                                IssmDouble F2 = 0;                                 
-
-                                if (dW[i] < 0){                         // excess 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 = fmin(dMax, maxF[i]-F1);      // maximum refreeze
-                                        F2 = fmin(-dW[i], maxF2);            // pore water refreeze
+                                        F2 = fmin(-1*dW[i], maxF2);            // pore water refreeze
                                         m[i] = m[i] + F2;                   // mass after refreeze
                                         d[i] = m[i]/dz_0;
@@ -1446 +1445 @@
 
                                 // check if an ice layer forms 
-                                if (d[i] == dIce){
+                                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;
+                                        flxDn[i+1] = 0.0;
                                 }
                         }
@@ -1458 +1457 @@
 
                 //// GRID CELL SPACING AND MODEL DEPTH
-                for(int i=0;i<n;i++)if (W[i] < 0) _error_("negative pore water generated in melt equations");
+                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
@@ -1465 +1464 @@
 
                 // delete all cells with zero mass
-                D_size=0; for(int i=0;i<n;i++)if(m[i]!=0)D_size++; D=xNew<int>(D_size);
+                D_size=0; for(int i=0;i<n;i++)if(m[i]!=0)D_size++; 
+                D=xNew<int>(D_size);
                 D_size=0; for(int i=0;i<n;i++)if(m[i]!=0){ D[D_size] = i; D_size++;}
                 
@@ -1504 +1504 @@
     
         for (int i=0;i<n;i++){
-                if (Zcum[i]<=zTop){
+                if (Zcum[i]<=zTop+Dtol){
                         dzMin2[i]=dzMin;
                 }
@@ -1515 +1515 @@
         X = 0;
         for(int i=n-1;i>=0;i--){
-                if(dz[i]<dzMin2[i]){
+                if(dz[i]<dzMin2[i]-Dtol){
                         X=i;
                         break;
@@ -1531 +1531 @@
 
         for (int i = 0; i<=X;i++){
-                if (dz [i] < dzMin2[i]){                              // merge top cells
+                if (dz[i] < dzMin2[i]-Dtol){                              // merge top cells
                         // _printf_("dz > dzMin * 2')
                         // adjust variables as a linearly weighted function of mass
@@ -1542 +1542 @@
             
                         // merge with underlying grid cell and delete old cell
-                        dz [i+1] = dz[i] + dz[i+1];                 // combine cell depths
+                        dz[i+1] = dz[i] + dz[i+1];                 // combine cell depths
                         d[i+1] = m_new / dz[i+1];                   // combine top densities
                         W[i+1] = W[i+1] + W[i];                     // combine liquid water
@@ -1572 +1572 @@
         
                 // merge with underlying grid cell and delete old cell
-                dz [X1] = dz[X2] + dz[X1];                 // combine cell depths
+                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
@@ -1582 +1582 @@
 
         // delete combined cells
-        D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999)D_size++; D=xNew<int>(D_size); 
+        D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999)D_size++; 
+        D=xNew<int>(D_size); 
         D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999){ D[D_size] = i; D_size++;}
 
@@ -1602 +1603 @@
         X=0;
         for(int i=9;i>=0;i--){
-                if(dz[i]> 2* dzMin){
+                if(dz[i]> 2.0*dzMin+Dtol){
                         X=i;
                         break;
@@ -1608 +1609 @@
         }
         
-        int i=0;
-        while(i<=X){
-                if (dz [i] > dzMin *2){
+        int j=0;
+        while(j<=X){
+                if (dz[j] > dzMin*2.0+Dtol){
 
                                 // _printf_("dz > dzMin * 2");
                                 // split in two
-                                cellsplit(&dz, n, i,.5);
-                                cellsplit(&W, n, i,.5);
-                                cellsplit(&m, n, i,.5);
-                                cellsplit(&T, n, i,1.0);
-                                cellsplit(&d, n, i,1.0);
-                                cellsplit(&a, n, i,1.0);
-                                cellsplit(&EI, n, i,1.0);
-                                cellsplit(&EW, n, i,1.0);
-                                cellsplit(&re, n, i,1.0);
-                                cellsplit(&gdn, n, i,1.0);
-                                cellsplit(&gsp, n, i,1.0);
+                                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(&EI, n, j,1.0);
+                                cellsplit(&EW, n, j,1.0);
+                                cellsplit(&re, n, j,1.0);
+                                cellsplit(&gdn, n, j,1.0);
+                                cellsplit(&gsp, n, j,1.0);
                                 n++;
                                 X=X+1;
                 }
-                else i++;
+                else j++;
         }
 
@@ -1637 +1638 @@
         Ztot = cellsum(dz,n);
     
-        if (Ztot < zMin){
+        if (Ztot < zMin-Dtol){
                 // printf("Total depth < zMin %f \n", Ztot);
                 // mass and energy to be added
@@ -1671 +1672 @@
                 n=n+1;
         }
-        else if (Ztot > zMax){
+        else if (Ztot > zMax+Dtol){
                 // printf("Total depth > zMax %f \n", Ztot);
                 // mass and energy loss
@@ -1710 +1711 @@
 
         /*checks: */
-        for(int i=0;i<n;i++) if (W[i]<0) _error_("negative pore water generated in melt equations\n");
+        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: */
@@ -1753 +1754 @@
 
 } /*}}}*/ 
-void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,IssmDouble dIce, int m, int sid){ /*{{{*/
+void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, 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 COMPATION IS COMPNSATED FOR BY TRACES OF SNOW???]
@@ -1795 +1796 @@
         //// MAIN FUNCTION
         // specify constants
-        dt      = dt / 86400;  // convert from [s] to [d]
+        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
@@ -1802 +1803 @@
 
         /*intermediary: */
-        IssmDouble c0,c1,H;
+        IssmDouble c0=0.0;
+        IssmDouble c1=0.0;
+        IssmDouble H=0.0;
         
         if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   densification module\n");
@@ -1815 +1818 @@
 
         IssmDouble* obp = xNew<IssmDouble>(m);
-        obp[0]=0;
+        obp[0]=0.0;
         for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];
         
@@ -1826 +1829 @@
                 switch (denIdx){
                         case 1: // Herron and Langway (1980)
-                                c0 = (11 * exp(-10160 / (T[i] * 8.314))) * C/1000;
-                                c1 = (575 * exp(-21400 / (T[i]* 8.314))) * pow(C/1000,.5);
+                                c0 = (11.0 * exp(-10160.0 / (T[i] * 8.314))) * C/1000.0;
+                                c1 = (575.0 * exp(-21400.0 / (T[i]* 8.314))) * 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./(T[i] * 8.314)) + (42400./(Tmean * 8.314))) * (C * 9.81);
+                                H = exp((-60000.0/(T[i] * 8.314)) + (42400.0/(Tmean * 8.314))) * (C * 9.81);
                                 c0 = 0.07 * H;
                                 c1 = 0.03 * H;
@@ -1840 +1843 @@
 
                                 // common variable
-                                H = exp((-60/(T[i] * 8.314))) * obp[i] / pow(re[i]/1000,2.0);
+                                H = exp((-60.0/(T[i] * 8.314))) * obp[i] / pow(re[i]/1000.0,2.0);
                                 c0 = 9.2e-9 * H;
                                 c1 = 3.7e-9 * H;
@@ -1846 +1849 @@
 
                         case 4: // Li and Zwally (2004)
-                                c0 = (C/dIce) * (139.21 - 0.542*Tmean)*8.36*pow(273.15 - T[i],-2.061);
+                                c0 = (C/dIce) * (139.21 - 0.542*Tmean)*8.36*pow(CtoK - T[i],-2.061);
                                 c1 = c0;
                                 break;
@@ -1852 +1855 @@
                         case 5: // Helsen et al. (2008)
                                 // common variable
-                                c0 = (C/dIce) * (76.138 - 0.28965*Tmean)*8.36*pow(273.15 - T[i],-2.061);
+                                c0 = (C/dIce) * (76.138 - 0.28965*Tmean)*8.36*pow(CtoK - T[i],-2.061);
                                 c1 = c0;
                                 break;
@@ -1858 +1861 @@
 
                 // new snow density
-                if(d[i] <= 550) d[i] = d[i] + (c0 * (dIce - d[i]) / 365 * dt);
-                else            d[i] = d[i] + (c1 * (dIce - d[i]) / 365 * dt);
+                if(d[i] <= 550.0+Dtol) d[i] = d[i] + (c0 * (dIce - d[i]) / 365.0 * dt);
+                else            d[i] = d[i] + (c1 * (dIce - d[i]) / 365.0 * dt);
 
                 //disp((num2str(nanmean(c0 .* (dIce - d(idx)) / 365 * dt))))
 
                 // do not allow densities to exceed the density of ice
-                if(d[i]>dIce)d[i]=dIce;
+                if(d[i] > dIce-Dtol) d[i]=dIce;
 
                 // calculate new grid cell length
@@ -1876 +1879 @@
 
 } /*}}}*/
-void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, int sid){ /*{{{*/
+void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid){ /*{{{*/
 
         //// TURBULENT HEAT FLUX
@@ -1899 +1902 @@
         //   Tz: height above ground at which temperature (T) was sampled [m]
 
-        //// FUNCTION INITILIZATION 
-
-        // CA = 1005;                    // heat capacity of air (J kg-1 k-1)
-        // LF = 0.3345E6;                // latent heat of fusion(J kg-1)
-        // LV = 2.495E6;                 // latent heat of vaporization(J kg-1)
-        // dIce = 910;                   // density of ice [kg m-3]
-        
         /*intermediary:*/
-        IssmDouble d_air;
-        IssmDouble Ri;
-        IssmDouble z0;
-        IssmDouble coef_M,coef_H;
-        IssmDouble An, C;
-        IssmDouble L, eS;
+        IssmDouble d_air=0.0;
+        IssmDouble Ri=0.0;
+        IssmDouble z0=0.0;
+        IssmDouble coef_M=0.0;
+        IssmDouble coef_H=0.0;
+        IssmDouble An=0.0;
+        IssmDouble C=0.0;
+        IssmDouble L=0.0;
+        IssmDouble eS=0.0;
 
         /*output: */
-        IssmDouble shf, lhf, EC;
+        IssmDouble shf=0.0;
+        IssmDouble lhf=0.0;
+        IssmDouble EC=0.0;
         
         if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   turbulentFlux module\n");
 
         // calculated air density [kg/m3]
-        d_air = 0.029 * pAir /(8.314 * Ta);
+        d_air = 0.029 * pAir /(R * Ta);
 
         //// Determine Surface Roughness
         // Bougamont, 2006
         // wind/temperature surface roughness height [m]
-        if (ds < 910 && Ws == 0) z0 = 0.00012;               // 0.12 mm for dry snow
-        else if (ds >= 910) z0 = 0.0032;                // 3.2 mm for ice 
+        if (ds < dIce-Dtol && fabs(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
 
-        //// Monin–Obukhov Stability Correction
+        //// Monin-Obukhov Stability Correction
         // Reference:
         // Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
@@ -1935 +1936 @@
 
         // if V = 0 goes to infinity therfore if V = 0 change
-        if(V< .01) V=.01;
+        if(V < 0.01) V=0.01;
 
         // calculate the Bulk Richardson Number (Ri)
-        Ri = (2*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2));
-
-        // calculate Monin–Obukhov stability factors 'coef_M' and 'coef_H'
+        Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2));
+
+        // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
 
         // do not allow Ri to exceed 0.19
-        if(Ri>.19)Ri= 0.19;
+        if(Ri>0.19-Ttol)Ri= 0.19;
 
         // calculate momentum 'coef_M' stability factor
-        if (Ri > 0) coef_M = pow(1-5.2*Ri,-1); // if stable
-        else coef_M = pow(1-18*Ri,-0.25);
+        if (Ri > 0.0) coef_M = pow(1.0-5.2*Ri,-1.0); // if stable
+        else coef_M = pow(1.0-18*Ri,-0.25);
 
         // calculate heat/wind 'coef_H' stability factor
@@ -1959 +1960 @@
         //// Sensible Heat
         // calculate the sensible heat flux [W m-2](Patterson, 1998)
-        shf = C * 1005 * (Ta - Ts);
-
-        // adjust using Monin–Obukhov stability theory
+        shf = C * CA * (Ta - Ts);
+
+        // adjust using Monin-Obukhov stability theory
         shf = shf / (coef_M * coef_H);
 
@@ -1967 +1968 @@
         // determine if snow pack is melting & calcualte surface vapour pressure
         // over ice or liquid water
-        if (Ts >= 273.15){
-                L = 2.495E6;
+        if (Ts >= CtoK-Ttol){
+                L = LV;
 
                 // for an ice surface Murphy and Koop, 2005 [Equation 7]
@@ -1974 +1975 @@
         }
         else{
-                L = 2.8295E6; // latent heat of sublimation
+                L = LS; // latent heat of sublimation
                 // for liquid surface (assume liquid on surface when Ts == 0 deg C)
                 // Wright (1997), US Meteorological Handbook from Murphy and Koop,
                 // 2005 Apendix A
-                eS = 611.21 * exp(17.502 * (Ts - 273.15) / (240.97 + Ts - 273.15));
+                eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
         }
 
@@ -1984 +1985 @@
         lhf = C * L * (eAir - eS) * 0.622 / pAir;
 
-        // adjust using Monin–Obukhov stability theory (if lhf '+' then there is
+        // 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. 
@@ -1990 +1991 @@
 
         // mass loss (-)/acreation(+) due to evaporation/condensation [kg]
-        EC = lhf * 86400 / L;
+        EC = lhf * dts / L;
 
         /*assign output poitners: */

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

@@ -271 +271 @@
         IssmDouble  h = -1.54e-5;
         IssmDouble  smb,smbref,anomaly,yts,z;
-    
-    /* Get constants */
-    femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
-    /*iomodel->FindConstant(&yts,"md.constants.yts");*/
-    /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
-    /*Mathieu original*/
-    /*IssmDouble  smb,smbref,z;*/
-    
+
+        /* Get constants */
+        femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
+        /*iomodel->FindConstant(&yts,"md.constants.yts");*/
+        /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
+        /*Mathieu original*/
+        /*IssmDouble  smb,smbref,z;*/
+
         /*Loop over all the elements of this partition*/
         for(int i=0;i<femmodel->elements->Size();i++){
@@ -298 +298 @@
                         anomaly = smblistref[v];
 
-            /* Henning edited acc. to Riannes equations*/
-            /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
-            z_critical = z_critical + dz;
-            
-            /* Calculate smb acc. to the surface elevation z */
-            if(z<z_critical){
+                        /* Henning edited acc. to Riannes equations*/
+                        /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
+                        z_critical = z_critical + dz;
+
+                        /* Calculate smb acc. to the surface elevation z */
+                        if(z<z_critical){
                                 smb = a + b*z + c;
                         }
                         else{
-                          smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
-                        }
-            
+                                smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
+                        }
+
                         /* Compute smb including anomaly,
                                 correct for number of seconds in a year [s/yr]*/

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

@@ -25 +25 @@
 IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT);
 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);
-void albedo(IssmDouble* a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt,int m, int sid);
-void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid);
-void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, int sid);
-void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, int sid); 
-void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, int sid);
-void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,IssmDouble rho_ice,int m, int sid);
-void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, int sid);
+void albedo(IssmDouble* a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
+void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid);
+void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid);
+void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid); 
+void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid);
+void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid);
+void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid);
 #endif  /* _SurfaceMassBalancex_H*/

issm/trunk-jpl/test/NightlyRun/test243.m

@@ -45 +45 @@
 %Fields and tolerances to track changes
 field_names      ={'SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance'};
-field_tolerances ={5e-4,5e-5,0.0006,0.0002,1e-5,0.0003,2e-5,2e-7,1e-7};
+field_tolerances ={1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12};
 
 field_values={...