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

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

@@ -1332 +1332 @@
         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 +1381 @@
                         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 +1688 @@
                 // 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 +1722 @@
                 // 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 +1743 @@
                 celldelete(&EW,n,D,D_size);
                 n=D_size;
                 xDelete<int>(D);
-        
         }
 
         //// CHECK FOR MASS AND ENERGY CONSERVATION
@@ -1762 +1760 @@
 
         /*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");