diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 5c003600f4..d3ce44f499 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1995,7 +1995,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out) ! Pass FATES Harvested C to bc_out. call UpdateHarvestC(sites(s),bc_out(s)) - + end do ! This call to RecruitWaterStorage() makes an accounting of @@ -2009,7 +2009,6 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out) call RecruitWaterStorage(nsites,sites,bc_out) end if - end subroutine update_hlm_dynamics ! ===================================================================================== diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index ab0b77b737..73ea7231e7 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -14,6 +14,7 @@ module EDLoggingMortalityMod ! ==================================================================================== use FatesConstantsMod , only : r8 => fates_r8 + use FatesConstantsMod , only : rsnbl_math_prec use EDTypesMod , only : ed_cohort_type use EDTypesMod , only : ed_patch_type use EDTypesMod , only : site_massbal_type @@ -86,11 +87,6 @@ module EDLoggingMortalityMod real(r8), parameter :: harvest_litter_localization = 0.0_r8 - ! ! transfer factor from kg biomass (dry matter) to kg carbon - ! ! now we applied a simple fraction of 50% based on the IPCC - ! ! guideline - ! real(r8), parameter :: carbon_per_kg_biomass = 0.5_r8 - character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -99,6 +95,9 @@ module EDLoggingMortalityMod public :: logging_time public :: IsItLoggingTime public :: get_harvest_rate_area + public :: get_harvestable_carbon + public :: get_harvest_rate_carbon + public :: get_harvest_debt public :: UpdateHarvestC contains @@ -199,7 +198,8 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & hlm_harvest_rates, hlm_harvest_catnames, & hlm_harvest_units, & patch_anthro_disturbance_label, secondary_age, & - frac_site_primary) + frac_site_primary, harvestable_forest_c, & + harvest_tag) ! Arguments integer, intent(in) :: pft_i ! pft index @@ -210,6 +210,9 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & 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 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 + real(r8), intent(in) :: frac_site_primary real(r8), intent(out) :: lmort_direct ! direct (harvestable) mortality fraction real(r8), intent(out) :: lmort_collateral ! collateral damage mortality fraction real(r8), intent(out) :: lmort_infra ! infrastructure mortality fraction @@ -217,9 +220,15 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! but suffer from forest degradation (i.e. they ! are moved to newly-anthro-disturbed secondary ! forest patch) - real(r8), intent(in) :: frac_site_primary + integer, intent(out) :: harvest_tag(:) ! tag to record the harvest status + ! for the calculation of harvest debt in C-based + ! harvest mode + ! 0 - successful; + ! 1 - unsuccessful since not enough carbon + ! 2 - not applicable ! Local variables + integer :: cur_harvest_tag ! the harvest tag of the cohort today real(r8) :: harvest_rate ! the final harvest rate to apply to this cohort today ! todo: probably lower the dbhmin default value to 30 cm @@ -257,42 +266,44 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & call get_harvest_rate_area (patch_anthro_disturbance_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). + harvest_tag = 2 + cur_harvest_tag = 2 + if (fates_global_verbose()) then - write(fates_log(), *) 'Successfully Read Harvest Rate from HLM.' + write(fates_log(), *) 'Successfully Read Harvest Rate from HLM.', hlm_harvest_rates(:), harvest_rate end if else if (hlm_use_lu_harvest == itrue .and. hlm_harvest_units == hlm_harvest_carbon) then ! 2=use carbon from hlm - ! Shijie: Shall call another function, which transfer biomass/carbon into fraction? - ! Is it the correct place to call the function? - ! Inputs: patch_area, patch_biomass, what else? + ! shall call another subroutine, which transfers biomass/carbon into fraction - ! call get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_catnames, & - ! hlm_harvest_rates, frac_site_primary, secondary_age, harvest_rate) + call get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_catnames, & + hlm_harvest_rates, secondary_age, harvestable_forest_c, & + harvest_rate, harvest_tag, cur_harvest_tag) - ! if (fates_global_verbose()) then - ! write(fates_log(), *) 'Successfully Read Harvest Rate from HLM.', hlm_harvest_rates(:), harvest_rate - ! end if - - write(fates_log(),*) 'HLM harvest carbon data not implemented yet. Exiting.' - call endrun(msg=errMsg(sourcefile, __LINE__)) + if (fates_global_verbose()) then + write(fates_log(), *) 'Successfully Read Harvest Rate from HLM.', hlm_harvest_rates(:), harvest_rate, harvestable_forest_c + end if endif ! transfer of area to secondary land is based on overall area affected, not just logged crown area ! l_degrad accounts for the affected area between logged crowns if(prt_params%woody(pft_i) == itrue)then ! only set logging rates for trees - - ! direct logging rates, based on dbh min and max criteria - if (dbh >= logging_dbhmin .and. .not. & - ((logging_dbhmax < fates_check_param_set) .and. (dbh >= logging_dbhmax )) ) then - ! the logic of the above line is a bit unintuitive but allows turning off the dbhmax comparison entirely. - ! since there is an .and. .not. after the first conditional, the dbh:dbhmax comparison needs to be - ! the opposite of what would otherwise be expected... - lmort_direct = harvest_rate * logging_direct_frac - + if (cur_harvest_tag == 0) then + ! direct logging rates, based on dbh min and max criteria + if (dbh >= logging_dbhmin .and. .not. & + ((logging_dbhmax < fates_check_param_set) .and. (dbh >= logging_dbhmax )) ) then + ! the logic of the above line is a bit unintuitive but allows turning off the dbhmax comparison entirely. + ! since there is an .and. .not. after the first conditional, the dbh:dbhmax comparison needs to be + ! the opposite of what would otherwise be expected... + lmort_direct = harvest_rate * logging_direct_frac + else + lmort_direct = 0.0_r8 + end if else - lmort_direct = 0.0_r8 + lmort_direct = 0.0_r8 end if ! infrastructure (roads, skid trails, etc) mortality rates @@ -356,11 +367,12 @@ subroutine get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_ca integer :: h_index ! for looping over harvest categories integer :: icode ! Integer equivalent of the event code (parameter file only allows reals) - ! Loop around harvest categories to determine the annual hlm harvest rate for the current cohort based on patch history info + ! Loop around harvest categories to determine the annual hlm harvest rate for the current cohort based on patch history info + ! 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(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. & + 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 @@ -418,6 +430,253 @@ end subroutine get_harvest_rate_area ! ============================================================================ + subroutine get_harvestable_carbon (csite, site_area, hlm_harvest_catnames, harvestable_forest_c ) + + !USES: + use SFParamsMod, only : SF_val_cwd_frac + use EDTypesMod, only : AREA_INV + + + ! ------------------------------------------------------------------------------------------- + ! + ! DESCRIPTION: + ! get the total carbon availale for harvest for three different harvest categories: + ! primary forest, secondary mature forest and secondary young forest + ! under two different scenarios: + ! harvestable carbon: aggregate all cohorts matching the dbhmin harvest criteria + ! + ! this subroutine shall be called outside the patch loop + ! output will be used to estimate the area-based harvest rate (get_harvest_rate_carbon) + ! for each cohort. + + ! Arguments + type(ed_site_type), intent(in), target :: csite + real(r8), intent(in) :: site_area ! temporary variable + character(len=64), intent(in) :: hlm_harvest_catnames(:) ! names of hlm harvest categories + + real(r8), intent(out) :: harvestable_forest_c(hlm_num_lu_harvest_cats) + + ! Local Variables + type(ed_patch_type), pointer :: currentPatch + type(ed_cohort_type), pointer :: currentCohort + real(r8) :: harvestable_patch_c ! patch level total carbon available for harvest, kgC site-1 + real(r8) :: harvestable_cohort_c ! cohort level total carbon available for harvest, kgC site-1 + real(r8) :: sapw_m ! Biomass of sap wood + real(r8) :: struct_m ! Biomass of structural organs + integer :: pft ! Index of plant functional type + integer :: h_index ! for looping over harvest categories + + ! Initialization + harvestable_forest_c = 0._r8 + + ! loop over patches + currentPatch => csite%oldest_patch + do while (associated(currentPatch)) + harvestable_patch_c = 0._r8 + currentCohort => currentPatch%tallest + + do while (associated(currentCohort)) + pft = currentCohort%pft + + ! only account for cohorts matching the following conditions + if(int(prt_params%woody(pft)) == 1)then ! only set logging rates for trees + sapw_m = currentCohort%prt%GetState(sapw_organ, carbon12_element) + struct_m = currentCohort%prt%GetState(struct_organ, carbon12_element) + ! logging_direct_frac shall be 1 for LUH2 driven simulation and global simulation + ! in site level study logging_direct_frac shall be surveyed + ! unit: [kgC ] = [kgC/plant] * [plant/ha] * [ha/ 10k m2] * [ m2 area ] + harvestable_cohort_c = logging_direct_frac * ( sapw_m + struct_m ) * & + prt_params%allom_agb_frac(currentCohort%pft) * & + SF_val_CWD_frac(ncwd) * logging_export_frac * & + currentCohort%n * AREA_INV * site_area + + ! No harvest for trees without canopy + if (currentCohort%canopy_layer>=1) then + ! logging amount are based on dbh min and max criteria + if (currentCohort%dbh >= logging_dbhmin .and. .not. & + ((logging_dbhmax < fates_check_param_set) .and. (currentCohort%dbh >= logging_dbhmax )) ) then + ! Harvestable C: aggregate cohorts fit the criteria + harvestable_patch_c = harvestable_patch_c + harvestable_cohort_c + end if + end if + end if + currentCohort => currentCohort%shorter + end do + + ! judge which category the current patch belong to + ! 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 + ! 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. & + 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. & + currentPatch%age_since_anthro_disturbance < secondary_age_threshold) then + ! Secondary young + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2") then + harvestable_forest_c(h_index) = harvestable_forest_c(h_index) + harvestable_patch_c + end if + end if + end do + currentPatch => currentPatch%younger + end do + + end subroutine get_harvestable_carbon + + ! ============================================================================ + + subroutine get_harvest_rate_carbon (patch_anthro_disturbance_label, hlm_harvest_catnames, & + hlm_harvest_rates, secondary_age, harvestable_forest_c, & + harvest_rate, harvest_tag, cur_harvest_tag) + + ! ------------------------------------------------------------------------------------------- + ! + ! DESCRIPTION: + ! get the carbon-based harvest rates based on info passed to FATES from the boundary conditions in. + ! assumes logging_time == true + + ! 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 + 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 + integer, intent(inout) :: harvest_tag(:) ! 0. normal harvest; 1. current site does not have enough C but + ! can perform harvest by ignoring criteria; 2. current site does + ! not have enough carbon + ! This harvest tag shall be a patch level variable but since all + ! logging functions happen within cohort loop we can only put the + ! calculation here. Can think about optimizing the logging calculation + ! in the future. + integer, intent(out), optional :: cur_harvest_tag ! harvest tag of the current cohort + + ! Local Variables + integer :: h_index ! for looping over harvest categories + integer :: icode ! Integer equivalent of the event code (parameter file only allows reals) + real(r8) :: harvest_rate_c ! Temporary variable, kgC site-1 + real(r8) :: harvest_rate_supply ! Temporary variable, kgC site-1 + + ! This subroutine follows the same logic of get_harvest_rate_area + ! Loop over harvest categories to determine the hlm harvest rate demand and actual harvest rate for the + ! current cohort based on patch history info + + ! Initialize local variables + harvest_rate = 0._r8 + harvest_rate_c = 0._r8 + harvest_rate_supply = 0._r8 + harvest_tag(:) = 2 + + ! Since we have five harvest categories from forcing data but in FATES non-forest harvest + ! is merged with forest harvest, we only have three logging type in FATES (primary, secondary + ! 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(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. & + 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. & + 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 + harvest_rate_c = harvest_rate_c + hlm_harvest_rates(h_index) + endif + endif + end do + + ! 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(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) + harvest_tag(h_index) = 0 + else + harvest_tag(h_index) = 1 + end if + end if + else if (patch_anthro_disturbance_label .eq. secondaryforest .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 + harvest_rate_supply = harvest_rate_supply + harvestable_forest_c(h_index) + harvest_tag(h_index) = 0 + else + harvest_tag(h_index) = 1 + end if + end if + else if (patch_anthro_disturbance_label .eq. secondaryforest .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 + harvest_rate_supply = harvest_rate_supply + harvestable_forest_c(h_index) + harvest_tag(h_index) = 0 + else + harvest_tag(h_index) = 1 + end if + end if + end if + end do + + ! If any harvest category available, assign to cur_harvest_tag and trigger logging event + if(present(cur_harvest_tag))then + cur_harvest_tag = minval(harvest_tag) + end if + + ! Transfer carbon-based harvest rate to area-based harvest rate + if (harvest_rate_supply > rsnbl_math_prec .and. harvest_rate_supply > harvest_rate_c) then + harvest_rate = harvest_rate_c / harvest_rate_supply + else + ! If we force harvest rate to 1 when we don't have enough C, we will produce + ! primary patch with no area, which cannot be terminated under nocomp mode. + ! So we still keep the harvest rate to 0 for now. + harvest_rate = 0._r8 + end if + + ! Prevent the generation of tiny secondary patches + if(harvest_rate < 1e-8) harvest_rate = 0._r8 + + ! For carbon-based harvest rate, normalizing by site-level primary or secondary forest fraction + ! is not needed + + ! calculate today's harvest rate + ! whether to harvest today has already been determined by IsItLoggingTime + ! for icode == 2, icode < 0, and icode > 10000 apply the annual rate one time (no calc) + ! Bad logging event flag is caught in IsItLoggingTime, so don't check it here + icode = int(logging_event_code) + if(icode .eq. 1) then + ! Logging is turned off - not sure why we need another switch + harvest_rate = 0._r8 + else if(icode .eq. 3) then + ! Logging event every day - this may not work due to the mortality exclusivity + harvest_rate = harvest_rate / hlm_days_per_year + else if(icode .eq. 4) then + ! logging event once a month + if(hlm_current_day.eq.1 ) then + harvest_rate = harvest_rate / months_per_year + end if + end if + + end subroutine get_harvest_rate_carbon + + ! ============================================================================ + subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site_areadis, bc_in) ! ------------------------------------------------------------------------------------------- @@ -858,6 +1117,76 @@ subroutine UpdateHarvestC(currentSite,bc_out) AREA_INV * (1._r8 - pprodharv10_forest_mean) * unit_trans_factor return + end subroutine UpdateHarvestC + subroutine get_harvest_debt(site_in, bc_in, harvest_tag) + + ! + ! !DESCRIPTION: + ! + ! Calculate if we have harvest debt for primary and secondary land + ! Harvest debt is the accumulated total carbon + ! deficiency once the carbon amount available for harvest + ! is smaller than the harvest rate of forcing data. + ! Harvest debt is calculated on site level + ! TODO: we can define harvest debt as a fraction of the + ! harvest rate in the future + ! Note: Non-forest harvest is accounted for under forest + ! harvest, thus the harvest tag for non-forest is not applicable (= 2) + ! + ! !ARGUMENTS: + type(ed_site_type) , intent(inout), target :: site_in + type(bc_in_type), intent(in) :: bc_in + integer :: harvest_tag(hlm_num_lu_harvest_cats) + + ! !LOCAL VARIABLES: + integer :: h_index + real(r8) :: harvest_debt_pri + real(r8) :: harvest_debt_sec_mature + real(r8) :: harvest_debt_sec_young + + if(logging_time) then + + ! Initialize the local variables + harvest_debt_pri = 0._r8 + harvest_debt_sec_mature = 0._r8 + harvest_debt_sec_young = 0._r8 + + ! First we need to get harvest rate for all three categories + do h_index = 1, hlm_num_lu_harvest_cats + ! Primary forest harvest rate + if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. & + bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2" ) then + harvest_debt_pri = harvest_debt_pri + bc_in%hlm_harvest_rates(h_index) + else if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then + harvest_debt_sec_mature = harvest_debt_sec_mature + bc_in%hlm_harvest_rates(h_index) + else if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. & + bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then + harvest_debt_sec_young = harvest_debt_sec_young + bc_in%hlm_harvest_rates(h_index) + end if + end do + ! Next we get the harvest debt through the harvest tag + do h_index = 1, hlm_num_lu_harvest_cats + if (harvest_tag(h_index) .eq. 1) then + if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1") then + site_in%resources_management%harvest_debt = site_in%resources_management%harvest_debt + & + harvest_debt_pri + else if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then + site_in%resources_management%harvest_debt = site_in%resources_management%harvest_debt + & + harvest_debt_sec_mature + site_in%resources_management%harvest_debt_sec = site_in%resources_management%harvest_debt_sec + & + harvest_debt_sec_mature + else if(bc_in%hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2") then + site_in%resources_management%harvest_debt = site_in%resources_management%harvest_debt + & + harvest_debt_sec_young + site_in%resources_management%harvest_debt_sec = site_in%resources_management%harvest_debt_sec + & + harvest_debt_sec_young + end if + end if + end do + end if + + end subroutine get_harvest_debt + end module EDLoggingMortalityMod diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 388d021003..db5122a8d9 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -225,7 +225,8 @@ end subroutine mortality_rates ! ============================================================================ - subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary) + subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary, & + harvestable_forest_c, harvest_tag) ! ! !DESCRIPTION: @@ -234,7 +235,6 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr ! elsewhere). ! ! !USES: - use FatesInterfaceTypesMod, only : hlm_freq_day ! ! !ARGUMENTS @@ -242,6 +242,14 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr type(ed_cohort_type),intent(inout), target :: currentCohort type(bc_in_type), intent(in) :: bc_in real(r8), intent(in) :: frac_site_primary + + real(r8), intent(in) :: harvestable_forest_c(:) ! total carbon available for logging, kgC site-1 + integer, intent(out) :: harvest_tag(:) ! tag to record the harvest status + ! for the calculation of harvest debt in C-based + ! harvest mode + ! 0 - successful; + ! 1 - unsuccessful since not enough carbon + ! 2 - not applicable ! ! !LOCAL VARIABLES: real(r8) :: cmort ! starvation mortality rate (fraction per year) @@ -271,10 +279,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr bc_in%hlm_harvest_units, & currentCohort%patchptr%anthro_disturbance_label, & currentCohort%patchptr%age_since_anthro_disturbance, & - frac_site_primary) - - - + frac_site_primary, harvestable_forest_c, harvest_tag) if (currentCohort%canopy_layer > 1)then ! Include understory logging mortality rates not associated with disturbance diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index a5c2202430..a939ca6579 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -48,6 +48,7 @@ module EDPatchDynamicsMod use FatesInterfaceTypesMod , only : hlm_use_sp use FatesInterfaceTypesMod , only : hlm_use_nocomp use FatesInterfaceTypesMod , only : hlm_use_fixed_biogeog + use FatesInterfaceTypesMod , only : hlm_num_lu_harvest_cats use FatesGlobals , only : endrun => fates_endrun use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue, ifalse @@ -58,6 +59,9 @@ module EDPatchDynamicsMod use EDLoggingMortalityMod, only : logging_litter_fluxes use EDLoggingMortalityMod, only : logging_time use EDLoggingMortalityMod, only : get_harvest_rate_area + use EDLoggingMortalityMod, only : get_harvest_rate_carbon + use EDLoggingMortalityMod, only : get_harvestable_carbon + use EDLoggingMortalityMod, only : get_harvest_debt use EDParamsMod , only : fates_mortality_disturbance_fraction use FatesAllometryMod , only : carea_allom use FatesAllometryMod , only : set_root_fraction @@ -70,6 +74,7 @@ module EDPatchDynamicsMod use FatesConstantsMod , only : n_anthro_disturbance_categories use FatesConstantsMod , only : fates_unset_r8 use FatesConstantsMod , only : fates_unset_int + use FatesConstantsMod , only : hlm_harvest_carbon use EDCohortDynamicsMod , only : InitPRTObject use EDCohortDynamicsMod , only : InitPRTBoundaryConditions use ChecksBalancesMod, only : SiteMassStock @@ -187,9 +192,12 @@ subroutine disturbance_rates( site_in, bc_in) real(r8) :: dist_rate_ldist_notharvested integer :: threshold_sizeclass integer :: i_dist + integer :: h_index real(r8) :: frac_site_primary real(r8) :: harvest_rate real(r8) :: tempsum + real(r8) :: harvestable_forest_c(hlm_num_lu_harvest_cats) + integer :: harvest_tag(hlm_num_lu_harvest_cats) !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) @@ -198,6 +206,9 @@ subroutine disturbance_rates( site_in, bc_in) ! first calculate the fractino of the site that is primary land call get_frac_site_primary(site_in, frac_site_primary) + + ! get available biomass for harvest for all patches + call get_harvestable_carbon(site_in, bc_in%site_area, bc_in%hlm_harvest_catnames, harvestable_forest_c) currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -228,7 +239,9 @@ subroutine disturbance_rates( site_in, bc_in) bc_in%hlm_harvest_units, & currentPatch%anthro_disturbance_label, & currentPatch%age_since_anthro_disturbance, & - frac_site_primary) + frac_site_primary, & + harvestable_forest_c, & + harvest_tag) currentCohort%lmort_direct = lmort_direct currentCohort%lmort_collateral = lmort_collateral @@ -237,9 +250,12 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort => currentCohort%taller end do + currentPatch => currentPatch%younger end do + call get_harvest_debt(site_in, bc_in, harvest_tag) + ! --------------------------------------------------------------------------------------------- ! Calculate Disturbance Rates based on the mortality rates just calculated ! --------------------------------------------------------------------------------------------- @@ -289,6 +305,11 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort%lmort_infra + & currentCohort%l_degrad ) * & currentCohort%c_area/currentPatch%area + + if(currentPatch%disturbance_rates(dtype_ilog)>1.0) then + write(fates_log(),*) 'See luc mortalities:', currentCohort%lmort_direct, & + currentCohort%lmort_collateral, currentCohort%lmort_infra, currentCohort%l_degrad + end if ! Non-harvested part of the logging disturbance rate dist_rate_ldist_notharvested = dist_rate_ldist_notharvested + currentCohort%l_degrad * & @@ -299,13 +320,19 @@ subroutine disturbance_rates( site_in, bc_in) enddo !currentCohort ! for non-closed-canopy areas subject to logging, add an additional increment of area disturbed - ! equivalent to the fradction logged to account for transfer of interstitial ground area to new secondary lands + ! equivalent to the fraction logged to account for transfer of interstitial ground area to new secondary lands if ( logging_time .and. & (currentPatch%area - currentPatch%total_canopy_area) .gt. fates_tiny ) then ! The canopy is NOT closed. - call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, & - bc_in%hlm_harvest_rates, frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate) + if(bc_in%hlm_harvest_units == hlm_harvest_carbon) then + call get_harvest_rate_carbon (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, & + bc_in%hlm_harvest_rates, currentPatch%age_since_anthro_disturbance, harvestable_forest_c, & + harvest_rate, harvest_tag) + else + call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, & + bc_in%hlm_harvest_rates, frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate) + end if currentPatch%disturbance_rates(dtype_ilog) = currentPatch%disturbance_rates(dtype_ilog) + & (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area @@ -315,6 +342,12 @@ subroutine disturbance_rates( site_in, bc_in) (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area endif + ! For nocomp mode, we need to prevent producing too small patches, which may produce small patches + if ((hlm_use_nocomp .eq. itrue) .and. & + (currentPatch%disturbance_rates(dtype_ilog)*currentPatch%area .lt. min_patch_area_forced)) then + currentPatch%disturbance_rates(dtype_ilog) = 0._r8 + end if + ! fraction of the logging disturbance rate that is non-harvested if (currentPatch%disturbance_rates(dtype_ilog) .gt. nearzero) then currentPatch%fract_ldist_not_harvested = dist_rate_ldist_notharvested / & diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 778a02a84d..0c6c296ad2 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -1753,7 +1753,6 @@ subroutine SeedIn( currentSite, bc_in ) ! !USES: use EDTypesMod, only : area use EDTypesMod, only : homogenize_seed_pfts - ! ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 1a2b16283e..dcd4f7bc06 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -2515,7 +2515,6 @@ subroutine ForceDBH( ipft, crowndamage, canopy_trim, d, h, bdead, bl ) end if call h_allom(d,ipft,h) - if(counter>20)then write(fates_log(),*) 'dbh counter: ',counter,' is woody: ',& (prt_params%woody(ipft) == itrue) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 57f64c08cb..a232ea3249 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -303,6 +303,8 @@ subroutine zero_site( site_in ) site_in%fmort_cflux_ustory_damage(:,:) = 0._r8 ! Resources management (logging/harvesting, etc) + site_in%resources_management%harvest_debt = 0.0_r8 + site_in%resources_management%harvest_debt_sec = 0.0_r8 site_in%resources_management%trunk_product_site = 0.0_r8 ! canopy spread diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index a3bf63c419..427a687f8c 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -87,6 +87,7 @@ module EDMainMod use FatesPlantHydraulicsMod , only : AccumulateMortalityWaterStorage use FatesAllometryMod , only : h_allom,tree_sai,tree_lai use EDLoggingMortalityMod , only : IsItLoggingTime + use EDLoggingMortalityMod , only : get_harvestable_carbon use DamageMainMod , only : IsItDamageTime use EDPatchDynamicsMod , only : get_frac_site_primary use FatesGlobals , only : endrun => fates_endrun @@ -322,6 +323,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) ! FIX(SPM,032414) refactor so everything goes through interface ! ! !USES: + use FatesInterfaceTypesMod, only : hlm_num_lu_harvest_cats use FatesInterfaceTypesMod, only : nlevdamage use FatesAllometryMod , only : bleaf use FatesAllometryMod , only : carea_allom @@ -381,6 +383,9 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) real(r8) :: target_leaf_c real(r8) :: frac_site_primary + real(r8) :: harvestable_forest_c(hlm_num_lu_harvest_cats) + integer :: harvest_tag(hlm_num_lu_harvest_cats) + real(r8) :: n_old real(r8) :: n_recover real(r8) :: sapw_c @@ -414,6 +419,13 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) call get_frac_site_primary(currentSite, frac_site_primary) + ! Clear site GPP and AR passing to HLM + bc_out%gpp_site = 0._r8 + bc_out%ar_site = 0._r8 + + ! Patch level biomass are required for C-based harvest + call get_harvestable_carbon(currentSite, bc_in%site_area, bc_in%hlm_harvest_catnames, harvestable_forest_c) + ! Set a pointer to this sites carbon12 mass balance site_cmass => currentSite%mass_balance(element_pos(carbon12_element)) @@ -456,7 +468,6 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) do while(associated(currentCohort)) ft = currentCohort%pft - ! Some cohorts are created and inserted to the list while ! the loop is going. These are pointed to the "taller" position ! of current, and then inherit properties of their donor (current) @@ -465,8 +476,10 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) if_not_newlyrecovered: if(.not.newly_recovered) then + ! Calculate the mortality derivatives - call Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary ) + call Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary, & + harvestable_forest_c, harvest_tag ) ! ----------------------------------------------------------------------------- ! Apply Plant Allocation and Reactive Transport @@ -477,6 +490,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) ! decrement the available carbon pool to zero. ! ----------------------------------------------------------------------------- + if (hlm_use_ed_prescribed_phys .eq. itrue) then if (currentCohort%canopy_layer .eq. 1) then currentCohort%npp_acc = EDPftvarcon_inst%prescribed_npp_canopy(ft) & @@ -485,14 +499,14 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) currentCohort%npp_acc = EDPftvarcon_inst%prescribed_npp_understory(ft) & * currentCohort%c_area / currentCohort%n / hlm_days_per_year endif - + ! We don't explicitly define a respiration rate for prescribe phys ! but we do need to pass mass balance. So we say it is zero respiration currentCohort%gpp_acc = currentCohort%npp_acc currentCohort%resp_acc = 0._r8 - + end if - + ! ----------------------------------------------------------------------------- ! Save NPP/GPP/R in these "hold" style variables. These variables ! persist after this routine is complete, and used in I/O diagnostics. @@ -504,11 +518,17 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) ! _acc will be reset soon and will be accumulated on the next leaf ! photosynthesis step ! ----------------------------------------------------------------------------- - + currentCohort%npp_acc_hold = currentCohort%npp_acc * real(hlm_days_per_year,r8) currentCohort%gpp_acc_hold = currentCohort%gpp_acc * real(hlm_days_per_year,r8) currentCohort%resp_acc_hold = currentCohort%resp_acc * real(hlm_days_per_year,r8) - + + ! Passing gpp_acc_hold to HLM + bc_out%gpp_site = bc_out%gpp_site + currentCohort%gpp_acc_hold * & + AREA_INV * currentCohort%n / hlm_days_per_year / sec_per_day + bc_out%ar_site = bc_out%ar_site + currentCohort%resp_acc_hold * & + AREA_INV * currentCohort%n / hlm_days_per_year / sec_per_day + ! Conduct Maintenance Turnover (parteh) if(debug) call currentCohort%prt%CheckMassConservation(ft,3) if(any(currentSite%dstatus == [phen_dstat_moiston,phen_dstat_timeon])) then @@ -516,9 +536,8 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) else is_drought = .true. end if - call PRTMaintTurnover(currentCohort%prt,ft,is_drought) - + ! ----------------------------------------------------------------------------------- ! Call the routine that advances leaves in age. ! This will move a portion of the leaf mass in each @@ -538,7 +557,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) currentCohort%resp_excess = 0._r8 end if if_not_newlyrecovered - + ! If the current diameter of a plant is somehow less than what is consistent ! with what is allometrically consistent with the stuctural biomass, then ! correct the dbh to match. diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 6c27c2aa72..78f3e88edc 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -613,6 +613,9 @@ module EDTypesMod type, public :: ed_resources_management_type real(r8) :: trunk_product_site ! Actual trunk product at site level KgC/site + real(r8) :: harvest_debt ! the amount of kgC per site that did not successfully harvested + real(r8) :: harvest_debt_sec ! the amount of kgC per site from secondary patches that did + ! not successfully harvested !debug variables real(r8) :: delta_litter_stock ! kgC/site = kgC/ha @@ -781,7 +784,6 @@ module EDTypesMod ! in runs that are restarted, regardless of ! the conditions of restart - real(r8) :: water_memory(numWaterMem) ! last 10 days of soil moisture memory... @@ -1166,6 +1168,10 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%dgmort = ', ccohort%dgmort write(fates_log(),*) 'co%hmort = ', ccohort%hmort write(fates_log(),*) 'co%frmort = ', ccohort%frmort + write(fates_log(),*) 'co%asmort = ', ccohort%asmort + write(fates_log(),*) 'co%lmort_direct = ', ccohort%lmort_direct + write(fates_log(),*) 'co%lmort_collateral = ', ccohort%lmort_collateral + write(fates_log(),*) 'co%lmort_infra = ', ccohort%lmort_infra write(fates_log(),*) 'co%isnew = ', ccohort%isnew write(fates_log(),*) 'co%dndt = ', ccohort%dndt write(fates_log(),*) 'co%dhdt = ', ccohort%dhdt @@ -1177,6 +1183,7 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%cambial_mort = ', ccohort%cambial_mort write(fates_log(),*) 'co%size_class = ', ccohort%size_class write(fates_log(),*) 'co%size_by_pft_class = ', ccohort%size_by_pft_class + if (associated(ccohort%co_hydr) ) then call dump_cohort_hydr(ccohort) endif diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 6c12a78a06..54be8b6e98 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -294,6 +294,12 @@ module FatesHistoryInterfaceMod integer :: ih_canopy_biomass_si integer :: ih_understory_biomass_si + integer :: ih_npp_secondary_si + integer :: ih_gpp_secondary_si + integer :: ih_aresp_secondary_si + integer :: ih_maint_resp_secondary_si + integer :: ih_growth_resp_secondary_si + integer :: ih_primaryland_fusion_error_si integer :: ih_disturbance_rate_p2p_si integer :: ih_disturbance_rate_p2s_si @@ -303,6 +309,8 @@ module FatesHistoryInterfaceMod integer :: ih_fall_disturbance_rate_si integer :: ih_potential_disturbance_rate_si integer :: ih_harvest_carbonflux_si + integer :: ih_harvest_debt_si + integer :: ih_harvest_debt_sec_si ! Indices to site by size-class by age variables integer :: ih_nplant_si_scag @@ -338,7 +346,9 @@ module FatesHistoryInterfaceMod integer :: ih_err_fates_si integer :: ih_npatches_si + integer :: ih_npatches_sec_si integer :: ih_ncohorts_si + integer :: ih_ncohorts_sec_si integer :: ih_demotion_carbonflux_si integer :: ih_promotion_carbonflux_si integer :: ih_canopy_mortality_carbonflux_si @@ -474,6 +484,7 @@ module FatesHistoryInterfaceMod integer :: ih_ddbh_understory_si_scls integer :: ih_agb_si_scls integer :: ih_biomass_si_scls + integer :: ih_mortality_canopy_secondary_si_scls ! mortality vars integer :: ih_m1_si_scls @@ -487,6 +498,14 @@ module FatesHistoryInterfaceMod integer :: ih_m9_si_scls integer :: ih_m10_si_scls + integer :: ih_m1_sec_si_scls + integer :: ih_m2_sec_si_scls + integer :: ih_m3_sec_si_scls + integer :: ih_m7_sec_si_scls + integer :: ih_m8_sec_si_scls + integer :: ih_m9_sec_si_scls + integer :: ih_m10_sec_si_scls + integer :: ih_m10_si_cacls integer :: ih_nplant_si_cacls @@ -536,9 +555,11 @@ module FatesHistoryInterfaceMod ! indices to (site x pft) variables integer :: ih_biomass_si_pft + integer :: ih_biomass_sec_si_pft integer :: ih_leafbiomass_si_pft integer :: ih_storebiomass_si_pft integer :: ih_nindivs_si_pft + integer :: ih_nindivs_sec_si_pft integer :: ih_recruitment_si_pft integer :: ih_mortality_si_pft integer :: ih_mortality_carbonflux_si_pft @@ -548,7 +569,9 @@ module FatesHistoryInterfaceMod integer :: ih_crownarea_si_pft integer :: ih_canopycrownarea_si_pft integer :: ih_gpp_si_pft + integer :: ih_gpp_sec_si_pft integer :: ih_npp_si_pft + integer :: ih_npp_sec_si_pft integer :: ih_nocomp_pftpatchfraction_si_pft integer :: ih_nocomp_pftnpatches_si_pft integer :: ih_nocomp_pftburnedarea_si_pft @@ -572,6 +595,8 @@ module FatesHistoryInterfaceMod integer :: ih_fire_intensity_si_age integer :: ih_fire_sum_fuel_si_age + integer :: ih_lai_secondary_si + ! indices to (site x height) variables integer :: ih_canopy_height_dist_si_height integer :: ih_leaf_height_dist_si_height @@ -2180,7 +2205,9 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) integer :: tmp associate( hio_npatches_si => this%hvars(ih_npatches_si)%r81d, & + hio_npatches_sec_si => this%hvars(ih_npatches_sec_si)%r81d, & hio_ncohorts_si => this%hvars(ih_ncohorts_si)%r81d, & + hio_ncohorts_sec_si => this%hvars(ih_ncohorts_sec_si)%r81d, & hio_trimming_si => this%hvars(ih_trimming_si)%r81d, & hio_area_plant_si => this%hvars(ih_area_plant_si)%r81d, & hio_area_trees_si => this%hvars(ih_area_trees_si)%r81d, & @@ -2189,9 +2216,11 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_ca_weighted_height_si => this%hvars(ih_ca_weighted_height_si)%r81d, & hio_canopy_spread_si => this%hvars(ih_canopy_spread_si)%r81d, & hio_biomass_si_pft => this%hvars(ih_biomass_si_pft)%r82d, & + hio_biomass_sec_si_pft => this%hvars(ih_biomass_sec_si_pft)%r82d, & hio_leafbiomass_si_pft => this%hvars(ih_leafbiomass_si_pft)%r82d, & hio_storebiomass_si_pft => this%hvars(ih_storebiomass_si_pft)%r82d, & hio_nindivs_si_pft => this%hvars(ih_nindivs_si_pft)%r82d, & + hio_nindivs_sec_si_pft => this%hvars(ih_nindivs_sec_si_pft)%r82d, & hio_recruitment_si_pft => this%hvars(ih_recruitment_si_pft)%r82d, & hio_mortality_si_pft => this%hvars(ih_mortality_si_pft)%r82d, & hio_mortality_carbonflux_si_pft => this%hvars(ih_mortality_carbonflux_si_pft)%r82d, & @@ -2201,7 +2230,9 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_crownarea_si_pft => this%hvars(ih_crownarea_si_pft)%r82d, & hio_canopycrownarea_si_pft => this%hvars(ih_canopycrownarea_si_pft)%r82d, & hio_gpp_si_pft => this%hvars(ih_gpp_si_pft)%r82d, & + hio_gpp_sec_si_pft => this%hvars(ih_gpp_sec_si_pft)%r82d, & hio_npp_si_pft => this%hvars(ih_npp_si_pft)%r82d, & + hio_npp_sec_si_pft => this%hvars(ih_npp_sec_si_pft)%r82d, & hio_nesterov_fire_danger_si => this%hvars(ih_nesterov_fire_danger_si)%r81d, & hio_fire_nignitions_si => this%hvars(ih_fire_nignitions_si)%r81d, & hio_fire_fdi_si => this%hvars(ih_fire_fdi_si)%r81d, & @@ -2242,6 +2273,8 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_fall_disturbance_rate_si => this%hvars(ih_fall_disturbance_rate_si)%r81d, & hio_potential_disturbance_rate_si => this%hvars(ih_potential_disturbance_rate_si)%r81d, & hio_harvest_carbonflux_si => this%hvars(ih_harvest_carbonflux_si)%r81d, & + hio_harvest_debt_si => this%hvars(ih_harvest_debt_si)%r81d, & + hio_harvest_debt_sec_si => this%hvars(ih_harvest_debt_sec_si)%r81d, & hio_gpp_si_scpf => this%hvars(ih_gpp_si_scpf)%r82d, & hio_npp_totl_si_scpf => this%hvars(ih_npp_totl_si_scpf)%r82d, & hio_npp_leaf_si_scpf => this%hvars(ih_npp_leaf_si_scpf)%r82d, & @@ -2265,6 +2298,7 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_lai_canopy_si_scpf => this%hvars(ih_lai_canopy_si_scpf)%r82d, & hio_lai_understory_si_scpf => this%hvars(ih_lai_understory_si_scpf)%r82d, & hio_mortality_canopy_si_scpf => this%hvars(ih_mortality_canopy_si_scpf)%r82d, & + hio_mortality_canopy_secondary_si_scls => this%hvars(ih_mortality_canopy_secondary_si_scls)%r82d, & hio_mortality_understory_si_scpf => this%hvars(ih_mortality_understory_si_scpf)%r82d, & hio_m3_mortality_canopy_si_scpf => this%hvars(ih_m3_mortality_canopy_si_scpf)%r82d, & hio_m3_mortality_understory_si_scpf => this%hvars(ih_m3_mortality_understory_si_scpf)%r82d, & @@ -2320,6 +2354,14 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_m10_si_scls => this%hvars(ih_m10_si_scls)%r82d, & hio_m10_si_cacls => this%hvars(ih_m10_si_cacls)%r82d, & + hio_m1_sec_si_scls => this%hvars(ih_m1_sec_si_scls)%r82d, & + hio_m2_sec_si_scls => this%hvars(ih_m2_sec_si_scls)%r82d, & + hio_m3_sec_si_scls => this%hvars(ih_m3_sec_si_scls)%r82d, & + hio_m7_sec_si_scls => this%hvars(ih_m7_sec_si_scls)%r82d, & + hio_m8_sec_si_scls => this%hvars(ih_m8_sec_si_scls)%r82d, & + hio_m9_sec_si_scls => this%hvars(ih_m9_sec_si_scls)%r82d, & + hio_m10_sec_si_scls => this%hvars(ih_m10_sec_si_scls)%r82d, & + hio_c13disc_si_scpf => this%hvars(ih_c13disc_si_scpf)%r82d, & hio_cwd_elcwd => this%hvars(ih_cwd_elcwd)%r82d, & @@ -2384,6 +2426,7 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_yesterdaycanopylevel_understory_si_scls => this%hvars(ih_yesterdaycanopylevel_understory_si_scls)%r82d, & hio_area_si_age => this%hvars(ih_area_si_age)%r82d, & hio_lai_si_age => this%hvars(ih_lai_si_age)%r82d, & + hio_lai_secondary_si => this%hvars(ih_lai_secondary_si)%r81d, & hio_canopy_area_si_age => this%hvars(ih_canopy_area_si_age)%r82d, & hio_ncl_si_age => this%hvars(ih_ncl_si_age)%r82d, & hio_npatches_si_age => this%hvars(ih_npatches_si_age)%r82d, & @@ -2548,6 +2591,8 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) this%hvars(ih_h2oveg_recruit_si)%r81d(io_si) = sites(s)%si_hydr%h2oveg_recruit this%hvars(ih_h2oveg_growturn_err_si)%r81d(io_si) = sites(s)%si_hydr%h2oveg_growturn_err end if + hio_harvest_debt_si(io_si) = sites(s)%resources_management%harvest_debt + hio_harvest_debt_sec_si(io_si) = sites(s)%resources_management%harvest_debt_sec ! error in primary lands from patch fusion [m2 m-2 day-1] -> [m2 m-2 yr-1] hio_primaryland_fusion_error_si(io_si) = sites(s)%primary_land_patchfusion_error * days_per_year @@ -2583,6 +2628,9 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) ! Increment the number of patches per site hio_npatches_si(io_si) = hio_npatches_si(io_si) + 1._r8 + if ( cpatch%anthro_disturbance_label .eq. secondaryforest ) then + hio_npatches_sec_si(io_si) = hio_npatches_sec_si(io_si) + 1._r8 + end if cpatch%age_class = get_age_class_index(cpatch%age) @@ -2595,9 +2643,9 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) (cpatch%tveg24%GetMean()- t_water_freeze_k_1atm)*cpatch%area*AREA_INV ! Increment some patch-age-resolved diagnostics - hio_lai_si_age(io_si,cpatch%age_class) = hio_lai_si_age(io_si,cpatch%age_class) & + sum(cpatch%tlai_profile(:,:,:)) * cpatch%area + hio_ncl_si_age(io_si,cpatch%age_class) = hio_ncl_si_age(io_si,cpatch%age_class) & + cpatch%ncl_p * cpatch%area hio_npatches_si_age(io_si,cpatch%age_class) = hio_npatches_si_age(io_si,cpatch%age_class) + 1._r8 @@ -2623,6 +2671,12 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) + cpatch%area * AREA_INV endif + ! Secondary forest mean LAI + if ( cpatch%anthro_disturbance_label .eq. secondaryforest ) then + hio_lai_secondary_si(io_si) = hio_lai_secondary_si(io_si) & + + sum(cpatch%tlai_profile(:,:,:)) * cpatch%total_canopy_area + end if + ! patch-age-resolved fire variables do i_pft = 1,numpft ! for scorch height, weight the value by patch area within any @@ -2688,6 +2742,10 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) ! Increment the number of cohorts per site hio_ncohorts_si(io_si) = hio_ncohorts_si(io_si) + 1._r8 + if ( cpatch%anthro_disturbance_label .eq. secondaryforest ) then + hio_ncohorts_sec_si(io_si) = hio_ncohorts_sec_si(io_si) + 1._r8 + end if + n_perm2 = ccohort%n * AREA_INV hio_canopy_area_si_age(io_si,cpatch%age_class) = hio_canopy_area_si_age(io_si,cpatch%age_class) & @@ -2800,9 +2858,19 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_nindivs_si_pft(io_si,ft) = hio_nindivs_si_pft(io_si,ft) + & ccohort%n * AREA_INV + if ( cpatch%anthro_disturbance_label .eq. secondaryforest ) then + hio_nindivs_sec_si_pft(io_si,ft) = hio_nindivs_sec_si_pft(io_si,ft) + & + ccohort%n * AREA_INV + end if + hio_biomass_si_pft(io_si, ft) = hio_biomass_si_pft(io_si, ft) + & (ccohort%n * AREA_INV) * total_m + if ( cpatch%anthro_disturbance_label .eq. secondaryforest ) then + hio_biomass_sec_si_pft(io_si, ft) = hio_biomass_sec_si_pft(io_si, ft) + & + (ccohort%n * AREA_INV) * total_m + end if + ! update total biomass per age bin hio_biomass_si_age(io_si,cpatch%age_class) = hio_biomass_si_age(io_si,cpatch%age_class) & + total_m * ccohort%n * AREA_INV @@ -2929,6 +2997,13 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_npp_si_pft(io_si, ft) = hio_npp_si_pft(io_si, ft) + & ccohort%npp_acc_hold * n_perm2 / days_per_year / sec_per_day + if ( cpatch%anthro_disturbance_label .eq. secondaryforest ) then + hio_gpp_sec_si_pft(io_si, ft) = hio_gpp_sec_si_pft(io_si, ft) + & + ccohort%gpp_acc_hold * n_perm2 / days_per_year / sec_per_day + hio_npp_sec_si_pft(io_si, ft) = hio_npp_sec_si_pft(io_si, ft) + & + ccohort%npp_acc_hold * n_perm2 / days_per_year / sec_per_day + end if + ! Site by Size-Class x PFT (SCPF) ! ------------------------------------------------------------------------ @@ -3067,6 +3142,23 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) ccohort%frmort*ccohort%n / m2_per_ha hio_m9_si_scls(io_si,scls) = hio_m9_si_scls(io_si,scls) + ccohort%smort*ccohort%n / m2_per_ha + ! Examine secondary forest mortality and mortality rates + if(cpatch%anthro_disturbance_label .eq. secondaryforest) then + + if (hlm_use_cohort_age_tracking .eq.itrue) then + hio_m10_sec_si_scls(io_si,scls) = hio_m10_sec_si_scls(io_si,scls) + & + ccohort%asmort*ccohort%n / m2_per_ha + end if + + hio_m1_sec_si_scls(io_si,scls) = hio_m1_sec_si_scls(io_si,scls) + ccohort%bmort*ccohort%n / m2_per_ha + hio_m2_sec_si_scls(io_si,scls) = hio_m2_sec_si_scls(io_si,scls) + ccohort%hmort*ccohort%n / m2_per_ha + hio_m3_sec_si_scls(io_si,scls) = hio_m3_sec_si_scls(io_si,scls) + ccohort%cmort*ccohort%n / m2_per_ha + hio_m7_sec_si_scls(io_si,scls) = hio_m7_sec_si_scls(io_si,scls) + & + (ccohort%lmort_direct+ccohort%lmort_collateral+ccohort%lmort_infra) * ccohort%n / m2_per_ha + hio_m8_sec_si_scls(io_si,scls) = hio_m8_sec_si_scls(io_si,scls) + & + ccohort%frmort*ccohort%n / m2_per_ha + hio_m9_sec_si_scls(io_si,scls) = hio_m9_sec_si_scls(io_si,scls) + ccohort%smort*ccohort%n / m2_per_ha + end if !C13 discrimination if(gpp_cached + ccohort%gpp_acc_hold > 0.0_r8)then @@ -3339,6 +3431,15 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & ccohort%n * sec_per_day * days_per_year / m2_per_ha + if ( cpatch%anthro_disturbance_label .eq. secondaryforest ) then + hio_mortality_canopy_secondary_si_scls(io_si,scls) = hio_mortality_canopy_secondary_si_scls(io_si,scls) + & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + & + ccohort%frmort + ccohort%smort + ccohort%asmort) * ccohort%n / m2_per_ha + & + (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & + ccohort%n * sec_per_day * days_per_year / m2_per_ha + end if + + hio_nplant_understory_si_scpf(io_si,scpf) = hio_nplant_understory_si_scpf(io_si,scpf) + ccohort%n / m2_per_ha hio_nplant_understory_si_scls(io_si,scls) = hio_nplant_understory_si_scls(io_si,scls) + ccohort%n / m2_per_ha hio_lai_understory_si_scls(io_si,scls) = hio_lai_understory_si_scls(io_si,scls) + & @@ -3577,6 +3678,13 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) endif end do + ! divide secondary plant leaf area by secondary forest area to get the secondary forest LAI + if (hio_fraction_secondary_forest_si(io_si) .gt. nearzero) then + hio_lai_secondary_si(io_si) = hio_lai_secondary_si(io_si) / (hio_fraction_secondary_forest_si(io_si)*AREA) + else + hio_lai_secondary_si(io_si) = 0._r8 + end if + ! pass the cohort termination mortality as a flux to the history, and then reset the termination mortality buffer ! note there are various ways of reporting the total mortality, so pass to these as well do i_pft = 1, numpft @@ -3639,6 +3747,12 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_mortality_canopy_si_scls(io_si,i_scls) = hio_mortality_canopy_si_scls(io_si,i_scls) + & sites(s)%fmort_rate_canopy(i_scls, i_pft) / m2_per_ha + ! Shijie: Think about how to add later? + !if ( cpatch%anthro_disturbance_label .eq. secondaryforest ) then + ! hio_mortality_canopy_secondary_si_scls(io_si,i_scls) = hio_mortality_canopy_secondary_si_scls(io_si,i_scls) + & + ! sites(s)%term_nindivs_canopy(i_scls,i_pft) * days_per_year / m2_per_ha + !end if + ! the fire mortality rates for each layer are total dead, since the usable ! output will then normalize by the counts, we are allowed to sum over layers hio_mortality_understory_si_scpf(io_si,i_scpf) = hio_mortality_understory_si_scpf(io_si,i_scpf) + & @@ -4183,10 +4297,15 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) real(r8) :: per_dt_tstep ! Time step in frequency units (/s) associate( hio_gpp_si => this%hvars(ih_gpp_si)%r81d, & + hio_gpp_secondary_si => this%hvars(ih_gpp_secondary_si)%r81d, & hio_npp_si => this%hvars(ih_npp_si)%r81d, & + hio_npp_secondary_si => this%hvars(ih_npp_secondary_si)%r81d, & hio_aresp_si => this%hvars(ih_aresp_si)%r81d, & + hio_aresp_secondary_si => this%hvars(ih_aresp_secondary_si)%r81d, & hio_maint_resp_si => this%hvars(ih_maint_resp_si)%r81d, & + hio_maint_resp_secondary_si => this%hvars(ih_maint_resp_secondary_si)%r81d, & hio_growth_resp_si => this%hvars(ih_growth_resp_si)%r81d, & + hio_growth_resp_secondary_si => this%hvars(ih_growth_resp_secondary_si)%r81d, & hio_c_stomata_si => this%hvars(ih_c_stomata_si)%r81d, & hio_c_lblayer_si => this%hvars(ih_c_lblayer_si)%r81d, & hio_rad_error_si => this%hvars(ih_rad_error_si)%r81d, & @@ -4306,12 +4425,12 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) hio_rad_error_si(io_si) = hio_rad_error_si(io_si) + & cpatch%radiation_error * cpatch%area * AREA_INV - ! Only accumulate the instantaneous vegetation temperature for vegetated patches - if (cpatch%patchno .ne. 0) then - hio_tveg(io_si) = hio_tveg(io_si) + & - (bc_in(s)%t_veg_pa(cpatch%patchno) - t_water_freeze_k_1atm) * & - cpatch%area / site_area_veg - end if + ! Only accumulate the instantaneous vegetation temperature for vegetated patches + if (cpatch%patchno .ne. 0) then + hio_tveg(io_si) = hio_tveg(io_si) + & + (bc_in(s)%t_veg_pa(cpatch%patchno) - t_water_freeze_k_1atm) * & + cpatch%area / site_area_veg + end if ccohort => cpatch%shortest do while(associated(ccohort)) @@ -4341,6 +4460,20 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep) hio_maint_resp_si(io_si) = hio_maint_resp_si(io_si) + & ccohort%resp_m * n_perm2 * per_dt_tstep + ! Secondary forest only + if ( cpatch%anthro_disturbance_label .eq. secondaryforest ) then + hio_npp_secondary_si(io_si) = hio_npp_secondary_si(io_si) + & + npp * n_perm2 * per_dt_tstep + hio_gpp_secondary_si(io_si) = hio_gpp_secondary_si(io_si) + & + ccohort%gpp_tstep * n_perm2 * per_dt_tstep + hio_aresp_secondary_si(io_si) = hio_aresp_secondary_si(io_si) + & + aresp * n_perm2 * per_dt_tstep + hio_growth_resp_secondary_si(io_si) = hio_growth_resp_secondary_si(io_si) + & + resp_g * n_perm2 * per_dt_tstep + hio_maint_resp_secondary_si(io_si) = hio_maint_resp_secondary_si(io_si) + & + ccohort%resp_m * n_perm2 * per_dt_tstep + end if + ! Add up the total Net Ecosystem Production ! for this timestep. [kgC/m2/s] hio_nep_si(io_si) = hio_nep_si(io_si) + & @@ -5068,6 +5201,19 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_ncohorts_si) + call this%set_history_var(vname='FATES_NPATCHES_SECONDARY', units='', & + long='total number of patches per site', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_npatches_sec_si) + + call this%set_history_var(vname='FATES_NCOHORTS_SECONDARY', units='', & + long='total number of cohorts per site', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_ncohorts_sec_si) + + ! Patch variables call this%set_history_var(vname='FATES_TRIMMING', units='1', & long='degree to which canopy expansion is limited by leaf economics (0-1)', & use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & @@ -5181,6 +5327,12 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_biomass_si_pft) + call this%set_history_var(vname='FATES_VEGC_SE_PF', units='kg m-2', & + long='total PFT-level biomass in kg of carbon per land area, secondary patches', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_biomass_sec_si_pft) + call this%set_history_var(vname='FATES_LEAFC_PF', units='kg m-2', & long='total PFT-level leaf biomass in kg carbon per m2 land area', & use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & @@ -5217,12 +5369,30 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_npp_si_pft) + call this%set_history_var(vname='FATES_GPP_SE_PF', units='kg m-2 s-1', & + long='total PFT-level GPP in kg carbon per m2 land area per second, secondary patches', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_gpp_sec_si_pft) + + call this%set_history_var(vname='FATES_NPP_SE_PF', units='kg m-2 yr-1', & + long='total PFT-level NPP in kg carbon per m2 land area per second, secondary patches', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_npp_sec_si_pft) + call this%set_history_var(vname='FATES_NPLANT_PF', units='m-2', & long='total PFT-level number of individuals per m2 land area', & use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_nindivs_si_pft) + call this%set_history_var(vname='FATES_NPLANT_SEC_PF', units='m-2', & + long='total PFT-level number of individuals per m2 land area, secondary patches', & + use_default='active', avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_nindivs_sec_si_pft) + call this%set_history_var(vname='FATES_RECRUITMENT_PF', & units='m-2 yr-1', & long='PFT-level recruitment rate in number of individuals per m2 land area per year', & @@ -5268,6 +5438,12 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_lai_si_age) + call this%set_history_var(vname='FATES_LAI_SECONDARY', units='m2 m-2', & + long='leaf area index per m2 land area, secondary patches', & + use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index=ih_lai_secondary_si) + call this%set_history_var(vname='FATES_CANOPYAREA_AP', units='m2 m-2', & long='canopy area by age bin per m2 land area', use_default='active', & avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', upfreq=1, ivar=ivar, & @@ -5319,20 +5495,20 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_SECONDARY_FOREST_FRACTION', & units='m2 m-2', long='secondary forest fraction', & - use_default='inactive', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_fraction_secondary_forest_si) call this%set_history_var(vname='FATES_WOOD_PRODUCT', units='kg m-2', & long='total wood product from logging in kg carbon per m2 land area', & - use_default='inactive', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_woodproduct_si) call this%set_history_var(vname='FATES_SECONDARY_FOREST_VEGC', & units='kg m-2', & long='biomass on secondary lands in kg carbon per m2 land area (mult by FATES_SECONDARY_FOREST_FRACTION to get per secondary forest area)', & - use_default='inactive', avgflag='A', vtype=site_r8, & + use_default='active', avgflag='A', vtype=site_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_biomass_secondary_forest_si) @@ -6053,13 +6229,24 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=2, & ivar=ivar, initialize=initialize_variables, index = ih_tveg_si ) - ! radiation error + ! radiation error - call this%set_history_var(vname='FATES_RAD_ERROR', units='W m-2 ', & + call this%set_history_var(vname='FATES_RAD_ERROR', units='W m-2 ', & long='radiation error in FATES RTM', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=2, & ivar=ivar, initialize=initialize_variables, index = ih_rad_error_si) + call this%set_history_var(vname='FATES_HARVEST_DEBT', units='kg C', & + long='Accumulated carbon failed to be harvested', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_harvest_debt_si ) + + call this%set_history_var(vname='FATES_HARVEST_DEBT_SEC', units='kg C', & + long='Accumulated carbon failed to be harvested from secondary patches', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_harvest_debt_sec_si ) + + ! Ecosystem Carbon Fluxes (updated rapidly, upfreq=2) @@ -6068,28 +6255,55 @@ subroutine define_history_vars(this, initialize_variables) use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & upfreq=2, ivar=ivar, initialize=initialize_variables, index = ih_npp_si) + call this%set_history_var(vname='FATES_NPP_SECONDARY', units='kg m-2 s-1', & + long='net primary production in kg carbon per m2 per second, secondary patches', & + use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=2, ivar=ivar, initialize=initialize_variables, index = ih_npp_secondary_si) + call this%set_history_var(vname='FATES_GPP', units='kg m-2 s-1', & long='gross primary production in kg carbon per m2 per second', & use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & upfreq=2, ivar=ivar, initialize=initialize_variables, index = ih_gpp_si) + call this%set_history_var(vname='FATES_GPP_SECONDARY', units='kg m-2 s-1', & + long='gross primary production in kg carbon per m2 per second, secondary patches', & + use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=2, ivar=ivar, initialize=initialize_variables, index = ih_gpp_secondary_si) + call this%set_history_var(vname='FATES_AUTORESP', units='kg m-2 s-1', & long='autotrophic respiration in kg carbon per m2 per second', & use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & upfreq=2, ivar=ivar, initialize=initialize_variables, index = ih_aresp_si) + call this%set_history_var(vname='FATES_AUTORESP_SECONDARY', units='kg m-2 s-1', & + long='autotrophic respiration in kg carbon per m2 per second, secondary patches', & + use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=2, ivar=ivar, initialize=initialize_variables, index = ih_aresp_secondary_si) + call this%set_history_var(vname='FATES_GROWTH_RESP', units='kg m-2 s-1', & long='growth respiration in kg carbon per m2 per second', & use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & upfreq=2, ivar=ivar, initialize=initialize_variables, & index = ih_growth_resp_si) + call this%set_history_var(vname='FATES_GROWTH_RESP_SECONDARY', units='kg m-2 s-1', & + long='growth respiration in kg carbon per m2 per second, secondary patches', & + use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=2, ivar=ivar, initialize=initialize_variables, & + index = ih_growth_resp_secondary_si) + call this%set_history_var(vname='FATES_MAINT_RESP', units='kg m-2 s-1', & - long='maintenance respiration in kg carbon per m2 land area per second', & + long='maintenance respiration in kg carbon per m2 land area per second, secondary patches', & use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & upfreq=2, ivar=ivar, initialize=initialize_variables, & index = ih_maint_resp_si) + call this%set_history_var(vname='FATES_MAINT_RESP_SECONDARY', units='kg m-2 s-1', & + long='maintenance respiration in kg carbon per m2 land area per second', & + use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=2, ivar=ivar, initialize=initialize_variables, & + index = ih_maint_resp_secondary_si) + call this%set_history_var(vname='FATES_EXCESS_RESP', units='kg m-2 s-1', & long='respiration of un-allocatable carbon gain', & use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & @@ -7004,6 +7218,13 @@ subroutine define_history_vars(this, initialize_variables) hlms='CLM:ALM', upfreq=1, ivar=ivar, & initialize=initialize_variables, index = ih_mortality_canopy_si_scls) + call this%set_history_var(vname='FATES_MORTALITY_CANOPY_SE_SZ', & + units = 'm-2 yr-1', & + long='total mortality of canopy trees by size class in number of plants per m2, secondary patches', & + use_default='active', avgflag='A', vtype=site_size_r8, & + hlms='CLM:ALM', upfreq=1, ivar=ivar, & + initialize=initialize_variables, index = ih_mortality_canopy_secondary_si_scls) + call this%set_history_var(vname='FATES_NPLANT_USTORY_SZ', & units = 'm-2', & long='number of understory plants per m2 by size class', & @@ -7071,6 +7292,27 @@ subroutine define_history_vars(this, initialize_variables) hlms='CLM:ALM', upfreq=1, ivar=ivar, & initialize=initialize_variables, index = ih_m3_si_scls) + call this%set_history_var(vname='FATES_MORTALITY_BACKGROUND_SE_SZ', & + units = 'm-2 yr-1', & + long='background mortality by size in number of plants per m2 per year, secondary patches', & + use_default='active', avgflag='A', vtype=site_size_r8, & + hlms='CLM:ALM', upfreq=1, ivar=ivar, & + initialize=initialize_variables, index = ih_m1_sec_si_scls) + + call this%set_history_var(vname='FATES_MORTALITY_HYDRAULIC_SE_SZ', & + units = 'm-2 yr-1', & + long='hydraulic mortality by size in number of plants per m2 per year, secondary patches', & + use_default='active', avgflag='A', vtype=site_size_r8, & + hlms='CLM:ALM', upfreq=1, ivar=ivar, & + initialize=initialize_variables, index = ih_m2_sec_si_scls) + + call this%set_history_var(vname='FATES_MORTALITY_CSTARV_SE_SZ', & + units = 'm-2 yr-1', & + long='carbon starvation mortality by size in number of plants per m2 per year, secondary patches', & + use_default='active', avgflag='A', vtype=site_size_r8, & + hlms='CLM:ALM', upfreq=1, ivar=ivar, & + initialize=initialize_variables, index = ih_m3_sec_si_scls) + call this%set_history_var(vname='FATES_MORTALITY_IMPACT_SZ', & units = 'm-2 yr-1', & long='impact mortality by size in number of plants per m2 per year', & @@ -7127,6 +7369,34 @@ subroutine define_history_vars(this, initialize_variables) hlms='CLM:ALM', upfreq=1, ivar=ivar, & initialize=initialize_variables, index = ih_m10_si_cacls) + call this%set_history_var(vname='FATES_MORTALITY_LOGGING_SE_SZ', & + units = 'm-2 yr-1', & + long='logging mortality by size in number of plants per m2 per event, secondary patches', & + use_default='active', avgflag='A', vtype=site_size_r8, & + hlms='CLM:ALM', upfreq=1, ivar=ivar, & + initialize=initialize_variables, index = ih_m7_sec_si_scls) + + call this%set_history_var(vname='FATES_MORTALITY_FREEZING_SE_SZ', & + units = 'm-2 event-1', & + long='freezing mortality by size in number of plants per m2 per event, secondary patches', & + use_default='active', avgflag='A', vtype=site_size_r8, & + hlms='CLM:ALM', upfreq=1, ivar=ivar, & + initialize=initialize_variables, index = ih_m8_sec_si_scls) + + call this%set_history_var(vname='FATES_MORTALITY_SENESCENCE_SE_SZ', & + units = 'm-2 yr-1', & + long='senescence mortality by size in number of plants per m2 per event, secondary patches', & + use_default='active', avgflag='A', vtype=site_size_r8, & + hlms='CLM:ALM', upfreq=1, ivar=ivar, & + initialize=initialize_variables, index = ih_m9_sec_si_scls) + + call this%set_history_var(vname='FATES_MORTALITY_AGESCEN_SE_SZ', & + units = 'm-2 yr-1', & + long='age senescence mortality by size in number of plants per m2 per year, secondary patches', & + use_default='active', avgflag='A', vtype=site_size_r8, & + hlms='CLM:ALM', upfreq=1, ivar=ivar, & + initialize=initialize_variables, index = ih_m10_sec_si_scls) + call this%set_history_var(vname='FATES_NPP_CANOPY_SZ', units = 'kg m-2 s-1', & long='NPP of canopy plants by size class in kg carbon per m2 per second', & use_default='inactive', avgflag='A', vtype=site_size_r8, & diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index c322fc5a6a..b314073e60 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -368,6 +368,9 @@ subroutine zero_bcs(fates,s) end if fates%bc_out(s)%plant_stored_h2o_si = 0.0_r8 + ! Land Use realated + fates%bc_out(s)%gpp_site = 0.0_r8 + fates%bc_out(s)%ar_site = 0.0_r8 fates%bc_out(s)%hrv_deadstemc_to_prod10c = 0.0_r8 fates%bc_out(s)%hrv_deadstemc_to_prod100c = 0.0_r8 diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index dd306e9734..e6aedc08f7 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -531,6 +531,8 @@ module FatesInterfaceTypesMod character(len=64), allocatable :: hlm_harvest_catnames(:) ! names of hlm_harvest d1 integer :: hlm_harvest_units ! what units are the harvest rates specified in? [area vs carbon] + + real(r8) :: site_area ! Actual area of current site [m2], only used in carbon-based harvest ! Fixed biogeography mode real(r8), allocatable :: pft_areafrac(:) ! Fractional area of the FATES column occupied by each PFT @@ -728,7 +730,9 @@ module FatesInterfaceTypesMod ! FATES LULCC real(r8) :: hrv_deadstemc_to_prod10c ! Harvested C flux to 10-yr wood product pool [Site-Level, gC m-2 s-1] real(r8) :: hrv_deadstemc_to_prod100c ! Harvested C flux to 100-yr wood product pool [Site-Level, gC m-2 s-1] - + real(r8) :: gpp_site ! Site level GPP, for NBP diagnosis in HLM [Site-Level, gC m-2 s-1] + real(r8) :: ar_site ! Site level Autotrophic Resp, for NBP diagnosis in HLM [Site-Level, gC m-2 s-1] + end type bc_out_type diff --git a/parteh/PRTGenericMod.F90 b/parteh/PRTGenericMod.F90 index 2075216da2..ca8d337b31 100644 --- a/parteh/PRTGenericMod.F90 +++ b/parteh/PRTGenericMod.F90 @@ -878,16 +878,17 @@ subroutine DeallocatePRTVartypes(this) ! Check to see if there is any value in these pools? ! SHould not deallocate if there is any carbon left - do i_var = 1, prt_global%num_vars - deallocate(this%variables(i_var)%val) - deallocate(this%variables(i_var)%val0) - deallocate(this%variables(i_var)%net_alloc) - deallocate(this%variables(i_var)%turnover) - deallocate(this%variables(i_var)%burned) - deallocate(this%variables(i_var)%damaged) - end do + if(allocated(this%variables)) then + do i_var = 1, prt_global%num_vars + deallocate(this%variables(i_var)%val) + deallocate(this%variables(i_var)%val0) + deallocate(this%variables(i_var)%net_alloc) + deallocate(this%variables(i_var)%turnover) + deallocate(this%variables(i_var)%burned) + end do - deallocate(this%variables) + deallocate(this%variables) + end if if(allocated(this%bc_in))then deallocate(this%bc_in)