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.

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

Modified crown allometric spread logic #279

Merged
merged 11 commits into from
Oct 6, 2017
46 changes: 19 additions & 27 deletions biogeochem/EDCanopyStructureMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -802,59 +802,51 @@ subroutine canopy_spread( currentSite )
! Calculates the spatial spread of tree canopies based on canopy closure.
!
! !USES:
use EDParamsMod , only : ED_val_maxspread, ED_val_minspread
use EDTypesMod , only : AREA
use EDParamsMod, only : ED_val_canopy_closure_thresh
!
! !ARGUMENTS
type (ed_site_type), intent(inout), target :: currentSite
!
! !LOCAL VARIABLES:
type (ed_cohort_type), pointer :: currentCohort
type (ed_patch_type) , pointer :: currentPatch
real(r8) :: arealayer(nlevleaf) ! Amount of canopy in each layer.
real(r8) :: arealayer ! Amount of canopy in each layer.
real(r8) :: inc ! Arbitrary daily incremental change in canopy area
integer :: z
!----------------------------------------------------------------------

inc = 0.005_r8
inc = 0.05_r8

currentPatch => currentSite%oldest_patch

arealayer = 0.0_r8
do while (associated(currentPatch))

!calculate canopy area in each canopy storey...
arealayer = 0.0_r8
!calculate canopy area in each patch...
currentCohort => currentPatch%tallest
do while (associated(currentCohort))
currentCohort%c_area = c_area(currentCohort)
if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then
arealayer(currentCohort%canopy_layer) = arealayer(currentCohort%canopy_layer) + currentCohort%c_area
if(EDPftvarcon_inst%woody(currentCohort%pft) .eq. 1 .and. currentCohort%canopy_layer .eq. 1 ) then
arealayer = arealayer + currentCohort%c_area
endif
currentCohort => currentCohort%shorter
enddo

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So this is not really arealayer, but area of all upper canopy layers in the site.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

right, i guess i ought to rename it

!If the canopy area is approaching closure, squash the tree canopies and make them taller and thinner
do z = 1,nclmax

if(arealayer(z)/currentPatch%area > 0.9_r8)then
currentPatch%spread(z) = currentPatch%spread(z) - inc
else
currentPatch%spread(z) = currentPatch%spread(z) + inc
endif
if(currentPatch%spread(z) >= ED_val_maxspread)then
currentPatch%spread(z) = ED_val_maxspread
endif
if(currentPatch%spread(z) <= ED_val_minspread)then
currentPatch%spread(z) = ED_val_minspread
endif
enddo !z
!write(fates_log(),*) 'spread',currentPatch%spread(1:2)
!currentPatch%spread(:) = ED_val_maxspread
!FIX(RF,033114) spread is off
!write(fates_log(),*) 'canopy_spread',currentPatch%area,currentPatch%spread(1:2)
currentPatch => currentPatch%younger
currentPatch => currentPatch%younger

enddo !currentPatch

!If the canopy area is approaching closure, squash the tree canopies and make them taller and thinner
if( arealayer/AREA .gt. ED_val_canopy_closure_thresh ) then
currentSite%spread = currentSite%spread - inc
else
currentSite%spread = currentSite%spread + inc
endif

! put within bounds to make sure it stays between 0 and 1
currentSite%spread = max(min(currentSite%spread, 1._r8), 0._r8)

end subroutine canopy_spread


Expand Down
21 changes: 8 additions & 13 deletions biogeochem/EDGrowthFunctionsMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -225,14 +225,13 @@ real(r8) function c_area( cohort_in )
! Function of DBH (cm) canopy spread (m/cm) and number of individuals.
! ============================================================================

use EDParamsMod , only : ED_val_grass_spread
use EDTypesMod , only : nclmax

type(ed_cohort_type), intent(in) :: cohort_in

real(r8) :: dbh ! Tree diameter at breat height. cm.
real(r8) :: crown_area_to_dbh_exponent
integer :: can_layer_index
real(r8) :: spreadterm

! default is to use the same exponent as the dbh to bleaf exponent so that per-plant canopy depth remains invariant during growth,
! but allowed to vary via the allom_blca_expnt_diff term (which has default value of zero)
Expand All @@ -244,9 +243,8 @@ real(r8) function c_area( cohort_in )
write(fates_log(),*) 'z_area 2',EDPftvarcon_inst%max_dbh
write(fates_log(),*) 'z_area 3',EDPftvarcon_inst%woody
write(fates_log(),*) 'z_area 4',cohort_in%n
write(fates_log(),*) 'z_area 5',cohort_in%patchptr%spread
write(fates_log(),*) 'z_area 5',cohort_in%siteptr%spread
write(fates_log(),*) 'z_area 6',cohort_in%canopy_layer
write(fates_log(),*) 'z_area 7',ED_val_grass_spread
end if

dbh = min(cohort_in%dbh,EDPftvarcon_inst%max_dbh(cohort_in%pft))
Expand All @@ -260,15 +258,12 @@ real(r8) function c_area( cohort_in )
! So we allow layer index exceedence here and force it down to max.
! (rgk/cdk 05/2017)
! ----------------------------------------------------------------------------------

can_layer_index = min(cohort_in%canopy_layer,nclmax)

if(EDPftvarcon_inst%woody(cohort_in%pft) == 1)then
c_area = 3.142_r8 * cohort_in%n * &
(cohort_in%patchptr%spread(can_layer_index)*dbh)**crown_area_to_dbh_exponent
else
c_area = 3.142_r8 * cohort_in%n * (ED_val_grass_spread*dbh)**crown_area_to_dbh_exponent
end if

! apply site-level spread elasticity to the cohort crown allometry term
spreadterm = cohort_in%siteptr%spread * EDPftvarcon_inst%allom_d2ca_coefficient_max(cohort_in%pft) + &
(1._r8 - cohort_in%siteptr%spread) * EDPftvarcon_inst%allom_d2ca_coefficient_min(cohort_in%pft)
!
c_area = 3.142_r8 * cohort_in%n * (spreadterm * dbh)**crown_area_to_dbh_exponent
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what is this constant?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

its kinda like pi but only has 4 significant figures. happy to switch to actual pi now that i've confirmed that it's bit for bit as it currently stands.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

interesting, that's the circle thing

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

which is itself slightly weird. I feel like that made more sense physically when the crown_area_to_dbh_exponent was 2 (as in the original PPA papers) but at this point its sort of an artifact. that term could easily be pulled into the allom_d2ca_coefficient but for simplicity I left as-is since doing so would change the dependency between allom_d2ca_coefficient and crown_area_to_dbh_exponent.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok, so then would people prefer I subsume it into the default value of the allom_d2ca_coefficient given whatever default crown_area_to_dbh_exponent we are currently using? or leave as-is and make it a named pi?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

well, honestly my comment started off as a joke, but now I see where you are going. There is no math reflecting a circular form anymore, so the pi is somewhat misleading. imho, the precision on that pi specifically isn't important and everyone knows what 3.14 is. If you want to remove it and subsume in the parameter, up to you, I would see it as a bonus prize.

Copy link
Contributor Author

@ckoven ckoven Oct 4, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

actually, while doing this, i realize that its weird that this particular power law is written as (ax)^b whereas all the others are just ax^b. does anyone mind if I pull the a out of the parentheses to conform to what the other pieces of code are doing (and adjust the value accordingly)?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

that comment got mangled by formatting, look at edited version

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks. that makes it easier to understand the crown plasticity as just the ratio of the max:min terms rather than that ratio raised to some power.


end function c_area

Expand Down
13 changes: 3 additions & 10 deletions biogeochem/EDPatchDynamicsMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ subroutine spawn_patches( currentSite, bc_in)
!
! !USES:

use EDParamsMod , only : ED_val_maxspread, ED_val_understorey_death
use EDParamsMod , only : ED_val_understorey_death
use EDCohortDynamicsMod , only : zero_cohort, copy_cohort, terminate_cohorts

!
Expand All @@ -300,7 +300,6 @@ subroutine spawn_patches( currentSite, bc_in)
real(r8) :: leaf_litter_local(maxpft) ! initial value of leaf litter. KgC/m2
real(r8) :: cwd_ag_local(ncwd) ! initial value of above ground coarse woody debris. KgC/m2
real(r8) :: cwd_bg_local(ncwd) ! initial value of below ground coarse woody debris. KgC/m2
real(r8) :: spread_local(nclmax) ! initial value of canopy spread parameter.no units
!---------------------------------------------------------------------

storesmallcohort => null() ! storage of the smallest cohort for insertion routine
Expand Down Expand Up @@ -328,12 +327,11 @@ subroutine spawn_patches( currentSite, bc_in)
cwd_bg_local = 0.0_r8
leaf_litter_local = 0.0_r8
root_litter_local = 0.0_r8
spread_local(1:nclmax) = ED_val_maxspread
age = 0.0_r8

allocate(new_patch)
call create_patch(currentSite, new_patch, age, site_areadis, &
spread_local, cwd_ag_local, cwd_bg_local, leaf_litter_local, &
cwd_ag_local, cwd_bg_local, leaf_litter_local, &
root_litter_local)

new_patch%tallest => null()
Expand Down Expand Up @@ -748,8 +746,6 @@ subroutine average_patch_properties( currentPatch, newPatch, patch_site_areadis

enddo

newPatch%spread = newPatch%spread + currentPatch%spread * patch_site_areadis/newPatch%area

end subroutine average_patch_properties

! ============================================================================
Expand Down Expand Up @@ -1079,7 +1075,7 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat
end subroutine mortality_litter_fluxes

! ============================================================================
subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_local,cwd_bg_local, &
subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_local, &
leaf_litter_local,root_litter_local)
!
! !DESCRIPTION:
Expand All @@ -1096,7 +1092,6 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_
real(r8), intent(in) :: cwd_bg_local(:) ! initial value of below ground coarse woody debris. KgC/m2
real(r8), intent(in) :: root_litter_local(:)! initial value of root litter. KgC/m2
real(r8), intent(in) :: leaf_litter_local(:)! initial value of leaf litter. KgC/m2
real(r8), intent(in) :: spread_local(:) ! initial value of canopy spread parameter.no units
!
! !LOCAL VARIABLES:
!---------------------------------------------------------------------
Expand Down Expand Up @@ -1126,7 +1121,6 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_
new_patch%age = age
new_patch%age_class = 1
new_patch%area = areap
new_patch%spread = spread_local
new_patch%cwd_ag = cwd_ag_local
new_patch%cwd_bg = cwd_bg_local
new_patch%leaf_litter = leaf_litter_local
Expand Down Expand Up @@ -1221,7 +1215,6 @@ subroutine zero_patch(cp_p)
currentPatch%nrad(:,:) = 999 ! number of exposed leaf layers for each canopy layer and pft
currentPatch%ncan(:,:) = 999 ! number of total leaf layers for each canopy layer and pft
currentPatch%lai = nan ! leaf area index of patch
currentPatch%spread(:) = nan ! dynamic ratio of dbh to canopy area.
currentPatch%pft_agb_profile(:,:) = nan

! DISTURBANCE
Expand Down
10 changes: 5 additions & 5 deletions main/EDInitMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,9 @@ subroutine zero_site( site_in )
! Resources management (logging/harvesting, etc)
site_in%resources_management%trunk_product_site = 0.0_r8

! canopy spread
site_in%spread = 0._r8

end subroutine zero_site

! ============================================================================
Expand Down Expand Up @@ -187,7 +190,7 @@ subroutine set_site_properties( nsites, sites)
sites(s)%frac_burnt = 0.0_r8
sites(s)%old_stock = 0.0_r8


sites(s)%spread = 1.0_r8
end do

return
Expand All @@ -204,7 +207,6 @@ subroutine init_patches( nsites, sites, bc_in)
!


use EDParamsMod , only : ED_val_maxspread
use FatesPlantHydraulicsMod, only : updateSizeDepRhizHydProps
use FatesInventoryInitMod, only : initialize_sites_by_inventory

Expand All @@ -218,7 +220,6 @@ subroutine init_patches( nsites, sites, bc_in)
integer :: s
real(r8) :: cwd_ag_local(ncwd)
real(r8) :: cwd_bg_local(ncwd)
real(r8) :: spread_local(nclmax)
real(r8) :: leaf_litter_local(maxpft)
real(r8) :: root_litter_local(maxpft)
real(r8) :: age !notional age of this patch
Expand All @@ -237,7 +238,6 @@ subroutine init_patches( nsites, sites, bc_in)
cwd_bg_local(:) = 0.0_r8 !ED_val_init_litter
leaf_litter_local(:) = 0.0_r8
root_litter_local(:) = 0.0_r8
spread_local(:) = ED_val_maxspread
age = 0.0_r8
! ---------------------------------------------------------------------------------------------

Expand Down Expand Up @@ -276,7 +276,7 @@ subroutine init_patches( nsites, sites, bc_in)

! make new patch...
call create_patch(sites(s), newp, age, AREA, &
spread_local, cwd_ag_local, cwd_bg_local, leaf_litter_local, &
cwd_ag_local, cwd_bg_local, leaf_litter_local, &
root_litter_local)

call init_cohorts(newp, bc_in(s))
Expand Down
44 changes: 12 additions & 32 deletions main/EDParamsMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -23,11 +23,8 @@ module EDParamsMod
real(r8),protected :: ED_size_diagnostic_scale ! Flag to switch between a linear and exponential
! scale on the plant size axis in diagnostics (NOT USED YET)
real(r8),protected :: fates_mortality_disturbance_fraction ! the fraction of canopy mortality that results in disturbance
real(r8),protected :: ED_val_grass_spread
real(r8),protected :: ED_val_comp_excln
real(r8),protected :: ED_val_stress_mort
real(r8),protected :: ED_val_maxspread ! maximum ratio of dbh to canopy area (cm/m2)
real(r8),protected :: ED_val_minspread ! minimum ratio of dbh to canopy area (cm/m2)
real(r8),protected :: ED_val_init_litter
real(r8),protected :: ED_val_nignitions
real(r8),protected :: ED_val_understorey_death
Expand All @@ -47,14 +44,12 @@ module EDParamsMod
real(r8),protected :: ED_val_phen_coldtemp
real(r8),protected :: ED_val_cohort_fusion_tol
real(r8),protected :: ED_val_patch_fusion_tol
real(r8),protected :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry

character(len=param_string_length),parameter :: ED_name_size_diagnostic_scale = "fates_size_diagnostic_scale"
character(len=param_string_length),parameter :: ED_name_mort_disturb_frac = "fates_mort_disturb_frac"
character(len=param_string_length),parameter :: ED_name_grass_spread = "fates_grass_spread"
character(len=param_string_length),parameter :: ED_name_comp_excln = "fates_comp_excln"
character(len=param_string_length),parameter :: ED_name_stress_mort = "fates_stress_mort"
character(len=param_string_length),parameter :: ED_name_maxspread = "fates_maxspread"
character(len=param_string_length),parameter :: ED_name_minspread = "fates_minspread"
character(len=param_string_length),parameter :: ED_name_init_litter = "fates_init_litter"
character(len=param_string_length),parameter :: ED_name_nignitions = "fates_nignitions"
character(len=param_string_length),parameter :: ED_name_understorey_death = "fates_understorey_death"
Expand All @@ -73,7 +68,8 @@ module EDParamsMod
character(len=param_string_length),parameter :: ED_name_phen_ncolddayslim= "fates_phen_ncolddayslim"
character(len=param_string_length),parameter :: ED_name_phen_coldtemp= "fates_phen_coldtemp"
character(len=param_string_length),parameter :: ED_name_cohort_fusion_tol= "fates_cohort_fusion_tol"
character(len=param_string_length),parameter :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol"
character(len=param_string_length),parameter :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol"
character(len=param_string_length),parameter :: ED_name_canopy_closure_thresh= "fates_canopy_closure_thresh"

! Hydraulics Control Parameters (ONLY RELEVANT WHEN USE_FATES_HYDR = TRUE)
! ----------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -121,11 +117,8 @@ subroutine FatesParamsInit()

ED_size_diagnostic_scale = nan
fates_mortality_disturbance_fraction = nan
ED_val_grass_spread = nan
ED_val_comp_excln = nan
ED_val_stress_mort = nan
ED_val_maxspread = nan
ED_val_minspread = nan
ED_val_init_litter = nan
ED_val_nignitions = nan
ED_val_understorey_death = nan
Expand All @@ -145,6 +138,7 @@ subroutine FatesParamsInit()
ED_val_phen_coldtemp = nan
ED_val_cohort_fusion_tol = nan
ED_val_patch_fusion_tol = nan
ED_val_canopy_closure_thresh = nan

hydr_psi0 = nan
hydr_psicap = nan
Expand Down Expand Up @@ -179,21 +173,12 @@ subroutine FatesRegisterParams(fates_params)
call fates_params%RegisterParameter(name=ED_name_mort_disturb_frac, dimension_shape=dimension_shape_1d, &
dimension_names=dim_names)

call fates_params%RegisterParameter(name=ED_name_grass_spread, dimension_shape=dimension_shape_1d, &
dimension_names=dim_names)

call fates_params%RegisterParameter(name=ED_name_comp_excln, dimension_shape=dimension_shape_1d, &
dimension_names=dim_names)

call fates_params%RegisterParameter(name=ED_name_stress_mort, dimension_shape=dimension_shape_1d, &
dimension_names=dim_names)

call fates_params%RegisterParameter(name=ED_name_maxspread, dimension_shape=dimension_shape_1d, &
dimension_names=dim_names)

call fates_params%RegisterParameter(name=ED_name_minspread, dimension_shape=dimension_shape_1d, &
dimension_names=dim_names)

call fates_params%RegisterParameter(name=ED_name_init_litter, dimension_shape=dimension_shape_1d, &
dimension_names=dim_names)

Expand Down Expand Up @@ -251,6 +236,9 @@ subroutine FatesRegisterParams(fates_params)
call fates_params%RegisterParameter(name=ED_name_patch_fusion_tol, dimension_shape=dimension_shape_1d, &
dimension_names=dim_names)

call fates_params%RegisterParameter(name=ED_name_canopy_closure_thresh, dimension_shape=dimension_shape_1d, &
dimension_names=dim_names)

call fates_params%RegisterParameter(name=hydr_name_psi0, dimension_shape=dimension_shape_1d, &
dimension_names=dim_names)

Expand Down Expand Up @@ -291,21 +279,12 @@ subroutine FatesReceiveParams(fates_params)
call fates_params%RetreiveParameter(name=ED_name_mort_disturb_frac, &
data=fates_mortality_disturbance_fraction)

call fates_params%RetreiveParameter(name=ED_name_grass_spread, &
data=ED_val_grass_spread)

call fates_params%RetreiveParameter(name=ED_name_comp_excln, &
data=ED_val_comp_excln)

call fates_params%RetreiveParameter(name=ED_name_stress_mort, &
data=ED_val_stress_mort)

call fates_params%RetreiveParameter(name=ED_name_maxspread, &
data=ED_val_maxspread)

call fates_params%RetreiveParameter(name=ED_name_minspread, &
data=ED_val_minspread)

call fates_params%RetreiveParameter(name=ED_name_init_litter, &
data=ED_val_init_litter)

Expand Down Expand Up @@ -363,6 +342,9 @@ subroutine FatesReceiveParams(fates_params)
call fates_params%RetreiveParameter(name=ED_name_patch_fusion_tol, &
data=ED_val_patch_fusion_tol)

call fates_params%RetreiveParameter(name=ED_name_canopy_closure_thresh, &
data=ED_val_canopy_closure_thresh)

call fates_params%RetreiveParameter(name=hydr_name_psi0, &
data=hydr_psi0)

Expand Down Expand Up @@ -400,11 +382,8 @@ subroutine FatesReportParams(is_master)
write(fates_log(),*) '----------- FATES Scalar Parameters -----------------'
write(fates_log(),fmt0) 'ED_size_diagnostic_scale = ',ED_size_diagnostic_scale
write(fates_log(),fmt0) 'fates_mortality_disturbance_fraction = ',fates_mortality_disturbance_fraction
write(fates_log(),fmt0) 'ED_val_grass_spread = ',ED_val_grass_spread
write(fates_log(),fmt0) 'ED_val_comp_excln = ', ED_val_comp_excln
write(fates_log(),fmt0) 'ED_val_comp_excln = ',ED_val_comp_excln
write(fates_log(),fmt0) 'ED_val_stress_mort = ',ED_val_stress_mort
write(fates_log(),fmt0) 'ED_val_maxspread = ',ED_val_maxspread
write(fates_log(),fmt0) 'ED_val_minspread = ',ED_val_minspread
write(fates_log(),fmt0) 'ED_val_init_litter = ',ED_val_init_litter
write(fates_log(),fmt0) 'ED_val_nignitions = ',ED_val_nignitions
write(fates_log(),fmt0) 'ED_val_understorey_death = ',ED_val_understorey_death
Expand All @@ -424,6 +403,7 @@ subroutine FatesReportParams(is_master)
write(fates_log(),fmt0) 'ED_val_phen_coldtemp = ',ED_val_phen_coldtemp
write(fates_log(),fmt0) 'ED_val_cohort_fusion_tol = ',ED_val_cohort_fusion_tol
write(fates_log(),fmt0) 'ED_val_patch_fusion_tol = ',ED_val_patch_fusion_tol
write(fates_log(),fmt0) 'ED_val_canopy_closure_thresh = ',ED_val_canopy_closure_thresh
write(fates_log(),fmt0) 'hydr_psi0 = ',hydr_psi0
write(fates_log(),fmt0) 'hydr_psicap = ',hydr_psicap
write(fates_log(),fmt0) 'logging_dbhmin = ',logging_dbhmin
Expand Down
Loading