Changeset 23189 for issm/trunk/src/c/modules/SurfaceMassBalancex/Gembx.cpp

Timestamp:

08/28/18 09:45:51 (7 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 23187

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/c/modules/SurfaceMassBalancex/Gembx.cpp (modified) (66 diffs)

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

@@ -24 +24 @@
 const double LS = 2.8295E6;             // latent heat of sublimation (J kg-1)
 const double SB = 5.67E-8;                // Stefan-Boltzmann constant (W m-2 K-4)
-const double CA = 1005.0;                    // heat capacity of air (J kg-1 k-1)
-const double R = 8.314;                      // gas constant (mol-1 K-1)
+const double CA = 1005.0;                    // heat capacity of air (J kg-1 K-1)
+const double R = 8.314;                      // gas constant (J mol-1 K-1)
 
 void Gembx(FemModel* femmodel){  /*{{{*/
@@ -107 +107 @@
         xDelete(dzT);
         xDelete(dzB);
-        
+
         //---------NEED TO IMPLEMENT A PROPER GRID STRECHING ALGORITHM------------
 
@@ -202 +202 @@
         gdn=*pgdn;
         gsp=*pgsp;
-        
+
         /*only when aIdx = 1 or 2 do we run grainGrowth: */
         if(aIdx!=1 & aIdx!=2){
@@ -220 +220 @@
         //set maximum water content by mass to 9 percent (Brun, 1980)
         for(int i=0;i<m;i++)if(lwc[i]>9.0-Wtol) lwc[i]=9.0;
-
 
         /* Calculate temperature gradiant across grid cells 
@@ -244 +243 @@
         // take absolute value of temperature gradient
         for(int i=0;i<m;i++)dT[i]=fabs(dT[i]);
-        
+
         /*Snow metamorphism. Depends on value of snow dendricity and wetness of the snowpack: */
         for(int i=0;i<m;i++){
@@ -324 +323 @@
                 re[i]=gsz[i]/2.0;
         }
-        
+
         /*Free ressources:*/
         xDelete<IssmDouble>(gsz);
@@ -348 +347 @@
         //// Inputs
         // aIdx      = albedo method to use
-        
+
         // Method 0
         //  aValue   = direct input value for albedo, override all changes to albedo
@@ -397 +396 @@
         /*Recover pointers: */
         a=*pa;
-        
+
         //some constants:
         const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
@@ -513 +512 @@
 
         /* ENGLACIAL THERMODYNAMICS*/
-         
+
         /* Description: 
            computes new temperature profile accounting for energy absorption and 
@@ -537 +536 @@
         // T: grid cell temperature [k]
         // EC: evaporation/condensation [kg]
-        
+
         /*intermediary: */
         IssmDouble* K = NULL;
@@ -580 +579 @@
         IssmDouble dlw=0.0;
         IssmDouble dT_dlw=0.0;
-        
+
         /*outputs:*/
         IssmDouble EC=0.0;
@@ -597 +596 @@
         //initialize Evaporation - Condenstation 
         EC = 0.0;
-        
+
         // check if all SW applied to surface or distributed throught subsurface
         // swIdx = length(swf) > 1
@@ -609 +608 @@
         // if V = 0 goes to infinity therfore if V = 0 change
         if(V<0.01-Dtol)V=0.01;
-        
+
         // Bulk-transfer coefficient for turbulent fluxes
         An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
@@ -627 +626 @@
 
         // THERMAL DIFFUSION COEFFICIENTS
- 
+
         // A discretization scheme which truncates the Taylor-Series expansion
         // after the 3rd term is used. See Patankar 1980, Ch. 3&4
- 
+
         // discretized heat equation:
- 
+
         //                 Tp = (Au*Tu° + Ad*Td° + (Ap-Au-Ad)Tp° + S) / Ap
- 
+
         // where neighbor coefficients Au, Ap, & Ad are
- 
+
         //                   Au = [dz_u/2KP + dz_p/2KE]^-1
         //                   Ad = [dz_d/2KP + dz_d/2KD]^-1 
@@ -649 +648 @@
         KD = xNew<IssmDouble>(m);
         KP = xNew<IssmDouble>(m);
-        
+
         KU[0] = UNDEF;
         KD[m-1] = UNDEF;
@@ -661 +660 @@
         dzU[0]=UNDEF;
         dzD[m-1]=UNDEF;
-        
+
         for(int i=1;i<m;i++) dzU[i]= dz[i-1];
         for(int i=0;i<m-1;i++) dzD[i] = dz[i+1];
@@ -693 +692 @@
                 Ap[i] = (d[i]*dz[i]*CI)/dt;
         }
-        
+
         // create "neighbor" coefficient matrix
         Nu = xNew<IssmDouble>(m);
@@ -703 +702 @@
                 Np[i]= 1.0 - Nu[i] - Nd[i];
         }
-        
+
         // specify boundary conditions: constant flux at bottom
         Nu[m-1] = 0.0;
         Np[m-1] = 1.0;
-        
+
         // zero flux at surface
         Np[0] = 1.0 - Nd[0];
-        
+
         // Create neighbor arrays for diffusion calculations instead of a tridiagonal matrix
         Nu[0] = 0.0;
         Nd[m-1] = 0.0;
-        
+
         /* RADIATIVE FLUXES*/
 
@@ -720 +719 @@
         sw = xNew<IssmDouble>(m);
         for(int i=0;i<m;i++) sw[i]= swf[i] * dt;
-        
+
         // temperature change due to SW
         dT_sw = xNew<IssmDouble>(m);
@@ -746 +745 @@
                 // store initial temperature
                 //T_init = T;
-    
+
                 // calculate temperature of snow surface (Ts)
                 // when incoming SW radition is allowed to penetrate the surface,
@@ -754 +753 @@
                 Ts = (T[0] + T[1])/2.0;
                 Ts = fmin(CtoK,Ts);    // don't allow Ts to exceed 273.15 K (0 degC)
-                
+
                 //TURBULENT HEAT FLUX
-    
+
                 // Monin-Obukhov Stability Correction
                 // Reference:
@@ -764 +763 @@
                 // calculate the Bulk Richardson Number (Ri)
                 Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
-                
+
                 // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
-    
+
                 // do not allow Ri to exceed 0.19
                 Ri = fmin(Ri, 0.19);
@@ -778 +777 @@
                         coefM =pow (1.0-18*Ri,-0.25);
                 }
-                
+
                 // calculate heat/wind 'coef_H' stability factor
                 if (Ri <= -0.03+Ttol) coefH = coefM/1.3;
                 else coefH = coefM;
-                
+
                 //// Sensible Heat
                 // calculate the sensible heat flux [W m-2](Patterson, 1998)
@@ -801 +800 @@
                 else{
                         L = LS; // latent heat of sublimation 
-                        
+
                         // for an ice surface Murphy and Koop, 2005 [Equation 7]
                         eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
@@ -823 +822 @@
                 ulw = - (SB * pow(Ts,4.0)* teValue) * dt ;
                 dT_ulw = ulw / TCs;
-                
+
                 // new grid point temperature
-    
+
                 //SW penetrates surface
                 for(int j=0;j<m;j++) T[j] = T[j] + dT_sw[j];
                 T[0] = T[0] + dT_dlw + dT_ulw + dT_turb;
-                
+
                 // temperature diffusion
                 for(int j=0;j<m;j++) T0[1+j]=T[j];
@@ -835 +834 @@
                 for(int j=0;j<m;j++) Td[j] = T0[2+j];
                 for(int j=0;j<m;j++) T[j] = (Np[j] * T[j]) + (Nu[j] * Tu[j]) + (Nd[j] * Td[j]);
-                
+
                 // calculate cumulative evaporation (+)/condensation(-)
                 EC = EC + (EC_day/dts)*dt;
@@ -873 +872 @@
         xDelete<IssmDouble>(Td);
 
-
         /*Assign output pointers:*/
         *pEC=EC;
@@ -901 +899 @@
         //   swf     = absorbed shortwave radiation [W m-2]
         //
-        
+
         /*outputs: */
         IssmDouble* swf=NULL;
@@ -910 +908 @@
         swf=xNewZeroInit<IssmDouble>(m);
 
-
         // SHORTWAVE FUNCTION
         if (swIdx == 0) {// all sw radation is absorbed in by the top grid cell
-        
+
                 // calculate surface shortwave radiation fluxes [W m-2]
                 swf[0] = (1.0 - as) * dsw;
@@ -948 +945 @@
                         }
 
-
                         // spectral albedos:
                         // 0.3 - 0.8um
@@ -988 +984 @@
                                 B2_cum[i+1]=cum2;
                         }
-
 
                         // flux across grid cell boundaries
@@ -1015 +1010 @@
                         xDelete<IssmDouble>(Qs1);
                         xDelete<IssmDouble>(Qs2);
-                        
-                        
+
                 }
                 else{  //function of grid cell density
@@ -1144 +1138 @@
 
         if (P > 0.0+Dtol){
-                        
 
                 if (T_air <= CtoK+Ttol){ // if snow
@@ -1210 +1203 @@
 
                 mass_diff = mass - massinit - P;
-                
+
                 #ifndef _HAVE_ADOLC_  //avoid round operation. only check in forward mode.
                 mass_diff = round(mass_diff * 100.0)/100.0;
@@ -1253 +1246 @@
         IssmDouble* M=NULL;
         int*       D=NULL;
-        
+
         IssmDouble sumM=0.0;
         IssmDouble sumER=0.0;
@@ -1265 +1258 @@
         IssmDouble X=0.0;
         IssmDouble Wi=0.0;
-    
+
         IssmDouble Ztot=0.0;
         IssmDouble T_bot=0.0;
@@ -1279 +1272 @@
         IssmDouble EW_bot=0.0;
         bool        top=false;
-    
+
         IssmDouble* Zcum=NULL;
         IssmDouble* dzMin2=NULL;
@@ -1286 +1279 @@
         int X1=0;
         int X2=0;
-    
+
         int        D_size = 0;
 
@@ -1303 +1296 @@
         int         n=*pn;
         IssmDouble* R=NULL;
-        
+
         if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   melt module\n");
 
@@ -1335 +1328 @@
         // 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]
         const IssmDouble Swi = 0.07;                     // assumed constant after Colbeck, 1974
@@ -1493 +1486 @@
                                 T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
 
-
                                 // check if an ice layer forms 
                                 if (fabs(d[i] - dIce) < Dtol){
@@ -1505 +1497 @@
                 }
 
-
                 //// GRID CELL SPACING AND MODEL DEPTH
                 for(int i=0;i<n;i++)if (W[i] < 0.0 - Wtol) _error_("negative pore water generated in melt equations");
-                
+
                 // delete all cells with zero mass
                 // adjust pore water
@@ -1533 +1524 @@
                 n=D_size;
                 xDelete<int>(D);
-        
+
                 // calculate new grid lengths
                 for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
@@ -1545 +1536 @@
         Zcum=xNew<IssmDouble>(n);
         dzMin2=xNew<IssmDouble>(n);
-    
+
         Zcum[0]=dz[0]; // Compute a cumulative depth vector
-    
+
         for (int i=1;i<n;i++){
                 Zcum[i]=Zcum[i-1]+dz[i];
         }
-    
+
         for (int i=0;i<n;i++){
                 if (Zcum[i]<=zTop+Dtol){
@@ -1569 +1560 @@
                 }
         }
-    
+
         //Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell)
         if(X==n-1){
@@ -1589 +1580 @@
                         gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
                         gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
-            
+
                         // merge with underlying grid cell and delete old cell
                         dz[i+1] = dz[i] + dz[i+1];                 // combine cell depths
@@ -1595 +1586 @@
                         W[i+1] = W[i+1] + W[i];                     // combine liquid water
                         m[i+1] = m_new;                             // combine top masses
-            
+
                         // set cell to 99999 for deletion
                         m[i] = -99999;
@@ -1611 +1602 @@
                         }
                 }
-        
+
                 // adjust variables as a linearly weighted function of mass
                 IssmDouble m_new = m[X2] + m[X1];
@@ -1619 +1610 @@
                 gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
                 gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
-        
+
                 // merge with underlying grid cell and delete old cell
                 dz[X1] = dz[X2] + dz[X1];                 // combine cell depths
@@ -1625 +1616 @@
                 W[X1] = W[X1] + W[X2];                     // combine liquid water
                 m[X1] = m_new;                             // combine top masses
-        
+
                 // set cell to 99999 for deletion
                 m[X2] = -99999;
@@ -1648 +1639 @@
         n=D_size;
         xDelete<int>(D);
-    
+
         // check if any of the top 10 cell depths are too large
         X=0;
@@ -1657 +1648 @@
                 }
         }
-        
+
         int j=0;
         while(j<=X){
@@ -1686 +1677 @@
         // calculate total model depth
         Ztot = cellsum(dz,n);
-    
+
         if (Ztot < zMin-Dtol){
                 // printf("Total depth < zMin %f \n", Ztot);
@@ -1692 +1683 @@
                 mAdd = m[n-1]+W[n-1];
                 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];
@@ -1705 +1696 @@
                 EI_bot=EI[n-1];
                 EW_bot=EW[n-1];
-        
+
                 dz_add=dz_bot;
-        
+
                 newcell(&dz,dz_bot,top,n);
                 newcell(&T,T_bot,top,n);
@@ -1727 +1718 @@
                 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
                 D_size=n-1;
                 D=xNew<int>(D_size);
-        
+
                 for(int i=0;i<n-1;i++) D[i]=i;
                 celldelete(&dz,n,D,D_size);
@@ -1768 +1759 @@
                         << "dm: " << dm << " dE: " << dE << "\n");
         #endif
-        
+
         /*Free ressources:*/
         if(m)xDelete<IssmDouble>(m);
@@ -1784 +1775 @@
         xDelete<IssmDouble>(Zcum);
         xDelete<IssmDouble>(dzMin2);
-    
+
         /*Assign output pointers:*/
         *pM=sumM;
@@ -1854 +1845 @@
         IssmDouble c1=0.0;
         IssmDouble H=0.0;
+        IssmDouble M0=0.0;
+        IssmDouble M1=0.0;
+        IssmDouble c0arth=0.0;
+        IssmDouble c1arth=0.0;
 
         /*output: */
         IssmDouble* dz=NULL;
         IssmDouble* d=NULL;
-        
+
         if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   densification module\n");
 
@@ -1867 +1862 @@
         // initial mass
         IssmDouble* mass_init = xNew<IssmDouble>(m);for(int i=0;i<m;i++) mass_init[i]=d[i] * dz[i];
-        
+
         /*allocations and initialization of overburden pressure and factor H: */
         IssmDouble* cumdz = xNew<IssmDouble>(m-1);
@@ -1876 +1871 @@
         obp[0]=0.0;
         for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];
-        
+
         // calculate new snow/firn density for:
         //   snow with densities <= 550 [kg m-3]
         //   snow with densities > 550 [kg m-3]
-                
-        
+
         for(int i=0;i<m;i++){
                 switch (denIdx){
                         case 1: // Herron and Langway (1980)
-                                c0 = (11.0 * exp(-10160.0 / (T[i] * 8.314))) * C/1000.0;
-                                c1 = (575.0 * exp(-21400.0 / (T[i]* 8.314))) * pow(C/1000.0,.5);
+                                c0 = (11.0 * exp(-10160.0 / (T[i] * R))) * C/1000.0;
+                                c1 = (575.0 * exp(-21400.0 / (T[i]* R))) * pow(C/1000.0,.5);
                                 break;
+
                         case 2: // Arthern et al. (2010) [semi-emperical]
                                 // common variable
                                 // NOTE: Ec=60000, Eg=42400 (i.e. should be in J not kJ)
-                                H = exp((-60000.0/(T[i] * 8.314)) + (42400.0/(Tmean * 8.314))) * (C * 9.81);
+                                H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
                                 c0 = 0.07 * H;
                                 c1 = 0.03 * H;
@@ -1899 +1894 @@
 
                                 // common variable
-                                H = exp((-60.0/(T[i] * 8.314))) * obp[i] / pow(re[i]/1000.0,2.0);
+                                H = exp((-60000.0/(T[i] * R))) * obp[i] / pow(re[i]/1000.0,2.0);
                                 c0 = 9.2e-9 * H;
                                 c1 = 3.7e-9 * H;
@@ -1913 +1908 @@
                                 c0 = (C/dIce) * (76.138 - 0.28965*Tmean)*8.36*pow(CtoK - T[i],-2.061);
                                 c1 = c0;
+                                break;
+
+                        case 6: // Ligtenberg and others (2011) [semi-emperical], Antarctica 
+                                // common variable
+                                H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
+                                c0arth = 0.07 * H;
+                                c1arth = 0.03 * H;
+                                M0 = fmax(1.435 - (0.151 * log(C)),0.25);
+                                M1 = fmax(2.366 - (0.293 * log(C)),0.25);
+                                c0 = M0*c0arth;
+                                c1 = M1*c1arth;
+                                break;
+
+                        case 7: // Kuipers Munneke and others (2015) [semi-emperical], Greenland
+                                // common variable
+                                H = exp((-60000.0/(T[i] * R)) + (42400.0/(T[i] * R))) * (C * 9.81);
+                                c0arth = 0.07 * H;
+                                c1arth = 0.03 * H;
+                                M0 = fmax(1.042 - (0.0916 * log(C)),0.25);
+                                M1 = fmax(1.734 - (0.2039 * log(C)),0.25);
+                                c0 = M0*c0arth;
+                                c1 = M1*c1arth;
                                 break;
                 }
@@ -1977 +1994 @@
         IssmDouble lhf=0.0;
         IssmDouble EC=0.0;
-        
+
         if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   turbulentFlux module\n");
 
@@ -2002 +2019 @@
 
         // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
-        
+
         // do not allow Ri to exceed 0.19
         Ri = fmin(Ri, 0.19);
@@ -2046 +2063 @@
         }
 
-
         // Latent heat flux [W m-2]
         lhf = C * L * (eAir - eS) * 0.622 / pAir;