Changeset 22736

Timestamp:

05/02/18 18:25:10 (7 years ago)

Author:

schlegel

Message:

CHG: conserve energy when adding or removing bottom layers

File:

• issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (modified) (6 diffs)

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

@@ -1333 +1333 @@
 
         // 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);
         
@@ -1382 +1382 @@
             
                         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);  
                 }
 
@@ -1690 +1689 @@
                 // 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
@@ -1724 +1723 @@
                 // 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]);
         
@@ -1745 +1744 @@
                 n=D_size;
                 xDelete<int>(D);
-        
         }
 
@@ -1763 +1761 @@
         /*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"