Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

merge fates hooks for biomass heat storage #3

Open
wants to merge 1 commit into
base: heat_storage_biomass
Choose a base branch
from
Open

merge fates hooks for biomass heat storage #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
159 changes: 90 additions & 69 deletions src/biogeophys/CanopyFluxesMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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 &
Expand Down
93 changes: 93 additions & 0 deletions src/utils/clmfates_interfaceMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down