Changeset 24220

Timestamp:

10/14/19 18:31:47 (5 years ago)

Author:

schlegel

Message:

CHG: minor, add in some checks for GEMB

File:

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

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

@@ -239 +239 @@
         // Take forward differences on left and right edges
         dT[0] = (T[1] - T[0])/(zGPC[1]-zGPC[0]);
-        dT[m-1] = (T[m-1] - T[m-2])/(zGPC[m-1]-zGPC[m-2]);
+        if(m>1) dT[m-1] = (T[m-1] - T[m-2])/(zGPC[m-1]-zGPC[m-2]);
 
         //Take centered differences on interior points
@@ -667 +667 @@
         KP = xNew<IssmDouble>(m);
 
-        KU[0] = UNDEF;
-        KD[m-1] = UNDEF;
+        KU[0] = Delflag;
+        KD[m-1] = Delflag;
         for(int i=1;i<m;i++) KU[i]= K[i-1];
         for(int i=0;i<m-1;i++) KD[i] = K[i+1];
@@ -676 +676 @@
         dzU = xNew<IssmDouble>(m);
         dzD = xNew<IssmDouble>(m);
-        dzU[0]=UNDEF;
-        dzD[m-1]=UNDEF;
+        dzU[0]=Delflag;
+        dzD[m-1]=Delflag;
 
         for(int i=1;i<m;i++) dzU[i]= dz[i-1];
@@ -1401 +1401 @@
                 for(int i=0;i<n;i++) m[i] += dW[i];                                            // new mass [kg]
                 for(int i=0;i<n;i++) d[i] = m[i] / dz[i];                                    // density [kg m-3]   
-                for(int i=0;i<n;i++) T[i] = T[i] + (dW[i]*(LF+(CtoK - T[i])*CI)/(m[i]*CI));      // temperature [K]
+                for(int i=0;i<n;i++) if(m[i]>Wtol) T[i] = T[i] + (dW[i]*(LF+(CtoK - T[i])*CI)/(m[i]*CI));      // temperature [K]
 
                 // if pore water froze in ice then adjust d and dz thickness
@@ -1433 +1433 @@
 
                         int i = 0;
-                        while (cellsum(surpE,n) > 0.0+Ttol){
+                        while (cellsum(surpE,n) > 0.0+Ttol && i<n){
                                 // use surplus energy to increase the temperature of lower cell
                                 T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];
@@ -1451 +1451 @@
 
                 // convert temperature excess to melt [kg]
-                for(int i=0;i<n;i++) M[i] = min(exsT[i] * d[i] * dz[i] * CI / LF, m[i]);  // melt
+                double Mmax=0.0;
+                for(int i=0;i<n;i++){
+                        Mmax=exsT[i] * d[i] * dz[i] * CI / LF;
+                        M[i] = min(Mmax, m[i]);  // melt
+                }
                 sumM = cellsum(M,n);                                                       // total melt [kg]
 
@@ -1511 +1515 @@
                         // some or all meltwater refreezes
                         else{
-                                // change in density density and temperature
+                                // change in density and temperature
                                 // _printf_("MELT REFREEZE");
                                 //-----------------------melt water-----------------------------
@@ -1541 +1545 @@
 
                                 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
+                                if (m[i]>Wtol){
+                                        T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
+                                }
 
                                 // check if an ice layer forms 
@@ -2003 +2009 @@
 
                 // do not allow densities to exceed the density of ice
-                if(d[i] > dIce-Dtol) d[i]=dIce;
+                if(d[i] > dIce-Ptol) d[i]=dIce;
 
                 // calculate new grid cell length