diff --git a/Externals.cfg b/Externals.cfg index 7fbf9a86e3..3abc4b91fc 100644 --- a/Externals.cfg +++ b/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 diff --git a/src/biogeochem/CNVegStructUpdateMod.F90 b/src/biogeochem/CNVegStructUpdateMod.F90 index 6f5b19476e..22bbea8614 100644 --- a/src/biogeochem/CNVegStructUpdateMod.F90 +++ b/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 diff --git a/src/biogeochem/SatellitePhenologyMod.F90 b/src/biogeochem/SatellitePhenologyMod.F90 index cd03a742a8..469c7ba37a 100644 --- a/src/biogeochem/SatellitePhenologyMod.F90 +++ b/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 diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index 946a8b25cf..31289b0ca3 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/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