more work on CanopyFluxes; removed some unnecessary variables

src/BareGroundFluxes.cc

@@ -167,8 +167,6 @@ forc_rho                   [double] density (kg/m**3)
 soilbeta                   [double] soil wetness relative to field capacity
 dqgdT                      [double] temperature derivative of "qg"
 htvp                       [double] latent heat of vapor of water (or sublimation) [j/kg]
-forc_u                     [double] atmospheric wind speed in east direction (m/s)
-forc_v                     [double] atmospheric wind speed in north direction (m/s)
 forc_q                     [double] atmospheric specific humidity (kg/kg)
 thm                        [double] intermediate variable (forc_t+0.0098*forc_hgt_t_patch)
 t_h2osfc                   [double] surface water temperature
@@ -178,12 +176,9 @@ qg_h2osfc                  [double] specific humidity at h2osfc surface [kg/kg]
 t_soisno[nlevgrnd+nlevsno] [double] col soil temperature (Kelvin)
 
 OUTPUTS:
-ram1                       [double] aerodynamical resistance (s/m)
 cgrnds                     [double] deriv, of soil sensible heat flux wrt soil temp [w/m2/k]
 cgrndl                     [double] deriv of soil latent heat flux wrt soil temp [w/m**2/k]
 cgrnd                      [double] deriv. of soil energy flux wrt to soil temp [w/m2/k]
-taux                       [double] wind (shear) stress: e-w (kg/m/s**2)
-tauy                       [double] wind (shear) stress: n-s (kg/m/s**2)
 eflx_sh_grnd               [double] sensible heat flux from ground (W/m**2) [+ to atm]
 eflx_sh_tot                [double] total sensible heat flux (W/m**2) [+ to atm]
 eflx_sh_snow               [double] sensible heat flux from snow (W/m**2) [+ to atm]
@@ -203,8 +198,6 @@ qflx_ev_h2osfc             [double] evaporation flux from h2osfc (W/m**2) [+ to
     const double& soilbeta,
     const double& dqgdT,
     const double& htvp,
-    const double& forc_u,
-    const double& forc_v,
     const double& forc_q,
     const double& thm,
     const double& t_h2osfc,
@@ -212,12 +205,9 @@ qflx_ev_h2osfc             [double] evaporation flux from h2osfc (W/m**2) [+ to
     const double& qg_soil,
     const double& qg_h2osfc,
     const double *t_soisno,
-    double& ram1,
     double& cgrnds,  
     double& cgrndl,
     double& cgrnd,
-    double& taux,
-    double& tauy,
     double& eflx_sh_grnd,
     double& eflx_sh_tot,
     double& eflx_sh_snow,
@@ -231,9 +221,8 @@ qflx_ev_h2osfc             [double] evaporation flux from h2osfc (W/m**2) [+ to
     
     if (!Land.lakpoi && !Land.urbpoi && frac_veg_nosno == 0) {
       
-      double ram, rah, raw, raih, raiw;
+      double rah, raw, raih, raiw;
       // Determine aerodynamic resistances
-      ram  = 1.0 / (ustar * ustar / um);
       rah  = 1.0 / (temp1 * ustar);
       raw  = 1.0 / (temp2 * ustar);
       raih = forc_rho * cpair / rah;
@@ -245,7 +234,6 @@ qflx_ev_h2osfc             [double] evaporation flux from h2osfc (W/m**2) [+ to
         // Lee and Pielke 1992 beta is applied
         raiw = soilbeta * forc_rho / raw;
       }
-      ram1 = ram;  // pass value to global variable (driver)
   
       // Output to pft-level data structures
       // Derivative of fluxes with respect to ground temperature
@@ -255,8 +243,6 @@ qflx_ev_h2osfc             [double] evaporation flux from h2osfc (W/m**2) [+ to
   
       // Surface fluxes of momentum, sensible and latent heat
       // using ground temperatures from previous time step
-      taux          = -forc_rho * forc_u / ram;
-      tauy          = -forc_rho * forc_v / ram;
       eflx_sh_grnd  = -raih * dth;
       eflx_sh_tot   = eflx_sh_grnd;
   

src/CMakeLists.txt

-add_library (physics physics_dummy.cc CanopyHydrology.cc coldstart_topography.cc coldstart_snow.cc SurfaceRadiation.cc SurfaceResistance.cc QSat.cc CanopyTemperature.cc FrictionVelocity.cc BareGroundFluxes.cc SoilMoistStress.cc)
+add_library (physics physics_dummy.cc CanopyHydrology.cc coldstart_topography.cc coldstart_snow.cc SurfaceRadiation.cc SurfaceResistance.cc QSat.cc CanopyTemperature.cc FrictionVelocity.cc BareGroundFluxes.cc CanopyFluxes.cc)
 target_include_directories (physics PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})

src/CanopyHydrology.cc

@@ -276,8 +276,8 @@ OUTPUTS:
 snow_depth                   [double] snow height (m)
 h2osno                       [double] snow water (mm H2O)
 int_snow                     [double] integrated snowfall [mm]
-h2osoi_liq[nlevgrnd+nlevsno] [double] ice lens (kg/m2)                       
-h2osoi_ice[nlevgrnd+nlevsno] [double] liquid water (kg/m2)
+h2osoi_liq[nlevgrnd+nlevsno] [double] liquid water (kg/m2)
+h2osoi_ice[nlevgrnd+nlevsno] [double] ice lens (kg/m2)
 t_soisno[nlevgrnd+nlevsno]   [double] soil temperature (Kelvin)
 frac_iceoldnlevgrnd+nlevsno] [double] fraction of ice relative to the tot water
 snl                          [int] number of snow layers

src/CanopyTemperature.cc

@@ -303,7 +303,6 @@ potential temp and wind speed.
 
 INPUTS:
 Land                         [LandType] struct containing information about landtype
-vtype                        [int] vegetation type
 snl                          [int] number of snow layers
 frac_sno                     [double] fraction of ground covered by snow (0 to 1)
 forc_th                      [double] atmospheric potential temperature (Kelvin)
@@ -337,7 +336,6 @@ cgrndl                       [double] deriv. of soil latent heat flux wrt soil t
 */
   void GroundProperties(
     const LandType& Land,
-    const bool& vtype,
     const int& snl,
     const double& frac_sno,
     const double& forc_th,
@@ -399,8 +397,8 @@ cgrndl                       [double] deriv. of soil latent heat flux wrt soil t
       }
       z0hg = z0mg;            // initial set only
       z0qg = z0mg;            // initial set only
-      z0m = z0mr[vtype] * htop;
-      displa = displar[vtype] * htop;
+      z0m = z0mr[Land.vtype] * htop;
+      displa = displar[Land.vtype] * htop;
   
       // vegetation roughness lengths
       z0mv = z0m;

src/RootBioPhys.cc

+#include <cmath>
+#include <assert.h>
+#include "clm_constants.h"
+
+
+
+/*
+DESCRIPTION:
+compute root profile for soil water uptake using equation from Zeng 2001, J. Hydrometeorology
+
+INPUTS:
+vtype [int] vegetation type
+roota_par [double] CLM rooting distribution parameter [1/m]
+rootb_par [double] CLM rooting distribution parameter [1/m]
+
+
+OUTPUT:
+rootfr_road_perv[nlevgrnd]   [double] fraction of roots in each soil layer
+*/
+void init_vegrootfr(const int& vtype, double rootfr[nlevgrnd]) {
+  // (computing from surface, d is depth in meter): Y = 1 -1/2 (exp(-ad)+exp(-bd) under the constraint that
+  // Y(d =0.1m) = 1-beta^(10 cm) and Y(d=d_obs)=0.99 with beta & d_obs given in Zeng et al. (1998).
+  if (vtype != noveg) {
+    double totrootfr = 0.0;
+    for (int i = 0; i < nlevsoi - 1; i++) {
+      rootfr[i] = 0.5 * (exp(-roota_par[vtype] * zi[i+nlevsno]) + exp(-rootb_par[vtype] * zi[i+nlevsno]) -
+                  exp(-roota_par[vtype] * zi[i+1+nlevsno]) - exp(-rootb_par[vtype] * zi[i+1+nlevsno]));
+      if (i < nlev2bed) { totrootfr += rootfr[i]; }
+    }
+    rootfr[nlevsoi-1] = 0.5 * (exp(-roota_par[vtype] * zi[nlevsoi-1+nlevsno]) +
+                   exp(-rootb_par[vtype] * zi[nlevsoi-1+nlevsno]));
+    // Adjust layer root fractions if nlev2bed < nlevsoi
+    if (use_var_soil_thick && nlev2bed < nlevsoi) {
+      for (int i = 0; i < nlev2bed; i++) { rootfr[i] /= totrootfr; }
+      for (int i = nlev2bed; i < nlevsoi; i++) { rootfr[i] = 0.0; }
+    }
+  } else {
+    for (int i = 0; i < nlevsoi; i++) { rootfr[i] = 0.0; }
+  }
+  for (int i = nlevsoi; i < nlevgrnd; i++) { rootfr[i] = 0.0; }
+}

src/SoilMoistStress.cc → src/SoilMoistStress.hh

@@ -22,7 +22,7 @@ void array_normalization (double *arr_inout)
 
   for (int i = 0; i < nlevgrnd; i++) { arr_sum += arr_inout[i]; }
   for (int i = 0; i < nlevgrnd; i++) {
-    if (arr_sum > 0.0) { arr_inout[i] = arr_inout[i] / arr_sum; }
+    if (arr_sum > 0.0) { arr_inout[i] /= arr_sum; }
   }
 }
 
@@ -88,7 +88,7 @@ DESCRIPTION: compute the effective soil porosity
 
 INPUTS:
 watsat[nlevgrnd]             [double] volumetric soil water at saturation (porosity)
-h2osoi_ice[nlevgrnd+nlevsno] [double] liquid water (kg/m2)
+h2osoi_ice[nlevgrnd+nlevsno] [double] ice lens (kg/m2)
 dz[nlevgrnd+nlevsno]         [double] layer thickness (m)
 
 OUTPUTS:
@@ -115,7 +115,7 @@ DESCRIPTION: compute the volumetric liquid water content
 
 INPUTS:
 eff_porosity[nlevgrnd]       [double] effective soil porosity
-h2osoi_ice[nlevgrnd+nlevsno] [double] ice lens (kg/m2)
+h2osoi_liq[nlevgrnd+nlevsno] [double] liquid water (kg/m2)
 dz[nlevgrnd+nlevsno]         [double] layer thickness (m)
 
 OUTPUTS:
@@ -136,45 +136,47 @@ void calc_volumetric_h2oliq(
 
 /*
 INPUTS:
-h2osoi_liqvol[nlevgrnd+nlevsno] [double] liquid volumetric moisture, will be used for BeTR
+h2osoi_liqvol[nlevgrnd+nlevsno] [double] liquid volumetric moisture
 rootfr[nlevgrnd]                [double] fraction of roots in each soil layer
-rootfr_unf[nlevgrnd]            [double] root fraction defined for unfrozen layers only
 t_soisno[nlevgrnd+nlevsno]      [double] col soil temperature (Kelvin)
 tc_stress                       [double] critical soil temperature for soil water stress (C)
 sucsat[nlevgrnd]                [double] minimum soil suction (mm)
 watsat[nlevgrnd]                [double] volumetric soil water at saturation (porosity)
-h2osoi_vol[nlevgrnd]            [double] volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3]
 bsw[nlevgrnd]                   [double] Clapp and Hornberger "b
 smpso[numpft]                   [double] soil water potential at full stomatal opening (mm)                    
 smpsc[numpft]                   [double] soil water potential at full stomatal closure (mm)
 eff_porosity[nlevgrnd]          [double] effective soil porosity
+altmax_indx                     [int] index corresponding to maximum active layer depth from current year
+altmax_lastyear_indx            [int] index corresponding to maximum active layer depth from prior year
 
 OUTPUTS:
 rootr[nlevgrnd]                 [double] effective fraction of roots in each soil layer
-rresis[nlevgrnd]                [double] root soil water stress (resistance) by layer (0-1)
 btran                           [double] transpiration wetness factor (0 to 1) (integrated soil water stress)
 */
 void calc_root_moist_stress(
   const int& vtype,
   const double *h2osoi_liqvol,
   const double *rootfr,
-  const double *rootfr_unf,
   const double *t_soisno,
   const double& tc_stress,
   const double *sucsat,
   const double *watsat,
-  const double *h2osoi_vol,
   const double *bsw,
   const double *smpso,
   const double *smpsc,
   const double *eff_porosity,
+  const int& altmax_indx,
+  const int& altmax_lastyear_indx,
 
   double *rootr,
-  double *rresis,
   double& btran)
 {
-  double s_node, smp_node, smp_node_lf;
+  double s_node, smp_node;
   const double btran0 = 0.0;
+  double rootfr_unf[nlevgrnd] = {0.0}; // unfrozen root fraction
+  double rresis[nlevgrnd] = {0.0}; // root soil water stress (resistance) by layer (0-1)
+
+  normalize_unfrozen_rootfr(t_soisno, rootfr, altmax_indx, altmax_lastyear_indx, rootfr_unf);
 
   for (int i = 0; i < nlevgrnd; i++) {
     // Root resistance factors
@@ -185,7 +187,7 @@ void calc_root_moist_stress(
       s_node = std::max(h2osoi_liqvol[nlevsno + i] / eff_porosity[i], 0.01);
       soil_suction(sucsat[i], s_node, bsw[i], smp_node);
       smp_node = std::max(smpsc[vtype], smp_node);
-      rresis[i] = std::min( (eff_porosity[i]/watsat[i]) * (smp_node - smpsc[vtype]) / (smpso[vtype] - smpsc[vtype]), 1.0);
+      rresis[i] = std::min((eff_porosity[i]/watsat[i]) * (smp_node - smpsc[vtype]) / (smpso[vtype] - smpsc[vtype]), 1.0);
 
       if (!perchroot && !perchroot_alt) {
         rootr[i] = rootfr[i] * rresis[i];
@@ -194,9 +196,6 @@ void calc_root_moist_stress(
       }
 
       btran += std::max(rootr[i], 0.0);
-      s_node = h2osoi_vol[i] / watsat[i];
-      soil_suction(sucsat[i], s_node, bsw[i], smp_node_lf);
-      smp_node_lf = std::max(smpsc[vtype], smp_node_lf);
     }
   }
   // Normalize root resistances to get layer contribution to ET

src/clm_constants.h

@@ -71,7 +71,7 @@ static const double ELM_RWV     = ELM_RGAS/ELM_MWWV;          // Water vapor gas
 
 
 //from clm_varcon.F90
-static const double snw_rds_min = 54.526;               //minimum allowed snow effective radius (also "fresh snow" value) [microns]
+static const double snw_rds_min = 54.526;               // minimum allowed snow effective radius (also "fresh snow" value) [microns]
 static const double zlnd = 0.01;                        // Roughness length for soil [m]
 static const double zsno = 0.0024;                      // Roughness length for snow [m]
 static const double spval = 1.0e36;                     // special value for real data
@@ -83,10 +83,11 @@ static const double hvap   = 2.501e6;                   // Latent heat of evap f
 static const double hfus   = 3.337e5;                   // Latent heat of fusion for ice [J/kg]
 static const double hsub   = hvap + hfus;               // Latent heat of sublimation    [J/kg]
 static const double grav   = 9.80616;                   // gravity constant [m/s2]
-static const double roverg = ELM_RWV/grav*1000.0; // Rw/g constant = (8.3144/0.018)/(9.80616)*1000. mm/K
+static const double roverg = ELM_RWV/grav*1000.0;       // Rw/g constant = (8.3144/0.018)/(9.80616)*1000. mm/K
 static const double vkc    = 0.4;                       // von Karman constant [-]
 static const double cpair  = 1.00464e3;                 // specific heat of dry air   ~ J/kg/K
-
+static const double tlsai_crit = 2.0;                   // critical value of elai+esai for which aerodynamic parameters are maximum
+static const double sb = 5.67e-8;                       // Stefan-Boltzmann constant ~ W/m^2/K^4
 
 
 static const bool perchroot = false;     // calculate root function only in unfrozen soil, based on instantaneous temp

src/landtype.h

@@ -2,10 +2,10 @@
 #define LANDTYPE_H_
 
 struct LandType {
-  int ltype;
-  int ctype;
-  bool urbpoi;
-  bool lakpoi;
+  // land unit, urban unit, vegetation unit
+  int ltype, ctype, vtype;
+  // urban point, lake point
+  bool urbpoi, lakpoi;
 };
 
 #endif