Snow burial fix

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
 

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