Changeset 22614

Timestamp:

03/22/18 19:44:48 (7 years ago)

Author:

schlegel

Message:

CHG: add runoff before delete of no mass layer

File:

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

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

@@ -1247 +1247 @@
         IssmDouble* surpT=NULL;
         IssmDouble* surpE=NULL;
-        IssmDouble* F=NULL;
         IssmDouble* flxDn=NULL;
         IssmDouble  ER=0.0;
@@ -1311 +1310 @@
         /*Allocations: */
         M=xNewZeroInit<IssmDouble>(n); 
-        maxF=xNew<IssmDouble>(n); 
-        dW=xNew<IssmDouble>(n); 
+        maxF=xNewZeroInit<IssmDouble>(n); 
+        dW=xNewZeroInit<IssmDouble>(n); 
 
         // store initial mass [kg] and energy [J]
@@ -1355 +1354 @@
 
                 // if pore water froze in ice then adjust d and dz thickness
-                for(int i=0;i<n;i++)if(d[i]>dIce)d[i]=dIce;
+                for(int i=0;i<n;i++)if(d[i]>dIce-Dtol)d[i]=dIce;
                 for(int i=0;i<n;i++) dz[i]= m[i]/d[i];
         }
@@ -1375 +1374 @@
                 // (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.0, exsT[i]- 159.1342);
+                surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0.0, exsT[i]- LF/CI);
         
-                if (cellsum(surpT,n) > 0.0 + Ttol ){
+                if (cellsum(surpT,n) > 0.0 + Ttol){
                         // _printf_("T Surplus");
                         // calculate surplus energy
@@ -1390 +1389 @@
                                 T[i+1] = fmin(CtoK, T[i+1]);
                 
-                                surpT[i+1] = fmax(0.0, exsT[i+1] - 159.1342);
+                                surpT[i+1] = fmax(0.0, exsT[i+1] - 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] = 159.1342;
+                                exsT[i] = LF/CI;
                                 surpE[i] = 0.0;
                                 i = i + 1;
@@ -1401 +1400 @@
 
                 // convert temperature excess to melt [kg]
-                for(int i=0;i<n;i++) M[i] = exsT[i] * d[i] * dz[i] * CI / LF;      // melt
-                sumM = cellsum(M,n);                                               // total melt [kg]
+                for(int i=0;i<n;i++) M[i] = fmin(exsT[i] * d[i] * dz[i] * CI / LF, m[i]);  // melt
+                sumM = cellsum(M,n);                                                       // total melt [kg]
 
                 // calculate maximum refreeze amount, maxF [kg]
@@ -1409 +1408 @@
                 // initialize refreeze, runoff, flxDn and dW vectors [kg]
                 IssmDouble* F = xNewZeroInit<IssmDouble>(n);
-                IssmDouble* R=xNewZeroInit<IssmDouble>(n);
+                IssmDouble* R = xNewZeroInit<IssmDouble>(n);
 
                 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];
+                flxDn=xNewZeroInit<IssmDouble>(n+1);
 
                 // determine the deepest grid cell where melt/pore water is generated
@@ -1442 +1441 @@
                                 m[i] = m[i] - M[i];                     // mass after melt
                                 Wi = (dIce-d[i]) * Swi * (m[i]/d[i]);    // irreducible water 
-                                dW[i] = fmin(inM, Wi - W[i]);            // change in pore water
+                                dW[i] = fmax(fmin(inM, Wi - W[i]),-1*W[i]);            // change in pore water
                                 R[i] = fmax(0.0, inM - dW[i]);             // runoff
+                                F[i] = 0.0;
                         }
                         // check if no energy to refreeze meltwater
@@ -1453 +1453 @@
                                 m[i] = m[i] - M[i];                     // mass after melt
                                 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.0, inM-dW[i]);         // meltwater out
+                                dW[i] = fmax(fmin(inM, Wi - W[i]),-1*W[i]);              // change in pore water
+                                flxDn[i+1] = fmax(0.0, inM - dW[i]);         // meltwater out
+                                R[i] = 0.0;
                                 F[i] = 0.0;                               // no freeze 
                         }
@@ -1462 +1463 @@
                                 // _printf_("MELT REFREEZE");
                                 //-----------------------melt water-----------------------------
+                                m[i] = m[i] - M[i];
                                 IssmDouble dz_0 = m[i]/d[i];          
                                 IssmDouble dMax = (dIce - d[i])*dz_0;              // d max = dIce
@@ -1471 +1473 @@
                                 Wi = (dIce-d[i])* Swi * dz_0;            // irreducible water 
                                 dW[i] = fmin(inM - F1, Wi-W[i]);         // change in pore water
-                                if (-1*dW[i]>W[i]-Wtol ){
+                                if (dW[i] < 0.0-Wtol && -1*dW[i]>W[i]-Wtol ){
                                         dW[i]= -1*W[i];
                                 }
@@ -1482 +1484 @@
                                         m[i] = m[i] + F2;                   // mass after refreeze
                                         d[i] = m[i]/dz_0;
+                                        dW[i] = dW[i] - F2;
                                 }
 
-                                flxDn[i+1] = inM - F1 - dW[i] - F2;     // meltwater out        
+                                F[i] = F1 + F2;
+
+                                flxDn[i+1] = inM - F1 - dW[i];     // meltwater out        
                                 T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
 
@@ -1507 +1512 @@
                 for(int i=0;i<n;i++)W[i] += dW[i];
 
+                //calculate Rsum:
+                Rsum=cellsum(R,n);
+
                 // delete all cells with zero mass
-                D_size=0; for(int i=0;i<n;i++)if(m[i]!=0)D_size++; 
+                D_size=0; for(int i=0;i<n;i++)if(m[i]>0.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++;}
-                
+                D_size=0; for(int i=0;i<n;i++)if(m[i]>0.0){ D[D_size] = i; D_size++;}
+
                 celldelete(&m,n,D,D_size);
                 celldelete(&W,n,D,D_size);
@@ -1522 +1530 @@
                 celldelete(&EI,n,D,D_size);
                 celldelete(&EW,n,D,D_size);
-                celldelete(&R,n,D,D_size);
                 n=D_size;
                 xDelete<int>(D);
@@ -1528 +1535 @@
                 // calculate new grid lengths
                 for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
-
-                //calculate Rsum:
-                Rsum=cellsum(R,n);
 
                 /*Free ressources:*/
@@ -1592 +1596 @@
             
                         // set cell to 99999 for deletion
-                        m[i] = 99999;
+                        m[i] = -99999;
                 }
         }
@@ -1600 +1604 @@
          //find closest cell to merge with
                 for(int i=n-2;i>=0;i--){
-                        if(m[i]!=99999){
+                        if(m[i]!=-99999){
                                 X2=i;
                                 X1=n-1;
@@ -1622 +1626 @@
         
                 // set cell to 99999 for deletion
-                m[X2] = 99999;
+                m[X2] = -99999;
         }
 
         // delete combined cells
-        D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999)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++;}
+        D_size=0; for(int i=0;i<n;i++)if(m[i]!=-99999){ D[D_size] = i; D_size++;}
 
         celldelete(&m,n,D,D_size);
@@ -1723 +1727 @@
                 dz_add=-(dz[n-1]);
         
-                // add a grid cell of the same size and temperature to the bottom
+                // remove a grid cell from the bottom
                 D_size=n-1;
                 D=xNew<int>(D_size);
@@ -1756 +1760 @@
         /*checks: */
         for(int i=0;i<n;i++) if (W[i]<0.0-Wtol) _error_("negative pore water generated in melt equations\n");
-        
+
         /*only in forward mode! avoid round in AD mode as it is not differentiable: */
         #ifndef _HAVE_ADOLC_