Changeset 24689

Timestamp:

04/03/20 11:57:14 (5 years ago)

Author:

schlegel

Message:

CHG: clean up some GEMB code and climatology for last time step, update Archives

Location:

issm/trunk-jpl

Files:

• src/c/modules/SurfaceMassBalancex/Gembx.cpp (modified) (5 diffs)
• test/Archives/Archive243.arch (modified) ( previous)
• test/Archives/Archive244.arch (modified) ( previous)
• test/Archives/Archive252.arch (modified) ( previous)

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

@@ -70 +70 @@
                                 IssmDouble time0     = Ta_input_tr->GetTimeByOffset(-1);
                                 IssmDouble timeend   = Ta_input_tr->GetTimeByOffset(offsetend);
+                                timeend   = ceil(timeend/dt)*dt;
                                 if (time>time0 & timeend>time0){
                                         delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
@@ -111 +112 @@
         #ifndef _HAVE_AD_  //avoid the round operation check!
         if (dgpTop != round(dgpTop)){ 
-                _error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop");
+                _error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop\n");
         }
         #endif
         if(dzTop < 0.05){
-                _printf_("initial top grid cell length (dzTop) is < 0.05 m");
+                _printf_("initial top grid cell length (dzTop) is < 0.05 m\n");
         }
         gpTop=reCast<int,IssmDouble>(dgpTop);
@@ -282 +283 @@
 
         // Take forward differences on left and right edges
-        dT[0] = (T[1] - T[0])/(zGPC[1]-zGPC[0]);
-        if(m>1) dT[m-1] = (T[m-1] - T[m-2])/(zGPC[m-1]-zGPC[m-2]);
+        if(m>1){
+                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]);
+        }
 
         //Take centered differences on interior points
@@ -343 +346 @@
                         // higher than 15 degC m-1 (Brun et al., 1992)
                         //When wet-snow grains (class 6) are submitted to a
-                        // temperature gradient higher than degC m-1, their sphericity
+                        // temperature gradient higher than 5 degC m-1, their sphericity
                         // decreases according to Equations (4). When sphericity
                         // reaches 0, their size increases according to the functions
                         // determined by Marbouty. (Brun et al., 1992)
-                        if(gsp[i]>0.0+Gdntol && gsp[i]<1.0-Gdntol && (dT[i]>15.0+Ttol || (T[i]>5.0+Ttol && W[i]>0.0+Wtol)) ){
-                                //determine coeff0icients
+                        if(gsp[i]>0.0+Gdntol && gsp[i]<1.0-Gdntol && (dT[i]>15.0+Ttol || (dT[i]>5.0+Ttol && W[i]>0.0+Wtol)) ){
+                                //determine coefficients
                                 IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4);
                                 gsp[i] +=C;
                                 if (gsp[i]<0.0+Gdntol)gsp[i]=0.0;
+                                //determine new grain size (mm)
+                                gsz[i] = 0.35 + (0.5-gsp[i])*0.1;
                         }
                         /*Dry snow metamorphism (Marbouty, 1980) grouped model coefficents
                          *from Marbouty, 1980: Figure 9*/
-                        else if(W[i]<=0.0+Wtol || (gsp[i]<=0.0+Gdntol && dT[i]>15.0+Ttol) || (gsp[i]<=0.0+Gdntol && T[i]>5.0+Ttol && W[i]>0.0+Wtol)){
+                        else if(W[i]<=0.0+Wtol || (gsp[i]<=0.0+Gdntol && dT[i]>15.0+Ttol) || (gsp[i]<=0.0+Gdntol && dT[i]>5.0+Ttol && W[i]>0.0+Wtol)){
                                 //_printf_("Nondendritic snow metamorphism\n");
                                 Q = Marbouty(T[i],d[i],dT[i]);
@@ -366 +371 @@
                                 //_printf_("Nondendritic wet snow metamorphism\n");
                                 //wet rate of change coefficient
-                                IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *dts);   // [mm^3 s^-1]
+                                IssmDouble E = (1.28e-8 + 4.22e-10 * pow(lwc[i],3.0))* (dt *dts);   // [mm^3 s^-1]
 
                                 // calculate change in grain volume and convert to grain size