diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 12385a8f9d..f403533748 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -1124,6 +1124,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) nextc%n*nextc%lmort_infra)/newn currentCohort%l_degrad = (currentCohort%n*currentCohort%l_degrad + & nextc%n*nextc%l_degrad)/newn + currentCohort%harv_c = (currentCohort%n*currentCohort%harv_c + & + nextc%n*nextc%harv_c)/newn ! biomass and dbh tendencies currentCohort%ddbhdt = (currentCohort%n*currentCohort%ddbhdt + & diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 old mode 100644 new mode 100755 index bf6ab7443c..c68d89d9e9 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -68,7 +68,13 @@ module EDLoggingMortalityMod use FatesConstantsMod , only : months_per_year, days_per_sec, years_per_day, g_per_kg use FatesConstantsMod , only : hlm_harvest_area_fraction use FatesConstantsMod , only : hlm_harvest_carbon - use FatesConstantsMod, only : fates_check_param_set + use FatesConstantsMod , only : fates_check_param_set + use FatesConstantsMod , only : logging_uniform_size + use FatesConstantsMod , only : logging_double_rotation + use FatesConstantsMod , only : logging_quadruple_rotation + use FatesConstantsMod , only : logging_logistic_size + use FatesConstantsMod , only : logging_inv_logistic_size + use FatesConstantsMod , only : logging_gaussian_size implicit none private @@ -96,11 +102,11 @@ 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 - + private :: get_cohort_harvest_rate_scale + contains subroutine IsItLoggingTime(is_master,currentSite) @@ -200,8 +206,9 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & hlm_harvest_rates, hlm_harvest_catnames, & hlm_harvest_units, & patch_land_use_label, secondary_age, & - frac_site_primary, harvestable_forest_c, & - harvest_tag) + frac_site_primary, frac_site_secondary_mature, & + harvestable_forest_c, harvest_rate_scale, harvest_tag, & + harvest_wp_scale, logging_preference_options ) ! Arguments integer, intent(in) :: pft_i ! pft index @@ -215,6 +222,13 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & 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(in) :: frac_site_secondary_mature + real(r8), intent(in) :: harvest_rate_scale ! scaling factor which increase/decrease + ! the disturbance rate proportionally considering patch age + real(r8), intent(in) :: harvest_wp_scale ! For IFM: scaling factor which increase/decrease + ! the harvest wood product, not apply to the disturbance area/rate + integer, intent(in) :: logging_preference_options ! Options for IFM management strategy + 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 @@ -232,6 +246,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! 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 + real(r8) :: harvest_rate_scale_cohort ! scaling factor after considering dbh size for each cohort ! todo: probably lower the dbhmin default value to 30 cm ! todo: change the default logging_event_code to 1 september (-244) @@ -240,15 +255,19 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! todo: eventually set up distinct harvest practices, each with a set of input paramaeters ! todo: implement harvested carbon inputs + ! For area-based harvest, harvest_tag shall always be 2 (not applicable). + ! For carbon-based harvest, initialize harvest_tag from 2 + harvest_tag = 2 + cur_harvest_tag = 2 + if (logging_time) then ! Pass logging rates to cohort level - if (hlm_use_lu_harvest == ifalse) then ! 0=use fates logging parameters directly when logging_time == .true. ! this means harvest the whole cohort area harvest_rate = 1._r8 - + else if (hlm_use_lu_harvest == itrue .and. hlm_harvest_units == hlm_harvest_area_fraction) then ! We are harvesting based on areal fraction, not carbon/biomass terms. ! 1=use area fraction from hlm @@ -266,15 +285,11 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! Get the area-based harvest rates based on info passed to FATES from the boundary condition call get_harvest_rate_area (patch_land_use_label, hlm_harvest_catnames, & - hlm_harvest_rates, frac_site_primary, secondary_age, harvest_rate) + hlm_harvest_rates, frac_site_primary, frac_site_secondary_mature, 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 + if (fates_global_verbose()) then write(fates_log(), *) 'Successfully Read Harvest Rate from HLM.', hlm_harvest_rates(:), harvest_rate - end if + end if else if (hlm_use_lu_harvest == itrue .and. hlm_harvest_units == hlm_harvest_carbon) then ! 2=use carbon from hlm @@ -290,17 +305,27 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & endif + + ! Logging size prefrence or change rotation length + call get_cohort_harvest_rate_scale (logging_preference_options, dbh, harvest_wp_scale, harvest_rate_scale, & + harvest_rate_scale_cohort) + ! 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 - if (cur_harvest_tag == 0) then + if (cur_harvest_tag == 0 .or. cur_harvest_tag == 2) 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 + ! For IFM calculation we may encounter issue with inappropriately high target wood product which exceeds supply + ! in which case, harvest_rate_scale_cohort >> harvest_rate_scale. We added the min() function to just + ! truncate any lmort_direct beyond supply and maximize the lmort_direct first, so under this circumstance + ! lmort_collateral and lmort_infra might not be the same value as settled in parameter file + ! Need further improvement for a more general case + lmort_direct = min(harvest_rate_scale_cohort * harvest_rate * logging_direct_frac, harvest_rate_scale * harvest_rate) else lmort_direct = 0.0_r8 end if @@ -312,13 +337,15 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & if (dbh >= logging_dbhmax_infra) then lmort_infra = 0.0_r8 else - lmort_infra = harvest_rate * logging_mechanical_frac + lmort_infra = min(harvest_rate_scale_cohort * harvest_rate * logging_mechanical_frac, & + harvest_rate_scale*harvest_rate - lmort_direct) end if ! Collateral damage to smaller plants below the direct logging size threshold ! will be applied via "understory_death" via the disturbance algorithm if (canopy_layer .eq. 1) then - lmort_collateral = harvest_rate * logging_collateral_frac + lmort_collateral = min(harvest_rate_scale_cohort * harvest_rate * logging_collateral_frac, & + harvest_rate_scale*harvest_rate - lmort_direct - lmort_infra) else lmort_collateral = 0._r8 endif @@ -326,12 +353,13 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & else ! non-woody plants still killed by infrastructure lmort_direct = 0.0_r8 lmort_collateral = 0.0_r8 - lmort_infra = harvest_rate * logging_mechanical_frac + lmort_infra = min(harvest_rate_scale_cohort * harvest_rate * logging_mechanical_frac, & + harvest_rate_scale*harvest_rate) end if ! the area occupied by all plants in the canopy that aren't killed is still disturbed at the harvest rate if (canopy_layer .eq. 1) then - l_degrad = harvest_rate - (lmort_direct + lmort_infra + lmort_collateral) ! fraction passed to 'degraded' forest. + l_degrad = harvest_rate_scale * harvest_rate - (lmort_direct + lmort_infra + lmort_collateral) ! fraction passed to 'degraded' forest. else l_degrad = 0._r8 endif @@ -345,11 +373,10 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & end subroutine LoggingMortality_frac - ! ============================================================================ subroutine get_harvest_rate_area (patch_land_use_label, hlm_harvest_catnames, hlm_harvest_rates, & - frac_site_primary, secondary_age, harvest_rate) + frac_site_primary, frac_site_secondary_mature, secondary_age, harvest_rate) ! ------------------------------------------------------------------------------------------- @@ -364,6 +391,7 @@ subroutine get_harvest_rate_area (patch_land_use_label, hlm_harvest_catnames, hl 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(in) :: frac_site_secondary_mature real(r8), intent(out) :: harvest_rate ! Local Variables @@ -393,21 +421,37 @@ subroutine get_harvest_rate_area (patch_land_use_label, hlm_harvest_catnames, hl endif end do - ! Normalize by site-level primary or secondary forest fraction + ! Normalize by site-level primary, secondary mature and secondary young 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_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) + harvest_rate = min((harvest_rate / frac_site_primary), 0.99_r8) else harvest_rate = 0._r8 endif else - if ((1._r8-frac_site_primary) .gt. fates_tiny) then - harvest_rate = min((harvest_rate / (1._r8-frac_site_primary)),& - (1._r8-frac_site_primary)) + ! if ((1._r8-frac_site_primary) .gt. fates_tiny) then + ! harvest_rate = min((harvest_rate / (1._r8-frac_site_primary)),& + ! (1._r8-frac_site_primary)) + ! else + ! harvest_rate = 0._r8 + ! endif + if(secondary_age >= secondary_age_threshold) then + ! Secondary mature + if (frac_site_secondary_mature .gt. fates_tiny) then + harvest_rate = min((harvest_rate / frac_site_secondary_mature), 0.99_r8) + else + harvest_rate = 0._r8 + endif else - harvest_rate = 0._r8 + ! Secondary young + if ((1._r8-frac_site_primary-frac_site_secondary_mature) .gt. fates_tiny) then + harvest_rate = min((harvest_rate / (1._r8 - frac_site_primary -& + frac_site_secondary_mature)) , 0.99_r8) + else + harvest_rate = 0._r8 + endif endif endif @@ -431,110 +475,6 @@ subroutine get_harvest_rate_area (patch_land_use_label, hlm_harvest_catnames, hl 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(fates_patch_type), pointer :: currentPatch - type(fates_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%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%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%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 - 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_land_use_label, hlm_harvest_catnames, & @@ -577,7 +517,6 @@ subroutine get_harvest_rate_carbon (patch_land_use_label, hlm_harvest_catnames, 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 @@ -790,6 +729,7 @@ subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site end if + !write(fates_log(), *) 'Wood product before.', currentSite%mass_balance(1)%wood_product do el = 1,num_elements element_id = element_list(el) @@ -874,7 +814,7 @@ subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site !adjust how wood is partitioned between the cwd classes based on cohort dbh call adjust_SF_CWD_frac(currentCohort%dbh,ncwd,SF_val_CWD_frac,SF_val_CWD_frac_adj) - do c = 1,ncwd-1 + do c = 1,ncwd-1 new_litt%ag_cwd(c) = new_litt%ag_cwd(c) + & ag_wood * SF_val_CWD_frac_adj(c) * donate_m2 @@ -986,6 +926,15 @@ subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site site_mass%wood_product = site_mass%wood_product + & ag_wood * logging_export_frac + site_mass%wood_product_sz(currentCohort%size_class) = site_mass%wood_product_sz(currentCohort%size_class) + & + ag_wood * logging_export_frac + ! Transfer to kg per indiv + !if(logging_time) then + currentCohort%harv_c = ag_wood / currentCohort%n * logging_export_frac + !else + ! currentCohort%harv_c = 0._r8 + !end if + new_litt%ag_cwd(ncwd) = new_litt%ag_cwd(ncwd) + ag_wood * & (1._r8-logging_export_frac)*donate_m2 @@ -1077,6 +1026,7 @@ subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site end if end do + write(fates_log(),*) 'How large is wood product pool?', site_mass%wood_product ! Not sure why this is called here, but I suppose it can't hurt ! (rgk 06-2019) @@ -1134,6 +1084,8 @@ subroutine UpdateHarvestC(currentSite,bc_out) end subroutine UpdateHarvestC + ! ===================================================================================== + subroutine get_harvest_debt(site_in, bc_in, harvest_tag) ! @@ -1152,7 +1104,7 @@ subroutine get_harvest_debt(site_in, bc_in, harvest_tag) ! !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) + integer, intent(in) :: harvest_tag(hlm_num_lu_harvest_cats) ! !LOCAL VARIABLES: integer :: h_index @@ -1160,6 +1112,11 @@ subroutine get_harvest_debt(site_in, bc_in, harvest_tag) real(r8) :: harvest_debt_sec_mature real(r8) :: harvest_debt_sec_young + ! Flush the older value + site_in%resources_management%harvest_debt = 0._r8 + site_in%resources_management%harvest_debt_sec_m = 0._r8 + site_in%resources_management%harvest_debt_sec_y = 0._r8 + if(logging_time) then ! Initialize the local variables @@ -1182,19 +1139,19 @@ subroutine get_harvest_debt(site_in, bc_in, harvest_tag) 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 (harvest_tag(h_index) .ne. 0 .and. bc_in%hlm_harvest_units == hlm_harvest_carbon) 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 + & + site_in%resources_management%harvest_debt_sec_m = site_in%resources_management%harvest_debt_sec_m + & 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 + & + site_in%resources_management%harvest_debt_sec_y = site_in%resources_management%harvest_debt_sec_y + & harvest_debt_sec_young end if end if @@ -1203,4 +1160,57 @@ subroutine get_harvest_debt(site_in, bc_in, harvest_tag) end subroutine get_harvest_debt + ! ============================================================================ + + subroutine get_cohort_harvest_rate_scale (logging_preference_options, dbh, & + harvest_wp_scale, harvest_rate_scale_patch, harvest_rate_scale_cohort) + + ! !USES: + + ! !ARGUMENTS: + integer, intent(in) :: logging_preference_options + real(r8), intent(in) :: dbh + real(r8), intent(in) :: harvest_wp_scale + real(r8), intent(in) :: harvest_rate_scale_patch + real(r8), intent(out) :: harvest_rate_scale_cohort + + ! [Cohort level] Logging size prefrence or change rotation length + select case (logging_preference_options) + + case (logging_uniform_size) + + harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale_patch + + case (logging_double_rotation) + + harvest_rate_scale_cohort = harvest_wp_scale * 0.5_r8 * harvest_rate_scale_patch + + case (logging_quadruple_rotation) + + harvest_rate_scale_cohort = harvest_wp_scale * 0.25_r8 * harvest_rate_scale_patch + + case (logging_logistic_size) + + harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale_patch * & + 1._r8/(1._r8 + exp(-0.1*(dbh-(logging_dbhmin+logging_dbhmax)/2._r8))) + + case (logging_inv_logistic_size) + + harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale_patch * & + 1._r8/(1._r8 + exp(0.1*(dbh-(logging_dbhmin+logging_dbhmax)/2._r8))) + + case (logging_gaussian_size) + + harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale_patch * & + (1._r8/(5._r8*2.5_r8)) * exp(-(dbh - (logging_dbhmin+logging_dbhmax)/2._r8) ** & + 2._r8 / (2._r8 * 5._r8 ** 2._r8)) + + case DEFAULT + + end select + + end subroutine get_cohort_harvest_rate_scale + + ! ============================================================================ + end module EDLoggingMortalityMod diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 old mode 100644 new mode 100755 index fc941d0371..e7cf37fbfc --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -22,6 +22,7 @@ module EDMortalityFunctionsMod use FatesInterfaceTypesMod , only : hlm_use_tree_damage use EDLoggingMortalityMod , only : LoggingMortality_frac use EDParamsMod , only : fates_mortality_disturbance_fraction + use EDParamsMod , only : logging_preference_options use PRTGenericMod, only : carbon12_element use PRTGenericMod, only : store_organ @@ -235,7 +236,8 @@ end subroutine mortality_rates subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, & mean_temp, land_use_label, age_since_anthro_disturbance, & - frac_site_primary, harvestable_forest_c, harvest_tag) + frac_site_primary, frac_site_secondary_mature, harvestable_forest_c, & + harvest_rate_scale, harvest_tag) ! ! !DESCRIPTION: @@ -255,8 +257,10 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, & integer, intent(in) :: land_use_label real(r8), intent(in) :: age_since_anthro_disturbance real(r8), intent(in) :: frac_site_primary + real(r8), intent(in) :: frac_site_secondary_mature real(r8), intent(in) :: harvestable_forest_c(:) ! total carbon available for logging, kgC site-1 + real(r8), intent(in) :: harvest_rate_scale ! scale factor for patch-specific harvest rate integer, intent(out) :: harvest_tag(:) ! tag to record the harvest status ! for the calculation of harvest debt in C-based ! harvest mode @@ -293,7 +297,10 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, & bc_in%hlm_harvest_units, & land_use_label, & age_since_anthro_disturbance, & - frac_site_primary, harvestable_forest_c, harvest_tag) + frac_site_primary, frac_site_secondary_mature, & + harvestable_forest_c, & + harvest_rate_scale, harvest_tag, & + currentSite%resources_management%harvest_wp_scale, logging_preference_options) 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 old mode 100644 new mode 100755 index 140c108d66..45f3308968 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -69,7 +69,6 @@ module EDPatchDynamicsMod 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 @@ -80,11 +79,16 @@ module EDPatchDynamicsMod use FatesConstantsMod , only : years_per_day use FatesConstantsMod , only : nearzero use FatesConstantsMod , only : primaryland, secondaryland, pastureland, rangeland, cropland + use FatesConstantsMod , only : secondary_age_threshold use FatesConstantsMod , only : n_landuse_cats + use FatesConstantsMod , only : min_harvest_rate, min_nocomp_pftfrac_perlanduse use FatesLandUseChangeMod, only : get_landuse_transition_rates use FatesConstantsMod , only : fates_unset_r8 use FatesConstantsMod , only : fates_unset_int use FatesConstantsMod , only : hlm_harvest_carbon + use FatesConstantsMod , only : logging_no_age_preference, logging_oldest_first + use FatesConstantsMod , only : logging_uniform_size, logging_double_rotation, logging_quadruple_rotation + use FatesConstantsMod , only : logging_logistic_size, logging_inv_logistic_size, logging_gaussian_size use EDCohortDynamicsMod , only : InitPRTObject use ChecksBalancesMod, only : SiteMassStock use PRTGenericMod, only : carbon12_element @@ -101,6 +105,8 @@ module EDPatchDynamicsMod use SFParamsMod, only : SF_VAL_CWD_FRAC use EDParamsMod, only : logging_event_code use EDParamsMod, only : logging_export_frac + use EDParamsMod, only : logging_age_preference, logging_preference_options + use EDParamsMod, only : logging_ifm_harvest_scale use EDParamsMod, only : maxpatches_by_landuse use FatesRunningMeanMod, only : ema_sdlng_mdd use FatesRunningMeanMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par, ema_sdlng2sap_par @@ -124,6 +130,7 @@ module EDPatchDynamicsMod public :: set_patchno private:: fuse_2_patches public :: get_frac_site_primary + public :: get_frac_site_secondary_mature character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -176,8 +183,10 @@ subroutine disturbance_rates( site_in, bc_in) use EDMortalityFunctionsMod , only : ExemptTreefallDist ! loging flux use EDLoggingMortalityMod , only : LoggingMortality_frac + ! IFM related + use FatesForestManagementMod, only : get_site_harvest_rate_scale, get_patch_harvest_rate_scale + use FatesForestManagementMod, only : get_harvestable_carbon - ! !ARGUMENTS: type(ed_site_type) , intent(inout) :: site_in type(bc_in_type) , intent(in) :: bc_in @@ -204,25 +213,42 @@ subroutine disturbance_rates( site_in, bc_in) integer :: threshold_sizeclass integer :: i_dist integer :: h_index + integer :: npatches + real(r8) :: frac_site real(r8) :: frac_site_primary + real(r8) :: frac_site_secondary_mature real(r8) :: harvest_rate real(r8) :: tempsum real(r8) :: mean_temp real(r8) :: harvestable_forest_c(hlm_num_lu_harvest_cats) + real(r8) :: harvestable_patch_c integer :: harvest_tag(hlm_num_lu_harvest_cats) + integer :: harvest_tag_csite(hlm_num_lu_harvest_cats) ! harvest tag of current site + real(r8) :: remain_harvest_rate(hlm_num_lu_harvest_cats) ! temporary variable holding the remining harvest demand real(r8) :: landuse_transition_matrix(n_landuse_cats, n_landuse_cats) ! [m2/m2/day] real(r8) :: current_fates_landuse_state_vector(n_landuse_cats) ! [m2/m2] + + real(r8) :: harvested_cohort_c ! used in IFM + !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) ! And the same rates in understory plants have already been applied to %dndt !---------------------------------------------------------------------------------------------- - - ! first calculate the fraction of the site that is primary land + + ! first calculate the fraction of the site that is primary land and secondary mature land call get_frac_site_primary(site_in, frac_site_primary) + call get_frac_site_secondary_mature(site_in, frac_site_secondary_mature) - ! get available biomass for harvest for all patches + ! get available biomass carbon for harvest for all patches call get_harvestable_carbon(site_in, bc_in%site_area, bc_in%hlm_harvest_catnames, harvestable_forest_c) + ! get patch level scaling factor for harvest rate (currentPatch%harvest_rate_scale) + call get_patch_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, frac_site_primary, frac_site_secondary_mature) + + ! get site level scaling factor for harvest rate (site_in%resources_management%harvest_wp_scale) + call get_site_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, frac_site_primary, frac_site_secondary_mature) + + !Loop through cohorts to get disturbance rates currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -254,21 +280,30 @@ subroutine disturbance_rates( site_in, bc_in) currentPatch%land_use_label, & currentPatch%age_since_anthro_disturbance, & frac_site_primary, & + frac_site_secondary_mature, & harvestable_forest_c, & - harvest_tag) - + currentPatch%harvest_rate_scale, & + harvest_tag, site_in%resources_management%harvest_wp_scale, logging_preference_options) + currentCohort%lmort_direct = lmort_direct currentCohort%lmort_collateral = lmort_collateral currentCohort%lmort_infra = lmort_infra currentCohort%l_degrad = l_degrad + if (logging_time) then + do h_index = 1,hlm_num_lu_harvest_cats + harvest_tag_csite(h_index) = min(harvest_tag_csite(h_index), harvest_tag(h_index)) + end do + end if + + currentCohort => currentCohort%taller end do currentPatch => currentPatch%younger end do - call get_harvest_debt(site_in, bc_in, harvest_tag) + call get_harvest_debt(site_in, bc_in, harvest_tag_csite) if ( hlm_use_luh .eq. itrue ) then call get_landuse_transition_rates(bc_in, landuse_transition_matrix) @@ -367,15 +402,18 @@ subroutine disturbance_rates( site_in, bc_in) harvest_rate, harvest_tag) else call get_harvest_rate_area (currentPatch%land_use_label, bc_in%hlm_harvest_catnames, & - bc_in%hlm_harvest_rates, frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate) + bc_in%hlm_harvest_rates, frac_site_primary, frac_site_secondary_mature, & + 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 + (currentPatch%area - currentPatch%total_canopy_area) * & + currentPatch%harvest_rate_scale * harvest_rate / currentPatch%area ! Non-harvested part of the logging disturbance rate dist_rate_ldist_notharvested = dist_rate_ldist_notharvested + & - (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area + (currentPatch%area - currentPatch%total_canopy_area) * & + currentPatch%harvest_rate_scale * harvest_rate / currentPatch%area endif ! For nocomp mode, we need to prevent producing too small patches, which may produce small patches @@ -569,7 +607,7 @@ subroutine spawn_patches( currentSite, bc_in) else disturbance_rate = currentPatch%landuse_transition_rates(i_landusechange_receiverpatchlabel) endif - + if(disturbance_rate > (1.0_r8 + rsnbl_math_prec)) then write(fates_log(),*) 'patch disturbance rate > 1 ?',disturbance_rate call currentPatch%Dump() @@ -797,6 +835,7 @@ subroutine spawn_patches( currentSite, bc_in) nc%lmort_collateral = nan nc%lmort_infra = nan nc%l_degrad = nan + nc%harv_c = nan else ! understory trees @@ -861,6 +900,7 @@ subroutine spawn_patches( currentSite, bc_in) nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral nc%lmort_infra = currentCohort%lmort_infra + nc%harv_c = currentCohort%harv_c ! understory trees that might potentially be knocked over in the disturbance. ! The existing (donor) patch should not have any impact mortality, it should @@ -885,9 +925,10 @@ subroutine spawn_patches( currentSite, bc_in) nc%asmort = currentCohort%asmort nc%dgmort = currentCohort%dgmort nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral nc%lmort_infra = currentCohort%lmort_infra + nc%harv_c = currentCohort%harv_c endif woody_if_falldtype endif in_canopy_if_falldtype @@ -957,6 +998,7 @@ subroutine spawn_patches( currentSite, bc_in) nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral nc%lmort_infra = currentCohort%lmort_infra + nc%harv_c = currentCohort%harv_c ! Some of of the leaf mass from living plants has been @@ -1055,6 +1097,7 @@ subroutine spawn_patches( currentSite, bc_in) nc%lmort_direct = 0._r8 nc%lmort_collateral = 0._r8 nc%lmort_infra = 0._r8 + nc%harv_c = 0._r8 else @@ -1121,6 +1164,7 @@ subroutine spawn_patches( currentSite, bc_in) nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral nc%lmort_infra = currentCohort%lmort_infra + nc%harv_c = currentCohort%harv_c else @@ -1144,6 +1188,7 @@ subroutine spawn_patches( currentSite, bc_in) nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral nc%lmort_infra = currentCohort%lmort_infra + nc%harv_c = currentCohort%harv_c endif woody_if_logdtype ! is/is-not woody @@ -2361,6 +2406,11 @@ subroutine landusechange_litter_fluxes(currentSite, currentPatch, & site_mass%wood_product = site_mass%wood_product + & woodproduct_mass + + site_mass%wood_product_sz(currentCohort%size_class) = & + site_mass%wood_product_sz(currentCohort%size_class) + woodproduct_mass + + currentCohort%harv_c = currentCohort%harv_c + woodproduct_mass / currentCohort%n endif new_litt%ag_cwd(c) = new_litt%ag_cwd(c) + donatable_mass * donate_m2 curr_litt%ag_cwd(c) = curr_litt%ag_cwd(c) + donatable_mass * retain_m2 @@ -2413,12 +2463,17 @@ subroutine fuse_patches( csite, bc_in ) integer :: i_pftlabel !nocomp pft iterator real(r8) :: primary_land_fraction_beforefusion,primary_land_fraction_afterfusion integer :: pftlabelmin, pftlabelmax + + integer, parameter :: merge_strategy = 1 ! 0 - biomass profile based only + ! 1 - age + biomass profile based ! + real(r8) :: agetol !tolerance of patch fusion routine. Starts off high and is reduced to 0.0 if there are too many patches. !--------------------------------------------------------------------- currentSite => csite profiletol = ED_val_patch_fusion_tol + agetol = 0.1_r8 primary_land_fraction_beforefusion = 0._r8 primary_land_fraction_afterfusion = 0._r8 @@ -2560,12 +2615,23 @@ subroutine fuse_patches( csite, bc_in ) !---------------------------------------------------------------------! ! Look for differences in profile biomass, above the minimum biomass ! !---------------------------------------------------------------------! - - if(norm > profiletol)then + if(currentPatch%land_use_label .eq. primaryland .and. & + norm > profiletol)then fuse_flag = 0 !do not fuse - keep apart. endif + + ! Conditions come later have higher priorities + if(merge_strategy == 1) then + if(currentPatch%land_use_label .eq. secondaryland .and. & + abs(currentPatch%age_since_anthro_disturbance - & + tpp%age_since_anthro_disturbance) >= agetol) then + + fuse_flag = 0 !do not fuse - keep apart. + + endif + endif endif agbprof_gt_zero_if enddo hgt_bin_loop enddo pft_loop @@ -2602,6 +2668,7 @@ subroutine fuse_patches( csite, bc_in ) !------------------------------------------------------------------------! profiletol = ED_val_patch_fusion_tol + agetol = 0.1_r8 endif fuseflagset_if endif different_patches_if @@ -2637,6 +2704,7 @@ subroutine fuse_patches( csite, bc_in ) if(nopatches(i_lulabel) > maxpatches_by_landuse(i_lulabel))then iterate = 1 profiletol = profiletol * patch_fusion_tolerance_relaxation_increment + agetol = agetol * patch_fusion_tolerance_relaxation_increment !(1._r8/patch_fusion_tolerance_relaxation_increment) !---------------------------------------------------------------------! ! Making profile tolerance larger means that more fusion will happen ! @@ -3027,6 +3095,7 @@ subroutine terminate_patches(currentSite) count_cycles = 0 else + write(fates_log(),*) 'currentPatch%area?', currentPatch%area, 'currentPatch%pft?', currentPatch%nocomp_pft_label count_cycles = count_cycles + 1 end if @@ -3038,7 +3107,7 @@ subroutine terminate_patches(currentSite) write(fates_log(),*) 'disabling the endrun statement following this message.' write(fates_log(),*) 'FATES may or may not continue to operate within error' write(fates_log(),*) 'tolerances, but will generate another fail if it does not.' - call endrun(msg=errMsg(sourcefile, __LINE__)) + !call endrun(msg=errMsg(sourcefile, __LINE__)) ! Note to user. If you DO decide to remove the end-run above this line ! Make sure that you keep the pointer below this line, or you will get @@ -3204,168 +3273,152 @@ end subroutine get_frac_site_primary ! ===================================================================================== - subroutine InsertPatch(currentSite, newPatch) + subroutine get_frac_site_secondary_mature(site_in, frac_site_secondary_mature) + ! ! !DESCRIPTION: - ! Insert patch into linked list + ! Calculate how much of a site is secondary mature land ! ! !USES: + use EDTypesMod , only : ed_site_type ! ! !ARGUMENTS: - type (ed_site_type), intent(inout) :: currentSite - type (fates_patch_type), intent(inout), pointer :: newPatch + type(ed_site_type) , intent(in), target :: site_in + real(r8) , intent(out) :: frac_site_secondary_mature ! !LOCAL VARIABLES: type (fates_patch_type), pointer :: currentPatch - integer :: insert_method ! Temporary dev - logical :: found_landuselabel_match - integer, parameter :: unordered_lul_groups= 1 - integer, parameter :: primaryland_oldest_group = 2 - integer, parameter :: numerical_order_lul_groups = 3 - integer, parameter :: age_order_only = 4 - - ! Insert new patch case options: - ! Option 1: Group the landuse types together, but the group order doesn't matter - ! Option 2: Option 1, but primarylands are forced to be the oldest group - ! Option 3: Option 1, but groups are in numerical order according to land use label index integer - ! (i.e. primarylands=1, secondarylands=2, ..., croplands=5) - ! Option 4: Don't group the patches by land use label. Simply add new patches to the youngest end. - - ! Hardcode the default insertion method. The options developed during FATES V1 land use are - ! currently being held for potential future usage. - insert_method = primaryland_oldest_group - - ! Start from the youngest patch and work to oldest, regarless of insertion_method - currentPatch => currentSite%youngest_patch - ! For the three grouped cases, if the land use label of the youngest patch on the site - ! is a match to the new patch land use label, simply insert it as the new youngest. - ! This is applicable to the non-grouped option 4 method as well. - if (currentPatch%land_use_label .eq. newPatch%land_use_label ) then - newPatch%older => currentPatch - newPatch%younger => null() - currentPatch%younger => newPatch - currentSite%youngest_patch => newPatch - else - - ! If the current site youngest patch land use label doesn't match the new patch - ! land use label then work through the list until you find the matching type. - ! Since we've just checked the youngest patch, move to the next patch and - ! initialize the match flag to false. - found_landuselabel_match = .false. - currentPatch => currentPatch%older - select case(insert_method) + frac_site_secondary_mature = 0._r8 + currentPatch => site_in%oldest_patch + do while (associated(currentPatch)) + ! Here we need to be careful since age_since_anthro_disturbance is not + ! initialized for primary land, so is more robust to be a second level + ! if condition instead of using and statement. + if (currentPatch%land_use_label .eq. secondaryland) then + if(currentPatch%age_since_anthro_disturbance .ge. secondary_age_threshold) then + frac_site_secondary_mature = frac_site_secondary_mature + currentPatch%area * AREA_INV + endif + endif + currentPatch => currentPatch%younger + end do - ! Option 1 - order of land use label groups does not matter - case (unordered_lul_groups) + end subroutine get_frac_site_secondary_mature - do while(associated(currentPatch) .and. .not. found_landuselabel_match) - if (currentPatch%land_use_label .eq. newPatch%land_use_label) then - found_landuselabel_match = .true. - else - currentPatch => currentPatch%older - end if - end do + ! ===================================================================================== - ! In the case where we've found a land use label matching the new patch label, - ! insert the newPatch will as the youngest patch for that land use type. - if (associated(currentPatch)) then - newPatch%older => currentPatch - newPatch%younger => currentPatch%younger - currentPatch%younger%older => newPatch - currentPatch%younger => newPatch - else - ! In the case in which we get to the end of the list and haven't found - ! a landuse label match simply add the new patch to the youngest end. - newPatch%older => currentSite%youngest_patch - newPatch%younger => null() - currentSite%youngest_patch%younger => newPatch - currentSite%youngest_patch => newPatch - endif + subroutine InsertPatch(currentSite, newPatch) - ! Option 2 - primaryland group must be on the oldest end - case (primaryland_oldest_group) + ! !DESCRIPTION: + ! Insert patch into linked list + ! + ! !USES: + ! + ! !ARGUMENTS: + type (ed_site_type), intent(inout) :: currentSite + type (fates_patch_type), intent(inout), pointer :: newPatch - do while(associated(currentPatch) .and. .not. found_landuselabel_match) - if (currentPatch%land_use_label .eq. newPatch%land_use_label) then - found_landuselabel_match = .true. - else - currentPatch => currentPatch%older - end if - end do + ! !LOCAL VARIABLES: + type (fates_patch_type), pointer :: currentPatch + logical :: patch_inserted - ! In the case where we've found a land use label matching the new patch label, - ! insert the newPatch will as the youngest patch for that land use type. - if (associated(currentPatch)) then - newPatch%older => currentPatch - newPatch%younger => currentPatch%younger - currentPatch%younger%older => newPatch - currentPatch%younger => newPatch - else - ! In the case in which we get to the end of the list and haven't found - ! a landuse label match. - - ! If the new patch is primarylands add it to the oldest end of the list - if (newPatch%land_use_label .eq. primaryland) then - newPatch%older => null() - newPatch%younger => currentSite%oldest_patch - currentSite%oldest_patch%older => newPatch - currentSite%oldest_patch => newPatch - else - ! If the new patch land use type is not primaryland and we are at the - ! oldest end of the list, add it to the youngest end - newPatch%older => currentSite%youngest_patch - newPatch%younger => null() - currentSite%youngest_patch%younger => newPatch - currentSite%youngest_patch => newPatch - endif - endif + ! The goal here is to have patches ordered in a specific way. That way is to have them + ! looped as the following, where LU refers to the land use label, PFT refers to the + ! nocomp PFT label, and Y and O refer to continuous patch ages. + ! + ! LU1 ---- LU2 ---- LU3 -- etc + ! / \ / \ / \ + ! PFT1 --- PFT2 | PFT1 --- PFT2 | PFT1 --- PFT2 -- etc + ! / \ / \ / \ / \ / \ / \ + ! O - Y O - Y O - Y O - Y O - Y O - Y -- etc + + ! I.e. to treat land use as the outermost loop element, then nocomp PFT as next loop element, + ! and then age as the innermost loop element. Visualizing the above as a linked list patches: + + ! LU1/PFT1/O <-> LU1/PFT1/Y <-> LU1/PFT2/O <- ... -> LU3/PFT2/O <-> LU3/PFT2/Y + + ! Mapping this setup onto the existing "older/younger" scheme means that lower number + ! land use and pft labels are considered "older". Note that this matches the current + ! initialization scheme in which patches are created and linked in increasing pft + ! numerical order starting from 1. This also aligns with the current set_patchno scheme + ! in which patches are given an indexable number for the API iteration loops. + + ! The way to accomplsh this most simply is to define a pseudo-age that includes all of the + ! above info and sort the patches based on the pseudo-age. i.e. take some number larger + ! than any patch will ever reach in actual age. Then take the LU, multiply it by the big + ! number squared, add it to the pft number multiplied by the big number, and add to the age. + ! And lastly to sort using that instead of the actual age. - ! Option 3 - groups are numerically ordered with primaryland group starting at oldest end. - case (numerical_order_lul_groups) + ! If land use is turned off or nocomp is turned off, then this should devolve to the prior + ! behavior of just age sorting. - ! If the youngest patch landuse label number is greater than the new - ! patch land use label number, the new patch must be inserted somewhere - ! in between oldest and youngest - do while(associated(currentPatch) .and. .not. found_landuselabel_match) - if (currentPatch%land_use_label .eq. newPatch%land_use_label .or. & - currentPatch%land_use_label .lt. newPatch%land_use_label) then - found_landuselabel_match = .true. - else - currentPatch => currentPatch%older - endif - end do + patch_inserted = .false. + + if (GetPseudoPatchAge(newPatch) .le. GetPseudoPatchAge(currentSite%youngest_patch)) then - ! In the case where we've found a landuse label matching the new patch label - ! insert the newPatch will as the youngest patch for that land use type. - if (associated(currentPatch)) then + ! insert new patch at the head of the linked list + newPatch%older => currentSite%youngest_patch + newPatch%younger => null() + currentSite%youngest_patch%younger => newPatch + currentSite%youngest_patch => newPatch - newPatch%older => currentPatch - newPatch%younger => currentPatch%younger - currentPatch%younger%older => newPatch - currentPatch%younger => newPatch + patch_inserted = .true. + else if (GetPseudoPatchAge(newPatch) .ge. GetPseudoPatchAge(currentSite%oldest_patch)) then - else + ! insert new patch at the end of the linked list + newPatch%younger => currentSite%oldest_patch + newPatch%older => null() + currentSite%oldest_patch%older => newPatch + currentSite%oldest_patch => newPatch - ! In the case were we get to the end, the new patch - ! must be numerically the smallest, so put it at the oldest position - newPatch%older => null() - newPatch%younger => currentSite%oldest_patch - currentSite%oldest_patch%older => newPatch - currentSite%oldest_patch => newPatch + patch_inserted = .true. + else + ! new patch has a pseudo-age somewhere within the linked list. find the first patch which + ! has a pseudo age older than it, and put it ahead of that patch + currentPatch => currentSite%youngest_patch + do while (associated(currentPatch) .and. ( .not. patch_inserted) ) + if (GetPseudoPatchAge(newPatch) .lt. GetPseudoPatchAge(currentPatch)) then + newPatch%older => currentPatch + newPatch%younger => currentPatch%younger + currentPatch%younger%older => newPatch + currentPatch%younger => newPatch + patch_inserted = .true. endif + currentPatch => currentPatch%older + end do + end if - ! Option 4 - always add the new patch as the youngest regardless of land use label - case (age_order_only) - ! Set the current patch to the youngest patch - newPatch%older => currentSite%youngest_patch - newPatch%younger => null() - currentSite%youngest_patch%younger => newPatch - currentSite%youngest_patch => newPatch - end select + if ( .not. patch_inserted) then + ! something has gone wrong. abort. + write(fates_log(),*) 'something has gone wrong in the patch insertion, because no place to put the new patch was found' + call endrun(msg=errMsg(sourcefile, __LINE__)) end if end subroutine InsertPatch + ! ===================================================================================== + + function GetPseudoPatchAge(CurrentPatch) result(pseudo_age) + + ! Purpose: we want to sort the patches in a way that takes into account both their + ! continuous and categorical variables. Calculate a pseudo age that does this, by taking + ! the integer labels, multiplying these by large numbers, and adding to the continuous age. + ! Note to ensure that lower integer land use label and pft label numbers are considered + ! "younger" (i.e higher index patchno) in the linked list, they are summed and multiplied by + ! negative one. The patch age is still added normally to this negative pseudoage calculation + ! as a higher age will result in a less negative number correlating with an "older" patch. + + type (fates_patch_type), intent(in), pointer :: CurrentPatch + real(r8) :: pseudo_age + real(r8), parameter :: max_actual_age = 1.e4 ! hard to imagine a patch older than 10,000 years + real(r8), parameter :: max_actual_age_squared = 1.e8 + + pseudo_age = -1.0_r8 * (real(CurrentPatch%land_use_label,r8) * max_actual_age_squared + & + real(CurrentPatch%nocomp_pft_label,r8) * max_actual_age) + CurrentPatch%age + + end function GetPseudoPatchAge + + ! ===================================================================================== + end module EDPatchDynamicsMod diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 727aaa7370..1c66d9dd0b 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -488,6 +488,10 @@ subroutine PreDisturbanceLitterFluxes( currentSite, currentPatch, bc_in ) site_mass => currentSite%mass_balance(el) + site_mass%frag_in = site_mass%frag_in + currentPatch%area * & + ( sum(litt%ag_cwd_in) + sum(litt%bg_cwd_in) + & + sum(litt%leaf_fines_in) + sum(litt%root_fines_in)) + ! Fragmentation flux to soil decomposition model [kg/site/day] site_mass%frag_out = site_mass%frag_out + currentPatch%area * & ( sum(litt%ag_cwd_frag) + sum(litt%bg_cwd_frag) + & @@ -2877,7 +2881,7 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in) flux_diags%leaf_litter_input(pft) = & flux_diags%leaf_litter_input(pft) + & - leaf_m_turnover * currentCohort%n + (leaf_m_turnover+repro_m_turnover) * currentCohort%n root_fines_tot = (fnrt_m_turnover + store_m_turnover ) * & plant_dens @@ -2898,7 +2902,6 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in) flux_diags%root_litter_input(pft) + & (fnrt_m_turnover + store_m_turnover ) * currentCohort%n - ! Assumption: turnover from deadwood and sapwood are lumped together in CWD pool !update partitioning of stem wood (struct + sapw) to cwd based on cohort dbh @@ -2926,10 +2929,8 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in) flux_diags%cwd_bg_input(c) = flux_diags%cwd_bg_input(c) + & bg_cwd_tot*currentPatch%area - enddo - ! --------------------------------------------------------------------------------- ! PART 2 Litter fluxes from non-disturbance inducing mortality. Kg/m2/day ! --------------------------------------------------------------------------------- @@ -2963,8 +2964,7 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in) flux_diags%leaf_litter_input(pft) = & flux_diags%leaf_litter_input(pft) + & - leaf_m * dead_n*currentPatch%area - + (leaf_m+repro_m) * dead_n*currentPatch%area ! %n has not been updated due to mortality yet, thus ! the litter flux has already been counted since it captured @@ -3022,6 +3022,13 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in) site_mass%wood_product = site_mass%wood_product + & trunk_wood * currentPatch%area * logging_export_frac + site_mass%wood_product_sz(currentCohort%size_class) = & + site_mass%wood_product_sz(currentCohort%size_class) + & + trunk_wood * currentPatch%area * logging_export_frac + + currentCohort%harv_c = currentCohort%harv_c + & + (trunk_wood * currentPatch%area * logging_export_frac) / currentCohort%n + ! Add AG wood to litter from the non-exported fraction of wood ! from direct anthro sources @@ -3037,7 +3044,7 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in) SF_val_CWD_frac_adj(c) * (dead_n_natural+dead_n_ilogging) * & prt_params%allom_agb_frac(pft) - flux_diags%cwd_ag_input(c) = flux_diags%cwd_ag_input(c) + & + flux_diags%cwd_ag_input(c) = flux_diags%cwd_ag_input(c) + (struct_m + sapw_m) * & SF_val_CWD_frac_adj(c) * (dead_n_natural+dead_n_ilogging) * & currentPatch%area * prt_params%allom_agb_frac(pft) @@ -3050,12 +3057,10 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in) flux_diags%cwd_ag_input(c) = flux_diags%cwd_ag_input(c) + & SF_val_CWD_frac_adj(c) * dead_n * (struct_m + sapw_m) * & currentPatch%area * prt_params%allom_agb_frac(pft) - end if end do - ! Update diagnostics that track resource management if( element_id .eq. carbon12_element ) then @@ -3097,7 +3102,6 @@ subroutine CWDInput( currentSite, currentPatch, litt, bc_in) currentCohort => currentCohort%taller enddo ! end loop over cohorts - return end subroutine CWDInput diff --git a/biogeochem/FatesCohortMod.F90 b/biogeochem/FatesCohortMod.F90 index e98c4a345a..c709c69316 100644 --- a/biogeochem/FatesCohortMod.F90 +++ b/biogeochem/FatesCohortMod.F90 @@ -239,6 +239,7 @@ module FatesCohortMod real(r8) :: l_degrad ! rate of trees that are not killed but suffer from forest degradation ! (i.e. they are moved to newly-anthro-disturbed secondary ! forest patch) [fraction/logging activity] + real(r8) :: harv_c ! harvested carbon [kgC/indiv] !--------------------------------------------------------------------------- @@ -434,6 +435,7 @@ subroutine NanValues(this) this%lmort_collateral = nan this%lmort_infra = nan this%l_degrad = nan + this%harv_c = nan ! GROWTH DERIVATIVES this%dndt = nan @@ -530,6 +532,7 @@ subroutine ZeroValues(this) this%lmort_collateral = 0._r8 this%lmort_infra = 0._r8 this%l_degrad = 0._r8 + this%harv_c = 0._r8 this%fraction_crown_burned = 0._r8 this%cambial_mort = 0._r8 this%crownfire_mort = 0._r8 @@ -769,6 +772,7 @@ subroutine Copy(this, copyCohort) copyCohort%lmort_collateral = this%lmort_collateral copyCohort%lmort_infra = this%lmort_infra copyCohort%l_degrad = this%l_degrad + copyCohort%harv_c = this%harv_c ! GROWTH DERIVATIVES copyCohort%dndt = this%dndt @@ -1074,6 +1078,7 @@ subroutine Dump(this) write(fates_log(),*) 'cohort%lmort_direct = ', this%lmort_direct write(fates_log(),*) 'cohort%lmort_collateral = ', this%lmort_collateral write(fates_log(),*) 'cohort%lmort_infra = ', this%lmort_infra + write(fates_log(),*) 'cohort%harv_c = ', this%harv_c write(fates_log(),*) 'cohort%isnew = ', this%isnew write(fates_log(),*) 'cohort%dndt = ', this%dndt write(fates_log(),*) 'cohort%dhdt = ', this%dhdt diff --git a/biogeochem/FatesForestManagementMod.F90 b/biogeochem/FatesForestManagementMod.F90 new file mode 100644 index 0000000000..885bd98a03 --- /dev/null +++ b/biogeochem/FatesForestManagementMod.F90 @@ -0,0 +1,684 @@ +module FatesForestManagementMod + + ! ============================================================================ + ! Created by Shijie Shu on Sep 8. 2025 + ! This module contains all subroutines and functions related to planned + ! wood harvest and improved forest management. + !----------------------------------------------------------------------------------- + ! Improved forest management practies can be quantified as modifiers to change the + ! logging disturbance rates at cohort-patch-site levels, thus quantify strategies as + ! different priorties regarding size/age for harvest activity. + ! Site: site_in%resources_management%harvest_wp_scale. Adjust the overall harvest rate + ! to match the demand. + ! Patch: currentPatch%harvest_rate_scale. Modifier considering age priority. + ! Cohort [LOCAL]: harvest_rate_scale_cohort. Modifier considering size priority. + ! ============================================================================ + + ! Uses + use EDParamsMod , only : logging_age_preference, logging_preference_options + use EDParamsMod , only : logging_ifm_harvest_scale + use EDTypesMod , only : AREA, AREA_INV + use EDTypesMod , only : ed_site_type + use EDLoggingMortalityMod , only : get_harvest_rate_area, get_harvest_rate_carbon + use FatesPatchMod , only : fates_patch_type + use FatesCohortMod , only : fates_cohort_type + use FatesGlobals , only : fates_log + use FatesConstantsMod , only : r8 => fates_r8 + use FatesConstantsMod , only : logging_no_age_preference, logging_oldest_first + use FatesConstantsMod , only : logging_uniform_size, logging_double_rotation, logging_quadruple_rotation + use FatesConstantsMod , only : logging_logistic_size, logging_inv_logistic_size, logging_gaussian_size + use FatesConstantsMod , only : fates_unset_r8, fates_tiny + use FatesConstantsMod , only : primaryland, secondaryland + use FatesConstantsMod , only : min_harvest_rate, min_nocomp_pftfrac_perlanduse + use FatesConstantsMod , only : hlm_harvest_carbon + use FatesLitterMod , only : ncwd + use FatesInterfaceTypesMod , only : bc_in_type + use FatesInterfaceTypesMod , only : hlm_num_lu_harvest_cats, numpft + use PRTParametersMod , only : prt_params + use PRTGenericMod , only : carbon12_element + use PRTGenericMod , only : leaf_organ + use PRTGenericMod , only : fnrt_organ + use PRTGenericMod , only : sapw_organ + use PRTGenericMod , only : store_organ + use PRTGenericMod , only : repro_organ + use PRTGenericMod , only : struct_organ + + ! Subroutines + public :: get_site_harvest_rate_scale + public :: get_patch_harvest_rate_scale + public :: get_harvestable_carbon + private :: get_harvestable_patch_carbon + private :: get_harvested_cohort_carbon + + contains + + ! ============================================================================ + + subroutine get_site_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, frac_site_primary, frac_site_secondary_mature) + + ! !USES: + use EDLoggingMortalityMod , only : logging_time, LoggingMortality_frac + + ! !ARGUMENTS: + type(ed_site_type) , intent(inout) :: site_in + type(bc_in_type) , intent(in) :: bc_in + real(r8), intent(in) :: harvestable_forest_c(hlm_num_lu_harvest_cats) + real(r8), intent(in) :: frac_site_primary + real(r8), intent(in) :: frac_site_secondary_mature + ! + ! !LOCAL VARIABLES: + type (fates_patch_type) , pointer :: currentPatch + type (fates_cohort_type), pointer :: currentCohort + + real(r8) :: lmort_direct + real(r8) :: lmort_collateral + real(r8) :: lmort_infra + real(r8) :: l_degrad + real(r8) :: harvested_cohort_c ! Variable to store cohort level harvested C + integer :: harvest_tag(hlm_num_lu_harvest_cats) + + real(r8) :: target_wp + real(r8) :: actual_wp + integer :: it + + ! !PARAMETERS + integer, parameter :: max_iteration = 5 ! Maximum iterations to obtain target harvest wood product. Currently I don't + ! plan to move it into parameter files + + ! ================================================================================================== + ! DESCRIPTION: + ! + ! For size-dependet harvest, cetain size class might have fewer or no harvest if the available + ! inventory is less than the demand. + ! We design the following algorithm to force the other nearby size class to fullfill these demands: + ! Loop through cohorts to preview the summed harvested product C for both no harvest size priority, + ! i.e., target wood product C + ! (target_wp), and with size priority, i.e., actual harvested C (actual_wp) + ! Then calculate the ratio of target_wp to actual_wp as adjustment ratio + ! The adjustment ratio will be applied as additional scaling factor to modify cohort-level harvest mortality + ! One time adjustment might still not able to match the demand, we iterate this process no more than 5 times. + ! ================================================================================================== + + !Initialize + site_in%resources_management%harvest_wp_scale = 1._r8 + write(fates_log(),*) 'See flag2' + + !Only run this part when considering size dependent priority + if (logging_time .and. logging_preference_options >= logging_logistic_size) then + + iteration_loop: do it = 1, max_iteration + + target_wp = 0._r8 + actual_wp = 0._r8 + + currentPatch => site_in%oldest_patch + do while (associated(currentPatch)) + + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + + ! Target harvest product + call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_layer, & + lmort_direct,lmort_collateral,lmort_infra,l_degrad,& + bc_in%hlm_harvest_rates, & + bc_in%hlm_harvest_catnames, & + bc_in%hlm_harvest_units, & + currentPatch%land_use_label, & + currentPatch%age_since_anthro_disturbance, & + frac_site_primary, & + frac_site_secondary_mature, & + harvestable_forest_c, & + currentPatch%harvest_rate_scale, & + harvest_tag, logging_ifm_harvest_scale, logging_uniform_size) + + call get_harvested_cohort_carbon(currentCohort, lmort_direct, harvested_cohort_c) + + target_wp = target_wp + harvested_cohort_c + + ! Actual harvest product + call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_layer, & + lmort_direct,lmort_collateral,lmort_infra,l_degrad,& + bc_in%hlm_harvest_rates, & + bc_in%hlm_harvest_catnames, & + bc_in%hlm_harvest_units, & + currentPatch%land_use_label, & + currentPatch%age_since_anthro_disturbance, & + frac_site_primary, & + frac_site_secondary_mature, & + harvestable_forest_c, & + currentPatch%harvest_rate_scale, & + harvest_tag, site_in%resources_management%harvest_wp_scale, logging_preference_options) + + call get_harvested_cohort_carbon(currentCohort, lmort_direct, harvested_cohort_c) + + actual_wp = actual_wp + harvested_cohort_c + + currentCohort => currentCohort%taller + + end do + + currentPatch => currentPatch%younger + + end do + + ! adjustment ratio + if( target_wp / (actual_wp + 1e-7_r8) < 1.01 .and. target_wp / (actual_wp + 1e-7_r8) > 0.99 ) exit iteration_loop + + site_in%resources_management%harvest_wp_scale = site_in%resources_management%harvest_wp_scale * target_wp/(actual_wp+1e-7_r8) + + write(fates_log(),*) 'See adjustment factor for harvest scale:', site_in%resources_management%harvest_wp_scale + + end do iteration_loop + + end if ! logging_time + + end subroutine get_site_harvest_rate_scale + + ! ============================================================================ + + subroutine get_patch_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, frac_site_primary, frac_site_secondary_mature) + + ! !USES: + use EDLoggingMortalityMod , only : logging_time + + ! !ARGUMENTS: + type(ed_site_type) , intent(inout) :: site_in + type(bc_in_type) , intent(in) :: bc_in + real(r8), intent(in) :: harvestable_forest_c(hlm_num_lu_harvest_cats) + real(r8), intent(in) :: frac_site_primary + real(r8), intent(in) :: frac_site_secondary_mature + + ! !LOCAL VARIABLES: + type (fates_patch_type) , pointer :: currentPatch + type (fates_cohort_type), pointer :: currentCohort + + integer :: npatches ! Number of patches + integer :: ipatch, curpft ! index, used in min-heap merge + integer, allocatable :: jpatch(:) ! index array, used in min-heap merge + integer, allocatable :: pftlen(:) ! array length, used in min-heap merge + real(r8), allocatable :: minheap(:) ! minimum heap array, used in min-heap merge + integer, allocatable :: order_of_patches(:) ! Record a copy of patch order + integer, allocatable :: order_of_patches_per_pft(:,:) ! Record a copy of patch order for each pft + real(r8), allocatable :: age_patches(:) ! Record a copy of patch age + real(r8), allocatable :: age_patches_per_pft(:,:) ! Record a copy of patch age for each pft + + real(r8) :: frac_site ! Temporary variable + real(r8) :: harvest_rate ! Temporary variable + real(r8) :: harvestable_patch_c + real(r8) :: remain_harvest_rate(hlm_num_lu_harvest_cats) ! Temporary variable + real(r8) :: temp_totc ! for debug test, can be removed in released version + integer :: h_index ! index for looping over harvest categories + integer :: harvest_tag(hlm_num_lu_harvest_cats) + + ! !PARAMETERS + + ! Initialize the rest of local variables + npatches = 0 + frac_site = 1._r8 + remain_harvest_rate = fates_unset_r8 ! set to negative value to indicate uninitialized + + ! ===================================================================== + ! DESCRIPTION: + ! + ! This modifier is for age_prioritized forest management + ! The following code will calculate patch level harvest_rate_scale + ! for each secondary patch + ! There're currently 2 options: + ! 1) uniformly harvest each patch (harvest_rate_scale = 1._r8) + ! 2) harvest the oldest patch first + + ! Initialize harvest_rate_scale + write(fates_log(),*) 'See flag1' + currentPatch => site_in%oldest_patch + + do while (associated(currentPatch)) + ! Initialize harvest rate scale to 1.0 first + currentPatch%harvest_rate_scale = 1._r8 + + currentPatch => currentPatch%younger + end do + + ! ===================================================================== + ! "Sort" patch based on patch age since anthropegenic disturbance + ! ===================================================================== + ! Patch sequence in the data list (patchno) follows pseudo-age defined + ! in function GetPseudoPatchAge to simplify LUC related calculation. + ! Here we don't modify the actual patch sequence but use another age + ! tag currentPatch%order_age_since_anthro to record patch age for + ! age-prioritized wood harvest strategy. + ! + ! In FATES, PFT-patch relation can be different under diffrent modes. + ! This can be separated into 2 cases: + ! Case 1: For most of FATES modes, we directly use pseudo-age since + ! there's no paired 1 patch - 1 PFT relation + ! Case 2: For the more general case in global simulation that harvest + ! all PFTs equally under no competition mode, we need to sort again + ! using, actual patch age. For each PFT the patch sequence is already + ! sorted, here we need to loop over PFTs and merge these sorted + ! sequences into one sorted array through min-heap merge + write(fates_log(),*) 'logging_preference_options:', logging_preference_options + write(fates_log(),*) 'logging_age_preference:', logging_age_preference + write(fates_log(),*) 'logging_time:', logging_time + write(fates_log(),*) 'logging_oldest_first:', logging_oldest_first + write(fates_log(),*) 'hlm_use_nocomp:', hlm_use_nocomp + + if (logging_time .and. logging_age_preference == logging_oldest_first) then + write(fates_log(),*) 'See flag3' + + ! First loop to count total number of patches and initialize order + ! using patchno + ! For most of the cases the process is over after this step + currentPatch => site_in%oldest_patch + npatches = 0 + + do while (associated(currentPatch)) + npatches = npatches + 1 + currentPatch%order_age_since_anthro = currentPatch%patchno + + currentPatch => currentPatch%younger + end do + + ! For case 2, i.e., nocomp mode only + if (hlm_use_nocomp.eq.itrue) then + write(fates_log(),*) 'See flag4' + + allocate(order_of_patches(1:npatches)) + allocate(order_of_patches_per_pft(1:numpft,1:npatches)) + allocate(age_patches(1:npatches)) + allocate(age_patches_per_pft(1:numpft,1:npatches)) + allocate(minheap(1:numpft)) + allocate(jpatch(1:numpft)) + allocate(pftlen(1:numpft)) + + ! Second loop to assign patch sequence, age and number of pacthes in each PFT to array + ! thus can avoid for loop over data list + currentPatch => site_in%oldest_patch + ipatch = 0 + pftlen(:) = 0 + do while (associated(currentPatch)) + ipatch = ipatch + 1 + order_of_patches(ipatch) = currentPatch%order_age_since_anthro + age_patches(ipatch) = currentPatch%age + + do curpft=1,numpft + if(currentPatch%nocomp_pft_label .eq. curpft) then + pftlen(curpft) = pftlen(curpft) + 1 + order_of_patches_per_pft(curpft, pftlen(curpft)) = currentPatch%order_age_since_anthro + age_patches_per_pft(curpft, pftlen(curpft)) = currentPatch%age + end if + end do + + currentPatch => currentPatch%younger + end do + + ! Third loop to initailize min-heap, first element of each array + minheap(:) = 0._r8 + do curpft=1,numpft + if(pftlen(curpft) > 0) then + minheap(curpft) = age_patches_per_pft(curpft, 1) + end if + end do + + ! Fourth loop to get min from min-heap then insert 'age_patches_per_pft' and record order in 'order_of_patches' + ! no need to loop over data list here + jpatch(:) = 1 + do ipatch = 1, npacthes + ! Find the minimum and assign the order + curpft = minloc(minheap, dim=1) + order_of_patches(ipatch) = order_of_patches_per_pft(curpft, jpatch(curpft)) + ! Move to the next element and insert into minheap array + jpatch(curpft) = jpatch(curpft) + 1 + if(jpatch(curpft) <= pftlen(curpft)) then + minheap(curpft) = age_patches_per_pft(ipft, jpatch(curpft)) + else + minheap(curpft) = 0._r8 + end if + end do + + ! Fifthth loop to assign order back to patches + do ipatch = 1, npatches + currentPatch => site_in%oldest_patch + + ! Look for the patch by usig patchno + inner_loop: do while (associated(currentPatch)) + + if(order_of_patches(ipatch) .eq. currentPatch%patchno) then + currentPatch%order_age_since_anthro = ipatch + exit inner_loop + end if + + currentPatch => currentPatch%younger + end do inner_loop + end do + + ! Clean temperary variables + deallocate(order_of_patches) + deallocate(order_of_patches_per_pft) + deallocate(age_patches) + deallocate(age_patches_per_pft) + deallocate(minheap) + deallocate(jpatch) + deallocate(pftlen) + + end if ! nocomp + + ! Loop through all patches and calculate the harvest rate scale based + ! on oldest first strategy + target_loop: do ipatch = 1, npatches + + currentPatch => site_in%oldest_patch + search_loop: do while (associated(currentPatch)) + + found_patch: if(currentPatch%order_age_since_anthro .eq. ipatch) then + + ! Here age priority is based on age since anthropogenic disturbance + ! thus we skip the primary land since they all have this age = 0 + if_secondary: if(currentPatch%land_use_label .eq. secondaryland) then + if(currentPatch%age_since_anthro_disturbance >= secondary_age_threshold) then + h_index = 3 + frac_site = frac_site_secondary_mature + else + h_index = 4 + frac_site = 1._r8 - frac_site_primary - frac_site_secondary_mature + end if + ! Obtain uniform harvest rate (area fraction) of the corresponding patch + if(bc_in%hlm_harvest_units == hlm_harvest_carbon) then + call get_harvest_rate_carbon (currentPatch%land_use_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%land_use_label, bc_in%hlm_harvest_catnames, & + bc_in%hlm_harvest_rates, frac_site_primary, frac_site_secondary_mature, & + currentPatch%age_since_anthro_disturbance, harvest_rate) + end if + ! For the first time, initialize remain_harvest_rate + if(remain_harvest_rate(h_index) < 0) then + if(bc_in%hlm_harvest_units == hlm_harvest_carbon) then + ! In kgC ha-1, no need to scale + remain_harvest_rate(h_index) = bc_in%hlm_harvest_rates(h_index) + else + ! In area fraction (0 - 1) after scaling to the area of per patch land use type + remain_harvest_rate(h_index) = harvest_rate + end if + end if + ! Only calculate harvest rate scale for non-zero harvest rate + if (harvest_rate > min_harvest_rate) then + if(bc_in%hlm_harvest_units == hlm_harvest_carbon) then + ! Compare the patch level harvestable carbon and see if larger than remain_harvest_rate + ! Note the scaling factor applied on area fraction as harvest unit + call get_harvestable_patch_carbon(currentPatch, harvestable_patch_c) + if(harvestable_patch_c >= remain_harvest_rate(h_index)) then + currentPatch%harvest_rate_scale = remain_harvest_rate(h_index) / & + (harvestable_patch_c + fates_tiny) / harvest_rate + remain_harvest_rate(h_index) = 0._r8 + else + ! harvest almost 100%, leave 0.01% to prevent removing the patch completely, which may cause + ! model crashing under nocomp mode + currentPatch%harvest_rate_scale = (1._r8 - min_nocomp_pftfrac_perlanduse) / harvest_rate + remain_harvest_rate(h_index) = remain_harvest_rate(h_index) - (1._r8 - min_nocomp_pftfrac_perlanduse) * harvestable_patch_c + end if + else + ! Compare the patch area and see if larger than remain_harvest_rate + if((currentPatch%area/(AREA*frac_site)) >= remain_harvest_rate(h_index)) then + currentPatch%harvest_rate_scale = remain_harvest_rate(h_index) / & + (currentPatch%area/(AREA*frac_site) + fates_tiny) / harvest_rate + remain_harvest_rate(h_index) = 0._r8 + else + ! leave 1% to prevent removing the patch completely + currentPatch%harvest_rate_scale = (1._r8 - min_nocomp_pftfrac_perlanduse) / harvest_rate + remain_harvest_rate(h_index) = remain_harvest_rate(h_index) - (1._r8 - min_nocomp_pftfrac_perlanduse) * (currentPatch%area/(AREA*frac_site)) + end if + end if + end if + + ! For diagnosing the harvestable carbon + ! These codes can be removed in released version + temp_totc = 0._r8 + currentCohort => currentPatch%tallest + + do while (associated(currentCohort)) + + if (currentCohort%canopy_layer>=1) then + temp_totc = temp_totc + (currentCohort%prt%GetState(sapw_organ, carbon12_element) + & + currentCohort%prt%GetState(struct_organ, carbon12_element)) * currentCohort%n * AREA_INV + end if + currentCohort => currentCohort%shorter + + end do + + write(fates_log(),*) 'See patch number:', currentPatch%patchno + write(fates_log(),*) 'See patch age:', currentPatch%age + write(fates_log(),*) 'See patch total carbon (kgc m-2):', temp_totc + write(fates_log(),*) 'See patch age sec:', currentPatch%age_since_anthro_disturbance + write(fates_log(),*) 'See patch order sec:', currentPatch%order_age_since_anthro + write(fates_log(),*) 'See harvest index:', h_index + write(fates_log(),*) 'See harvest rate scale:', currentPatch%harvest_rate_scale + write(fates_log(),*) 'See original harvest rate:', bc_in%hlm_harvest_rates(h_index) + write(fates_log(),*) 'See harvest rate:', harvest_rate + write(fates_log(),*) 'See harvestable_forest_c:', harvestable_forest_c + write(fates_log(),*) 'See site fraction:', frac_site + write(fates_log(),*) 'See sec mature fraction:', frac_site_secondary_mature + write(fates_log(),*) 'See sec young fraction:', 1 - frac_site_primary - frac_site_secondary_mature + write(fates_log(),*) 'See area:', currentPatch%area/AREA + write(fates_log(),*) '==================== Next patch ================' + + end if if_secondary + + ! Exit current inner loop to find the next old ipatch + exit search_loop + + end if found_patch + + currentPatch => currentPatch%younger + end do search_loop + + end do target_loop + + end if ! logging_time + + end subroutine get_patch_harvest_rate_scale + + ! ============================================================================ + + subroutine get_harvestable_carbon (csite, site_area, hlm_harvest_catnames, harvestable_forest_c ) + + !USES: + use SFParamsMod, only : SF_val_cwd_frac + + ! ------------------------------------------------------------------------------------------- + ! + ! DESCRIPTION: + ! get the total site level 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 ! The total site area + 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(fates_patch_type), pointer :: currentPatch + type(fates_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 (exclude 0) + 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%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%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%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 + 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_harvestable_patch_carbon ( currentPatch, harvestable_patch_c ) + + !USES: + use SFParamsMod, only : SF_val_cwd_frac + + ! ------------------------------------------------------------------------------------------- + ! + ! DESCRIPTION: + ! get the total carbon (ha-1) availale for harvest in current Patch + ! harvestable carbon: aggregate all cohorts matching the dbhmin and dbhmax harvest criteria + ! + ! this subroutine shall be called inside the patch loop + ! output will be used for diagnosis or for scaling patch level harvest rate + + ! Arguments + type(fates_patch_type) , intent(in), target :: currentPatch + + real(r8), intent(out) :: harvestable_patch_c + + ! Local Variables + type(fates_cohort_type), pointer :: currentCohort + 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 + + ! loop over cohorts + 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 + + end subroutine get_harvestable_patch_carbon + + + ! ============================================================================ + + subroutine get_harvested_cohort_carbon ( currentCohort, lmort_direct, harvested_cohort_c ) + + !USES: + use SFParamsMod, only : SF_val_cwd_frac + + + ! ------------------------------------------------------------------------------------------- + ! + ! DESCRIPTION: + ! get the harvested carbon in current Cohort through lmort_direct + ! + ! this subroutine shall be called inside the cohort loop + ! output will be used for adjusting IFM wood harvest with size priority + + ! Arguments + type(fates_cohort_type) , intent(in), target :: currentCohort + real(r8), intent(in) :: lmort_direct + + real(r8), intent(out) :: harvested_cohort_c + + ! Local Variables + real(r8) :: sapw_m ! Biomass of sap wood + real(r8) :: struct_m ! Biomass of structural organs + + ! only account for cohorts matching the following conditions + sapw_m = currentCohort%prt%GetState(sapw_organ, carbon12_element) + struct_m = currentCohort%prt%GetState(struct_organ, carbon12_element) + + harvested_cohort_c = lmort_direct * ( sapw_m + struct_m ) * & + prt_params%allom_agb_frac(currentCohort%pft) * & + SF_val_CWD_frac(ncwd) * logging_export_frac * & + currentCohort%n + + end subroutine get_harvested_cohort_carbon + + ! ============================================================================ + +end module FatesForestManagementMod diff --git a/biogeochem/FatesLitterMod.F90 b/biogeochem/FatesLitterMod.F90 index 7d55dc9aab..e0733fb94c 100644 --- a/biogeochem/FatesLitterMod.F90 +++ b/biogeochem/FatesLitterMod.F90 @@ -447,7 +447,7 @@ subroutine adjust_SF_CWD_frac(dbh,ncwd,SF_val_CWD_frac,SF_val_CWD_frac_adj) !ARGUMENTS real(r8), intent(in) :: dbh !dbh of cohort [cm] - integer, intent(in) :: ncwd !number of cwd pools + integer, intent(in) :: ncwd !number of cwd pools real(r8), intent(in) :: SF_val_CWD_frac(:) !fates parameter specifying the !fraction of struct + sapw going !to each CWD class diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 old mode 100644 new mode 100755 index d86e5c5d51..0eb2b4d6a1 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -178,6 +178,10 @@ module FatesPatchMod ! 4) land use change real(r8) :: landuse_transition_rates(n_landuse_cats) ! land use tranision rate real(r8) :: fract_ldist_not_harvested ! fraction of logged area that is canopy trees that weren't harvested [0-1] + real(r8) :: harvest_rate_scale ! scaling factor applied to logging disturbance rate (primaryland, secondaryland, etc) + ! purpose is to assign patch-specific harvest priority based on certain IFM + ! strategy + integer :: order_age_since_anthro ! order of the patch based on the age since anthropogenic disturbance !--------------------------------------------------------------------------- @@ -363,6 +367,8 @@ subroutine NanValues(this) ! DISTURBANCE this%disturbance_rates(:) = nan this%fract_ldist_not_harvested = nan + this%harvest_rate_scale = nan + this%order_age_since_anthro = fates_unset_int ! LAND USE this%landuse_transition_rates(:) = nan @@ -440,6 +446,7 @@ subroutine ZeroValues(this) ! DISTURBANCE this%disturbance_rates(:) = 0.0_r8 this%fract_ldist_not_harvested = 0.0_r8 + this%harvest_rate_scale = 0.0_r8 ! LAND USE this%landuse_transition_rates(:) = 0.0_r8 @@ -745,7 +752,9 @@ subroutine Dump(this) write(fates_log(),*) 'pa%c_stomata = ',this%c_stomata write(fates_log(),*) 'pa%c_lblayer = ',this%c_lblayer write(fates_log(),*) 'pa%disturbance_rates = ',this%disturbance_rates(:) + write(fates_log(),*) 'pa%harvest_rate_scale = ',this%harvest_rate_scale write(fates_log(),*) 'pa%land_use_label = ',this%land_use_label + write(fates_log(),*) 'pa%order_age_since_anthro = ',this%order_age_since_anthro write(fates_log(),*) '----------------------------------------' do el = 1, num_elements diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index fd44f07bbe..72c93301c8 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -209,6 +209,11 @@ subroutine init_site_vars( site_in, bc_in, bc_out ) allocate(site_in%flux_diags(el)%root_litter_input(1:numpft)) end do + ! Mass balance + do el=1,num_elements + allocate(site_in%mass_balance(el)%wood_product_sz(1:nlevsclass)) + end do + ! Initialize the static soil ! arrays from the boundary (initial) condition site_in%zi_soil(:) = bc_in%zi_sisl(:) @@ -334,7 +339,8 @@ subroutine zero_site( site_in ) ! 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%harvest_debt_sec_y = 0.0_r8 + site_in%resources_management%harvest_debt_sec_m = 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 159e942c6a..d6d013d88d 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -67,7 +67,7 @@ module EDMainMod use FatesPatchMod , only : fates_patch_type use FatesCohortMod , only : fates_cohort_type use EDTypesMod , only : AREA - use EDTypesMod , only : site_massbal_type + use EDTypesMod , only : site_massbal_type, site_fluxdiags_type use PRTGenericMod , only : num_elements use PRTGenericMod , only : element_list use PRTGenericMod , only : element_pos @@ -86,9 +86,10 @@ 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 FatesForestManagementMod , only : get_harvestable_carbon use DamageMainMod , only : IsItDamageTime use EDPatchDynamicsMod , only : get_frac_site_primary + use EDPatchDynamicsMod , only : get_frac_site_secondary_mature use FatesGlobals , only : endrun => fates_endrun use ChecksBalancesMod , only : SiteMassStock use EDMortalityFunctionsMod , only : Mortality_Derivative @@ -373,6 +374,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) ! because it inherited them (such as daily carbon balance) real(r8) :: target_leaf_c real(r8) :: frac_site_primary + real(r8) :: frac_site_secondary_mature real(r8) :: harvestable_forest_c(hlm_num_lu_harvest_cats) integer :: harvest_tag(hlm_num_lu_harvest_cats) @@ -409,6 +411,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) !----------------------------------------------------------------------- call get_frac_site_primary(currentSite, frac_site_primary) + call get_frac_site_secondary_mature(currentSite, frac_site_secondary_mature) ! Clear site GPP and AR passing to HLM bc_out%gpp_site = 0._r8 @@ -474,7 +477,8 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) currentPatch%btran_ft, mean_temp, & currentPatch%land_use_label, & currentPatch%age_since_anthro_disturbance, frac_site_primary, & - harvestable_forest_c, harvest_tag) + frac_site_secondary_mature, harvestable_forest_c, & + currentPatch%harvest_rate_scale, harvest_tag) ! ----------------------------------------------------------------------------- ! Apply Plant Allocation and Reactive Transport @@ -858,6 +862,7 @@ subroutine TotalBalanceCheck (currentSite, call_index ) ! ! !LOCAL VARIABLES: type(site_massbal_type),pointer :: site_mass + type(site_fluxdiags_type), pointer :: flux_diags real(r8) :: biomass_stock ! total biomass in Kg/site real(r8) :: litter_stock ! total litter in Kg/site real(r8) :: seed_stock ! total seed mass in Kg/site @@ -901,6 +906,7 @@ subroutine TotalBalanceCheck (currentSite, call_index ) do el = 1, num_elements site_mass => currentSite%mass_balance(el) + flux_diags => currentSite%flux_diags(element_pos(carbon12_element)) call SiteMassStock(currentSite,el,total_stock,biomass_stock,litter_stock,seed_stock) @@ -929,7 +935,10 @@ subroutine TotalBalanceCheck (currentSite, call_index ) error_frac = 0.0_r8 end if - if ( error_frac > 10e-6_r8 .or. (error /= error) ) then + ! For debug + write(fates_log(),*) 'call index: ',call_index + write(fates_log(),*) 'wood_product: ',site_mass%wood_product + if ( error_frac > 10e-3_r8 .or. (error /= error) ) then write(fates_log(),*) 'mass balance error detected' write(fates_log(),*) 'element type (see PRTGenericMod.F90): ',element_list(el) write(fates_log(),*) 'error fraction relative to biomass stock: ',error_frac diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 old mode 100644 new mode 100755 index beaf11e140..e6d4064576 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -235,9 +235,18 @@ module EDParamsMod character(len=param_string_length), parameter, public :: maxcohort_name = "fates_maxcohort" - ! Logging Control Parameters (ONLY RELEVANT WHEN USE_FATES_LOGGING = TRUE) ! ---------------------------------------------------------------------------------------------- + integer,protected, public :: logging_age_preference ! switch for choosing patch age preference for logging + ! 1 - uniform, 2 - from oldest + character(len=param_string_length),parameter,public :: logging_name_age_preference = "fates_landuse_logging_age_preference" + + integer,protected, public :: logging_preference_options ! switch for choosing size or rotation preference for logging + ! 1 for uniform, 2 for double rotation, 3 for quadruple rotation, 4 for logistic, 5 for inverse logistic, 6 for gaussian + character(len=param_string_length),parameter,public :: logging_name_preference_options = "fates_landuse_logging_preference_options" + + real(r8),protected, public :: logging_ifm_harvest_scale ! scaling factor of harvest rate under improved forest management + character(len=param_string_length),parameter,public :: logging_name_ifm_harvest_scale = "fates_landuse_logging_ifm_harvest_scale" real(r8),protected,public :: logging_dbhmin ! Minimum dbh at which logging is applied (cm) ! Typically associated with harvesting @@ -340,6 +349,9 @@ subroutine FatesParamsInit() hydr_psicap = nan hydr_solver = -9 bgc_soil_salinity = nan + logging_age_preference = -9 + logging_preference_options = -9 + logging_ifm_harvest_scale = nan logging_dbhmin = nan logging_dbhmax = nan logging_collateral_frac = nan @@ -510,6 +522,15 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=bgc_name_soil_salinity, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) + call fates_params%RegisterParameter(name=logging_name_age_preference, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + + call fates_params%RegisterParameter(name=logging_name_preference_options, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + + call fates_params%RegisterParameter(name=logging_name_ifm_harvest_scale, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + call fates_params%RegisterParameter(name=logging_name_dbhmin, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -729,6 +750,17 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetrieveParameter(name=bgc_name_soil_salinity, & data=bgc_soil_salinity) + call fates_params%RetrieveParameter(name=logging_name_age_preference, & + data=tmpreal) + logging_age_preference = nint(tmpreal) + + call fates_params%RetrieveParameter(name=logging_name_preference_options, & + data=tmpreal) + logging_preference_options = nint(tmpreal) + + call fates_params%RetrieveParameter(name=logging_name_ifm_harvest_scale, & + data=logging_ifm_harvest_scale) + call fates_params%RetrieveParameter(name=logging_name_dbhmin, & data=logging_dbhmin) @@ -874,6 +906,9 @@ subroutine FatesReportParams(is_master) write(fates_log(),fmt0) 'hydro_psicap = ',hydr_psicap write(fates_log(),fmt0) 'hydro_solver = ',hydr_solver write(fates_log(),fmt0) 'bgc_soil_salinity = ', bgc_soil_salinity + write(fates_log(),fmt0) 'logging_age_preference = ',logging_age_preference + write(fates_log(),fmt0) 'logging_preference_options = ',logging_preference_options + write(fates_log(),fmt0) 'logging_ifm_harvest_scale = ',logging_ifm_harvest_scale write(fates_log(),fmt0) 'logging_dbhmin = ',logging_dbhmin write(fates_log(),fmt0) 'logging_dbhmax = ',logging_dbhmax write(fates_log(),fmt0) 'logging_collateral_frac = ',logging_collateral_frac diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 5745141c70..5e672b637d 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -1797,7 +1797,7 @@ subroutine FatesCheckParams(is_master) write(fates_log(),*) 'fates_rad_model = 1 or 2 ...' write(fates_log(),*) 'You specified fates_rad_model = ',radiation_model write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + ! call endrun(msg=errMsg(sourcefile, __LINE__)) end if if(.not.any(regeneration_model == [default_regeneration, & diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 old mode 100644 new mode 100755 index 34e5f319d7..6243c28fde --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -134,8 +134,12 @@ module EDTypesMod 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 + real(r8) :: harvest_debt_sec_y ! the amount of kgC per site from young secondary patches that did ! not successfully harvested + real(r8) :: harvest_debt_sec_m ! the amount of kgC per site from mature secondary patches that did + ! not successfully harvested + real(r8) :: harvest_wp_scale ! For Size-dependent IFM. The adjustment factor revising harvest scale to match a target + ! harvest wood product amount !debug variables real(r8) :: delta_litter_stock ! kgC/site = kgC/ha @@ -202,6 +206,7 @@ module EDTypesMod real(r8) :: seed_out ! Total mass of seeds exported outside of fates site [kg/site/day] ! (this is not used currently, placeholder, rgk feb-2019) + real(r8) :: frag_in ! Litter and coarse woody debris fragmentation flux [kg/site/day] real(r8) :: frag_out ! Litter and coarse woody debris fragmentation flux [kg/site/day] real(r8) :: wood_product ! Total mass exported as wood product [kg/site/day] @@ -214,6 +219,7 @@ module EDTypesMod real(r8) :: patch_resize_err ! This is the amount of mass gained (or loss when negative) ! due to re-sizing patches when area math starts to lose ! precision + real(r8), allocatable :: wood_product_sz(:) ! Mass exported as wood product from different size bins [kg/site/day] contains @@ -488,8 +494,10 @@ subroutine ZeroMassBalFlux(this) this%net_root_uptake = 0._r8 this%seed_in = 0._r8 this%seed_out = 0._r8 + this%frag_in = 0._r8 this%frag_out = 0._r8 this%wood_product = 0._r8 + this%wood_product_sz(:)= 0._r8 this%burn_flux_to_atm = 0._r8 this%flux_generic_in = 0._r8 this%flux_generic_out = 0._r8 diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 old mode 100644 new mode 100755 index df3a719204..90c8dcdb1b --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -125,9 +125,11 @@ module FatesConstantsMod integer, public :: fates_np_comp_scaling = fates_unset_int - real(fates_r8), parameter, public :: secondary_age_threshold = 94._fates_r8 ! less than this value is young secondary land + real(fates_r8), parameter, public :: secondary_age_threshold = 1._fates_r8 ! less than this value is young secondary land ! based on average age of global ! secondary 1900s land in hurtt-2011 + ! here revise to the age threshold for + ! logging cycle of IFM/ARR site simulations ! integer labels for specifying harvest units integer, parameter, public :: photosynth_acclim_model_none = 1 @@ -141,6 +143,18 @@ module FatesConstantsMod integer, parameter, public :: lmrmodel_ryan_1991 = 1 integer, parameter, public :: lmrmodel_atkin_etal_2017 = 2 + ! integer labels for logging age preference + integer, parameter, public :: logging_no_age_preference = 1 + integer, parameter, public :: logging_oldest_first = 2 + + ! integer labels for logging size preference and rotation length options + integer, parameter, public :: logging_uniform_size = 1 + integer, parameter, public :: logging_double_rotation = 2 + integer, parameter, public :: logging_quadruple_rotation = 3 + integer, parameter, public :: logging_logistic_size = 4 + integer, parameter, public :: logging_inv_logistic_size = 5 + integer, parameter, public :: logging_gaussian_size = 6 + ! Error Tolerances ! Allowable error in carbon allocations, should be applied to estimates @@ -173,6 +187,12 @@ module FatesConstantsMod ! precisions are preventing perfect zero in comparison real(fates_r8), parameter, public :: nearzero = 1.0e-30_fates_r8 + ! in nocomp simulations, what is the minimum PFT fraction for any given land use type? + real(fates_r8), parameter, public :: min_nocomp_pftfrac_perlanduse = 0.01_fates_r8 + + ! The smallest harvest rate considered in calculation + real(fates_r8), parameter, public :: min_harvest_rate = 1.0e-7_fates_r8 + ! Unit conversion constants: ! Conversion factor umols of Carbon -> kg of Carbon (1 mol = 12g) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 old mode 100644 new mode 100755 index 41ab2ca1b6..fd30b5b104 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -58,7 +58,9 @@ module FatesHistoryInterfaceMod use FatesAllometryMod , only : CrownDepth use FatesAllometryMod , only : bstore_allom use FatesAllometryMod , only : set_root_fraction - + use FatesLitterMod , only : adjust_SF_CWD_frac + use SFParamsMod , only : SF_val_cwd_frac + use EDPftvarcon , only : EDPftvarcon_inst use PRTParametersMod , only : prt_params @@ -316,7 +318,8 @@ module FatesHistoryInterfaceMod integer :: ih_fall_disturbance_rate_si integer :: ih_harvest_carbonflux_si integer :: ih_harvest_debt_si - integer :: ih_harvest_debt_sec_si + integer :: ih_harvest_debt_sec_m_si + integer :: ih_harvest_debt_sec_y_si ! Indices to site by size-class by age variables integer :: ih_nplant_si_scag @@ -467,6 +470,7 @@ module FatesHistoryInterfaceMod integer :: ih_ar_frootm_si_scpf integer :: ih_c13disc_si_scpf + integer :: ih_harvest_carbonflux_si_scls ! indices to (site x scls [size class bins]) variables integer :: ih_ba_si_scls @@ -2311,6 +2315,10 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) real(r8) :: storep_understory_scpf(numpft*nlevsclass) real(r8) :: storec_canopy_scpf(numpft*nlevsclass) real(r8) :: storec_understory_scpf(numpft*nlevsclass) + + ! Updated wood partitioning to CWD based on dbh + real(r8) :: SF_val_CWD_frac_adj(ncwd) + integer :: return_code integer :: i_dist, j_dist @@ -2393,8 +2401,10 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_logging_disturbance_rate_si => this%hvars(ih_logging_disturbance_rate_si)%r81d, & hio_fall_disturbance_rate_si => this%hvars(ih_fall_disturbance_rate_si)%r81d, & hio_harvest_carbonflux_si => this%hvars(ih_harvest_carbonflux_si)%r81d, & + hio_harvest_carbonflux_si_scls => this%hvars(ih_harvest_carbonflux_si_scls)%r82d, & 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_harvest_debt_sec_m_si => this%hvars(ih_harvest_debt_sec_m_si)%r81d, & + hio_harvest_debt_sec_y_si => this%hvars(ih_harvest_debt_sec_y_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, & @@ -2737,8 +2747,11 @@ 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 + + ! output site-level harvest debt [kgC day-1] -> [kgC yr-1] + hio_harvest_debt_si(io_si) = sites(s)%resources_management%harvest_debt * days_per_year + hio_harvest_debt_sec_m_si(io_si) = sites(s)%resources_management%harvest_debt_sec_m * days_per_year + hio_harvest_debt_sec_y_si(io_si) = sites(s)%resources_management%harvest_debt_sec_y * days_per_year ! 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 @@ -2762,7 +2775,7 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_fall_disturbance_rate_si(io_si) = sum(sites(s)%disturbance_rates(dtype_ifall,1:n_landuse_cats,1:n_landuse_cats)) * & days_per_year - hio_harvest_carbonflux_si(io_si) = sites(s)%mass_balance(element_pos(carbon12_element))%wood_product * AREA_INV + hio_harvest_carbonflux_si(io_si) = sites(s)%mass_balance(element_pos(carbon12_element))%wood_product * AREA_INV * days_per_year ! Loop through patches to sum up diagonistics ipa = 0 @@ -3326,6 +3339,16 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_c13disc_si_scpf(io_si,scpf) = 0.0_r8 endif + !Harvest C flux (must be calculated before since cohort/patch has been modified) + !call adjust_SF_CWD_frac(ccohort%dbh, ncwd, SF_val_CWD_frac, SF_val_CWD_frac_adj) + !hio_harvest_carbonflux_si_scpf(io_si,scpf) = hio_harvest_carbonflux_si_scpf(io_si,scpf) + & + ! ccohort%n * ccohort%lmort_direct * ( struct_m + sapw_m ) * & + ! prt_params%allom_agb_frac(ccohort%pft) * SF_val_CWD_frac_adj(ncwd) * & + ! AREA_INV * days_per_year + hio_harvest_carbonflux_si_scls(io_si,scls) = sites(s)%mass_balance(element_pos(carbon12_element))%wood_product_sz(scls) * AREA_INV * days_per_year + ! Reset to zero, Not sure why I need this but it will not be automatically reset which seems strange to me + ccohort%harv_c = 0._r8 + ! number density [/m2] hio_nplant_si_scpf(io_si,scpf) = hio_nplant_si_scpf(io_si,scpf) + ccohort%n / m2_per_ha @@ -3507,7 +3530,7 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) ccohort%frmort + ccohort%smort + ccohort%asmort + ccohort%dgmort) * & total_m * ccohort%n * days_per_sec * years_per_day * ha_per_m2 + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * total_m * & - ccohort%n * ha_per_m2 + ccohort%n * ha_per_m2 * days_per_sec * years_per_day hio_canopy_mortality_crownarea_si(io_si) = hio_canopy_mortality_crownarea_si(io_si) + & (ccohort%bmort + ccohort%hmort + ccohort%cmort + & @@ -3657,7 +3680,7 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) ccohort%frmort + ccohort%smort + ccohort%asmort + ccohort%dgmort) * & total_m * ccohort%n * days_per_sec * years_per_day * ha_per_m2 + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * total_m * & - ccohort%n * ha_per_m2 + ccohort%n * ha_per_m2 * days_per_sec * years_per_day hio_understory_mortality_crownarea_si(io_si) = hio_understory_mortality_crownarea_si(io_si) + & (ccohort%bmort + ccohort%hmort + ccohort%cmort + & @@ -3964,14 +3987,14 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) ! ! carbon flux associated with mortality of trees dying by fire hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & - sum(sites(s)%fmort_carbonflux_canopy(:)) / g_per_kg + sum(sites(s)%fmort_carbonflux_canopy(:)) * years_per_day / g_per_kg hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - sum(sites(s)%fmort_carbonflux_ustory(:)) / g_per_kg + sum(sites(s)%fmort_carbonflux_ustory(:)) * years_per_day / g_per_kg ! treat carbon flux from imort the same way hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - sum(sites(s)%imort_carbonflux(:)) + sum(sites(s)%imort_carbonflux(:)) * years_per_day do i_pft = 1, numpft hio_mortality_carbonflux_si_pft(io_si,i_pft) = hio_mortality_carbonflux_si_pft(io_si,i_pft) + & @@ -4444,10 +4467,10 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) ! mortality-associated carbon fluxes hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & - sum(sites(s)%term_carbonflux_canopy(:)) * days_per_sec * ha_per_m2 + sum(sites(s)%term_carbonflux_canopy(:)) * days_per_sec * years_per_day * ha_per_m2 hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - sum(sites(s)%term_carbonflux_ustory(:)) * days_per_sec * ha_per_m2 + sum(sites(s)%term_carbonflux_ustory(:)) * days_per_sec * years_per_day * ha_per_m2 ! add site level mortality counting to crownarea diagnostic hio_canopy_mortality_crownarea_si(io_si) = hio_canopy_mortality_crownarea_si(io_si) + & @@ -5777,14 +5800,14 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_SECONDAREA_ANTHRODIST_AP', & units='m2 m-2', & long='secondary forest patch area age distribution since anthropgenic disturbance', & - use_default='inactive', avgflag='A', vtype=site_age_r8, & + use_default='active', avgflag='A', vtype=site_age_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_agesince_anthrodist_si_age) call this%set_history_var(vname='FATES_SECONDAREA_DIST_AP', & units='m2 m-2', & long='secondary forest patch area age distribution since any kind of disturbance', & - use_default='inactive', avgflag='A', vtype=site_age_r8, & + use_default='active', avgflag='A', vtype=site_age_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, initialize=initialize_variables, & index=ih_secondarylands_area_si_age) @@ -5941,7 +5964,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_SEED_BANK', units='kg m-2', & long='total seed mass of all PFTs in kg carbon per m2 land area', & - use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + use_default='active', avgflag='I', vtype=site_r8, hlms='CLM:ALM', & upfreq=1, ivar=ivar, initialize=initialize_variables, & index = ih_seed_bank_si) @@ -6040,7 +6063,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_VEGC', units='kg m-2', & long='total biomass in live plants in kg carbon per m2 land area', & - use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + use_default='active', avgflag='I', vtype=site_r8, hlms='CLM:ALM', & upfreq=1, ivar=ivar, initialize=initialize_variables, & index = ih_totvegc_si) @@ -6454,6 +6477,14 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index = ih_harvest_carbonflux_si) + call this%set_history_var(vname='FATES_HARVEST_CARBON_FLUX_SZ', & + units='kg m-2 s-1', & + long='harvest carbon flux in kg carbon per m2 per second per size class', & + use_default='active', avgflag='A', & + vtype=site_size_r8, hlms='CLM:ALM', & + upfreq=1, ivar=ivar, initialize=initialize_variables, & + index = ih_harvest_carbonflux_si_scls) + ! Canopy Resistance call this%set_history_var(vname='FATES_STOMATAL_COND', & @@ -6499,22 +6530,25 @@ subroutine define_history_vars(this, initialize_variables) 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_AR', units='gC/m^2/s', & - long='autotrophic respiration', 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_AR', units='gC/m^2/s', & +! long='autotrophic respiration', 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_HARVEST_DEBT', units='kg C', & + call this%set_history_var(vname='FATES_HARVEST_DEBT', units='kg C yr-1', & 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', & + call this%set_history_var(vname='FATES_HARVEST_DEBT_SEC_MATURE', units='kg C yr-1', & + long='Accumulated carbon failed to be harvested from mature 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 ) - - + ivar=ivar, initialize=initialize_variables, index = ih_harvest_debt_sec_m_si ) + + call this%set_history_var(vname='FATES_HARVEST_DEBT_SEC_YOUNG', units='kg C yr-1', & + long='Accumulated carbon failed to be harvested from young and non-forest 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_y_si ) ! Ecosystem Carbon Fluxes (updated rapidly, upfreq=2) @@ -6612,7 +6646,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_GPP_AP', units='kg m-2 s-1', & long='gross primary productivity by age bin in kg carbon per m2 per second', & - use_default='inactive', avgflag='A', vtype=site_age_r8, & + use_default='active', avgflag='A', vtype=site_age_r8, & hlms='CLM:ALM', upfreq=2, ivar=ivar, initialize=initialize_variables, & index = ih_gpp_si_age) @@ -6973,7 +7007,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_VEGC_APPF',units = 'kg m-2', & long='biomass per PFT in each age bin in kg carbon per m2', & - use_default='inactive', avgflag='A', vtype=site_agepft_r8, & + use_default='active', avgflag='A', vtype=site_agepft_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, & initialize=initialize_variables, index = ih_biomass_si_agepft) @@ -6990,7 +7024,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_GPP_SZPF', units='kg m-2 s-1', & long='gross primary production by pft/size in kg carbon per m2 per second', & - use_default='inactive', avgflag='A', vtype=site_size_pft_r8, & + use_default='active', avgflag='A', vtype=site_size_pft_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, & initialize=initialize_variables, index = ih_gpp_si_scpf) @@ -7024,7 +7058,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_NPP_SZPF', units='kg m-2 s-1', & long='total net primary production by pft/size in kg carbon per m2 per second', & - use_default='inactive', avgflag='A', vtype=site_size_pft_r8, & + use_default='active', avgflag='A', vtype=site_size_pft_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, & initialize=initialize_variables, index = ih_npp_totl_si_scpf) @@ -7123,7 +7157,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_VEGC_ABOVEGROUND_SZPF', & units = 'kg m-2', & long='aboveground biomass by pft/size in kg carbon per m2', & - use_default='inactive', avgflag='A', vtype=site_size_pft_r8, & + use_default='active', avgflag='A', vtype=site_size_pft_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, initialize=initialize_variables, & index = ih_agb_si_scpf) @@ -7367,7 +7401,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_AUTORESP_SZPF', & units = 'kg m-2 s-1', & long='total autotrophic respiration in kg carbon per m2 per second by pft/size', & - use_default='inactive', avgflag='A', vtype=site_size_pft_r8, & + use_default='active', avgflag='A', vtype=site_size_pft_r8, & hlms='CLM:ALM', upfreq=2, ivar=ivar, & initialize=initialize_variables, index = ih_ar_si_scpf) @@ -7464,14 +7498,14 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_DEMOTION_RATE_SZ', & units = 'm-2 yr-1', & long='demotion rate from canopy to understory by size class in number of plants per m2 per year', & - use_default='inactive', avgflag='A', vtype=site_size_r8, & + use_default='active', avgflag='A', vtype=site_size_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, & initialize=initialize_variables, index = ih_demotion_rate_si_scls) call this%set_history_var(vname='FATES_PROMOTION_RATE_SZ', & units = 'm-2 yr-1', & long='promotion rate from understory to canopy by size class', & - use_default='inactive', avgflag='A', vtype=site_size_r8, & + use_default='active', avgflag='A', vtype=site_size_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, & initialize=initialize_variables, index = ih_promotion_rate_si_scls) @@ -7787,7 +7821,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_SAPWOOD_ALLOC_CANOPY_SZ', & units = 'kg m-2 s-1', & long='allocation to sapwood C for canopy plants by size class in kg carbon per m2 per second', & - use_default='inactive', avgflag='A', vtype=site_size_r8, & + use_default='active', avgflag='A', vtype=site_size_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, initialize=initialize_variables, & index = ih_npp_sapw_canopy_si_scls) @@ -7948,7 +7982,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_SAPWOOD_ALLOC_USTORY_SZ', & units = 'kg m-2 s-1', & long='allocation to sapwood C for understory plants by size class in kg carbon per m2 per second', & - use_default='inactive', avgflag='A', vtype=site_size_r8, & + use_default='active', avgflag='A', vtype=site_size_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, initialize=initialize_variables, & index = ih_npp_sapw_understory_si_scls) @@ -8187,7 +8221,7 @@ subroutine define_history_vars(this, initialize_variables) ! CARBON call this%set_history_var(vname='FATES_VEGC_SZPF', units='kg m-2', & long='total vegetation biomass in live plants by size-class x pft in kg carbon per m2', & - use_default='inactive', avgflag='A', vtype=site_size_pft_r8, & + use_default='active', avgflag='A', vtype=site_size_pft_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, initialize=initialize_variables, & index = ih_totvegc_scpf) @@ -8259,7 +8293,6 @@ subroutine define_history_vars(this, initialize_variables) upfreq=1, ivar=ivar, initialize=initialize_variables, & index = ih_npp_stor_si) - ! PLANT HYDRAULICS hydro_active_if: if(hlm_use_planthydro.eq.itrue) then @@ -8376,7 +8409,7 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='FATES_BTRAN_SZPF', units='1', & long='mean individual level BTRAN by size class x pft', & - use_default='inactive', avgflag='A', vtype=site_size_pft_r8, & + use_default='active', avgflag='A', vtype=site_size_pft_r8, & hlms='CLM:ALM', upfreq=4, ivar=ivar, & initialize=initialize_variables, index = ih_btran_scpf) diff --git a/main/FatesParametersInterface.F90 b/main/FatesParametersInterface.F90 old mode 100644 new mode 100755 index aa0ef85287..7acbc24045 --- a/main/FatesParametersInterface.F90 +++ b/main/FatesParametersInterface.F90 @@ -12,7 +12,7 @@ module FatesParametersInterface integer, parameter, public :: max_params = 250 integer, parameter, public :: max_dimensions = 2 integer, parameter, public :: max_used_dimensions = 25 - integer, parameter, public :: param_string_length = 40 + integer, parameter, public :: param_string_length = 50 ! NOTE(bja, 2017-02) these are the values returned from netcdf after ! inquiring about the number of dimensions integer, parameter, public :: dimension_shape_scalar = 0 diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index e288664751..39a6c95982 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -762,6 +762,15 @@ variables: double fates_hydro_solver ; fates_hydro_solver:units = "unitless" ; fates_hydro_solver:long_name = "switch designating which numerical solver for plant hydraulics, 1 = 1D taylor, 2 = 2D Picard, 3 = 2D Newton (deprecated)" ; + double fates_landuse_logging_age_preference ; + fates_landuse_logging_age_preference:units = "index" ; + fates_landuse_logging_age_preference:long_name = "Integer code of logging age preference options, 1 = no preference, 2 = oldest first" ; + double fates_landuse_logging_preference_options ; + fates_landuse_logging_preference_options:units = "index" ; + fates_landuse_logging_preference_options:long_name = "Integer code of logging size and rotation length preference options, see codes for detail (1-6)" ; + double fates_landuse_logging_ifm_harvest_scale ; + fates_landuse_logging_ifm_harvest_scale:units = "fraction" ; + fates_landuse_logging_ifm_harvest_scale:long_name = "Scaling factor of target harvest rate under improved forest management, normally shall be 0 - 1" ; double fates_landuse_logging_coll_under_frac ; fates_landuse_logging_coll_under_frac:units = "fraction" ; fates_landuse_logging_coll_under_frac:long_name = "Fraction of stems killed in the understory when logging generates disturbance" ; @@ -1657,6 +1666,12 @@ data: fates_hydro_solver = 1 ; + fates_landuse_logging_age_preference = 1 ; + + fates_landuse_logging_preference_options = 1 ; + + fates_landuse_logging_ifm_harvest_scale = 1.0 ; + fates_landuse_logging_coll_under_frac = 0.55983 ; fates_landuse_logging_collateral_frac = 0.05 ;