diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index 01141fa7c3..6fdae37435 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -685,86 +685,107 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! Calculate biomass heat capacities ! if(use_biomass_heat_storage) then -bioms: do f = 1, fn - p = filterp(f) - ! fraction of stem receiving incoming radiation - frac_rad_abs_by_stem(p) = (esai(p))/(elai(p)+esai(p)) + fates_if: if(use_fates) then + + call clm_fates%BiomassHeatIndicesFromFates(filterp(1:fn), & ! in + k_internal, & ! in + canopystate_inst, & ! in + frac_rad_abs_by_stem(bounds%begp:bounds%endp), & ! in/out + sa_leaf(bounds%begp:bounds%endp), & ! in/out + sa_stem(bounds%begp:bounds%endp), & ! in/out + sa_internal(bounds%begp:bounds%endp), & ! in/out + cp_leaf(bounds%begp:bounds%endp), & ! in/out + cp_stem(bounds%begp:bounds%endp), & ! in/out + rstem(bounds%begp:bounds%endp)) ! in/out + + else + + bioms: do f = 1, fn + p = filterp(f) + + ! fraction of stem receiving incoming radiation + frac_rad_abs_by_stem(p) = (esai(p))/(elai(p)+esai(p)) - ! when elai = 0, do not multiply by k_vert (i.e. frac_rad_abs_by_stem = 1) - if(elai(p) > 0._r8) frac_rad_abs_by_stem(p) = k_vert * frac_rad_abs_by_stem(p) + ! when elai = 0, do not multiply by k_vert (i.e. frac_rad_abs_by_stem = 1) + if(elai(p) > 0._r8) frac_rad_abs_by_stem(p) = k_vert * frac_rad_abs_by_stem(p) + + ! if using Satellite Phenology mode, use values in parameter file + ! otherwise calculate dbh from stem biomass + + if(use_cn) then + if(stem_biomass(p) > 0._r8) then + dbh(p) = 2._r8 * sqrt(stem_biomass(p) * (1._r8 - fbw(patch%itype(p))) & + / ( shr_const_pi * htop(p) * k_cyl_vol & + * nstem(patch%itype(p)) * wood_density(patch%itype(p)))) + else + dbh(p) = 0._r8 + endif - ! if using Satellite Phenology mode, use values in parameter file - ! otherwise calculate dbh from stem biomass - if(use_cn) then - if(stem_biomass(p) > 0._r8) then - dbh(p) = 2._r8 * sqrt(stem_biomass(p) * (1._r8 - fbw(patch%itype(p))) & - / ( shr_const_pi * htop(p) * k_cyl_vol & - * nstem(patch%itype(p)) * wood_density(patch%itype(p)))) else - dbh(p) = 0._r8 + dbh(p) = dbh_param(patch%itype(p)) endif - else - dbh(p) = dbh_param(patch%itype(p)) - endif - ! leaf and stem surface area - sa_leaf(p) = elai(p) - ! double in spirit of full surface area for sensible heat - sa_leaf(p) = 2.*sa_leaf(p) - - ! Surface area for stem - sa_stem(p) = nstem(patch%itype(p))*(htop(p)*shr_const_pi*dbh(p)) - ! adjust for departure of cylindrical stem model - sa_stem(p) = k_cyl_area * sa_stem(p) - - ! - ! only calculate separate leaf/stem heat capacity for trees - ! and shrubs if dbh is greater than some minimum value - ! (set surface area for stem, and fraction absorbed by stem to zero) - if(.not.(is_tree(patch%itype(p)) .or. is_shrub(patch%itype(p))) & - .or. dbh(p) < min_stem_diameter) then - frac_rad_abs_by_stem(p) = 0.0 - sa_stem(p) = 0.0 - endif + ! leaf and stem surface area + sa_leaf(p) = elai(p) + ! double in spirit of full surface area for sensible heat + sa_leaf(p) = 2.*sa_leaf(p) + + ! Surface area for stem + sa_stem(p) = nstem(patch%itype(p))*(htop(p)*shr_const_pi*dbh(p)) + ! adjust for departure of cylindrical stem model + sa_stem(p) = k_cyl_area * sa_stem(p) + + ! + ! only calculate separate leaf/stem heat capacity for trees + ! and shrubs if dbh is greater than some minimum value + ! (set surface area for stem, and fraction absorbed by stem to zero) + if(.not.(is_tree(patch%itype(p)) .or. is_shrub(patch%itype(p))) & + .or. dbh(p) < min_stem_diameter) then + frac_rad_abs_by_stem(p) = 0.0 + sa_stem(p) = 0.0 + endif - ! cross-sectional area of stems - carea_stem = shr_const_pi * (dbh(p)*0.5)**2 - - ! if using Satellite Phenology mode, calculate leaf and stem biomass - if(.not. use_cn) then - ! boreal needleleaf lma*c2b ~ 0.25 kg dry mass/m2(leaf) - leaf_biomass(p) = 0.25_r8 * max(0.01_r8, sa_leaf(p)) & - / (1.-fbw(patch%itype(p))) - stem_biomass(p) = carea_stem * htop(p) * k_cyl_vol & - * nstem(patch%itype(p)) * wood_density(patch%itype(p)) & - /(1.-fbw(patch%itype(p))) - endif + ! cross-sectional area of stems + carea_stem = shr_const_pi * (dbh(p)*0.5)**2 + + ! if using Satellite Phenology mode, calculate leaf and stem biomass + if( .not. use_cn ) then + ! boreal needleleaf lma*c2b ~ 0.25 kg dry mass/m2(leaf) + leaf_biomass(p) = 0.25_r8 * max(0.01_r8, sa_leaf(p)) & + / (1.-fbw(patch%itype(p))) + stem_biomass(p) = carea_stem * htop(p) * k_cyl_vol & + * nstem(patch%itype(p)) * wood_density(patch%itype(p)) & + /(1.-fbw(patch%itype(p))) + endif - ! internal longwave fluxes between leaf and stem - ! (use same area of interaction i.e. ignore leaf <-> leaf) - sa_internal(p) = min(sa_leaf(p),sa_stem(p)) - sa_internal(p) = k_internal * sa_internal(p) + ! calculate specify heat capacity of vegetation + ! as weighted averaged of dry biomass and water + ! lma_dry has units of kg dry mass/m2 here + ! (Appendix B of Bonan et al., GMD, 2018) + + cp_leaf(p) = leaf_biomass(p) * (c_dry_biomass*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water) - ! calculate specify heat capacity of vegetation - ! as weighted averaged of dry biomass and water - ! lma_dry has units of kg dry mass/m2 here - ! (Appendix B of Bonan et al., GMD, 2018) + ! cp-stem will have units J/k/ground_area + cp_stem(p) = stem_biomass(p) * (c_dry_biomass*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water) + ! adjust for departure from cylindrical stem model + cp_stem(p) = k_cyl_vol * cp_stem(p) - cp_leaf(p) = leaf_biomass(p) * (c_dry_biomass*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water) + ! resistance between internal stem temperature and canopy air + rstem(p) = rstem_per_dbh(patch%itype(p))*dbh(p) - ! cp-stem will have units J/k/ground_area - cp_stem(p) = stem_biomass(p) * (c_dry_biomass*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water) - ! adjust for departure from cylindrical stem model - cp_stem(p) = k_cyl_vol * cp_stem(p) - ! resistance between internal stem temperature and canopy air - rstem(p) = rstem_per_dbh(patch%itype(p))*dbh(p) + ! internal longwave fluxes between leaf and stem + ! (use same area of interaction i.e. ignore leaf <-> leaf) + sa_internal(p) = min(sa_leaf(p),sa_stem(p)) + sa_internal(p) = k_internal * sa_internal(p) - enddo bioms + + enddo bioms + end if fates_if else - ! Otherwise set biomass heat storage terms to zero - do f = 1, fn + ! Otherwise set biomass heat storage terms to zero + do f = 1, fn p = filterp(f) sa_leaf(p) = (elai(p)+esai(p)) frac_rad_abs_by_stem(p) = 0._r8 @@ -773,9 +794,9 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, cp_leaf(p) = 0._r8 cp_stem(p) = 0._r8 rstem(p) = 0._r8 - end do + end do end if - + ! calculate daylength control for Vcmax do f = 1, fn @@ -1390,7 +1411,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! does not account for changes in SH or internal LW, ! as that would change result for t_veg above if (use_biomass_heat_storage) then - if (stem_biomass(p) > 0._r8) then + if (cp_stem(p) > 0._r8) then dt_stem(p) = (frac_rad_abs_by_stem(p)*(sabv(p) + air(p) + bir(p)*ts_ini(p)**4 & + cir(p)*lw_grnd) - eflx_sh_stem(p) & + lw_leaf(p)- lw_stem(p))/(cp_stem(p)/dtime & diff --git a/src/utils/clmfates_interfaceMod.F90 b/src/utils/clmfates_interfaceMod.F90 index 580340de01..c5d1bdd207 100644 --- a/src/utils/clmfates_interfaceMod.F90 +++ b/src/utils/clmfates_interfaceMod.F90 @@ -187,6 +187,7 @@ module CLMFatesInterfaceMod procedure, public :: wrap_canopy_radiation procedure, public :: wrap_bgc_summary procedure, public :: TransferZ0mDisp + procedure, public :: BiomassHeatIndicesFromFates procedure, private :: init_history_io procedure, private :: wrap_update_hlmfates_dyn procedure, private :: init_soil_depths @@ -1722,6 +1723,98 @@ subroutine wrap_photosynthesis(this, nc, bounds, fn, filterp, & end subroutine wrap_photosynthesis + ! ====================================================================================== + + subroutine BiomassHeatIndicesFromFates(this, & + filterp, & ! in + k_internal, & ! in + canopystate_inst, & ! in + frac_rad_abs_by_stem, & ! in/out + sa_leaf, & ! in/out + sa_stem, & ! in/out + sa_internal, & ! in/out + cp_leaf, & ! in/out + cp_stem, & ! in/out + rstem ) ! in/out + + + use clm_varcon, only : c_dry_biomass ! specific heat of dry biomass [J/kg/K] + use clm_varcon, only : c_water ! specific heat of water [J/kg/K] + + class(hlm_fates_interface_type), intent(inout) :: this + integer, intent(in) :: filterp(:) ! patch/pft filter + real(r8), intent(in) :: k_internal ! self-absorbtion parameter of leaf/stem longwave + type(canopystate_type), intent(in) :: canopystate_inst + real(r8), intent(inout) :: frac_rad_abs_by_stem(:) ! fraction of incoming radiation absorbed by stems + real(r8), intent(inout) :: sa_leaf(:) ! surface area of leaves [m2/m2_ground] + real(r8), intent(inout) :: sa_stem(:) ! surface area of stems [m2/m2_ground] + real(r8), intent(inout) :: sa_internal(:) ! min(sa_stem,sa_leaf) + real(r8), intent(inout) :: cp_leaf(:) ! heat capacity of leaves + real(r8), intent(inout) :: cp_stem(:) ! heat capacity of stems + real(r8), intent(inout) :: rstem(:) + + integer :: fn ! number of patches in filter + integer :: p ! Patch/pft index (CLM global position, ie begp:endp) + integer :: icp ! filter loop counter + integer :: ifp ! fates patch index associated with p, on this site/column + integer :: c ! column index (CLM global position, ie begc:endc) + + associate(elai => canopystate_inst%elai_patch , & + esai => canopystate_inst%esai_patch) + + write(iulog,*) 'biomass heat storage is only partially implemented' + write(iulog,*) 'this routine BiomassHeatingIndicesFromFates() should ' + write(iulog,*) 'not be callable, and is only a placeholder' + write(iulog,*) 'at the moment' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + fn = size(filterp) + + do icp = 1,fn + + p = filterp(icp) + c = patch%column(p) + ifp = p-col%patchi(c) + + ! fraction of stem receiving incoming radiation + frac_rad_abs_by_stem(p) = (esai(p))/(elai(p)+esai(p)) + + ! leaf and stem surface area + ! double in spirit of full surface area for sensible heat + sa_leaf(p) = 2.*elai(p) + + ! Surface area for stem +!! sa_stem(p) = this%fates(nc)%bc_out(s)%sa_stem(ifp) + + ! calculate specify heat capacity of vegetation + ! as weighted averaged of dry biomass and water + ! lma_dry has units of kg dry mass/m2 here + ! (Appendix B of Bonan et al., GMD, 2018) + +!! cp_leaf(p) = this%fates(nc)%bc_out(s)%patch_leaf_biomass(ifp) * & +!! (c_dry_biomass*(1. - this%fates(nc)%bc_out(s)%patch_leaf_fwb(ifp)) + & +!! c_water*this%fates(nc)%bc_out(s)%patch_leaf_fwb(ifp)) + + ! cp-stem will have units J/k/ground_area +!! cp_stem(p) = this%fates(nc)%bc_out(s)%patch_stem_biomass(ifp) * & +!! (c_dry_biomass*(1. - this%fates(nc)%bc_out(s)%patch_stem_fwb(ifp)) + & +!! c_water*this%fates(nc)%bc_out(s)%patch_stem_fwb(ifp)) + + ! resistance between internal stem temperature and canopy air + ! we pass in an effective dbh from the patch. This is volume + ! weighted mean dbh of all cohorts in the patch +!! rstem(p) = rstem_per_dbh * this%fates(nc)%bc_out(s)%patch_eff_dbh(ifp) + + ! internal longwave fluxes between leaf and stem + ! (use same area of interaction i.e. ignore leaf <-> leaf) + sa_internal(p) = k_internal * min(sa_leaf(p),sa_stem(p)) + + end do + end associate + return + end subroutine BiomassHeatIndicesFromFates + + ! ====================================================================================== subroutine wrap_accumulatefluxes(this, nc, fn, filterp)