Skip to content

Commit

Permalink
Merge pull request #9 from rgknox/mpaiao-pr-allom-modes-cstarve-merge
Browse files Browse the repository at this point in the history 
merge in cstarvation params and main, update xml patch file
  • Loading branch information
@mpaiao
mpaiao authored Feb 5, 2024
2 parents d29e182 + f450fb7 commit 92d0dff
Show file tree
Hide file tree
Showing 54 changed files with 13,004 additions and 3,706 deletions.
303 changes: 175 additions & 128 deletions biogeochem/EDCanopyStructureMod.F90
View file

Large diffs are not rendered by default.

36 changes: 27 additions & 9 deletions biogeochem/EDCohortDynamicsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,10 @@ Module EDCohortDynamicsMod
use PRTAllometricCNPMod, only : acnp_bc_out_id_pefflux, acnp_bc_out_id_limiter
use PRTAllometricCNPMod, only : acnp_bc_in_id_cdamage
use DamageMainMod, only : undamaged_class
use FatesConstantsMod, only : n_term_mort_types
use FatesConstantsMod, only : i_term_mort_type_cstarv
use FatesConstantsMod, only : i_term_mort_type_canlev
use FatesConstantsMod, only : i_term_mort_type_numdens

use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=)
use shr_log_mod, only : errMsg => shr_log_errMsg
Expand Down Expand Up @@ -382,12 +386,14 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
integer :: terminate ! do we terminate (itrue) or not (ifalse)
integer :: istat ! return status code
character(len=255) :: smsg
integer :: termination_type
!----------------------------------------------------------------------

currentCohort => currentPatch%shortest
do while (associated(currentCohort))

terminate = ifalse
termination_type = 0
tallerCohort => currentCohort%taller

leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element)
Expand All @@ -400,6 +406,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
! Check if number density is so low is breaks math (level 1)
if (currentcohort%n < min_n_safemath .and. level == 1) then
terminate = itrue
termination_type = i_term_mort_type_numdens
if ( debug ) then
write(fates_log(),*) 'terminating cohorts 0',currentCohort%n/currentPatch%area,currentCohort%dbh,currentCohort%pft,call_index
endif
Expand All @@ -413,6 +420,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
currentCohort%n <= min_nppatch .or. &
(currentCohort%dbh < 0.00001_r8 .and. store_c < 0._r8) ) then
terminate = itrue
termination_type = i_term_mort_type_numdens
if ( debug ) then
write(fates_log(),*) 'terminating cohorts 1',currentCohort%n/currentPatch%area,currentCohort%dbh,currentCohort%pft,call_index
endif
Expand All @@ -421,6 +429,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
! Outside the maximum canopy layer
if (currentCohort%canopy_layer > nclmax ) then
terminate = itrue
termination_type = i_term_mort_type_canlev
if ( debug ) then
write(fates_log(),*) 'terminating cohorts 2', currentCohort%canopy_layer,currentCohort%pft,call_index
endif
Expand All @@ -430,6 +439,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
if ( ( sapw_c+leaf_c+fnrt_c ) < 1e-10_r8 .or. &
store_c < 1e-10_r8) then
terminate = itrue
termination_type = i_term_mort_type_cstarv
if ( debug ) then
write(fates_log(),*) 'terminating cohorts 3', &
sapw_c,leaf_c,fnrt_c,store_c,currentCohort%pft,call_index
Expand All @@ -439,6 +449,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
! Total cohort biomass is negative
if ( ( struct_c+sapw_c+leaf_c+fnrt_c+store_c ) < 0._r8) then
terminate = itrue
termination_type = i_term_mort_type_cstarv
if ( debug ) then
write(fates_log(),*) 'terminating cohorts 4', &
struct_c,sapw_c,leaf_c,fnrt_c,store_c,currentCohort%pft,call_index
Expand All @@ -448,7 +459,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
endif ! if (.not.currentCohort%isnew .and. level == 2) then

if (terminate == itrue) then
call terminate_cohort(currentSite, currentPatch, currentCohort, bc_in)
call terminate_cohort(currentSite, currentPatch, currentCohort, bc_in, termination_type)
deallocate(currentCohort, stat=istat, errmsg=smsg)
if (istat/=0) then
write(fates_log(),*) 'dealloc001: fail on terminate_cohorts:deallocate(currentCohort):'//trim(smsg)
Expand All @@ -461,7 +472,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
end subroutine terminate_cohorts

!-------------------------------------------------------------------------------------!
subroutine terminate_cohort(currentSite, currentPatch, currentCohort, bc_in)
subroutine terminate_cohort(currentSite, currentPatch, currentCohort, bc_in, termination_type)
!
! !DESCRIPTION:
! Terminates an individual cohort and updates the site-level
Expand All @@ -474,6 +485,7 @@ subroutine terminate_cohort(currentSite, currentPatch, currentCohort, bc_in)
type (fates_patch_type) , intent(inout), target :: currentPatch
type (fates_cohort_type), intent(inout), target :: currentCohort
type(bc_in_type), intent(in) :: bc_in
integer, intent(in) :: termination_type

! !LOCAL VARIABLES:
type (fates_cohort_type) , pointer :: shorterCohort
Expand All @@ -491,6 +503,12 @@ subroutine terminate_cohort(currentSite, currentPatch, currentCohort, bc_in)

!----------------------------------------------------------------------

! check termination_type; it should not be 0
if (termination_type == 0) then
write(fates_log(),*) 'termination_type=0'
call endrun(msg=errMsg(sourcefile, __LINE__))
endif

leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element)
store_c = currentCohort%prt%GetState(store_organ, carbon12_element)
sapw_c = currentCohort%prt%GetState(sapw_organ, carbon12_element)
Expand All @@ -506,16 +524,16 @@ subroutine terminate_cohort(currentSite, currentPatch, currentCohort, bc_in)

! Update the site-level carbon flux and individuals count for the appropriate canopy layer
if(levcan==ican_upper) then
currentSite%term_nindivs_canopy(currentCohort%size_class,currentCohort%pft) = &
currentSite%term_nindivs_canopy(currentCohort%size_class,currentCohort%pft) + currentCohort%n
currentSite%term_nindivs_canopy(termination_type,currentCohort%size_class,currentCohort%pft) = &
currentSite%term_nindivs_canopy(termination_type,currentCohort%size_class,currentCohort%pft) + currentCohort%n

currentSite%term_carbonflux_canopy(currentCohort%pft) = currentSite%term_carbonflux_canopy(currentCohort%pft) + &
currentSite%term_carbonflux_canopy(termination_type,currentCohort%pft) = currentSite%term_carbonflux_canopy(termination_type,currentCohort%pft) + &
currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c)
else
currentSite%term_nindivs_ustory(currentCohort%size_class,currentCohort%pft) = &
currentSite%term_nindivs_ustory(currentCohort%size_class,currentCohort%pft) + currentCohort%n
currentSite%term_nindivs_ustory(termination_type,currentCohort%size_class,currentCohort%pft) = &
currentSite%term_nindivs_ustory(termination_type,currentCohort%size_class,currentCohort%pft) + currentCohort%n

currentSite%term_carbonflux_ustory(currentCohort%pft) = currentSite%term_carbonflux_ustory(currentCohort%pft) + &
currentSite%term_carbonflux_ustory(termination_type,currentCohort%pft) = currentSite%term_carbonflux_ustory(termination_type,currentCohort%pft) + &
currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c)
end if

Expand Down Expand Up @@ -553,7 +571,7 @@ subroutine terminate_cohort(currentSite, currentPatch, currentCohort, bc_in)

call currentCohort%FreeMemory()

end subroutine terminate_cohort
end subroutine terminate_cohort

! =====================================================================================

Expand Down
44 changes: 22 additions & 22 deletions biogeochem/EDLoggingMortalityMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ module EDLoggingMortalityMod
use PRTGenericMod , only : sapw_organ, struct_organ, leaf_organ
use PRTGenericMod , only : fnrt_organ, store_organ, repro_organ
use FatesAllometryMod , only : set_root_fraction
use FatesConstantsMod , only : primaryforest, secondaryforest, secondary_age_threshold
use FatesConstantsMod , only : primaryland, secondaryland, secondary_age_threshold
use FatesConstantsMod , only : fates_tiny
use FatesConstantsMod , only : months_per_year, days_per_sec, years_per_day, g_per_kg
use FatesConstantsMod , only : hlm_harvest_area_fraction
Expand Down Expand Up @@ -199,7 +199,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, &
lmort_collateral,lmort_infra, l_degrad, &
hlm_harvest_rates, hlm_harvest_catnames, &
hlm_harvest_units, &
patch_anthro_disturbance_label, secondary_age, &
patch_land_use_label, secondary_age, &
frac_site_primary, harvestable_forest_c, &
harvest_tag)

Expand All @@ -210,7 +210,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, &
real(r8), intent(in) :: hlm_harvest_rates(:) ! annual harvest rate per hlm category
character(len=64), intent(in) :: hlm_harvest_catnames(:) ! names of hlm harvest categories
integer, intent(in) :: hlm_harvest_units ! unit type of hlm harvest rates: [area vs. mass]
integer, intent(in) :: patch_anthro_disturbance_label ! patch level anthro_disturbance_label
integer, intent(in) :: patch_land_use_label ! patch level land_use_label
real(r8), intent(in) :: secondary_age ! patch level age_since_anthro_disturbance
real(r8), intent(in) :: harvestable_forest_c(:) ! total harvestable forest carbon
! of all hlm harvest categories
Expand Down Expand Up @@ -265,7 +265,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, &
! HARVEST_SH3 = harvest from secondary non-forest (assume this is young for biomass)

! Get the area-based harvest rates based on info passed to FATES from the boundary condition
call get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_catnames, &
call get_harvest_rate_area (patch_land_use_label, hlm_harvest_catnames, &
hlm_harvest_rates, frac_site_primary, secondary_age, harvest_rate)

! For area-based harvest, harvest_tag shall always be 2 (not applicable).
Expand All @@ -280,7 +280,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, &
! 2=use carbon from hlm
! shall call another subroutine, which transfers biomass/carbon into fraction

call get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_catnames, &
call get_harvest_rate_carbon (patch_land_use_label, hlm_harvest_catnames, &
hlm_harvest_rates, secondary_age, harvestable_forest_c, &
harvest_rate, harvest_tag, cur_harvest_tag)

Expand Down Expand Up @@ -348,7 +348,7 @@ end subroutine LoggingMortality_frac

! ============================================================================

subroutine get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_catnames, hlm_harvest_rates, &
subroutine get_harvest_rate_area (patch_land_use_label, hlm_harvest_catnames, hlm_harvest_rates, &
frac_site_primary, secondary_age, harvest_rate)


Expand All @@ -361,7 +361,7 @@ subroutine get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_ca
! Arguments
real(r8), intent(in) :: hlm_harvest_rates(:) ! annual harvest rate per hlm category
character(len=64), intent(in) :: hlm_harvest_catnames(:) ! names of hlm harvest categories
integer, intent(in) :: patch_anthro_disturbance_label ! patch level anthro_disturbance_label
integer, intent(in) :: patch_land_use_label ! patch level land_use_label
real(r8), intent(in) :: secondary_age ! patch level age_since_anthro_disturbance
real(r8), intent(in) :: frac_site_primary
real(r8), intent(out) :: harvest_rate
Expand All @@ -374,17 +374,17 @@ subroutine get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_ca
! We do account forest only since non-forest harvest has geographical mismatch to LUH2 dataset
harvest_rate = 0._r8
do h_index = 1,hlm_num_lu_harvest_cats
if (patch_anthro_disturbance_label .eq. primaryforest) then
if (patch_land_use_label .eq. primaryland) then
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. &
hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") then
harvest_rate = harvest_rate + hlm_harvest_rates(h_index)
endif
else if (patch_anthro_disturbance_label .eq. secondaryforest .and. &
else if (patch_land_use_label .eq. secondaryland .and. &
secondary_age >= secondary_age_threshold) then
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then
harvest_rate = harvest_rate + hlm_harvest_rates(h_index)
endif
else if (patch_anthro_disturbance_label .eq. secondaryforest .and. &
else if (patch_land_use_label .eq. secondaryland .and. &
secondary_age < secondary_age_threshold) then
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. &
hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then
Expand All @@ -396,7 +396,7 @@ subroutine get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_ca
! Normalize by site-level primary or secondary forest fraction
! since harvest_rate is specified as a fraction of the gridcell
! also need to put a cap so as not to harvest more primary or secondary area than there is in a gridcell
if (patch_anthro_disturbance_label .eq. primaryforest) then
if (patch_land_use_label .eq. primaryland) then
if (frac_site_primary .gt. fates_tiny) then
harvest_rate = min((harvest_rate / frac_site_primary),frac_site_primary)
else
Expand Down Expand Up @@ -511,18 +511,18 @@ subroutine get_harvestable_carbon (csite, site_area, hlm_harvest_catnames, harve
! since we have not separated forest vs. non-forest
! all carbon belongs to the forest categories
do h_index = 1,hlm_num_lu_harvest_cats
if (currentPatch%anthro_disturbance_label .eq. primaryforest) then
if (currentPatch%land_use_label .eq. primaryland) then
! Primary
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1") then
harvestable_forest_c(h_index) = harvestable_forest_c(h_index) + harvestable_patch_c
end if
else if (currentPatch%anthro_disturbance_label .eq. secondaryforest .and. &
else if (currentPatch%land_use_label .eq. secondaryland .and. &
currentPatch%age_since_anthro_disturbance >= secondary_age_threshold) then
! Secondary mature
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then
harvestable_forest_c(h_index) = harvestable_forest_c(h_index) + harvestable_patch_c
end if
else if (currentPatch%anthro_disturbance_label .eq. secondaryforest .and. &
else if (currentPatch%land_use_label .eq. secondaryland .and. &
currentPatch%age_since_anthro_disturbance < secondary_age_threshold) then
! Secondary young
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2") then
Expand All @@ -537,7 +537,7 @@ end subroutine get_harvestable_carbon

! ============================================================================

subroutine get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_catnames, &
subroutine get_harvest_rate_carbon (patch_land_use_label, hlm_harvest_catnames, &
hlm_harvest_rates, secondary_age, harvestable_forest_c, &
harvest_rate, harvest_tag, cur_harvest_tag)

Expand All @@ -550,7 +550,7 @@ subroutine get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_
! Arguments
real(r8), intent(in) :: hlm_harvest_rates(:) ! annual harvest rate per hlm category
character(len=64), intent(in) :: hlm_harvest_catnames(:) ! names of hlm harvest categories
integer, intent(in) :: patch_anthro_disturbance_label ! patch level anthro_disturbance_label
integer, intent(in) :: patch_land_use_label ! patch level land_use_label
real(r8), intent(in) :: secondary_age ! patch level age_since_anthro_disturbance
real(r8), intent(in) :: harvestable_forest_c(:) ! site level forest c matching criteria available for harvest, kgC site-1
real(r8), intent(out) :: harvest_rate ! area fraction
Expand Down Expand Up @@ -584,17 +584,17 @@ subroutine get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_
! mature and secondary young).
! Get the harvest rate from HLM
do h_index = 1,hlm_num_lu_harvest_cats
if (patch_anthro_disturbance_label .eq. primaryforest) then
if (patch_land_use_label .eq. primaryland) then
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. &
hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") then
harvest_rate_c = harvest_rate_c + hlm_harvest_rates(h_index)
endif
else if (patch_anthro_disturbance_label .eq. secondaryforest .and. &
else if (patch_land_use_label .eq. secondaryland .and. &
secondary_age >= secondary_age_threshold) then
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then
harvest_rate_c = harvest_rate_c + hlm_harvest_rates(h_index)
endif
else if (patch_anthro_disturbance_label .eq. secondaryforest .and. &
else if (patch_land_use_label .eq. secondaryland .and. &
secondary_age < secondary_age_threshold) then
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. &
hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then
Expand All @@ -606,7 +606,7 @@ subroutine get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_
! Determine harvest status (succesful or not)
! Here only three categories are used
do h_index = 1,hlm_num_lu_harvest_cats
if (patch_anthro_disturbance_label .eq. primaryforest) then
if (patch_land_use_label .eq. primaryland) then
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" ) then
if(harvestable_forest_c(h_index) >= harvest_rate_c) then
harvest_rate_supply = harvest_rate_supply + harvestable_forest_c(h_index)
Expand All @@ -615,7 +615,7 @@ subroutine get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_
harvest_tag(h_index) = 1
end if
end if
else if (patch_anthro_disturbance_label .eq. secondaryforest .and. &
else if (patch_land_use_label .eq. secondaryland .and. &
secondary_age >= secondary_age_threshold) then
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1" ) then
if(harvestable_forest_c(h_index) >= harvest_rate_c) then
Expand All @@ -625,7 +625,7 @@ subroutine get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_
harvest_tag(h_index) = 1
end if
end if
else if (patch_anthro_disturbance_label .eq. secondaryforest .and. &
else if (patch_land_use_label .eq. secondaryland .and. &
secondary_age < secondary_age_threshold) then
if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" ) then
if(harvestable_forest_c(h_index) >= harvest_rate_c) then
Expand Down
Loading

0 comments on commit 92d0dff

Please sign in to comment.