Changeset 25993

Timestamp:

02/15/21 13:49:48 (4 years ago)

Author:

schlegel

Message:

CHG: add options for Alex's newest Albedo scheme for GEMB

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Element.cpp (modified) (4 diffs)
• modules/SurfaceMassBalancex/Gembx.cpp (modified) (3 diffs)
• modules/SurfaceMassBalancex/SurfaceMassBalancex.h (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -3551 +3551 @@
         IssmDouble teValue=1.0;
         IssmDouble aValue=0.0;
+        IssmDouble szaValue=0.0;
+        IssmDouble cotValue=0.0;
+        IssmDouble ccsnowValue=0.0;
+        IssmDouble cciceValue=0.0;
         IssmDouble dt,time,smb_dt;
         int        aIdx=0;
@@ -3833 +3837 @@
         teValue_input->GetInputValue(&teValue,gauss);  // Emissivity [0-1]
         aValue_input->GetInputValue(&aValue,gauss);  // Albedo [0 1]
+        //Not implemented as input
+        //szaValue;  // Solar Zenith Angle [degree]
+        //cotValue;  // Cloud Optical Thickness
+        //ccsnowValue;//concentration of light absorbing carbon for snow [ppm1]
+        //cciceValue;//concentration of light absorbing carbon for ice [ppm1]
         //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
         /*}}}*/
@@ -3840 +3849 @@
 
         /*Snow, firn and ice albedo:*/
-        if(isalbedo)albedo(&a,aIdx,re,dz,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,Msurf,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
+        if(isalbedo)albedo(&a,aIdx,re,dz,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,Msurf,ccsnowValue,cciceValue,szaValue,cotValue,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
 
         /*Distribution of absorbed short wave radation with depth:*/
@@ -3927 +3936 @@
         /*Check mass conservation:*/
         if (dMass != 0.0){
-                _printf_("total system mass not conserved in MB function \n");
+                _printf_("Time: " << setprecision(8) << timeinputs/365.0/24.0/3600.0 << ",total system mass not conserved in MB function \n");
+                _printf_("sumMass: " << sumMass << " sumR: " << sumR << " sumW: " << sumW << " sumP: " << sumP << " sumEC: " << sumEC << " initMass: " << initMass << " sumMassAdd: " << sumMassAdd << " \n");
         }
         #endif

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

@@ -410 +410 @@
 
 }  /*}}}*/
-void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
+void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
 
         //// Calculates Snow, firn and ice albedo as a function of:
@@ -432 +432 @@
         // Methods 1 & 2
         //   re      = surface effective grain radius [mm]
+        // Method 1, optional
+        //  clabSnow = concentration of light absorbing carbon  [ppm1], default 0
+        //  SZA      = solar zenith angle of the incident radiation [deg], default 0
+        //  COT      = cloud optical thickness, default 0
+        //  For TWO LAYER
+        //  clabIce  = concentration of light absorbing carbon of first ice layer [ppm1], default 0
 
         // Methods 3
@@ -485 +491 @@
                         //function of effective grain radius
 
+                        //gardnerAlb(S1, c1, SZA, t, z1, S2, c2)
+                        //This is an implementation of the snow and ice broadband albedo
+                        //  parameterization developed by Alex Gardner.
+                        //Created By: Alex S. Gardner, Jet Propulsion Laboratory [alex.s.gardner@jpl.nasa.gov]
+                        //  Last Modified: June, 2014
+                        //Full Reference: Gardner, A. S., and Sharp, M. J.: A review of snow and
+                        //  ice albedo and the development of a new physically based broadband albedo
+                        //  parameterization, J. Geophys. Res., 115, F01009, 10.1029/2009jf001444,
+                        //  2010.
+
+                        //INPUTS
+                        // ONE LAYER
+                        //  - S1    : specific surface area of the snow or ice [cm^2 g-1]
+                        //  - c1    : concentration of light absorbing carbon  [ppm1]
+                        //  - SZA   : solar zenith angle of the incident radiation [deg]
+                        //  - t     : cloud optical thickness
+                        // TWO LAYER
+                        //  - z1    : depth of snow suface layer [mm w.e.]
+                        //  - S2    : specific surface area of bottom ice layer [cm^2 g-1]
+                        //  - c2    : concentration of light absorbing carbon of bottom ice
+                        //             layer [ppm1]
+                        IssmDouble c1=clabSnow;
+                        IssmDouble c2=clabIce;
+                        IssmDouble t=COT;
+
+                        //Single layer albedo parameterization
                         //convert effective radius to specific surface area [cm2 g-1]
-                        IssmDouble S = 3.0 / (0.091 * re[0]);
-
-                        //determine broadband albedo
-                        a[0]= 1.48 - pow(S,-0.07);
+                        IssmDouble S1 = 3.0 / (0.091 * re[0]);
+
+                        //effective solar zenith angle
+                        IssmDouble x = min(pow(t/(3*cos(Pi*SZA/180.0)),0.5), 1.0);
+                        IssmDouble u = 0.64*x + (1-x)*cos(Pi*SZA/180.0);
+
+                        // pure snow albedo
+                        IssmDouble as = 1.48 - pow(S1,-0.07);
+
+                        //change in pure snow albedo due to soot loading
+                        IssmDouble dac = max(0.04 - as, pow(-c1,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c1,0.6))*pow(S1,-0.25)));
+
+                        //Two layer albedo parameterization
+                        //  do two layer calculation if there is more than 1 layer
+                        IssmDouble z1=0.0;
+                        int lice=0;
+                        for(int l=0;(l<m & d[l]<dPHC-Dtol);l++){
+                                z1=z1+dz[l]*d[l]; //mm
+                                lice=l+1;
+                        }
+                        if (m>0 & lice<m & z1 > Dtol){
+                                // determine albedo values for bottom layer
+                                IssmDouble S2 = 3.0 / (0.091 * re[lice]);
+
+                                // pure snow albedo
+                                IssmDouble as2 = 1.48 - pow(S2,-0.07);
+
+                                // change in pure snow albedo due to soot loading
+                                IssmDouble dac2 = max(0.04 - as2, pow(-c2,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c2,0.6))*pow(S2,-0.25)));
+
+                                // determine the effective change due to finite depth and soot loading
+                                IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1-as) - 0.1*c1 - 0.13)));
+
+                                dac =  (as2 + dac2 - as) + A*((as + dac) - (as2 + dac2));
+                        }
+
+                        // change in albedo due to solar zenith angle
+                        IssmDouble dasz = 0.53*as*(1 - (as + dac))*pow(1 - u,1.2);
+
+                        // change in albedo due to cloud (apart from change in diffuse fraction)
+                        IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1 + 1.5*t,as));
+
+                        // Broadband albedo
+                        a[0] = as + dac + dasz + dat;
                 }
                 else if(aIdx==2){

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

@@ -29 +29 @@
 IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT);
 void grainGrowth(IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx, int sid);
-void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
+void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
 void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid);
 void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT);