Merge branch 'SnowBurialFix'

Grass/crop snow burial fraction fix, and some other fixes

(1) Crop/grass snow burial changed based on height and lodging (Danica
    Lombardozzi) - resolves #516. Updates grass and crop snow
    burial ('fb' in SatellitePhenologyMod.F90) from being 0.2m to change
    with PFT height and accounts for 20% lodging (falling over).

(2) In BGC code, make elai and esai depend on frac_sno (Danica
    Lombardozzi, Keith Oleson & Bill Sacks) - resolves #1116

(3) Fix to quadratic solution error bug caused by negative shaded
    photosynthesis (Danica Lombardozzi & Keith Oleson) - resolves
    #756

(4) Point to cime branch tag that avoids building CDEPS with LILAC, to
    avoid dependence on unreleased ESMF code (Bill Sacks) - resolves
    #1203

(5) Some tweaks to documentation and how the documentation is built
    (Bill Sacks)

Resolves #516 (crop and grass snow burial)

Resolves #1116 (elai and esai for non-SP runs should use
frac_sno)

Resolves #756 (Quadratic solution error)

Resolves #1203 (Rework LILAC build to avoid dependence on
unreleased ESMF code)

Externals.cfg

@@ -37,7 +37,7 @@ required = True
 local_path = cime
 protocol = git
 repo_url = https://github.com/ESMCI/cime
-tag = cime5.8.32
+tag = branch_tags/cime5.8.32_a01
 externals = ../Externals_cime.cfg
 required = True
 

doc/ChangeLog

@@ -1,4 +1,151 @@
 ===============================================================
+Tag name:  ctsm5.1.dev014
+Originator(s): dll (Danica Lombardozzi) / oleson (Keith Oleson) / sacks (Bill Sacks)
+Date:  Sat Nov  7 12:19:14 MST 2020
+One-line Summary: Grass/crop snow burial fraction fix, and some other fixes
+
+Purpose of changes
+------------------
+
+(1) Crop/grass snow burial changed based on height and lodging (Danica
+    Lombardozzi) - resolves ESCOMP/CTSM#516. Updates grass and crop snow
+    burial ('fb' in SatellitePhenologyMod.F90) from being 0.2m to change
+    with PFT height and accounts for 20% lodging (falling over).
+
+(2) In BGC code, make elai and esai depend on frac_sno (Danica
+    Lombardozzi, Keith Oleson & Bill Sacks) - resolves ESCOMP/CTSM#1116
+
+(3) Fix to quadratic solution error bug caused by negative shaded
+    photosynthesis (Danica Lombardozzi & Keith Oleson) - resolves
+    ESCOMP/CTSM#756
+
+(4) Point to cime branch tag that avoids building CDEPS with LILAC, to
+    avoid dependence on unreleased ESMF code (Bill Sacks) - resolves
+    ESCOMP/CTSM#1203
+
+(5) Some tweaks to documentation and how the documentation is built
+    (Bill Sacks)
+
+Bugs fixed or introduced
+------------------------
+
+Issues fixed (include CTSM Issue #):
+- Resolves ESCOMP/CTSM#516 (crop and grass snow burial)
+- Resolves ESCOMP/CTSM#1116 (elai and esai for non-SP runs should use
+  frac_sno)
+- Resolves ESCOMP/CTSM#756 (Quadratic solution error)
+- Resolves ESCOMP/CTSM#1203 (Rework LILAC build to avoid dependence on
+  unreleased ESMF code)
+
+CIME Issues fixed (include issue #):
+- https://github.com/ESMCI/cime/issues/3769 (Avoid building CDEPS with
+  CTSM's LILAC)
+
+Significant changes to scientifically-supported configurations
+--------------------------------------------------------------
+
+Does this tag change answers significantly for any of the following physics configurations?
+(Details of any changes will be given in the "Answer changes" section below.)
+
+    [Put an [X] in the box for any configuration with significant answer changes.]
+
+[X] clm5_1
+
+[X] clm5_0
+
+[X] ctsm5_0-nwp
+
+[X] clm4_5
+
+Notes of particular relevance for users
+---------------------------------------
+
+Caveats for users (e.g., need to interpolate initial conditions): none
+
+Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): none
+
+Changes made to namelist defaults (e.g., changed parameter values): none
+
+Changes to the datasets (e.g., parameter, surface or initial files): none
+
+Substantial timing or memory changes: none
+
+Notes of particular relevance for developers: (including Code reviews and testing)
+---------------------------------------------
+NOTE: Be sure to review the steps in README.CHECKLIST.master_tags as well as the coding style in the Developers Guide
+
+Caveats for developers (e.g., code that is duplicated that requires double maintenance): none
+
+Changes to tests or testing: none
+
+CTSM testing:
+
+ [PASS means all tests PASS and OK means tests PASS other than expected fails.]
+
+  build-namelist tests:
+
+    cheyenne - not run
+
+  tools-tests (test/tools):
+
+    cheyenne - not run
+
+  PTCLM testing (tools/shared/PTCLM/test):
+
+    cheyenne - not run
+
+  python testing (see instructions in python/README.md; document testing done):
+
+    (any machine) - not run
+
+  regular tests (aux_clm):
+
+    cheyenne ---- ok
+    izumi ------- ok
+
+    ok: tests pass, baselines fail as expected
+
+If the tag used for baseline comparisons was NOT the previous tag, note that here:
+
+
+Answer changes
+--------------
+
+Changes answers relative to baseline: YES
+
+  Summarize any changes to answers, i.e.,
+    - what code configurations: all (or at least most)
+    - what platforms/compilers: all
+    - nature of change (roundoff; larger than roundoff/same climate; new climate): 
+      new climate, at least to some extent
+
+      Changes arise from (1) & (2) (large changes) and from (3) (small
+      changes). The change in snow burial parameterizations results in
+      changes to albedo and water & energy fluxes.
+
+   If bitwise differences were observed, how did you show they were no worse
+   than roundoff? N/A
+
+   If this tag changes climate describe the run(s) done to evaluate the new
+   climate (put details of the simulations in the experiment database)
+       - casename: N/A
+
+   URL for LMWG diagnostics output used to validate new climate:
+       For analysis of diffs, see
+       https://github.com/danicalombardozzi/ctsm_py/blob/7543b0f3b413bae9974c11b467fdbc0413c3b7fa/notebooks/SnowBurial.ipynb
+
+Detailed list of changes
+------------------------
+
+List any externals directories updated (cime, rtm, mosart, cism, fates, etc.):
+- cime: cime5.8.32 -> branch_tags/cime5.8.32_a01
+- doc-builder: v1.0.2 -> v1.0.4
+
+Pull Requests that document the changes (include PR ids):
+https://github.com/ESCOMP/CTSM/pull/1112
+
+===============================================================
+===============================================================
 Tag name:  ctsm5.1.dev013
 Originator(s): jedwards (Jim Edwards) / sacks (Bill Sacks)
 Date:  Wed Nov  4 14:24:55 MST 2020

doc/ChangeSum

@@ -1,5 +1,6 @@
 Tag                   Who      Date  Summary
 ============================================================================================================================
+  ctsm5.1.dev014 dll/oles 11/07/2020 Grass/crop snow burial fraction fix, and some other fixes
   ctsm5.1.dev013 jedwards 11/04/2020 Allow pnetcdf for vector history files
   ctsm5.1.dev012 oleson/s 11/03/2020 Allow nlevgrnd < nlevurb
   ctsm5.1.dev011    sacks 11/02/2020 Change CISM2%NOEVOLVE compsets to use SGLC; documentation updates

src/biogeochem/CNVegStructUpdateMod.F90

@@ -74,8 +74,10 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, &
     real(r8) :: tsai_min   ! PATCH derived minimum tsai
     real(r8) :: tsai_alpha ! monthly decay rate of tsai
     real(r8) :: dt         ! radiation time step (sec)
+    real(r8) :: frac_sno_adjusted ! frac_sno adjusted per frac_sno_threshold
 
     real(r8), parameter :: dtsmonth = 2592000._r8 ! number of seconds in a 30 day month (60x60x24x30)
+    real(r8), parameter :: frac_sno_threshold = 0.999_r8  ! frac_sno values greater than this are treated as 1
     !-----------------------------------------------------------------------
     ! tsai formula from Zeng et. al. 2002, Journal of Climate, p1835
     !
@@ -106,7 +108,8 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, &
          nind               =>  dgvs_inst%nind_patch                    , & ! Input:  [real(r8) (:) ] number of individuals (#/m**2)                    
          fpcgrid            =>  dgvs_inst%fpcgrid_patch                 , & ! Input:  [real(r8) (:) ] fractional area of patch (pft area/nat veg area)    
 
-         snow_depth         =>  waterdiagnosticbulk_inst%snow_depth_col          , & ! Input:  [real(r8) (:) ] snow height (m)                                   
+         frac_sno           =>  waterdiagnosticbulk_inst%frac_sno_col   , & ! Input:  [real(r8) (:) ] fraction of ground covered by snow (0 to 1)
+         snow_depth         =>  waterdiagnosticbulk_inst%snow_depth_col , & ! Input:  [real(r8) (:) ] snow height (m)                                   
 
          forc_hgt_u_patch   =>  frictionvel_inst%forc_hgt_u_patch       , & ! Input:  [real(r8) (:) ] observational height of wind at patch-level [m]     
 
@@ -282,18 +285,29 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, &
          end if
 
          ! adjust lai and sai for burying by snow. 
-         ! snow burial fraction for short vegetation (e.g. grasses) as in
-         ! Wang and Zeng, 2007.
+         ! snow burial fraction for short vegetation (e.g. grasses, crops) changes with vegetation height 
+         ! accounts for a 20% bending factor, as used in Lombardozzi et al. (2018) GRL 45(18), 9889-9897
+
+         ! NOTE: The following snow burial code is duplicated in SatellitePhenologyMod.
+         ! Changes in one place should be accompanied by similar changes in the other.
+
          if (ivt(p) > noveg .and. ivt(p) <= nbrdlf_dcd_brl_shrub ) then
             ol = min( max(snow_depth(c)-hbot(p), 0._r8), htop(p)-hbot(p))
             fb = 1._r8 - ol / max(1.e-06_r8, htop(p)-hbot(p))
          else
-            fb = 1._r8 - max(min(snow_depth(c),0.2_r8),0._r8)/0.2_r8   ! 0.2m is assumed
+            fb = 1._r8 - (max(min(snow_depth(c),max(0.05,htop(p)*0.8_r8)),0._r8)/(max(0.05,htop(p)*0.8_r8)))
             !depth of snow required for complete burial of grasses
          endif
 
-         elai(p) = max(tlai(p)*fb, 0.0_r8)
-         esai(p) = max(tsai(p)*fb, 0.0_r8)
+         if (frac_sno(c) <= frac_sno_threshold) then
+            frac_sno_adjusted = frac_sno(c)
+         else
+            ! avoid tiny but non-zero elai and esai that can cause radiation and/or photosynthesis code to blow up
+            frac_sno_adjusted = 1._r8
+         end if
+
+         elai(p) = max(tlai(p)*(1.0_r8 - frac_sno_adjusted) + tlai(p)*fb*frac_sno_adjusted, 0.0_r8)
+         esai(p) = max(tsai(p)*(1.0_r8 - frac_sno_adjusted) + tsai(p)*fb*frac_sno_adjusted, 0.0_r8)
 
          ! Fraction of vegetation free of snow
          if ((elai(p) + esai(p)) > 0._r8) then

src/biogeochem/SatellitePhenologyMod.F90

@@ -403,19 +403,20 @@ subroutine SatellitePhenology(bounds, num_nolakep, filter_nolakep, &
          ! are less than 0.05, set equal to zero to prevent numerical
          ! problems associated with very small lai and sai.
 
-         ! snow burial fraction for short vegetation (e.g. grasses) as in
-         ! Wang and Zeng, 2007. 
+         ! snow burial fraction for short vegetation (e.g. grasses, crops) changes with vegetation height 
+         ! accounts for a 20% bending factor, as used in Lombardozzi et al. (2018) GRL 45(18), 9889-9897
+
+         ! NOTE: The following snow burial code is duplicated in CNVegStructUpdateMod.
+         ! Changes in one place should be accompanied by similar changes in the other.
 
          if (patch%itype(p) > noveg .and. patch%itype(p) <= nbrdlf_dcd_brl_shrub ) then
             ol = min( max(snow_depth(c)-hbot(p), 0._r8), htop(p)-hbot(p))
             fb = 1._r8 - ol / max(1.e-06_r8, htop(p)-hbot(p))
          else
-            fb = 1._r8 - max(min(snow_depth(c),0.2_r8),0._r8)/0.2_r8   ! 0.2m is assumed
-            !depth of snow required for complete burial of grasses
+            fb = 1._r8 - (max(min(snow_depth(c),max(0.05,htop(p)*0.8_r8)),0._r8)/(max(0.05,htop(p)*0.8_r8)))
          endif
 
          ! area weight by snow covered fraction
-
          elai(p) = max(tlai(p)*(1.0_r8 - frac_sno(c)) + tlai(p)*fb*frac_sno(c), 0.0_r8)
          esai(p) = max(tsai(p)*(1.0_r8 - frac_sno(c)) + tsai(p)*fb*frac_sno(c), 0.0_r8)
          if (elai(p) < 0.05_r8) elai(p) = 0._r8

src/biogeophys/PhotosynthesisMod.F90

@@ -4138,51 +4138,61 @@ subroutine ci_func_PHS(x,cisun, cisha, fvalsun, fvalsha, p, iv, c, bsun, bsha, b
     ! With an <= 0, then gs_mol = bbb
 
     ! Sunlit
-    cs_sun = cair - 1.4_r8/gb_mol * an_sun(p,iv) * forc_pbot(c)
-    cs_sun = max(cs_sun,10.e-06_r8)
+    if (an_sun(p,iv) >= 0._r8) then
+       cs_sun = cair - 1.4_r8/gb_mol * an_sun(p,iv) * forc_pbot(c)
+       cs_sun = max(cs_sun,10.e-06_r8)
+    end if
 
     if ( stomatalcond_mtd == stomatalcond_mtd_medlyn2011 )then
-       term = 1.6_r8 * an_sun(p,iv) / (cs_sun / forc_pbot(c) * 1.e06_r8)
-       aquad = 1.0_r8
-       bquad = -(2.0 * (medlynintercept(patch%itype(p))*1.e-06_r8 + term) + (medlynslope(patch%itype(p)) * term)**2 / &
+       if (an_sun(p,iv) >= 0._r8) then
+          term = 1.6_r8 * an_sun(p,iv) / (cs_sun / forc_pbot(c) * 1.e06_r8)
+          aquad = 1.0_r8
+          bquad = -(2.0 * (medlynintercept(patch%itype(p))*1.e-06_r8 + term) + (medlynslope(patch%itype(p)) * term)**2 / &
                (gb_mol*1.e-06_r8 * rh_can))
-       cquad = medlynintercept(patch%itype(p))*medlynintercept(patch%itype(p))*1.e-12_r8 + &
+          cquad = medlynintercept(patch%itype(p))*medlynintercept(patch%itype(p))*1.e-12_r8 + &
                (2.0*medlynintercept(patch%itype(p))*1.e-06_r8 + term * &
                (1.0 - medlynslope(patch%itype(p))* medlynslope(patch%itype(p)) / rh_can)) * term
 
-       call quadratic (aquad, bquad, cquad, r1, r2)
-       gs_mol_sun = max(r1,r2) * 1.e06_r8
-
+          call quadratic (aquad, bquad, cquad, r1, r2)
+          gs_mol_sun = max(r1,r2) * 1.e06_r8
+       end if
+
        ! Shaded
-       cs_sha = cair - 1.4_r8/gb_mol * an_sha(p,iv) * forc_pbot(c)
-       cs_sha = max(cs_sha,10.e-06_r8)
-
-       term = 1.6_r8 * an_sha(p,iv) / (cs_sha / forc_pbot(c) * 1.e06_r8)
-       aquad = 1.0_r8
-       bquad = -(2.0 * (medlynintercept(patch%itype(p))*1.e-06_r8 + term) + (medlynslope(patch%itype(p)) * term)**2 / &
+       if (an_sha(p,iv) >= 0._r8) then
+          cs_sha = cair - 1.4_r8/gb_mol * an_sha(p,iv) * forc_pbot(c)
+          cs_sha = max(cs_sha,10.e-06_r8)
+
+          term = 1.6_r8 * an_sha(p,iv) / (cs_sha / forc_pbot(c) * 1.e06_r8)
+          aquad = 1.0_r8
+          bquad = -(2.0 * (medlynintercept(patch%itype(p))*1.e-06_r8 + term) + (medlynslope(patch%itype(p)) * term)**2 / &
                (gb_mol*1.e-06_r8 * rh_can))
-       cquad = medlynintercept(patch%itype(p))*medlynintercept(patch%itype(p))*1.e-12_r8 + &
+          cquad = medlynintercept(patch%itype(p))*medlynintercept(patch%itype(p))*1.e-12_r8 + &
                (2.0*medlynintercept(patch%itype(p))*1.e-06_r8 + term * (1.0 - medlynslope(patch%itype(p))* &
                medlynslope(patch%itype(p)) / rh_can)) * term
 
-       call quadratic (aquad, bquad, cquad, r1, r2)
-       gs_mol_sha = max(r1,r2)* 1.e06_r8
+          call quadratic (aquad, bquad, cquad, r1, r2)
+          gs_mol_sha = max(r1,r2)* 1.e06_r8
+       end if
     else if ( stomatalcond_mtd == stomatalcond_mtd_bb1987 )then
-       aquad = cs_sun
-       bquad = cs_sun*(gb_mol - max(bsun*bbb(p),1._r8)) - mbb(p)*an_sun(p,iv)*forc_pbot(c)
-       cquad = -gb_mol*(cs_sun*max(bsun*bbb(p),1._r8) + mbb(p)*an_sun(p,iv)*forc_pbot(c)*rh_can)
-       call quadratic (aquad, bquad, cquad, r1, r2)
-       gs_mol_sun = max(r1,r2)
-
+       if (an_sun(p,iv) >= 0._r8) then
+          aquad = cs_sun
+          bquad = cs_sun*(gb_mol - max(bsun*bbb(p),1._r8)) - mbb(p)*an_sun(p,iv)*forc_pbot(c)
+          cquad = -gb_mol*(cs_sun*max(bsun*bbb(p),1._r8) + mbb(p)*an_sun(p,iv)*forc_pbot(c)*rh_can)
+          call quadratic (aquad, bquad, cquad, r1, r2)
+          gs_mol_sun = max(r1,r2)
+       end if
+
        ! Shaded
-       cs_sha = cair - 1.4_r8/gb_mol * an_sha(p,iv) * forc_pbot(c)
-       cs_sha = max(cs_sha,10.e-06_r8)
-
-       aquad = cs_sha
-       bquad = cs_sha*(gb_mol - max(bsha*bbb(p),1._r8)) - mbb(p)*an_sha(p,iv)*forc_pbot(c)
-       cquad = -gb_mol*(cs_sha*max(bsha*bbb(p),1._r8) + mbb(p)*an_sha(p,iv)*forc_pbot(c)*rh_can)
-       call quadratic (aquad, bquad, cquad, r1, r2)
-       gs_mol_sha = max(r1,r2)
+       if (an_sha(p,iv) >= 0._r8) then
+          cs_sha = cair - 1.4_r8/gb_mol * an_sha(p,iv) * forc_pbot(c)
+          cs_sha = max(cs_sha,10.e-06_r8)
+
+          aquad = cs_sha
+          bquad = cs_sha*(gb_mol - max(bsha*bbb(p),1._r8)) - mbb(p)*an_sha(p,iv)*forc_pbot(c)
+          cquad = -gb_mol*(cs_sha*max(bsha*bbb(p),1._r8) + mbb(p)*an_sha(p,iv)*forc_pbot(c)*rh_can)
+          call quadratic (aquad, bquad, cquad, r1, r2)
+          gs_mol_sha = max(r1,r2)
+       end if
     end if
 
     ! Derive new estimate for cisun,cisha