From 09aed9fa84a1d3032583cdd6372b763b35ba5928 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Thu, 18 Jan 2024 18:20:49 -0800 Subject: [PATCH 01/16] Bunch of minor modification for MRV initial test --- biogeochem/EDLoggingMortalityMod.F90 | 4 +++- biogeochem/EDPatchDynamicsMod.F90 | 3 ++- biogeochem/FatesLitterMod.F90 | 2 +- main/FatesConstantsMod.F90 | 4 +++- main/FatesHistoryInterfaceMod.F90 | 4 ++-- 5 files changed, 11 insertions(+), 6 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index bf6ab7443c..f8be95bda1 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -248,7 +248,9 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! 0=use fates logging parameters directly when logging_time == .true. ! this means harvest the whole cohort area harvest_rate = 1._r8 - + harvest_tag = 2 + cur_harvest_tag = 2 + 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 diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 140c108d66..54e96c561b 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -3027,6 +3027,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 +3039,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 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/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index df3a719204..c9a0920a53 100644 --- 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 = 10._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 diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 41ab2ca1b6..2bb97c96ad 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -5777,14 +5777,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) From 860229e7227d2a38464ded63a90af6554e97a401 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Thu, 18 Jan 2024 18:20:49 -0800 Subject: [PATCH 02/16] Bunch of minor modification for MRV initial test --- biogeochem/EDLoggingMortalityMod.F90 | 4 +++- biogeochem/EDPatchDynamicsMod.F90 | 3 ++- biogeochem/FatesLitterMod.F90 | 2 +- main/FatesConstantsMod.F90 | 4 +++- main/FatesHistoryInterfaceMod.F90 | 4 ++-- 5 files changed, 11 insertions(+), 6 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index bf6ab7443c..f8be95bda1 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -248,7 +248,9 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! 0=use fates logging parameters directly when logging_time == .true. ! this means harvest the whole cohort area harvest_rate = 1._r8 - + harvest_tag = 2 + cur_harvest_tag = 2 + 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 diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 140c108d66..54e96c561b 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -3027,6 +3027,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 +3039,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 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/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index df3a719204..c9a0920a53 100644 --- 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 = 10._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 diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 41ab2ca1b6..2bb97c96ad 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -5777,14 +5777,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) From 4780b1c5169727ba926dc60742f5d9733bc3e00c Mon Sep 17 00:00:00 2001 From: sshu88 Date: Fri, 19 Jan 2024 15:29:27 -0800 Subject: [PATCH 03/16] Some more minor revisions. --- biogeochem/EDLoggingMortalityMod.F90 | 28 +++++++++++++++------------- biogeochem/EDPatchDynamicsMod.F90 | 13 ++++++++++++- main/EDInitMod.F90 | 3 ++- main/EDPftvarcon.F90 | 2 +- main/EDTypesMod.F90 | 4 +++- main/FatesHistoryInterfaceMod.F90 | 28 ++++++++++++++++++---------- 6 files changed, 51 insertions(+), 27 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index f8be95bda1..90dfe85735 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -240,16 +240,18 @@ 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 - harvest_tag = 2 - cur_harvest_tag = 2 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. @@ -270,10 +272,6 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & call get_harvest_rate_area (patch_land_use_label, hlm_harvest_catnames, & hlm_harvest_rates, frac_site_primary, secondary_age, harvest_rate) - ! For area-based harvest, harvest_tag shall always be 2 (not applicable). - harvest_tag = 2 - cur_harvest_tag = 2 - if (fates_global_verbose()) then write(fates_log(), *) 'Successfully Read Harvest Rate from HLM.', hlm_harvest_rates(:), harvest_rate end if @@ -295,7 +293,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! 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 @@ -579,7 +577,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 @@ -1154,7 +1151,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 @@ -1162,6 +1159,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 @@ -1184,19 +1186,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 diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 54e96c561b..bb80f249ff 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -212,6 +212,7 @@ subroutine disturbance_rates( site_in, bc_in) integer :: harvest_tag(hlm_num_lu_harvest_cats) 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] + !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) ! And the same rates in understory plants have already been applied to %dndt @@ -223,6 +224,9 @@ subroutine disturbance_rates( site_in, bc_in) ! get available biomass for harvest for all patches call get_harvestable_carbon(site_in, bc_in%site_area, bc_in%hlm_harvest_catnames, harvestable_forest_c) + ! Initialize local variables + harvest_tag_csite = 2 + currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -262,13 +266,20 @@ subroutine disturbance_rates( site_in, bc_in) 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) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index fd44f07bbe..02a6d8999a 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -334,7 +334,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/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 index 34e5f319d7..16570c1a68 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -134,7 +134,9 @@ 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 !debug variables diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 2bb97c96ad..04b96ff8cc 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -316,7 +316,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 @@ -2394,7 +2395,8 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_fall_disturbance_rate_si => this%hvars(ih_fall_disturbance_rate_si)%r81d, & hio_harvest_carbonflux_si => this%hvars(ih_harvest_carbonflux_si)%r81d, & hio_harvest_debt_si => this%hvars(ih_harvest_debt_si)%r81d, & - hio_harvest_debt_sec_si => this%hvars(ih_harvest_debt_sec_si)%r81d, & + hio_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 +2739,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 @@ -6504,17 +6509,20 @@ subroutine define_history_vars(this, initialize_variables) 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) From 5d40c9e635fedc58d99d082975d3d318da103252 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Fri, 19 Jan 2024 17:42:44 -0800 Subject: [PATCH 04/16] Added lost variables when merging. --- biogeochem/EDPatchDynamicsMod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index bb80f249ff..3842d670e4 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -210,6 +210,7 @@ subroutine disturbance_rates( site_in, bc_in) real(r8) :: mean_temp real(r8) :: harvestable_forest_c(hlm_num_lu_harvest_cats) integer :: harvest_tag(hlm_num_lu_harvest_cats) + integer :: harvest_tag_csite(hlm_num_lu_harvest_cats) ! harvest tag of current site 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] From 371da0398a3465ac58d5e69bc2d86be05c736bbe Mon Sep 17 00:00:00 2001 From: sshu88 Date: Mon, 29 Apr 2024 14:24:47 -0700 Subject: [PATCH 05/16] Add harvest strategy to maximize harvest from the oldest patch first. --- biogeochem/EDLoggingMortalityMod.F90 | 55 +++++--- biogeochem/EDMortalityFunctionsMod.F90 | 9 +- biogeochem/EDPatchDynamicsMod.F90 | 171 ++++++++++++++++++++++++- biogeochem/FatesPatchMod.F90 | 6 + main/EDMainMod.F90 | 6 +- 5 files changed, 220 insertions(+), 27 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 90dfe85735..7595de7025 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -200,8 +200,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) ! Arguments integer, intent(in) :: pft_i ! pft index @@ -215,6 +216,9 @@ 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 harvest rate proportionally 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 @@ -270,11 +274,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) - 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 @@ -300,7 +304,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! 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 + lmort_direct = harvest_rate_scale * harvest_rate * logging_direct_frac else lmort_direct = 0.0_r8 end if @@ -312,13 +316,13 @@ 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 = harvest_rate_scale * harvest_rate * logging_mechanical_frac 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 = harvest_rate_scale * harvest_rate * logging_collateral_frac else lmort_collateral = 0._r8 endif @@ -326,12 +330,12 @@ 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 = harvest_rate_scale * harvest_rate * logging_mechanical_frac 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 @@ -349,7 +353,7 @@ 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 +368,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 +398,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 diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index fc941d0371..57b133f543 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -235,7 +235,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 +256,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 +296,9 @@ 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) if (currentCohort%canopy_layer > 1)then ! Include understory logging mortality rates not associated with disturbance diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 3842d670e4..4b89b42434 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -80,6 +80,7 @@ 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 FatesLandUseChangeMod, only : get_landuse_transition_rates use FatesConstantsMod , only : fates_unset_r8 @@ -124,6 +125,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__ @@ -204,30 +206,115 @@ subroutine disturbance_rates( site_in, bc_in) integer :: threshold_sizeclass integer :: i_dist integer :: h_index + 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) 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] + ! Control vars, will be moved into parameter file once the test is done + integer, parameter :: harvest_age_priority = 1 ! Tag to determine if we give priority to the oldest + ! patch + ! 0 - no; + ! 1 - yes, maximize the oldest; + ! 2 - yes, sigmoid (To be added). + !---------------------------------------------------------------------------------------------- ! 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 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 call get_harvestable_carbon(site_in, bc_in%site_area, bc_in%hlm_harvest_catnames, harvestable_forest_c) ! Initialize local variables + frac_site = 1._r8 harvest_tag_csite = 2 + remain_harvest_rate = fates_unset_r8 ! set to negative value to indicate uninitialized + ! Loop through all patches and calculate the harvest rate scale based + ! on pre-defined age priority strategy + currentPatch => site_in%oldest_patch + + do while (associated(currentPatch)) + ! Initialize harvest rate scale to 1.0 first + currentPatch%harvest_rate_scale = 1._r8 + + ! Only calculate harvest rate scale when we have logging event + ! to reduce computational cost + if (logging_time) then + ! Maximize harvest rate (99%) for older patch first until harvest demand is fullfilled + ! Right now we skip the calculation for the primary land + if(harvest_age_priority == 1 .and. currentPatch%land_use_label .eq. secondaryland) then + if(currentPatch%land_use_label .eq. primaryland) then + h_index = 1 + frac_site = frac_site_primary + else if (currentPatch%land_use_label .eq. secondaryland .and. & + 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 harvest rate 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 + ! Check if is the first time, if so initialize remain_harvest_rate + if(remain_harvest_rate(h_index) < 0) then + remain_harvest_rate(h_index) = harvest_rate + end if + ! Only calculate harvest rate scale for non-zero harvest rate + if (harvest_rate > 1e-7) then + ! 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)+1e-7) / harvest_rate + remain_harvest_rate(h_index) = 0._r8 + else + ! harvest almost 100%, leave 1% to prevent removing the patch completely, which may cause + ! model crash under nocomp mode + currentPatch%harvest_rate_scale = 0.99_r8 / harvest_rate + remain_harvest_rate(h_index) = remain_harvest_rate(h_index) - 0.99_r8 * (currentPatch%area/(AREA*frac_site)) + end if + end if + end if ! harvest age priority strategy + write(fates_log(),*) 'See patch number:', currentPatch%patchno + write(fates_log(),*) 'See patch age:', currentPatch%age + write(fates_log(),*) 'See patch age sec:', currentPatch%age_since_anthro_disturbance + write(fates_log(),*) 'See harvest index:', h_index + write(fates_log(),*) 'See harvest rate scale:', currentPatch%harvest_rate_scale + write(fates_log(),*) 'See harvest rate:', harvest_rate + 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 ! logging_time + + currentPatch => currentPatch%younger + end do + write(fates_log(),*) '<<<<<<<<<<<<<<<<<<<<<<< Next section >>>>>>>>>>>>>>>>>>>>>>>>>>' + + !Loop through cohorts to get disturbance rates currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -259,9 +346,18 @@ 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, & + currentPatch%harvest_rate_scale, & harvest_tag) - + if (logging_time) then + write(fates_log(),*) 'before logging, See patch number:', currentPatch%patchno + if(lmort_direct > 0._r8) write(fates_log(),*) 'lmort_direct:', lmort_direct + if(lmort_collateral > 0._r8) write(fates_log(),*) 'lmort_collateral:', lmort_collateral + if(lmort_infra > 0._r8) write(fates_log(),*) 'lmort_infra:', lmort_infra + if(l_degrad > 0._r8) write(fates_log(),*) 'l_degrad:', l_degrad + write(fates_log(),*) '================== Next Cohort ==================' + end if currentCohort%lmort_direct = lmort_direct currentCohort%lmort_collateral = lmort_collateral currentCohort%lmort_infra = lmort_infra @@ -379,15 +475,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 @@ -422,6 +521,13 @@ subroutine disturbance_rates( site_in, bc_in) end do endif + if (logging_time) then + write(fates_log(),*) 'See patch number:', currentPatch%patchno + write(fates_log(),*) 'See harvest rate scale:', currentPatch%harvest_rate_scale + write(fates_log(),*) 'See logging disturbance rate:', currentPatch%disturbance_rates(dtype_ilog) + write(fates_log(),*) 'See area:', currentPatch%area/AREA + write(fates_log(),*) '==================== Next patch ================' + end if currentPatch => currentPatch%younger enddo !patch loop @@ -581,7 +687,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() @@ -2425,12 +2531,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 = 0 ! 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 = 5._r8 primary_land_fraction_beforefusion = 0._r8 primary_land_fraction_afterfusion = 0._r8 @@ -2572,12 +2683,22 @@ subroutine fuse_patches( csite, bc_in ) !---------------------------------------------------------------------! ! Look for differences in profile biomass, above the minimum biomass ! !---------------------------------------------------------------------! - if(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 @@ -2614,6 +2735,7 @@ subroutine fuse_patches( csite, bc_in ) !------------------------------------------------------------------------! profiletol = ED_val_patch_fusion_tol + agetol = 5._r8 endif fuseflagset_if endif different_patches_if @@ -2649,6 +2771,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 = 0._r8 !---------------------------------------------------------------------! ! Making profile tolerance larger means that more fusion will happen ! @@ -3217,6 +3340,40 @@ end subroutine get_frac_site_primary ! ===================================================================================== + subroutine get_frac_site_secondary_mature(site_in, frac_site_secondary_mature) + + ! + ! !DESCRIPTION: + ! Calculate how much of a site is secondary mature land + ! + ! !USES: + use EDTypesMod , only : ed_site_type + ! + ! !ARGUMENTS: + 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 + + 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 + + end subroutine get_frac_site_secondary_mature + + ! ===================================================================================== + subroutine InsertPatch(currentSite, newPatch) ! !DESCRIPTION: diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index d86e5c5d51..d856dbc55f 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -178,6 +178,9 @@ 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 + ! strategy !--------------------------------------------------------------------------- @@ -363,6 +366,7 @@ subroutine NanValues(this) ! DISTURBANCE this%disturbance_rates(:) = nan this%fract_ldist_not_harvested = nan + this%harvest_rate_scale = nan ! LAND USE this%landuse_transition_rates(:) = nan @@ -440,6 +444,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,6 +750,7 @@ 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(),*) '----------------------------------------' diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 159e942c6a..66613a0bec 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -89,6 +89,7 @@ module EDMainMod use EDLoggingMortalityMod , 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 From 83d6d98ca1ee617d32c74320922d56222c004969 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Tue, 30 Apr 2024 11:14:29 -0700 Subject: [PATCH 06/16] Added bubble sort for patches. --- biogeochem/EDPatchDynamicsMod.F90 | 44 ++++++++++++++++++++++++++++++- biogeochem/FatesPatchMod.F90 | 3 +++ 2 files changed, 46 insertions(+), 1 deletion(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 4b89b42434..4a710c529e 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -206,6 +206,7 @@ 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 @@ -224,7 +225,10 @@ subroutine disturbance_rates( site_in, bc_in) ! patch ! 0 - no; ! 1 - yes, maximize the oldest; - ! 2 - yes, sigmoid (To be added). + ! 2 - yes, sigmoid (To be added). and more other? + integer :: exchanged ! used in bubble sorting + integer :: temp_order ! used in bubble sorting + integer :: ipatch ! used in bubble sorting !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) @@ -239,10 +243,47 @@ subroutine disturbance_rates( site_in, bc_in) call get_harvestable_carbon(site_in, bc_in%site_area, bc_in%hlm_harvest_catnames, harvestable_forest_c) ! Initialize local variables + exchanged = 1 + npatches = 0 frac_site = 1._r8 harvest_tag_csite = 2 remain_harvest_rate = fates_unset_r8 ! set to negative value to indicate uninitialized + ! We need to have an order based on age_since_anthro_disturbance + if (logging_time) then + + ! First loop to count number of patches and initialize order + ! with patchno + currentPatch => site_in%oldest_patch + + do while (associated(currentPatch)) + npatches = npatches + 1 + currentPatch%order_age_since_anthro = currentPatch%patchno + + currentPatch => currentPatch%younger + end do + + ! Second loop to perform bubble sorting + do ipatch = 1, npatches - 1 + if(exchanged .eq. 1) then + exchanged = 0 + currentPatch => site_in%oldest_patch + do while (associated(currentPatch)) + if(currentPatch%patchno <= (npatches - ipatch)) then + if(currentPatch%age_since_anthro_disturbance < & + currentPatch%younger%age_since_anthro_disturbance) then + temp_order = currentPatch%order_age_since_anthro + currentPatch%order_age_since_anthro = currentPatch%younger%order_age_since_anthro + currentPatch%younger%order_age_since_anthro = temp_order + exchanged = 1 + end if + end if + end do + end if + end do + + end if ! logging_time + ! Loop through all patches and calculate the harvest rate scale based ! on pre-defined age priority strategy currentPatch => site_in%oldest_patch @@ -300,6 +341,7 @@ subroutine disturbance_rates( site_in, bc_in) write(fates_log(),*) 'See patch number:', currentPatch%patchno write(fates_log(),*) 'See patch age:', currentPatch%age 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 harvest rate:', harvest_rate diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index d856dbc55f..919d619681 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -181,6 +181,7 @@ module FatesPatchMod 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 ! strategy + integer :: order_age_since_anthro ! order of the patch based on the age since anthropogenic disturbance !--------------------------------------------------------------------------- @@ -367,6 +368,7 @@ subroutine NanValues(this) 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 @@ -752,6 +754,7 @@ subroutine Dump(this) 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 From 21b9c164bc2b747ecd6cbc8163f88f440a12274f Mon Sep 17 00:00:00 2001 From: sshu88 Date: Thu, 2 May 2024 21:20:10 -0700 Subject: [PATCH 07/16] Revise the bubble sort algorithm. --- biogeochem/EDPatchDynamicsMod.F90 | 77 ++++++++++++++++++++++++++----- 1 file changed, 65 insertions(+), 12 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 4a710c529e..08d0cd98e5 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -228,7 +228,11 @@ subroutine disturbance_rates( site_in, bc_in) ! 2 - yes, sigmoid (To be added). and more other? integer :: exchanged ! used in bubble sorting integer :: temp_order ! used in bubble sorting - integer :: ipatch ! used in bubble sorting + real(r8) :: temp_age ! used in bubble sorting + integer :: ipatch, jpatch ! used in bubble sorting + integer, allocatable :: order_of_patches(:) ! Record a copy of patch order for bubble sort + real(r8), allocatable :: sec_age_patches(:) ! Record a copy of secondary forest age for bubble sort + !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) @@ -263,25 +267,74 @@ subroutine disturbance_rates( site_in, bc_in) currentPatch => currentPatch%younger end do - ! Second loop to perform bubble sorting + ! Second loop to assign order patches and age_since_anthro_disturbance + allocate(order_of_patches(1:npatches)) + allocate(sec_age_patches(1:npatches)) + + currentPatch => site_in%oldest_patch + ipatch = 0 + do while (associated(currentPatch)) + ipatch = ipatch + 1 + order_of_patches(ipatch) = currentPatch%patchno + sec_age_patches(ipatch) = currentPatch%age_since_anthro_disturbance + + currentPatch => currentPatch%younger + end do + + !write(fates_log(),*) 'See patch order:', order_of_patches(:) + !write(fates_log(),*) 'See patch age:', sec_age_patches(:) + ! Third loops to perform bubble sorting, since we have collected the order + ! and age and stored them into array, no need to have a patch loop do ipatch = 1, npatches - 1 if(exchanged .eq. 1) then exchanged = 0 - currentPatch => site_in%oldest_patch - do while (associated(currentPatch)) - if(currentPatch%patchno <= (npatches - ipatch)) then - if(currentPatch%age_since_anthro_disturbance < & - currentPatch%younger%age_since_anthro_disturbance) then - temp_order = currentPatch%order_age_since_anthro - currentPatch%order_age_since_anthro = currentPatch%younger%order_age_since_anthro - currentPatch%younger%order_age_since_anthro = temp_order - exchanged = 1 - end if + do jpatch = 1, npatches - ipatch + !write(fates_log(),*) 'ipatch:', ipatch + !write(fates_log(),*) 'jpatch:', jpatch + !write(fates_log(),*) 'nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn' + if(sec_age_patches(jpatch) < sec_age_patches(jpatch+1)) then + ! Swap order + temp_order = order_of_patches(jpatch) + order_of_patches(jpatch) = order_of_patches(jpatch + 1) + order_of_patches(jpatch + 1) = temp_order + ! Swap age + temp_age = sec_age_patches(jpatch) + sec_age_patches(jpatch) = sec_age_patches(jpatch + 1) + sec_age_patches(jpatch + 1) = temp_age + exchanged = 1 + !write(fates_log(),*) 'See patch 1 age :', currentPatch%age_since_anthro_disturbance + !write(fates_log(),*) 'See patch 2 age :', currentPatch%younger%age_since_anthro_disturbance + !write(fates_log(),*) 'See patch 1 order sec:', currentPatch%order_age_since_anthro + !write(fates_log(),*) 'See patch 2 order sec:', currentPatch%younger%order_age_since_anthro end if end do end if end do + write(fates_log(),*) 'After: See patch order:', order_of_patches(:) + write(fates_log(),*) 'After: See patch age:', sec_age_patches(:) + write(fates_log(),*) 'nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn' + + ! Fourth loop to assign order back to patches + do ipatch = 0, 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 + + ! Eliminate temperary variables + deallocate(order_of_patches) + deallocate(sec_age_patches) + end if ! logging_time ! Loop through all patches and calculate the harvest rate scale based From 664db78d3aa6abc8ffc979628d8b7b27221f9cf7 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Fri, 3 May 2024 22:55:31 -0700 Subject: [PATCH 08/16] Clear unnecessary write statements. --- biogeochem/EDPatchDynamicsMod.F90 | 188 +++++++++++++++--------------- 1 file changed, 92 insertions(+), 96 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 08d0cd98e5..29a4776ba6 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -239,7 +239,7 @@ subroutine disturbance_rates( site_in, bc_in) ! 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) @@ -253,6 +253,22 @@ subroutine disturbance_rates( site_in, bc_in) harvest_tag_csite = 2 remain_harvest_rate = fates_unset_r8 ! set to negative value to indicate uninitialized + ! The following code will calculate 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 + 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 + + ! We need to have an order based on age_since_anthro_disturbance if (logging_time) then @@ -281,17 +297,12 @@ subroutine disturbance_rates( site_in, bc_in) currentPatch => currentPatch%younger end do - !write(fates_log(),*) 'See patch order:', order_of_patches(:) - !write(fates_log(),*) 'See patch age:', sec_age_patches(:) ! Third loops to perform bubble sorting, since we have collected the order ! and age and stored them into array, no need to have a patch loop do ipatch = 1, npatches - 1 if(exchanged .eq. 1) then exchanged = 0 do jpatch = 1, npatches - ipatch - !write(fates_log(),*) 'ipatch:', ipatch - !write(fates_log(),*) 'jpatch:', jpatch - !write(fates_log(),*) 'nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn' if(sec_age_patches(jpatch) < sec_age_patches(jpatch+1)) then ! Swap order temp_order = order_of_patches(jpatch) @@ -302,21 +313,13 @@ subroutine disturbance_rates( site_in, bc_in) sec_age_patches(jpatch) = sec_age_patches(jpatch + 1) sec_age_patches(jpatch + 1) = temp_age exchanged = 1 - !write(fates_log(),*) 'See patch 1 age :', currentPatch%age_since_anthro_disturbance - !write(fates_log(),*) 'See patch 2 age :', currentPatch%younger%age_since_anthro_disturbance - !write(fates_log(),*) 'See patch 1 order sec:', currentPatch%order_age_since_anthro - !write(fates_log(),*) 'See patch 2 order sec:', currentPatch%younger%order_age_since_anthro end if end do end if end do - write(fates_log(),*) 'After: See patch order:', order_of_patches(:) - write(fates_log(),*) 'After: See patch age:', sec_age_patches(:) - write(fates_log(),*) 'nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn' - ! Fourth loop to assign order back to patches - do ipatch = 0, npatches + do ipatch = 1, npatches currentPatch => site_in%oldest_patch ! Look for the patch by usig patchno @@ -339,75 +342,82 @@ subroutine disturbance_rates( site_in, bc_in) ! Loop through all patches and calculate the harvest rate scale based ! on pre-defined age priority strategy - currentPatch => site_in%oldest_patch + ! Calculate the harvest rate scale based on pre-defined age priority strategy + ! Only calculate harvest rate scale when we have logging event + ! to reduce computational cost + if (logging_time .and. harvest_age_priority == 1) then + ! Maximize harvest rate (99%) for older patch first until harvest demand is fullfilled + ! ipatch is the rank of the secondary patch age + target_loop: do ipatch = 1, npatches - do while (associated(currentPatch)) - ! Initialize harvest rate scale to 1.0 first - currentPatch%harvest_rate_scale = 1._r8 + currentPatch => site_in%oldest_patch + search_loop: do while (associated(currentPatch)) + + found_patch: if(currentPatch%order_age_since_anthro .eq. ipatch) then + + ! Right now we skip the calculation for the primary land + 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 harvest rate 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 + ! Check if is the first time, if so initialize remain_harvest_rate + if(remain_harvest_rate(h_index) < 0) then + remain_harvest_rate(h_index) = harvest_rate + end if + ! Only calculate harvest rate scale for non-zero harvest rate + if (harvest_rate > 1e-7) then + ! 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)+1e-7) / harvest_rate + remain_harvest_rate(h_index) = 0._r8 + else + ! harvest almost 100%, leave 1% to prevent removing the patch completely, which may cause + ! model crash under nocomp mode + currentPatch%harvest_rate_scale = 0.99_r8 / harvest_rate + remain_harvest_rate(h_index) = remain_harvest_rate(h_index) - 0.99_r8 * (currentPatch%area/(AREA*frac_site)) + end if + end if + write(fates_log(),*) 'See patch number:', currentPatch%patchno + write(fates_log(),*) 'See patch age:', currentPatch%age + 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 harvest rate:', harvest_rate + 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 ================' - ! Only calculate harvest rate scale when we have logging event - ! to reduce computational cost - if (logging_time) then - ! Maximize harvest rate (99%) for older patch first until harvest demand is fullfilled - ! Right now we skip the calculation for the primary land - if(harvest_age_priority == 1 .and. currentPatch%land_use_label .eq. secondaryland) then - if(currentPatch%land_use_label .eq. primaryland) then - h_index = 1 - frac_site = frac_site_primary - else if (currentPatch%land_use_label .eq. secondaryland .and. & - 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 harvest rate 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 - ! Check if is the first time, if so initialize remain_harvest_rate - if(remain_harvest_rate(h_index) < 0) then - remain_harvest_rate(h_index) = harvest_rate - end if - ! Only calculate harvest rate scale for non-zero harvest rate - if (harvest_rate > 1e-7) then - ! 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)+1e-7) / harvest_rate - remain_harvest_rate(h_index) = 0._r8 - else - ! harvest almost 100%, leave 1% to prevent removing the patch completely, which may cause - ! model crash under nocomp mode - currentPatch%harvest_rate_scale = 0.99_r8 / harvest_rate - remain_harvest_rate(h_index) = remain_harvest_rate(h_index) - 0.99_r8 * (currentPatch%area/(AREA*frac_site)) - end if - end if - end if ! harvest age priority strategy - write(fates_log(),*) 'See patch number:', currentPatch%patchno - write(fates_log(),*) 'See patch age:', currentPatch%age - 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 harvest rate:', harvest_rate - 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 ! logging_time + end if if_secondary - currentPatch => currentPatch%younger - end do - write(fates_log(),*) '<<<<<<<<<<<<<<<<<<<<<<< Next section >>>>>>>>>>>>>>>>>>>>>>>>>>' + ! 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 !Loop through cohorts to get disturbance rates currentPatch => site_in%oldest_patch @@ -445,14 +455,7 @@ subroutine disturbance_rates( site_in, bc_in) harvestable_forest_c, & currentPatch%harvest_rate_scale, & harvest_tag) - if (logging_time) then - write(fates_log(),*) 'before logging, See patch number:', currentPatch%patchno - if(lmort_direct > 0._r8) write(fates_log(),*) 'lmort_direct:', lmort_direct - if(lmort_collateral > 0._r8) write(fates_log(),*) 'lmort_collateral:', lmort_collateral - if(lmort_infra > 0._r8) write(fates_log(),*) 'lmort_infra:', lmort_infra - if(l_degrad > 0._r8) write(fates_log(),*) 'l_degrad:', l_degrad - write(fates_log(),*) '================== Next Cohort ==================' - end if + currentCohort%lmort_direct = lmort_direct currentCohort%lmort_collateral = lmort_collateral currentCohort%lmort_infra = lmort_infra @@ -616,13 +619,6 @@ subroutine disturbance_rates( site_in, bc_in) end do endif - if (logging_time) then - write(fates_log(),*) 'See patch number:', currentPatch%patchno - write(fates_log(),*) 'See harvest rate scale:', currentPatch%harvest_rate_scale - write(fates_log(),*) 'See logging disturbance rate:', currentPatch%disturbance_rates(dtype_ilog) - write(fates_log(),*) 'See area:', currentPatch%area/AREA - write(fates_log(),*) '==================== Next patch ================' - end if currentPatch => currentPatch%younger enddo !patch loop From d8f0eb12435df8ca5b4cdf1a219296e55066e3e9 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Mon, 12 Aug 2024 15:06:17 -0700 Subject: [PATCH 09/16] Changes on calculating harvest rate for C-based harvest with maximizing age strategy. --- biogeochem/EDLoggingMortalityMod.F90 | 96 +++++++++++++++++++++++++--- biogeochem/EDPatchDynamicsMod.F90 | 94 +++++++++++++++++++-------- main/EDMainMod.F90 | 2 +- main/FatesConstantsMod.F90 | 2 +- main/FatesHistoryInterfaceMod.F90 | 8 +-- 5 files changed, 159 insertions(+), 43 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 7595de7025..398bd5e414 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -97,6 +97,7 @@ module EDLoggingMortalityMod public :: IsItLoggingTime public :: get_harvest_rate_area public :: get_harvestable_carbon + public :: get_harvestable_patch_carbon public :: get_harvest_rate_carbon public :: get_harvest_debt public :: UpdateHarvestC @@ -218,7 +219,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & 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 harvest rate proportionally + ! the harvest rate proportionally considering patch age 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 @@ -236,6 +237,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 ! todo: probably lower the dbhmin default value to 30 cm ! todo: change the default logging_event_code to 1 september (-244) @@ -276,9 +278,9 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & call get_harvest_rate_area (patch_land_use_label, hlm_harvest_catnames, & hlm_harvest_rates, frac_site_primary, frac_site_secondary_mature, secondary_age, harvest_rate) - ! 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 @@ -294,6 +296,15 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & endif + ! Logistic curve to account for preserving young trees + !harvest_rate_scale_cohort = harvest_rate_scale + ! Logistic + !harvest_rate_scale_cohort = harvest_rate_scale * 1._r8/(1._r8 + exp(-0.1*(dbh-(logging_dbhmin+logging_dbhmax)/2._r8))) + ! Inverse logistic + !harvest_rate_scale_cohort = harvest_rate_scale * 1._r8/(1._r8 + exp(0.1*(dbh-(logging_dbhmin+logging_dbhmax)/2._r8))) + ! Gaussian + harvest_rate_scale_cohort = harvest_rate_scale * 3._r8 * exp(-(dbh - (logging_dbhmin+logging_dbhmax)/2._r8) ** 2._r8 / (2._r8 * 5._r8 ** 2._r8)) + ! 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 @@ -304,7 +315,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! 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_scale * harvest_rate * logging_direct_frac + lmort_direct = harvest_rate_scale_cohort * harvest_rate * logging_direct_frac else lmort_direct = 0.0_r8 end if @@ -316,13 +327,13 @@ 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_scale * harvest_rate * logging_mechanical_frac + lmort_infra = harvest_rate_scale_cohort * harvest_rate * logging_mechanical_frac 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_scale * harvest_rate * logging_collateral_frac + lmort_collateral = harvest_rate_scale_cohort * harvest_rate * logging_collateral_frac else lmort_collateral = 0._r8 endif @@ -330,7 +341,7 @@ 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_scale * harvest_rate * logging_mechanical_frac + lmort_infra = harvest_rate_scale_cohort * harvest_rate * logging_mechanical_frac end if ! the area occupied by all plants in the canopy that aren't killed is still disturbed at the harvest rate @@ -476,7 +487,7 @@ subroutine get_harvestable_carbon (csite, site_area, hlm_harvest_catnames, harve ! Arguments type(ed_site_type), intent(in), target :: csite - real(r8), intent(in) :: site_area ! temporary variable + 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) @@ -504,7 +515,7 @@ subroutine get_harvestable_carbon (csite, site_area, hlm_harvest_catnames, harve 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 + 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 @@ -513,7 +524,7 @@ subroutine get_harvestable_carbon (csite, site_area, hlm_harvest_catnames, harve 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 + currentCohort%n! * AREA_INV * site_area ! No harvest for trees without canopy if (currentCohort%canopy_layer>=1) then @@ -558,6 +569,70 @@ end subroutine get_harvestable_carbon ! ============================================================================ + subroutine get_harvestable_patch_carbon ( currentPatch, harvestable_patch_c ) + + !USES: + use SFParamsMod, only : SF_val_cwd_frac + use EDTypesMod, only : AREA_INV + + + ! ------------------------------------------------------------------------------------------- + ! + ! 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_harvest_rate_carbon (patch_land_use_label, hlm_harvest_catnames, & hlm_harvest_rates, secondary_age, harvestable_forest_c, & harvest_rate, harvest_tag, cur_harvest_tag) @@ -1097,6 +1172,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) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 29a4776ba6..d356f677ad 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -69,7 +69,7 @@ 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_harvestable_carbon, get_harvestable_patch_carbon use EDLoggingMortalityMod, only : get_harvest_debt use EDParamsMod , only : fates_mortality_disturbance_fraction use FatesAllometryMod , only : carea_allom @@ -214,6 +214,7 @@ subroutine disturbance_rates( site_in, bc_in) 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 @@ -224,16 +225,15 @@ subroutine disturbance_rates( site_in, bc_in) integer, parameter :: harvest_age_priority = 1 ! Tag to determine if we give priority to the oldest ! patch ! 0 - no; - ! 1 - yes, maximize the oldest; - ! 2 - yes, sigmoid (To be added). and more other? + ! 1 - yes, maximize the oldest. integer :: exchanged ! used in bubble sorting integer :: temp_order ! used in bubble sorting real(r8) :: temp_age ! used in bubble sorting + real(r8) :: temp_totc ! for debug test integer :: ipatch, jpatch ! used in bubble sorting integer, allocatable :: order_of_patches(:) ! Record a copy of patch order for bubble sort real(r8), allocatable :: sec_age_patches(:) ! Record a copy of secondary forest age for bubble sort - !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) ! And the same rates in understory plants have already been applied to %dndt @@ -243,7 +243,7 @@ subroutine disturbance_rates( site_in, bc_in) 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) ! Initialize local variables @@ -273,7 +273,7 @@ subroutine disturbance_rates( site_in, bc_in) if (logging_time) then ! First loop to count number of patches and initialize order - ! with patchno + ! using patchno currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -283,7 +283,7 @@ subroutine disturbance_rates( site_in, bc_in) currentPatch => currentPatch%younger end do - ! Second loop to assign order patches and age_since_anthro_disturbance + ! Second loop to assign order of patches and age_since_anthro_disturbance allocate(order_of_patches(1:npatches)) allocate(sec_age_patches(1:npatches)) @@ -364,41 +364,80 @@ subroutine disturbance_rates( site_in, bc_in) h_index = 4 frac_site = 1._r8 - frac_site_primary - frac_site_secondary_mature end if - ! Obtain harvest rate of the corresponding patch + ! 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) + 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 - ! Check if is the first time, if so initialize remain_harvest_rate + ! For the first time, initialize remain_harvest_rate if(remain_harvest_rate(h_index) < 0) then - remain_harvest_rate(h_index) = harvest_rate + 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 > 1e-7) then - ! 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)+1e-7) / harvest_rate - remain_harvest_rate(h_index) = 0._r8 + 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 + 1e-7) / harvest_rate + remain_harvest_rate(h_index) = 0._r8 + else + ! harvest almost 100%, leave 1% to prevent removing the patch completely, which may cause + ! model crashing under nocomp mode + currentPatch%harvest_rate_scale = 0.99_r8 / harvest_rate + remain_harvest_rate(h_index) = remain_harvest_rate(h_index) - 0.99_r8 * harvestable_patch_c + end if else - ! harvest almost 100%, leave 1% to prevent removing the patch completely, which may cause - ! model crash under nocomp mode - currentPatch%harvest_rate_scale = 0.99_r8 / harvest_rate - remain_harvest_rate(h_index) = remain_harvest_rate(h_index) - 0.99_r8 * (currentPatch%area/(AREA*frac_site)) + ! 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)+1e-7) / harvest_rate + remain_harvest_rate(h_index) = 0._r8 + else + ! leave 1% to prevent removing the patch completely + currentPatch%harvest_rate_scale = 0.99_r8 / harvest_rate + remain_harvest_rate(h_index) = remain_harvest_rate(h_index) - 0.99_r8 * (currentPatch%area/(AREA*frac_site)) + end if end if end if + + ! Shijie: For diagnosing the harvestable carbon + 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 @@ -2623,7 +2662,7 @@ subroutine fuse_patches( csite, bc_in ) real(r8) :: primary_land_fraction_beforefusion,primary_land_fraction_afterfusion integer :: pftlabelmin, pftlabelmax - integer, parameter :: merge_strategy = 0 ! 0 - biomass profile based only + 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. @@ -2632,7 +2671,7 @@ subroutine fuse_patches( csite, bc_in ) currentSite => csite profiletol = ED_val_patch_fusion_tol - agetol = 5._r8 + agetol = 0.1_r8 primary_land_fraction_beforefusion = 0._r8 primary_land_fraction_afterfusion = 0._r8 @@ -2774,7 +2813,8 @@ 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. @@ -2784,7 +2824,7 @@ subroutine fuse_patches( csite, bc_in ) 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 + tpp%age_since_anthro_disturbance) >= agetol) then fuse_flag = 0 !do not fuse - keep apart. @@ -2826,7 +2866,7 @@ subroutine fuse_patches( csite, bc_in ) !------------------------------------------------------------------------! profiletol = ED_val_patch_fusion_tol - agetol = 5._r8 + agetol = 0.1_r8 endif fuseflagset_if endif different_patches_if @@ -2862,7 +2902,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 = 0._r8 + agetol = agetol * patch_fusion_tolerance_relaxation_increment !(1._r8/patch_fusion_tolerance_relaxation_increment) !---------------------------------------------------------------------! ! Making profile tolerance larger means that more fusion will happen ! diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 66613a0bec..fbe6821626 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -933,7 +933,7 @@ subroutine TotalBalanceCheck (currentSite, call_index ) error_frac = 0.0_r8 end if - if ( error_frac > 10e-6_r8 .or. (error /= error) ) then + 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/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index c9a0920a53..c70e8e1720 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -125,7 +125,7 @@ module FatesConstantsMod integer, public :: fates_np_comp_scaling = fates_unset_int - real(fates_r8), parameter, public :: secondary_age_threshold = 10._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 diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 04b96ff8cc..401fca3577 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -2767,7 +2767,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 @@ -6620,7 +6620,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) @@ -6981,7 +6981,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) @@ -7131,7 +7131,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) From c14dc2f4f599351dd3d757f0312a15f42fa7bb2b Mon Sep 17 00:00:00 2001 From: sshu88 Date: Wed, 23 Apr 2025 00:03:45 -0700 Subject: [PATCH 10/16] Added IFM options in Logging Module. --- biogeochem/EDLoggingMortalityMod.F90 | 61 +++++++++++++++++------- biogeochem/EDMortalityFunctionsMod.F90 | 4 +- biogeochem/EDPatchDynamicsMod.F90 | 64 +++++++++++++++++++++++++- biogeochem/FatesPatchMod.F90 | 2 +- main/EDParamsMod.F90 | 13 +++++- main/EDTypesMod.F90 | 2 + main/FatesConstantsMod.F90 | 8 ++++ main/FatesHistoryInterfaceMod.F90 | 8 ++-- main/FatesParametersInterface.F90 | 2 +- 9 files changed, 138 insertions(+), 26 deletions(-) mode change 100644 => 100755 biogeochem/EDLoggingMortalityMod.F90 mode change 100644 => 100755 biogeochem/EDMortalityFunctionsMod.F90 mode change 100644 => 100755 biogeochem/EDPatchDynamicsMod.F90 mode change 100644 => 100755 biogeochem/FatesPatchMod.F90 mode change 100644 => 100755 main/EDParamsMod.F90 mode change 100644 => 100755 main/EDTypesMod.F90 mode change 100644 => 100755 main/FatesConstantsMod.F90 mode change 100644 => 100755 main/FatesHistoryInterfaceMod.F90 mode change 100644 => 100755 main/FatesParametersInterface.F90 diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 old mode 100644 new mode 100755 index 398bd5e414..c815d0a698 --- 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 @@ -202,8 +208,8 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & hlm_harvest_units, & patch_land_use_label, secondary_age, & frac_site_primary, frac_site_secondary_mature, & - harvestable_forest_c, & - harvest_rate_scale, harvest_tag) + harvestable_forest_c, harvest_rate_scale, harvest_tag, & + harvest_wp_scale, logging_preference_options ) ! Arguments integer, intent(in) :: pft_i ! pft index @@ -219,7 +225,11 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & 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 harvest rate proportionally considering patch age + ! 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 @@ -237,7 +247,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 + 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) @@ -296,14 +306,23 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & endif - ! Logistic curve to account for preserving young trees - !harvest_rate_scale_cohort = harvest_rate_scale - ! Logistic - !harvest_rate_scale_cohort = harvest_rate_scale * 1._r8/(1._r8 + exp(-0.1*(dbh-(logging_dbhmin+logging_dbhmax)/2._r8))) - ! Inverse logistic - !harvest_rate_scale_cohort = harvest_rate_scale * 1._r8/(1._r8 + exp(0.1*(dbh-(logging_dbhmin+logging_dbhmax)/2._r8))) - ! Gaussian - harvest_rate_scale_cohort = harvest_rate_scale * 3._r8 * exp(-(dbh - (logging_dbhmin+logging_dbhmax)/2._r8) ** 2._r8 / (2._r8 * 5._r8 ** 2._r8)) + ! Logging size prefrence or change rotation length + if (logging_preference_options .eq. logging_uniform_size) then + harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale + else if(logging_preference_options .eq. logging_double_rotation) then + harvest_rate_scale_cohort = harvest_wp_scale * 0.5_r8 * harvest_rate_scale + else if(logging_preference_options .eq. logging_quadruple_rotation) then + harvest_rate_scale_cohort = harvest_wp_scale * 0.25_r8 * harvest_rate_scale + else if(logging_preference_options .eq. logging_logistic_size) then + harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale * & + 1._r8/(1._r8 + exp(-0.1*(dbh-(logging_dbhmin+logging_dbhmax)/2._r8))) + else if(logging_preference_options .eq. logging_inv_logistic_size) then + harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale * & + 1._r8/(1._r8 + exp(0.1*(dbh-(logging_dbhmin+logging_dbhmax)/2._r8))) + else if(logging_preference_options .eq. logging_gaussian_size) then + harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale * & + (1._r8/(5._r8*2.5_r8)) * exp(-(dbh - (logging_dbhmin+logging_dbhmax)/2._r8) ** 2._r8 / (2._r8 * 5._r8 ** 2._r8)) + end if ! 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 @@ -315,7 +334,12 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! 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_scale_cohort * 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 @@ -327,13 +351,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_scale_cohort * 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_scale_cohort * 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 @@ -341,7 +367,8 @@ 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_scale_cohort * 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 diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 old mode 100644 new mode 100755 index 57b133f543..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 @@ -298,7 +299,8 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, btran_ft, & age_since_anthro_disturbance, & frac_site_primary, frac_site_secondary_mature, & harvestable_forest_c, & - harvest_rate_scale, harvest_tag) + 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 d356f677ad..b385ca85c7 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -86,6 +86,8 @@ module EDPatchDynamicsMod use FatesConstantsMod , only : fates_unset_r8 use FatesConstantsMod , only : fates_unset_int use FatesConstantsMod , only : hlm_harvest_carbon + 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 @@ -102,6 +104,7 @@ 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_preference_options 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 @@ -202,6 +205,9 @@ subroutine disturbance_rates( site_in, bc_in) real(r8) :: l_degrad ! fraction of trees that are not killed but suffer from forest ! degradation (i.e. they are moved to newly-anthro-disturbed ! secondary forest patch) + real(r8) :: target_wp ! For size-dependent IFM: target harvest wood product for adjustment ratio calculation + real(r8) :: actual_wp ! For size-dependent IFM: actual harvest wood product after considering size priority + ! for adjustment ratio calculation real(r8) :: dist_rate_ldist_notharvested integer :: threshold_sizeclass integer :: i_dist @@ -226,6 +232,8 @@ subroutine disturbance_rates( site_in, bc_in) ! patch ! 0 - no; ! 1 - yes, maximize the oldest. + real(r8), parameter :: target_harvest_wp_scale = 0.5 ! For Size-dependent IFM. The scale is a target fraction of wood product + ! expected afer considering size-priority of IFM reduce the total harvest amount integer :: exchanged ! used in bubble sorting integer :: temp_order ! used in bubble sorting real(r8) :: temp_age ! used in bubble sorting @@ -251,7 +259,10 @@ subroutine disturbance_rates( site_in, bc_in) npatches = 0 frac_site = 1._r8 harvest_tag_csite = 2 + target_wp = 0._r8 + actual_wp = 0._r8 remain_harvest_rate = fates_unset_r8 ! set to negative value to indicate uninitialized + site_in%resources_management%harvest_wp_scale = 1._r8 ! The following code will calculate harvest_rate_scale for each secondary patch ! There're currently 2 options: @@ -458,6 +469,57 @@ subroutine disturbance_rates( site_in, bc_in) end if ! logging_time + !Loop through cohorts to preview the harvest product assuming no harvest size priority and with size priority + !Then calculate the ratio between these two cases to calculate adjustment ratio + if (logging_time .and. logging_preference_options >= logging_logistic_size) then + !Revise logging rate scale to match the target harvest demand after + !considering the logging size preference + 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, target_harvest_wp_scale, logging_uniform_size) + target_wp = target_wp + lmort_direct + + ! 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) + actual_wp = actual_wp + lmort_direct + + currentCohort => currentCohort%taller + end do + + currentPatch => currentPatch%younger + end do + ! adjustment ratio + 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 if ! logging_time + !Loop through cohorts to get disturbance rates currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -493,7 +555,7 @@ subroutine disturbance_rates( site_in, bc_in) frac_site_secondary_mature, & harvestable_forest_c, & currentPatch%harvest_rate_scale, & - harvest_tag) + harvest_tag, site_in%resources_management%harvest_wp_scale, logging_preference_options) currentCohort%lmort_direct = lmort_direct currentCohort%lmort_collateral = lmort_collateral diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 old mode 100644 new mode 100755 index 919d619681..0eb2b4d6a1 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -179,7 +179,7 @@ module FatesPatchMod 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 + ! 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 diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 old mode 100644 new mode 100755 index beaf11e140..513c20a59a --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -235,9 +235,11 @@ 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_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_dbhmin ! Minimum dbh at which logging is applied (cm) ! Typically associated with harvesting @@ -340,6 +342,7 @@ subroutine FatesParamsInit() hydr_psicap = nan hydr_solver = -9 bgc_soil_salinity = nan + logging_preference_options = -9 logging_dbhmin = nan logging_dbhmax = nan logging_collateral_frac = nan @@ -510,6 +513,9 @@ 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_preference_options, 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 +735,10 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetrieveParameter(name=bgc_name_soil_salinity, & data=bgc_soil_salinity) + call fates_params%RetrieveParameter(name=logging_name_preference_options, & + data=tmpreal) + logging_preference_options = nint(tmpreal) + call fates_params%RetrieveParameter(name=logging_name_dbhmin, & data=logging_dbhmin) @@ -874,6 +884,7 @@ 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_preference_options = ',logging_preference_options 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/EDTypesMod.F90 b/main/EDTypesMod.F90 old mode 100644 new mode 100755 index 16570c1a68..bcf11c987d --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -138,6 +138,8 @@ module EDTypesMod ! 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 diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 old mode 100644 new mode 100755 index c70e8e1720..cfc58fdbce --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -143,6 +143,14 @@ module FatesConstantsMod integer, parameter, public :: lmrmodel_ryan_1991 = 1 integer, parameter, public :: lmrmodel_atkin_etal_2017 = 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 diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 old mode 100644 new mode 100755 index 401fca3577..ad3e5c16eb --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -7472,14 +7472,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) @@ -7795,7 +7795,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) @@ -7956,7 +7956,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) 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 From 3bd23788acd9c22fdfb4e3ceb61a122158c1deb4 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Wed, 3 Sep 2025 13:45:33 -0700 Subject: [PATCH 11/16] Add function to obtain carbon from each cohort, add harvested carbon per size bin as history variable --- biogeochem/EDCohortDynamicsMod.F90 | 2 + biogeochem/EDLoggingMortalityMod.F90 | 52 +++++++++- biogeochem/EDPatchDynamicsMod.F90 | 136 ++++++++++++++++----------- biogeochem/EDPhysiologyMod.F90 | 24 +++-- biogeochem/FatesCohortMod.F90 | 5 + main/EDInitMod.F90 | 5 + main/EDMainMod.F90 | 7 +- main/EDTypesMod.F90 | 4 + main/FatesConstantsMod.F90 | 6 ++ main/FatesHistoryInterfaceMod.F90 | 65 +++++++++---- 10 files changed, 220 insertions(+), 86 deletions(-) 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 index c815d0a698..94c3870f11 100755 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -104,6 +104,7 @@ module EDLoggingMortalityMod public :: get_harvest_rate_area public :: get_harvestable_carbon public :: get_harvestable_patch_carbon + public :: get_harvested_cohort_carbon public :: get_harvest_rate_carbon public :: get_harvest_debt public :: UpdateHarvestC @@ -553,7 +554,7 @@ subroutine get_harvestable_carbon (csite, site_area, hlm_harvest_catnames, harve SF_val_CWD_frac(ncwd) * logging_export_frac * & currentCohort%n! * AREA_INV * site_area - ! No harvest for trees without canopy + ! 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. & @@ -660,6 +661,43 @@ 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 + + ! ============================================================================ + subroutine get_harvest_rate_carbon (patch_land_use_label, hlm_harvest_catnames, & hlm_harvest_rates, secondary_age, harvestable_forest_c, & harvest_rate, harvest_tag, cur_harvest_tag) @@ -912,6 +950,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) @@ -996,7 +1035,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 @@ -1108,6 +1147,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 diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index b385ca85c7..ae7efbbb59 100755 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -70,7 +70,7 @@ module EDPatchDynamicsMod use EDLoggingMortalityMod, only : get_harvest_rate_area use EDLoggingMortalityMod, only : get_harvest_rate_carbon use EDLoggingMortalityMod, only : get_harvestable_carbon, get_harvestable_patch_carbon - use EDLoggingMortalityMod, only : get_harvest_debt + use EDLoggingMortalityMod, only : get_harvested_cohort_carbon, get_harvest_debt use EDParamsMod , only : fates_mortality_disturbance_fraction use FatesAllometryMod , only : carea_allom use FatesAllometryMod , only : set_root_fraction @@ -82,6 +82,7 @@ module EDPatchDynamicsMod 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 @@ -234,13 +235,16 @@ subroutine disturbance_rates( site_in, bc_in) ! 1 - yes, maximize the oldest. real(r8), parameter :: target_harvest_wp_scale = 0.5 ! For Size-dependent IFM. The scale is a target fraction of wood product ! expected afer considering size-priority of IFM reduce the total harvest amount + real(r8), parameter :: max_iteration = 5 ! Maximum iterations to obtain target harvest wood product integer :: exchanged ! used in bubble sorting integer :: temp_order ! used in bubble sorting real(r8) :: temp_age ! used in bubble sorting real(r8) :: temp_totc ! for debug test + real(r8) :: harvested_cohort_c ! used in IFM integer :: ipatch, jpatch ! used in bubble sorting integer, allocatable :: order_of_patches(:) ! Record a copy of patch order for bubble sort real(r8), allocatable :: sec_age_patches(:) ! Record a copy of secondary forest age for bubble sort + integer :: it ! index of iteration !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) @@ -259,8 +263,6 @@ subroutine disturbance_rates( site_in, bc_in) npatches = 0 frac_site = 1._r8 harvest_tag_csite = 2 - target_wp = 0._r8 - actual_wp = 0._r8 remain_harvest_rate = fates_unset_r8 ! set to negative value to indicate uninitialized site_in%resources_management%harvest_wp_scale = 1._r8 @@ -357,7 +359,7 @@ subroutine disturbance_rates( site_in, bc_in) ! Only calculate harvest rate scale when we have logging event ! to reduce computational cost if (logging_time .and. harvest_age_priority == 1) then - ! Maximize harvest rate (99%) for older patch first until harvest demand is fullfilled + ! Maximize harvest rate (99.99%) for older patch first until harvest demand is fullfilled ! ipatch is the rank of the secondary patch age target_loop: do ipatch = 1, npatches @@ -396,31 +398,31 @@ subroutine disturbance_rates( site_in, bc_in) end if end if ! Only calculate harvest rate scale for non-zero harvest rate - if (harvest_rate > 1e-7) then + 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 + 1e-7) / harvest_rate + (harvestable_patch_c + fates_tiny) / harvest_rate remain_harvest_rate(h_index) = 0._r8 else - ! harvest almost 100%, leave 1% to prevent removing the patch completely, which may cause + ! harvest almost 100%, leave 0.01% to prevent removing the patch completely, which may cause ! model crashing under nocomp mode - currentPatch%harvest_rate_scale = 0.99_r8 / harvest_rate - remain_harvest_rate(h_index) = remain_harvest_rate(h_index) - 0.99_r8 * harvestable_patch_c + 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)+1e-7) / harvest_rate + (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 = 0.99_r8 / harvest_rate - remain_harvest_rate(h_index) = remain_harvest_rate(h_index) - 0.99_r8 * (currentPatch%area/(AREA*frac_site)) + 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 @@ -474,50 +476,66 @@ subroutine disturbance_rates( site_in, bc_in) if (logging_time .and. logging_preference_options >= logging_logistic_size) then !Revise logging rate scale to match the target harvest demand after !considering the logging size preference - currentPatch => site_in%oldest_patch - do while (associated(currentPatch)) + !Really need to add an iteration here to force the calculated harvested + !carbon match the forcing + !Let's test 3 times + iteration_loop: do it = 1, max_iteration - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) + target_wp = 0._r8 + actual_wp = 0._r8 - ! 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, target_harvest_wp_scale, logging_uniform_size) - target_wp = target_wp + lmort_direct - - ! 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) - actual_wp = actual_wp + lmort_direct + 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, target_harvest_wp_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 - currentCohort => currentCohort%taller + currentPatch => currentPatch%younger end do - - currentPatch => currentPatch%younger - end do - ! adjustment ratio - 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 + ! 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 !Loop through cohorts to get disturbance rates @@ -1107,6 +1125,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 @@ -1171,6 +1190,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 @@ -1195,9 +1215,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 @@ -1267,6 +1288,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 @@ -1365,6 +1387,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 @@ -1431,6 +1454,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 @@ -1454,6 +1478,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 @@ -2671,6 +2696,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 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/main/EDInitMod.F90 b/main/EDInitMod.F90 index 02a6d8999a..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(:) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index fbe6821626..35b3d92ced 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 @@ -862,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 @@ -905,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) @@ -933,6 +935,9 @@ subroutine TotalBalanceCheck (currentSite, call_index ) error_frac = 0.0_r8 end if + ! 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) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index bcf11c987d..6243c28fde 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -206,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] @@ -218,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 @@ -492,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 index cfc58fdbce..0c68df74c4 100755 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -183,6 +183,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 index ad3e5c16eb..fd30b5b104 100755 --- 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 @@ -468,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 @@ -2312,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 @@ -2394,6 +2401,7 @@ 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_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, & @@ -3331,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 @@ -3512,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 + & @@ -3662,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 + & @@ -3969,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) + & @@ -4449,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) + & @@ -5946,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) @@ -6045,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) @@ -6459,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', & @@ -6504,10 +6530,10 @@ 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 yr-1', & long='Accumulated carbon failed to be harvested', use_default='active', & @@ -6998,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) @@ -7032,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) @@ -7375,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) @@ -8195,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) @@ -8267,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 @@ -8384,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) From 01c322d47f19eb7a2f065a0821bf0a9582f34be3 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Mon, 8 Sep 2025 14:56:35 -0700 Subject: [PATCH 12/16] replace bubble sort with sort algorithm using pseudo-age --- biogeochem/EDPatchDynamicsMod.F90 | 404 ++++++++++++++++-------------- 1 file changed, 209 insertions(+), 195 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index ae7efbbb59..c31a7adb4e 100755 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -190,6 +190,7 @@ subroutine disturbance_rates( site_in, bc_in) ! ! !LOCAL VARIABLES: type (fates_patch_type) , pointer :: currentPatch + type (fates_patch_type) , pointer :: referredPatch type (fates_cohort_type), pointer :: currentCohort real(r8) :: cmort @@ -228,7 +229,7 @@ subroutine disturbance_rates( site_in, bc_in) 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] - ! Control vars, will be moved into parameter file once the test is done + ! Control vars, need to be moved into parameter file once the test is done integer, parameter :: harvest_age_priority = 1 ! Tag to determine if we give priority to the oldest ! patch ! 0 - no; @@ -236,12 +237,11 @@ subroutine disturbance_rates( site_in, bc_in) real(r8), parameter :: target_harvest_wp_scale = 0.5 ! For Size-dependent IFM. The scale is a target fraction of wood product ! expected afer considering size-priority of IFM reduce the total harvest amount real(r8), parameter :: max_iteration = 5 ! Maximum iterations to obtain target harvest wood product - integer :: exchanged ! used in bubble sorting - integer :: temp_order ! used in bubble sorting - real(r8) :: temp_age ! used in bubble sorting real(r8) :: temp_totc ! for debug test real(r8) :: harvested_cohort_c ! used in IFM - integer :: ipatch, jpatch ! used in bubble sorting + integer :: ipatch, , ipft, curpft ! used in min-heap merge + integer, allocatable :: jpatch(:), pftlen(:) ! used in min-heap merge + real(r8), allocatable :: minheap(:) ! array recording minimum heap integer, allocatable :: order_of_patches(:) ! Record a copy of patch order for bubble sort real(r8), allocatable :: sec_age_patches(:) ! Record a copy of secondary forest age for bubble sort integer :: it ! index of iteration @@ -250,6 +250,12 @@ subroutine disturbance_rates( site_in, bc_in) ! Calculate Mortality Rates (these were previously calculated during growth derivatives) ! And the same rates in understory plants have already been applied to %dndt !---------------------------------------------------------------------------------------------- + ! Sep 5. 2025 + ! New: improved forest management practies can be quantified as modifiers to change the + ! logging disturbance rates at cohort-patch-site levels: + ! site_in%resources_management%harvest_wp_scale + ! + !---------------------------------------------------------------------------------------------- ! 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) @@ -258,7 +264,7 @@ subroutine disturbance_rates( site_in, bc_in) ! 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) - ! Initialize local variables + ! Initialize the rest of local variables exchanged = 1 npatches = 0 frac_site = 1._r8 @@ -266,6 +272,9 @@ subroutine disturbance_rates( site_in, bc_in) remain_harvest_rate = fates_unset_r8 ! set to negative value to indicate uninitialized site_in%resources_management%harvest_wp_scale = 1._r8 + ! ===================================================================== + ! Age_prioritized forest management + ! ===================================================================== ! The following code will calculate harvest_rate_scale for each secondary patch ! There're currently 2 options: ! 1) uniformly harvest each patch (harvest_rate_scale = 1._r8) @@ -282,12 +291,25 @@ subroutine disturbance_rates( site_in, bc_in) end do - ! We need to have an order based on age_since_anthro_disturbance + ! ===================================================================== + ! "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 actually modify the patch sequence but use another age tag called 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 last 4 digits of pseudo-age, i.e., actual patch age. Since each PFT + ! the patch sequence is already sorted, here we need to loop through PFTs and merge these sorted sequences into one + ! sorted array using min-heap merge if (logging_time) then - ! First loop to count number of patches and initialize order + ! First loop to count total number of patches and initialize order ! using patchno currentPatch => site_in%oldest_patch + npatches = 0 do while (associated(currentPatch)) npatches = npatches + 1 @@ -296,60 +318,89 @@ subroutine disturbance_rates( site_in, bc_in) currentPatch => currentPatch%younger end do - ! Second loop to assign order of patches and age_since_anthro_disturbance - allocate(order_of_patches(1:npatches)) - allocate(sec_age_patches(1:npatches)) - - currentPatch => site_in%oldest_patch - ipatch = 0 - do while (associated(currentPatch)) - ipatch = ipatch + 1 - order_of_patches(ipatch) = currentPatch%patchno - sec_age_patches(ipatch) = currentPatch%age_since_anthro_disturbance + ! For case 2, i.e., nocomp mode only + if (hlm_use_nocomp.eq.itrue) then - currentPatch => currentPatch%younger - end do + allocate(order_of_patches(1:npatches)) + allocate(order_of_patches_per_pft(1:numpft,1:npatches)) + allocate(sec_age_patches(1:npatches)) + allocate(sec_age_patches_per_pft(1:numpft,1:npatches)) + allocate(minheap(1:numpft)) + allocate(jpatch(1:numpft)) + allocate(pftlen(1:numpft)) - ! Third loops to perform bubble sorting, since we have collected the order - ! and age and stored them into array, no need to have a patch loop - do ipatch = 1, npatches - 1 - if(exchanged .eq. 1) then - exchanged = 0 - do jpatch = 1, npatches - ipatch - if(sec_age_patches(jpatch) < sec_age_patches(jpatch+1)) then - ! Swap order - temp_order = order_of_patches(jpatch) - order_of_patches(jpatch) = order_of_patches(jpatch + 1) - order_of_patches(jpatch + 1) = temp_order - ! Swap age - temp_age = sec_age_patches(jpatch) - sec_age_patches(jpatch) = sec_age_patches(jpatch + 1) - sec_age_patches(jpatch + 1) = temp_age - exchanged = 1 + ! 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 + sec_age_patches(ipatch) = mod(GetPseudoPatchAge(currentPatch), 1.e4) + + do ipft=1,numpft + if(currentPatch%nocomp_pft_label .eq. ipft) + pftlen(ipft) = pftlen(ipft) + 1 + order_of_patches_per_pft(ipft,pftlen(ipft)) = currentPatch%order_age_since_anthro + sec_age_patches_per_pft(ipft, pftlen(ipft)) = mod(GetPseudoPatchAge(currentPatch), 1.e4) end if end do - end if - end do - ! Fourth loop to assign order back to patches - do ipatch = 1, npatches - currentPatch => site_in%oldest_patch + currentPatch => currentPatch%younger + end do - ! Look for the patch by usig patchno - inner_loop: do while (associated(currentPatch)) + ! Third loop to initailize min-heap, first element of each array + do ipft=1,numpft + if(pftlen(ipft) > 0) then + minheap(ipft) = sec_age_patches_per_pft(ipft, 0) + end if + end do - if(order_of_patches(ipatch) .eq. currentPatch%patchno) then - currentPatch%order_age_since_anthro = ipatch - exit inner_loop + ! Fourth loop to get min from min-heap then insert 'sec_age_patches_per_pft' and record order in 'order_of_patches' + ! no need to loop over data list here + ipatch = 0 + jpatch(:) = 1 + do ipatch = 1, npacthes + ! Find the minimum and assign the order + curpft = minloc(minheap) + 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) = sec_age_patches_per_pft(ipft, jpatch(curpft)) + else + minheap(curpft) = 0._r8 + end if end if + end do - currentPatch => currentPatch%younger - end do inner_loop - end do + ! Fifthth loop to assign order back to patches + do ipatch = 1, npatches + currentPatch => site_in%oldest_patch - ! Eliminate temperary variables - deallocate(order_of_patches) - deallocate(sec_age_patches) + ! 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(sec_age_patches) + deallocate(sec_age_patches_per_pft) + deallocate(minheap) + deallocate(jpatch) + deallocate(pftlen) + + end if ! nocomp end if ! logging_time @@ -368,7 +419,8 @@ subroutine disturbance_rates( site_in, bc_in) found_patch: if(currentPatch%order_age_since_anthro .eq. ipatch) then - ! Right now we skip the calculation for the primary land + ! 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 @@ -427,7 +479,8 @@ subroutine disturbance_rates( site_in, bc_in) end if end if - ! Shijie: For diagnosing the harvestable carbon + ! For diagnosing the harvestable carbon + ! These codes can be removed in released version temp_totc = 0._r8 currentCohort => currentPatch%tallest @@ -471,14 +524,19 @@ subroutine disturbance_rates( site_in, bc_in) end if ! logging_time - !Loop through cohorts to preview the harvest product assuming no harvest size priority and with size priority - !Then calculate the ratio between these two cases to calculate adjustment ratio + !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 to actual cases 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. + + !Only run this part when considering size dependent priority if (logging_time .and. logging_preference_options >= logging_logistic_size) then - !Revise logging rate scale to match the target harvest demand after - !considering the logging size preference - !Really need to add an iteration here to force the calculated harvested - !carbon match the forcing - !Let's test 3 times + iteration_loop: do it = 1, max_iteration target_wp = 0._r8 @@ -530,12 +588,18 @@ subroutine disturbance_rates( site_in, bc_in) 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 + 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 !Loop through cohorts to get disturbance rates @@ -3610,155 +3674,105 @@ subroutine InsertPatch(currentSite, newPatch) ! !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 + logical :: patch_inserted - ! 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) - - ! Option 1 - order of land use label groups does not matter - case (unordered_lul_groups) + ! 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 - 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 + ! 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: - ! 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 + ! LU1/PFT1/O <-> LU1/PFT1/Y <-> LU1/PFT2/O <- ... -> LU3/PFT2/O <-> LU3/PFT2/Y - ! Option 2 - primaryland group must be on the oldest end - case (primaryland_oldest_group) + ! 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. - 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 + ! If land use is turned off or nocomp is turned off, then this should devolve to the prior + ! behavior of just age sorting. - ! 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 + patch_inserted = .false. + + if (GetPseudoPatchAge(newPatch) .le. GetPseudoPatchAge(currentSite%youngest_patch)) then - ! Option 3 - groups are numerically ordered with primaryland group starting at oldest end. - case (numerical_order_lul_groups) + ! 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 - ! 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 = .true. + else if (GetPseudoPatchAge(newPatch) .ge. GetPseudoPatchAge(currentSite%oldest_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 end of the linked list + newPatch%younger => currentSite%oldest_patch + newPatch%older => null() + currentSite%oldest_patch%older => newPatch + currentSite%oldest_patch => newPatch - newPatch%older => currentPatch - newPatch%younger => currentPatch%younger + 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 - - else - - ! 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 + 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 From cc69cda21fbaaea0e3b8f1ba38858a74f0a19247 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Mon, 8 Sep 2025 17:00:10 -0700 Subject: [PATCH 13/16] Migrate options from code to parameter file. --- biogeochem/EDPatchDynamicsMod.F90 | 21 +++++++++------------ main/EDParamsMod.F90 | 18 ++++++++++++++++++ main/FatesConstantsMod.F90 | 4 ++++ parameter_files/fates_params_default.cdl | 15 +++++++++++++++ 4 files changed, 46 insertions(+), 12 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index c31a7adb4e..604400fedc 100755 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -87,6 +87,7 @@ module EDPatchDynamicsMod 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 @@ -105,7 +106,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_preference_options + 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 @@ -230,21 +232,16 @@ subroutine disturbance_rates( site_in, bc_in) real(r8) :: current_fates_landuse_state_vector(n_landuse_cats) ! [m2/m2] ! Control vars, need to be moved into parameter file once the test is done - integer, parameter :: harvest_age_priority = 1 ! Tag to determine if we give priority to the oldest - ! patch - ! 0 - no; - ! 1 - yes, maximize the oldest. - real(r8), parameter :: target_harvest_wp_scale = 0.5 ! For Size-dependent IFM. The scale is a target fraction of wood product - ! expected afer considering size-priority of IFM reduce the total harvest amount - real(r8), parameter :: max_iteration = 5 ! Maximum iterations to obtain target harvest wood product - real(r8) :: temp_totc ! for debug test + real(r8), parameter :: max_iteration = 5 ! Maximum iterations to obtain target harvest wood product. Currently I don't + ! plan to move it into parameter files + real(r8) :: temp_totc ! for debug test, can be removed in released version real(r8) :: harvested_cohort_c ! used in IFM integer :: ipatch, , ipft, curpft ! used in min-heap merge integer, allocatable :: jpatch(:), pftlen(:) ! used in min-heap merge real(r8), allocatable :: minheap(:) ! array recording minimum heap integer, allocatable :: order_of_patches(:) ! Record a copy of patch order for bubble sort real(r8), allocatable :: sec_age_patches(:) ! Record a copy of secondary forest age for bubble sort - integer :: it ! index of iteration + integer :: it ! iteration index !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) @@ -409,7 +406,7 @@ subroutine disturbance_rates( site_in, bc_in) ! Calculate the harvest rate scale based on pre-defined age priority strategy ! Only calculate harvest rate scale when we have logging event ! to reduce computational cost - if (logging_time .and. harvest_age_priority == 1) then + if (logging_time .and. logging_age_preference == logging_oldest_first) then ! Maximize harvest rate (99.99%) for older patch first until harvest demand is fullfilled ! ipatch is the rank of the secondary patch age target_loop: do ipatch = 1, npatches @@ -560,7 +557,7 @@ subroutine disturbance_rates( site_in, bc_in) frac_site_secondary_mature, & harvestable_forest_c, & currentPatch%harvest_rate_scale, & - harvest_tag, target_harvest_wp_scale, logging_uniform_size) + harvest_tag, logging_ifm_harvest_scale, logging_uniform_size) call get_harvested_cohort_carbon(currentCohort, lmort_direct, harvested_cohort_c) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 513c20a59a..661c1bda67 100755 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -237,10 +237,17 @@ module EDParamsMod ! 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 character(len=param_string_length),parameter,public :: logging_name_dbhmin = "fates_landuse_logging_dbhmin" @@ -342,7 +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 @@ -735,10 +744,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) @@ -884,7 +900,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/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 0c68df74c4..90c8dcdb1b 100755 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -143,6 +143,10 @@ 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 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 ; From 2956797d0808337934084504e80aabc48e284d82 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Tue, 9 Sep 2025 10:53:07 -0700 Subject: [PATCH 14/16] Refactor code. Put all IFM related subroutines in FatesForestManagementMod --- biogeochem/EDLoggingMortalityMod.F90 | 232 +------- biogeochem/EDPatchDynamicsMod.F90 | 369 +----------- biogeochem/FatesForestManagementMod.F90 | 724 ++++++++++++++++++++++++ 3 files changed, 739 insertions(+), 586 deletions(-) create mode 100644 biogeochem/FatesForestManagementMod.F90 diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 94c3870f11..4b74547ac6 100755 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -63,6 +63,7 @@ module EDLoggingMortalityMod use PRTGenericMod , only : sapw_organ, struct_organ, leaf_organ use PRTGenericMod , only : fnrt_organ, store_organ, repro_organ use FatesAllometryMod , only : set_root_fraction + use FatesForestManagementMod, only : get_cohort_harvest_rate_scale use FatesConstantsMod , only : primaryland, secondaryland, secondary_age_threshold use FatesConstantsMod , only : fates_tiny use FatesConstantsMod , only : months_per_year, days_per_sec, years_per_day, g_per_kg @@ -102,9 +103,6 @@ module EDLoggingMortalityMod public :: logging_time public :: IsItLoggingTime public :: get_harvest_rate_area - public :: get_harvestable_carbon - public :: get_harvestable_patch_carbon - public :: get_harvested_cohort_carbon public :: get_harvest_rate_carbon public :: get_harvest_debt public :: UpdateHarvestC @@ -307,23 +305,9 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & endif + ! Logging size prefrence or change rotation length - if (logging_preference_options .eq. logging_uniform_size) then - harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale - else if(logging_preference_options .eq. logging_double_rotation) then - harvest_rate_scale_cohort = harvest_wp_scale * 0.5_r8 * harvest_rate_scale - else if(logging_preference_options .eq. logging_quadruple_rotation) then - harvest_rate_scale_cohort = harvest_wp_scale * 0.25_r8 * harvest_rate_scale - else if(logging_preference_options .eq. logging_logistic_size) then - harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale * & - 1._r8/(1._r8 + exp(-0.1*(dbh-(logging_dbhmin+logging_dbhmax)/2._r8))) - else if(logging_preference_options .eq. logging_inv_logistic_size) then - harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale * & - 1._r8/(1._r8 + exp(0.1*(dbh-(logging_dbhmin+logging_dbhmax)/2._r8))) - else if(logging_preference_options .eq. logging_gaussian_size) then - harvest_rate_scale_cohort = harvest_wp_scale * harvest_rate_scale * & - (1._r8/(5._r8*2.5_r8)) * exp(-(dbh - (logging_dbhmin+logging_dbhmax)/2._r8) ** 2._r8 / (2._r8 * 5._r8 ** 2._r8)) - end if + call get_cohort_harvest_rate_scale (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 @@ -340,7 +324,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! 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) + 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 @@ -388,7 +372,6 @@ 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, & @@ -491,211 +474,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 ! 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 - use EDTypesMod, only : AREA_INV - - - ! ------------------------------------------------------------------------------------------- - ! - ! 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 - ! ============================================================================ subroutine get_harvest_rate_carbon (patch_land_use_label, hlm_harvest_catnames, & @@ -1305,6 +1083,8 @@ subroutine UpdateHarvestC(currentSite,bc_out) end subroutine UpdateHarvestC + ! ===================================================================================== + subroutine get_harvest_debt(site_in, bc_in, harvest_tag) ! diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 604400fedc..45f3308968 100755 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -69,8 +69,7 @@ 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, get_harvestable_patch_carbon - use EDLoggingMortalityMod, only : get_harvested_cohort_carbon, get_harvest_debt + use EDLoggingMortalityMod, only : get_harvest_debt use EDParamsMod , only : fates_mortality_disturbance_fraction use FatesAllometryMod , only : carea_allom use FatesAllometryMod , only : set_root_fraction @@ -184,15 +183,16 @@ 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 ! ! !LOCAL VARIABLES: type (fates_patch_type) , pointer :: currentPatch - type (fates_patch_type) , pointer :: referredPatch type (fates_cohort_type), pointer :: currentCohort real(r8) :: cmort @@ -209,9 +209,6 @@ subroutine disturbance_rates( site_in, bc_in) real(r8) :: l_degrad ! fraction of trees that are not killed but suffer from forest ! degradation (i.e. they are moved to newly-anthro-disturbed ! secondary forest patch) - real(r8) :: target_wp ! For size-dependent IFM: target harvest wood product for adjustment ratio calculation - real(r8) :: actual_wp ! For size-dependent IFM: actual harvest wood product after considering size priority - ! for adjustment ratio calculation real(r8) :: dist_rate_ldist_notharvested integer :: threshold_sizeclass integer :: i_dist @@ -231,29 +228,13 @@ subroutine disturbance_rates( site_in, bc_in) 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] - ! Control vars, need to be moved into parameter file once the test is done - real(r8), parameter :: max_iteration = 5 ! Maximum iterations to obtain target harvest wood product. Currently I don't - ! plan to move it into parameter files - real(r8) :: temp_totc ! for debug test, can be removed in released version real(r8) :: harvested_cohort_c ! used in IFM - integer :: ipatch, , ipft, curpft ! used in min-heap merge - integer, allocatable :: jpatch(:), pftlen(:) ! used in min-heap merge - real(r8), allocatable :: minheap(:) ! array recording minimum heap - integer, allocatable :: order_of_patches(:) ! Record a copy of patch order for bubble sort - real(r8), allocatable :: sec_age_patches(:) ! Record a copy of secondary forest age for bubble sort - integer :: it ! iteration index !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) ! And the same rates in understory plants have already been applied to %dndt !---------------------------------------------------------------------------------------------- - ! Sep 5. 2025 - ! New: improved forest management practies can be quantified as modifiers to change the - ! logging disturbance rates at cohort-patch-site levels: - ! site_in%resources_management%harvest_wp_scale - ! - !---------------------------------------------------------------------------------------------- - + ! 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) @@ -261,343 +242,11 @@ subroutine disturbance_rates( site_in, bc_in) ! 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) - ! Initialize the rest of local variables - exchanged = 1 - npatches = 0 - frac_site = 1._r8 - harvest_tag_csite = 2 - remain_harvest_rate = fates_unset_r8 ! set to negative value to indicate uninitialized - site_in%resources_management%harvest_wp_scale = 1._r8 - - ! ===================================================================== - ! Age_prioritized forest management - ! ===================================================================== - ! The following code will calculate 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 - 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 actually modify the patch sequence but use another age tag called 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 last 4 digits of pseudo-age, i.e., actual patch age. Since each PFT - ! the patch sequence is already sorted, here we need to loop through PFTs and merge these sorted sequences into one - ! sorted array using min-heap merge - if (logging_time) then - - ! First loop to count total number of patches and initialize order - ! using patchno - 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 - - allocate(order_of_patches(1:npatches)) - allocate(order_of_patches_per_pft(1:numpft,1:npatches)) - allocate(sec_age_patches(1:npatches)) - allocate(sec_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 - sec_age_patches(ipatch) = mod(GetPseudoPatchAge(currentPatch), 1.e4) - - do ipft=1,numpft - if(currentPatch%nocomp_pft_label .eq. ipft) - pftlen(ipft) = pftlen(ipft) + 1 - order_of_patches_per_pft(ipft,pftlen(ipft)) = currentPatch%order_age_since_anthro - sec_age_patches_per_pft(ipft, pftlen(ipft)) = mod(GetPseudoPatchAge(currentPatch), 1.e4) - end if - end do - - currentPatch => currentPatch%younger - end do - - ! Third loop to initailize min-heap, first element of each array - do ipft=1,numpft - if(pftlen(ipft) > 0) then - minheap(ipft) = sec_age_patches_per_pft(ipft, 0) - end if - end do - - ! Fourth loop to get min from min-heap then insert 'sec_age_patches_per_pft' and record order in 'order_of_patches' - ! no need to loop over data list here - ipatch = 0 - jpatch(:) = 1 - do ipatch = 1, npacthes - ! Find the minimum and assign the order - curpft = minloc(minheap) - 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) = sec_age_patches_per_pft(ipft, jpatch(curpft)) - else - minheap(curpft) = 0._r8 - end if - 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(sec_age_patches) - deallocate(sec_age_patches_per_pft) - deallocate(minheap) - deallocate(jpatch) - deallocate(pftlen) - - end if ! nocomp - - end if ! logging_time - - ! Loop through all patches and calculate the harvest rate scale based - ! on pre-defined age priority strategy - ! Calculate the harvest rate scale based on pre-defined age priority strategy - ! Only calculate harvest rate scale when we have logging event - ! to reduce computational cost - if (logging_time .and. logging_age_preference == logging_oldest_first) then - ! Maximize harvest rate (99.99%) for older patch first until harvest demand is fullfilled - ! ipatch is the rank of the secondary patch age - 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 - - !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 to actual cases 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. - - !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 + ! 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) - end if ! logging_time + ! 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 diff --git a/biogeochem/FatesForestManagementMod.F90 b/biogeochem/FatesForestManagementMod.F90 new file mode 100644 index 0000000000..df94cc3252 --- /dev/null +++ b/biogeochem/FatesForestManagementMod.F90 @@ -0,0 +1,724 @@ +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 EDLoggingMortalityMod , only : logging_time + use EDParamsMod , only : logging_age_preference, logging_preference_options + use EDParamsMod , only : logging_ifm_harvest_scale + use EDTypesMod , only : AREA, AREA_INV + use FatesGlobals , only : fates_log + 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 : 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_cohort_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 : 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 + + !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 EDPatchDynamicsMod , only : GetPseudoPatchAge, countPatches + + ! !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 + 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 last 4 digits of pseudo-age, i.e., actual patch age. Since 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 + + if (logging_time .and. logging_age_preference == logging_oldest_first) then + + ! 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 + + 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) = mod(GetPseudoPatchAge(currentPatch), 1.e4) + + do curpft=1,numpft + if(currentPatch%nocomp_pft_label .eq. curpft) + pftlen(curpft) = pftlen(curpft) + 1 + order_of_patches_per_pft(curpft, pftlen(curpftpft)) = currentPatch%order_age_since_anthro + age_patches_per_pft(curpft, pftlen(curpftpft)) = mod(GetPseudoPatchAge(currentPatch), 1.e4) + 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) + 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 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_cohort_harvest_rate_scale (dbh, harvest_wp_scale, harvest_rate_scale_patch, harvest_rate_scale_cohort) + + ! !USES: + use EDParamsMod , only : logging_dbhmin + use EDParamsMod , only : logging_dbhmax + + ! !ARGUMENTS: + 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 = csite%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 + + ! ============================================================================ + + 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 From 29d75952771fbcc31f0b30fe83ff6cbe7030cba5 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Wed, 10 Sep 2025 19:25:45 -0700 Subject: [PATCH 15/16] Some minor fixes --- biogeochem/EDLoggingMortalityMod.F90 | 60 +++++++++++++++- biogeochem/FatesForestManagementMod.F90 | 91 ++++++------------------- main/EDMainMod.F90 | 2 +- 3 files changed, 78 insertions(+), 75 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 4b74547ac6..c68d89d9e9 100755 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -63,7 +63,6 @@ module EDLoggingMortalityMod use PRTGenericMod , only : sapw_organ, struct_organ, leaf_organ use PRTGenericMod , only : fnrt_organ, store_organ, repro_organ use FatesAllometryMod , only : set_root_fraction - use FatesForestManagementMod, only : get_cohort_harvest_rate_scale use FatesConstantsMod , only : primaryland, secondaryland, secondary_age_threshold use FatesConstantsMod , only : fates_tiny use FatesConstantsMod , only : months_per_year, days_per_sec, years_per_day, g_per_kg @@ -106,7 +105,8 @@ module EDLoggingMortalityMod public :: get_harvest_rate_carbon public :: get_harvest_debt public :: UpdateHarvestC - + private :: get_cohort_harvest_rate_scale + contains subroutine IsItLoggingTime(is_master,currentSite) @@ -307,7 +307,8 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! Logging size prefrence or change rotation length - call get_cohort_harvest_rate_scale (dbh, harvest_wp_scale, harvest_rate_scale, harvest_rate_scale_cohort) + 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 @@ -1159,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/FatesForestManagementMod.F90 b/biogeochem/FatesForestManagementMod.F90 index df94cc3252..c7c1a0ce21 100644 --- a/biogeochem/FatesForestManagementMod.F90 +++ b/biogeochem/FatesForestManagementMod.F90 @@ -15,11 +15,14 @@ module FatesForestManagementMod ! ============================================================================ ! Uses - use EDLoggingMortalityMod , only : logging_time 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 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 @@ -28,6 +31,7 @@ module FatesForestManagementMod 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 @@ -41,7 +45,6 @@ module FatesForestManagementMod ! Subroutines public :: get_site_harvest_rate_scale public :: get_patch_harvest_rate_scale - public :: get_cohort_harvest_rate_scale public :: get_harvestable_carbon private :: get_harvestable_patch_carbon private :: get_harvested_cohort_carbon @@ -53,7 +56,7 @@ module FatesForestManagementMod subroutine get_site_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, frac_site_primary, frac_site_secondary_mature) ! !USES: - use EDLoggingMortalityMod , only : LoggingMortality_frac + use EDLoggingMortalityMod , only : logging_time, LoggingMortality_frac ! !ARGUMENTS: type(ed_site_type) , intent(inout) :: site_in @@ -174,7 +177,6 @@ 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 EDPatchDynamicsMod , only : GetPseudoPatchAge, countPatches ! !ARGUMENTS: type(ed_site_type) , intent(inout) :: site_in @@ -188,14 +190,14 @@ subroutine get_patch_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, fr type (fates_cohort_type), pointer :: currentCohort integer :: npatches ! Number of patches - integer :: ipatch, , curpft ! index, used in min-heap merge + 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 + 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), 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 @@ -247,10 +249,9 @@ subroutine get_patch_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, fr ! 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 last 4 digits of pseudo-age, i.e., actual patch age. Since 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 + ! 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 if (logging_time .and. logging_age_preference == logging_oldest_first) then @@ -286,13 +287,13 @@ subroutine get_patch_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, fr do while (associated(currentPatch)) ipatch = ipatch + 1 order_of_patches(ipatch) = currentPatch%order_age_since_anthro - age_patches(ipatch) = mod(GetPseudoPatchAge(currentPatch), 1.e4) + age_patches(ipatch) = currentPatch%age do curpft=1,numpft - if(currentPatch%nocomp_pft_label .eq. curpft) + if(currentPatch%nocomp_pft_label .eq. curpft) then pftlen(curpft) = pftlen(curpft) + 1 - order_of_patches_per_pft(curpft, pftlen(curpftpft)) = currentPatch%order_age_since_anthro - age_patches_per_pft(curpft, pftlen(curpftpft)) = mod(GetPseudoPatchAge(currentPatch), 1.e4) + 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 @@ -312,15 +313,14 @@ subroutine get_patch_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, fr jpatch(:) = 1 do ipatch = 1, npacthes ! Find the minimum and assign the order - curpft = minloc(minheap) + 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 + minheap(curpft) = age_patches_per_pft(ipft, jpatch(curpft)) + else + minheap(curpft) = 0._r8 end if end do @@ -469,57 +469,6 @@ end subroutine get_patch_harvest_rate_scale ! ============================================================================ - subroutine get_cohort_harvest_rate_scale (dbh, harvest_wp_scale, harvest_rate_scale_patch, harvest_rate_scale_cohort) - - ! !USES: - use EDParamsMod , only : logging_dbhmin - use EDParamsMod , only : logging_dbhmax - - ! !ARGUMENTS: - 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 = csite%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 - - ! ============================================================================ - subroutine get_harvestable_carbon (csite, site_area, hlm_harvest_catnames, harvestable_forest_c ) !USES: diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 35b3d92ced..d6d013d88d 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -86,7 +86,7 @@ module EDMainMod use FatesPlantHydraulicsMod , only : AccumulateMortalityWaterStorage use FatesAllometryMod , only : h_allom,tree_sai,tree_lai use EDLoggingMortalityMod , only : IsItLoggingTime - use EDLoggingMortalityMod , only : get_harvestable_carbon + use FatesForestManagementMod , only : get_harvestable_carbon use DamageMainMod , only : IsItDamageTime use EDPatchDynamicsMod , only : get_frac_site_primary use EDPatchDynamicsMod , only : get_frac_site_secondary_mature From 30bf63c712b419a8a0b236dcff008045574bb401 Mon Sep 17 00:00:00 2001 From: sshu88 Date: Tue, 16 Sep 2025 13:29:40 -0700 Subject: [PATCH 16/16] More minor fixes. Model version passed 2-years test run. --- biogeochem/FatesForestManagementMod.F90 | 13 ++++++++++++- main/EDParamsMod.F90 | 6 ++++++ 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/biogeochem/FatesForestManagementMod.F90 b/biogeochem/FatesForestManagementMod.F90 index c7c1a0ce21..885bd98a03 100644 --- a/biogeochem/FatesForestManagementMod.F90 +++ b/biogeochem/FatesForestManagementMod.F90 @@ -19,6 +19,7 @@ module FatesForestManagementMod 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 @@ -56,7 +57,7 @@ module FatesForestManagementMod 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 + use EDLoggingMortalityMod , only : logging_time, LoggingMortality_frac ! !ARGUMENTS: type(ed_site_type) , intent(inout) :: site_in @@ -100,6 +101,7 @@ subroutine get_site_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, fra !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 @@ -177,6 +179,7 @@ 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 @@ -225,6 +228,7 @@ subroutine get_patch_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, fr ! 2) harvest the oldest patch first ! Initialize harvest_rate_scale + write(fates_log(),*) 'See flag1' currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -252,8 +256,14 @@ subroutine get_patch_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, fr ! 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 @@ -270,6 +280,7 @@ subroutine get_patch_harvest_rate_scale(site_in, bc_in, harvestable_forest_c, fr ! 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)) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 661c1bda67..e6d4064576 100755 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -522,9 +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)