source: issm/oecreview/Archive/21724-22754/ISSM-22735-22736.diff@ 22755

Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp    (revision 22735)
+++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp    (revision 22736)
@@ -1332,7 +1332,7 @@
        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]
@@ -1381,21 +1381,20 @@
                        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){
+                       while (cellsum(surpE,n) > 0.0 + Ttol && 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] = fmax(0.0, T[i+1] - CtoK) + exsT[i+1];
-                               T[i+1] = fmin(CtoK, T[i+1]);
-                
-                               surpT[i+1] = fmax(0.0, exsT[i+1] - LF/CI);
+                               surpT[i+1] = fmax(0.0, (T[i+1] - CtoK - LF/CI));
                                surpE[i+1] = surpT[i+1] * CI * m[i+1];
-                
+
                                // adjust current cell properties (again 159.1342 is the max T)
-                               exsT[i] = LF/CI;
+                               T[i] = T[i]-surpE[i]/m[i+1]/CI;
                                surpE[i] = 0.0;
                                i = i + 1;
                        }
+                       // recalculate temperature excess above 0 deg C 
+                       for(int i=0;i<n;i++) exsT[i] = fmax(0.0, T[i] - CtoK);  
                }
 
                // convert temperature excess to melt [kg]
@@ -1689,7 +1688,7 @@
                // 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;
+               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];
@@ -1723,7 +1722,7 @@
                // 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);
+               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
@@ -1744,7 +1743,6 @@
                celldelete(&EW,n,D,D_size);
                n=D_size;
                xDelete<int>(D);
-        
        }
 
        //// CHECK FOR MASS AND ENERGY CONSERVATION
@@ -1762,7 +1760,7 @@
 
        /*only in forward mode! avoid round in AD mode as it is not differentiable: */
        #ifndef _HAVE_ADOLC_
-       dm = round(mSum0 - mSum1 + mAdd);
+       dm = round((mSum0 - mSum1 + mAdd)*1000);
        dE = round(sumE0 - sumE1 - sumER +  addE);
        if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
                        << "dm: " << dm << " dE: " << dE << "\n");