From 14e34f4d84bbd96209f93f935c3b564760960cd8 Mon Sep 17 00:00:00 2001 From: adivi Date: Thu, 7 May 2020 09:33:57 -0700 Subject: [PATCH 01/36] Initial setup for receiving harvest from hlm Updated fates interface to receive harvest data from hlm. This includes flags denoting whether to use hlm harvest or fates logging, and whether the hlm harvest is area or carbon. Have not implemented the actual harvest yet, and plan to use fates logging parameters to direct it. Changes to be committed: modified: biogeochem/EDLoggingMortalityMod.F90 modified: main/FatesInterfaceMod.F90 --- biogeochem/EDLoggingMortalityMod.F90 | 6 +++ main/FatesInterfaceMod.F90 | 76 ++++++++++++++++++++++++++-- 2 files changed, 78 insertions(+), 4 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index f8c9a4cef8..41b7bced61 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -44,6 +44,8 @@ module EDLoggingMortalityMod use FatesInterfaceMod , only : hlm_model_day use FatesInterfaceMod , only : hlm_day_of_year use FatesInterfaceMod , only : hlm_days_per_year + use FatesInterfaceMod , only : hlm_use_harvest_area + use FatesInterfaceMod , only : hlm_use_harvest_c use FatesInterfaceMod , only : hlm_use_logging use FatesInterfaceMod , only : hlm_use_planthydro use FatesConstantsMod , only : itrue,ifalse @@ -107,6 +109,8 @@ subroutine IsItLoggingTime(is_master,currentSite) if(hlm_use_logging.eq.ifalse) return +! todo: deal with hlm_use_harvest_area and the do_harvest in bc_in + if(icode .eq. 1) then ! Logging is turned off logging_time = .false. @@ -189,6 +193,8 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & if (logging_time) then if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees +! todo: deal with hlm_use_harvest_area vs. biomass here also + ! Pass logging rates to cohort level if (dbh >= logging_dbhmin ) then diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 0254c89123..60dc42bf35 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -44,7 +44,7 @@ module FatesInterfaceMod use PRTGenericMod , only : phosphorus_element use PRTAllometricCarbonMod, only : InitPRTGlobalAllometricCarbon ! use PRTAllometricCNPMod, only : InitPRTGlobalAllometricCNP - + use decompMod , only : bounds_type ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -130,8 +130,32 @@ module FatesInterfaceMod ! 1 = TRUE, 0 = FALSE + integer, public, protected :: hlm_use_harvest_area ! This flag signals whether or not to use + ! harvest area data from the hlm + ! If TRUE, it automatically sets + ! hlm_use_logging to TRUE + ! Only one of hlm_use_harvest_area + ! or hlm_use_harvest_c can be TRUE + + integer, public, protected :: hlm_use_harvest_c ! This flag signals whether or not to use + ! harvest carbon data from the hlm + ! If TRUE, it automatically sets + ! hlm_use_logging to TRUE + ! Only one of hlm_use_harvest_area + ! or hlm_use_harvest_c can be TRUE + ! These data are not available yet + integer, public, protected :: hlm_use_logging ! This flag signals whether or not to use ! the logging module + ! If hlm_use_harvest_area or + ! hlm_use_harvest_c are NOT true, + ! then logging is determined by + ! the fates parameter file + ! If hlm_use_harvest_area or + ! hlm_use_harvest_c are TRUE, + ! then this flag is automatically + ! set and logging is determined + ! by the harvest input from the hlm integer, public, protected :: hlm_use_planthydro ! This flag signals whether or not to use ! plant hydraulics (bchristo/xu methods) @@ -461,7 +485,14 @@ module FatesInterfaceMod real(r8),allocatable :: hksat_sisl(:) ! hydraulic conductivity at saturation (mm H2O /s) real(r8),allocatable :: h2o_liq_sisl(:) ! Liquid water mass in each layer (kg/m2) real(r8) :: smpmin_si ! restriction for min of soil potential (mm) - + + ! Land use + ! --------------------------------------------------------------------------------- + + logical :: hlm_do_harvest_today ! denotes whether hlm harvest data + ! are available for today + real(r8),allocatable :: hlm_harvest(:) ! harvest data from hlm + end type bc_in_type @@ -662,7 +693,7 @@ end subroutine fates_clean ! ==================================================================================== - subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in) + subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, harvest_bounds) ! --------------------------------------------------------------------------------- ! Allocate and Initialze the FATES boundary condition vectors @@ -672,10 +703,10 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in) type(bc_in_type), intent(inout) :: bc_in integer,intent(in) :: nlevsoil_in integer,intent(in) :: nlevdecomp_in + type(bounds_type),intent(in) :: harvest_bounds ! Allocate input boundaries - bc_in%nlevsoil = nlevsoil_in if(nlevsoil_in > numlevsoil_max) then @@ -785,6 +816,15 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in) allocate(bc_in%h2o_liq_sisl(nlevsoil_in)); bc_in%h2o_liq_sisl = nan end if + ! Land use + + ! harvest flags denote data from hlm, + ! while the logging flag signifies only that logging is occurring, + ! which could also just be FATES logging + if (hlm_use_harvest_area .or. hlm_use_harvest_c) then + allocate(bc_in%harvest(harvest_bounds%begg:harvest_bounds%endg)) + end if + return end subroutine allocate_bcin @@ -1429,6 +1469,8 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_parteh_mode = unset_int hlm_use_spitfire = unset_int hlm_use_planthydro = unset_int + hlm_use_harvest_area = unset_int + hlm_use_harvest_c = unset_int hlm_use_logging = unset_int hlm_use_ed_st3 = unset_int hlm_use_ed_prescribed_phys = unset_int @@ -1479,6 +1521,20 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(), *) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' end if + if ( .not.((hlm_use_harvest_area .eq.1).or.(hlm_use_harvest_area.eq.0)) ) then + if (fates_global_verbose()) then + write(fates_log(), *) 'The FATES namelist do_harvest flag must be 0 or 1, exiting' + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if ( .not.((hlm_use_harvest_c .eq.1).or.(hlm_use_harvest_c.eq.0)) ) then + if (fates_global_verbose()) then + write(fates_log(), *) 'The FATES namelist do_harvest flag must be 0 or 1, exiting' + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + if ( .not.((hlm_use_logging .eq.1).or.(hlm_use_logging.eq.0)) ) then if (fates_global_verbose()) then write(fates_log(), *) 'The FATES namelist use_logging flag must be 0 or 1, exiting' @@ -1691,6 +1747,18 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(),*) 'Transfering hlm_use_planthydro= ',ival,' to FATES' end if + case('use_harvest_area') + hlm_use_harvest_area = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering hlm_use_harvest_area= ',ival,' to FATES' + end if + + case('use_harvest_c') + hlm_use_harvest_c = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering hlm_use_harvest_c= ',ival,' to FATES' + end if + case('use_logging') hlm_use_logging = ival if (fates_global_verbose()) then From f0dc220edc43df31b6cf529565c1f0599524c3e3 Mon Sep 17 00:00:00 2001 From: adivi Date: Tue, 12 May 2020 09:53:35 -0700 Subject: [PATCH 02/36] Initial code for land use harvest passed to fates Completed the initial code for passing land use harvest from ELM to FATES. Still need to debug and test. Changes to be committed: modified: biogeochem/EDLoggingMortalityMod.F90 modified: main/FatesConstantsMod.F90 modified: main/FatesInterfaceMod.F90 --- biogeochem/EDLoggingMortalityMod.F90 | 117 +++++++++++++++++++++++---- main/FatesConstantsMod.F90 | 3 + main/FatesInterfaceMod.F90 | 77 +++++++++--------- 3 files changed, 143 insertions(+), 54 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 41b7bced61..d0ba1c4946 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -44,8 +44,8 @@ module EDLoggingMortalityMod use FatesInterfaceMod , only : hlm_model_day use FatesInterfaceMod , only : hlm_day_of_year use FatesInterfaceMod , only : hlm_days_per_year - use FatesInterfaceMod , only : hlm_use_harvest_area - use FatesInterfaceMod , only : hlm_use_harvest_c + use FatesInterfaceMod , only : hlm_use_lu_harvest + use FatesInterfaceMod , only : hlm_use_num_lu_harvest_cats use FatesInterfaceMod , only : hlm_use_logging use FatesInterfaceMod , only : hlm_use_planthydro use FatesConstantsMod , only : itrue,ifalse @@ -58,6 +58,7 @@ module EDLoggingMortalityMod use PRTGenericMod , only : fnrt_organ, store_organ, repro_organ use FatesAllometryMod , only : set_root_fraction use FatesAllometryMod , only : i_biomass_rootprof_context + use FatesConstantsMod , only : primaryforest, secondaryforest, secondary_age_threshold implicit none private @@ -107,12 +108,17 @@ subroutine IsItLoggingTime(is_master,currentSite) logging_time = .false. icode = int(logging_event_code) + ! this is true for either hlm harvest or fates logging if(hlm_use_logging.eq.ifalse) return -! todo: deal with hlm_use_harvest_area and the do_harvest in bc_in + ! all of these are valid for hlm harvest + ! so adjust annual harvest inputs accordingly in LoggingMortality_frac + ! note that the specific event will allow only one hlm harvest, regardless of input + ! code 3, every day, may not work properly because of the exclusive mortality selection + ! even less fequent low harvest rates may be excluded - may need to give harvest priority if(icode .eq. 1) then - ! Logging is turned off + ! Logging is turned off - not sure why we need another switch logging_time = .false. else if(icode .eq. 2) then @@ -173,12 +179,18 @@ end subroutine IsItLoggingTime ! ====================================================================================== subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & - lmort_collateral,lmort_infra, l_degrad ) + lmort_collateral,lmort_infra, l_degrad, & + hlm_harvest, hlm_harvest_catnames, & + use_history, secondary_age) ! Arguments integer, intent(in) :: pft_i ! pft index real(r8), intent(in) :: dbh ! diameter at breast height (cm) integer, intent(in) :: canopy_layer ! canopy layer of this cohort + real(r8), intent(in) :: hlm_harvest ! annual harvest rate per hlm category + character(len=64), intent(in) :: hlm_harvest_catnames ! names of hlm harvest categories + integer, intent(in) :: use_history ! patch level anthro_disturbance_label + real(r8), intent(in) :: secondary_age ! patch level age_since_anthro_disturbance 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 @@ -188,28 +200,103 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! forest patch) ! Parameters + integer :: icode ! Integer equivalent of the event code (parameter file only allows reals) + real(r89) :: harvest_rate ! the final harvest rate to apply to this cohort today real(r8), parameter :: adjustment = 1.0 ! adjustment for mortality rates - + real(r8), parameter :: months_per_year = 12.0 + + icode = int(logging_event_code) + if (logging_time) then if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees -! todo: deal with hlm_use_harvest_area vs. biomass here also +! todo: add a logging_dbhmax parameter, and probably lower the dbhmin one to 30 cm +! todo: change the default logging_event_code to 1 september (-244) +! todo: change the default logging_direct_frac to 0.7, which is closer to a clearcut event +! todo: check outputs against the LUH2 carbon data +! todo: eventually set up distinct harvest practices, each with a set of input paramaeters ! Pass logging rates to cohort level + if (hlm_use_lu_harvest == 0) then + ! 0=use fates logging when logging_time == .true. + ! this means harvest the whole cohort + harvest_rate = 1._r8 + + else if (hlm_use_lu_harvest == 1) then + ! 1=use area fraction from hlm + ! combine forest and non-forest fracs and then apply: + ! primary and secondary area fractions to the logging rates, which are fates parameters + + ! Definitions of the underlying harvest land category variables + ! these are hardcoded to match the LUH input data via landuse.timseries file (see dynHarvestMod) + ! these are fraction of vegetated area harvested, split into five land category variables + ! HARVEST_VH1 = harvest from primary forest + ! HARVEST_VH2 = harvest from primary non-forest + ! HARVEST_SH1 = harvest from secondary mature forest + ! HARVEST_SH2 = harvest from secondary young forest + ! HARVEST_SH3 = harvest from secondary non-forest (assume this is young for biomass) + + ! determine the annual hlm harvest rate for the current cohort based on patch info + harvest_rate = 0._r8 + do h_index = 1,hlm_use_num_lu_harvest_cats + if (use_history .eq. primaryforest) then + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. & + hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") + harvest_rate = harvest_rate + hlm_harvest(h_index) + endif + else if (use_history .eq. secondaryforest .and. + secondary_age >= secondary_age_threshold) then + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") + harvest_rate = harvest_rate + hlm_harvest(h_index) + endif + else if (use_history .eq. secondaryforest .and. + secondary_age < secondary_age_threshold) then + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. & + hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") + harvest_rate = harvest_rate + hlm_harvest(h_index) + endif + endif + end do + + ! calculate today's harvest rate + ! the timing has already been determined by IsItLoggingTime + ! for icode == 2, icode < 0, and icode > 10000 apply the annual rate one time (no calc) + ! Bad logging event flag is caught in IsItLoggingTime, so don't check it here + icode = int(logging_event_code) + !harvest_rate = harvest_rate / hlm_days_per_year + if(icode .eq. 1) then + ! Logging is turned off - not sure why we need another switch + harvest_rate = 0._r8 + else if(icode .eq. 3) then + ! Logging event every day - this may not work due to the mortality exclusivity + harvest_rate = harvest_rate / hlm_days_per_year + else if(icode .eq. 4) then + ! logging event once a month + if(hlm_current_day.eq.1 ) then + harvest_rate = harvest_rate / months_per_year + end if + + else if (hlm_use_lu_harvest == 2) then + ! 2=use carbon from hlm + ! not implemented yet + write(fates_log(),*) 'HLM harvest carbon data not implemented yet. Exiting.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + if (dbh >= logging_dbhmin ) then - lmort_direct = logging_direct_frac * adjustment + lmort_direct = harvest_rate * logging_direct_frac * adjustment l_degrad = 0._r8 else - lmort_direct = 0.0_r8 - l_degrad = logging_direct_frac * adjustment + lmort_direct = 0.0_r8 + l_degrad = harvest_rate * logging_direct_frac * adjustment end if - + if (dbh >= logging_dbhmax_infra) then lmort_infra = 0.0_r8 - l_degrad = l_degrad + logging_mechanical_frac * adjustment + l_degrad = l_degrad + harvest_rate * logging_mechanical_frac * adjustment else - lmort_infra = logging_mechanical_frac * adjustment + lmort_infra = harvest_rate * logging_mechanical_frac * adjustment end if !damage rates for size class < & > threshold_size need to be specified seperately @@ -220,9 +307,9 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! for collateral damage, even understory collateral damage. if (canopy_layer .eq. 1) then - lmort_collateral = logging_collateral_frac * adjustment + lmort_collateral = harvest_rate * logging_collateral_frac * adjustment else - lmort_collateral = 0._r8 + lmort_collateral = 0._r8 endif else diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 3f8aaafc57..f14e8fa4dd 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -30,6 +30,9 @@ module FatesConstantsMod integer, parameter, public :: n_anthro_disturbance_categories = 2 integer, parameter, public :: primaryforest = 1 integer, parameter, public :: secondaryforest = 2 + integer, parameter, public :: secondary_age_threshold = 94 ! less than this value is young + ! based on average age of global + ! secondary 1900s land in hurtt-2011 ! Error Tolerances diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 60dc42bf35..5ab2d262d0 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -130,32 +130,30 @@ module FatesInterfaceMod ! 1 = TRUE, 0 = FALSE - integer, public, protected :: hlm_use_harvest_area ! This flag signals whether or not to use + integer, public, protected :: hlm_use_lu_harvest ! This flag signals whether or not to use ! harvest area data from the hlm - ! If TRUE, it automatically sets - ! hlm_use_logging to TRUE - ! Only one of hlm_use_harvest_area - ! or hlm_use_harvest_c can be TRUE - - integer, public, protected :: hlm_use_harvest_c ! This flag signals whether or not to use - ! harvest carbon data from the hlm - ! If TRUE, it automatically sets - ! hlm_use_logging to TRUE - ! Only one of hlm_use_harvest_area - ! or hlm_use_harvest_c can be TRUE - ! These data are not available yet + ! 0 = do not use lu harvest from hlm + ! 1 = use area fraction of vegetated land + ! from from hlm + ! 2 = use carbon from hlm + ! If 1 or 2, it automatically sets + ! hlm_use_logging to 1 + + integer, public, protected :: hlm_num_lu_harvest_cats ! number of hlm harvest categories + ! this is the first dimension of: + ! harvest_rates in dynHarvestMod + ! bc_in%hlm_harvest and bc_in%hlm_harvest_catnames + integer, public, protected :: hlm_use_logging ! This flag signals whether or not to use ! the logging module - ! If hlm_use_harvest_area or - ! hlm_use_harvest_c are NOT true, + ! If hlm_use_lu_harvest is zero, ! then logging is determined by ! the fates parameter file - ! If hlm_use_harvest_area or - ! hlm_use_harvest_c are TRUE, + ! If hlm_use_lu_harvest is non-zero, ! then this flag is automatically - ! set and logging is determined - ! by the harvest input from the hlm + ! set to 1 and logging is determined + ! by the lu harvest input from the hlm integer, public, protected :: hlm_use_planthydro ! This flag signals whether or not to use ! plant hydraulics (bchristo/xu methods) @@ -491,7 +489,9 @@ module FatesInterfaceMod logical :: hlm_do_harvest_today ! denotes whether hlm harvest data ! are available for today - real(r8),allocatable :: hlm_harvest(:) ! harvest data from hlm + real(r8),allocatable :: hlm_harvest(:) ! annual harvest rate per cat from hlm for a site + character(len=64), allocatable :: hlm_harvest_catnames(:) ! names of hlm_harvest d1 + end type bc_in_type @@ -693,7 +693,7 @@ end subroutine fates_clean ! ==================================================================================== - subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, harvest_bounds) + subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats) ! --------------------------------------------------------------------------------- ! Allocate and Initialze the FATES boundary condition vectors @@ -703,7 +703,7 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, harvest_bounds) type(bc_in_type), intent(inout) :: bc_in integer,intent(in) :: nlevsoil_in integer,intent(in) :: nlevdecomp_in - type(bounds_type),intent(in) :: harvest_bounds + integer,intent(in) :: num_lu_harvest_cats ! Allocate input boundaries @@ -818,11 +818,10 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, harvest_bounds) ! Land use - ! harvest flags denote data from hlm, - ! while the logging flag signifies only that logging is occurring, - ! which could also just be FATES logging - if (hlm_use_harvest_area .or. hlm_use_harvest_c) then - allocate(bc_in%harvest(harvest_bounds%begg:harvest_bounds%endg)) + ! harvest flag denote data from hlm, + ! while the logging flag signifies only that logging is occurring (which could just be FATES logging) + if (hlm_use_lu_harvest .gt. 0) then + allocate(bc_in%hlm_harvest(num_lu_harvest_cats)) end if return @@ -1469,8 +1468,8 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_parteh_mode = unset_int hlm_use_spitfire = unset_int hlm_use_planthydro = unset_int - hlm_use_harvest_area = unset_int - hlm_use_harvest_c = unset_int + hlm_use_lu_harvest = unset_int + hlm_use_num_lu_harvest_cats = unset_int hlm_use_logging = unset_int hlm_use_ed_st3 = unset_int hlm_use_ed_prescribed_phys = unset_int @@ -1521,16 +1520,16 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(), *) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' end if - if ( .not.((hlm_use_harvest_area .eq.1).or.(hlm_use_harvest_area.eq.0)) ) then + if ( (hlm_use_lu_harvest .lt. 0).or.(hlm_use_lu_harvest .gt. 2) ) then if (fates_global_verbose()) then - write(fates_log(), *) 'The FATES namelist do_harvest flag must be 0 or 1, exiting' + write(fates_log(), *) 'The FATES lu_harvest flag must be 0 or 1 or 2, exiting' end if call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if ( .not.((hlm_use_harvest_c .eq.1).or.(hlm_use_harvest_c.eq.0)) ) then + if ( (hlm_use_num_lu_harvest_cats .lt. 0) ) then if (fates_global_verbose()) then - write(fates_log(), *) 'The FATES namelist do_harvest flag must be 0 or 1, exiting' + write(fates_log(), *) 'The FATES number of hlm harvest cats must be >= 0, exiting' end if call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -1747,16 +1746,16 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(),*) 'Transfering hlm_use_planthydro= ',ival,' to FATES' end if - case('use_harvest_area') - hlm_use_harvest_area = ival + case('use_lu_harvest') + hlm_use_lu_harvest = ival if (fates_global_verbose()) then - write(fates_log(),*) 'Transfering hlm_use_harvest_area= ',ival,' to FATES' + write(fates_log(),*) 'Transfering hlm_use_lu_harvest= ',ival,' to FATES' end if - case('use_harvest_c') - hlm_use_harvest_c = ival + case('use_num_lu_harvest_cats') + hlm_use_num_lu_harvest_cats = ival if (fates_global_verbose()) then - write(fates_log(),*) 'Transfering hlm_use_harvest_c= ',ival,' to FATES' + write(fates_log(),*) 'Transfering hlm_use_num_lu_harvest_cats= ',ival,' to FATES' end if case('use_logging') From 22adbed29ec6cea7f56d36be64cc8dd2382768be Mon Sep 17 00:00:00 2001 From: adivi Date: Tue, 12 May 2020 10:56:23 -0700 Subject: [PATCH 03/36] Made sure that LoggingMortality_frac had arguments Changes to be committed: modified: biogeochem/EDLoggingMortalityMod.F90 modified: biogeochem/EDMortalityFunctionsMod.F90 modified: biogeochem/EDPatchDynamicsMod.F90 --- biogeochem/EDLoggingMortalityMod.F90 | 4 ++-- biogeochem/EDMortalityFunctionsMod.F90 | 6 +++++- biogeochem/EDPatchDynamicsMod.F90 | 6 +++++- 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index d0ba1c4946..fc2ad0d6a9 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -187,8 +187,8 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & integer, intent(in) :: pft_i ! pft index real(r8), intent(in) :: dbh ! diameter at breast height (cm) integer, intent(in) :: canopy_layer ! canopy layer of this cohort - real(r8), intent(in) :: hlm_harvest ! annual harvest rate per hlm category - character(len=64), intent(in) :: hlm_harvest_catnames ! names of hlm harvest categories + real(r8), intent(in) :: hlm_harvest(:) ! annual harvest rate per hlm category + character(len=64), intent(in) :: hlm_harvest_catnames(:) ! names of hlm harvest categories integer, intent(in) :: use_history ! patch level anthro_disturbance_label real(r8), intent(in) :: secondary_age ! patch level age_since_anthro_disturbance real(r8), intent(out) :: lmort_direct ! direct (harvestable) mortality fraction diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index bb39e6cf39..d7310cefb8 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -198,7 +198,11 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) currentCohort%lmort_direct, & currentCohort%lmort_collateral, & currentCohort%lmort_infra, & - currentCohort%l_degrad) + currentCohort%l_degrad, & + bc_in%hlm_harvest, & + bc_in%hlm_harvest_catnames, & + currentCohort%patchptr%anthro_disturbance_label, & + currentCohort%patchptr%age_since_anthro_disturbance) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index a68a47a65e..2e8fee92b5 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -193,7 +193,11 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort%frmort = frmort call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_layer, & - lmort_direct,lmort_collateral,lmort_infra,l_degrad ) + lmort_direct,lmort_collateral,lmort_infra,l_degrad,& + bc_in%hlm_harvest, & + bc_in%hlm_harvest_catnames, & + currentPatch%anthro_disturbance_label, & + currentPatch%age_since_anthro_disturbance) currentCohort%lmort_direct = lmort_direct currentCohort%lmort_collateral = lmort_collateral From d29262908124a496afadb3ca1282d3a5f7e9dcd7 Mon Sep 17 00:00:00 2001 From: adivi Date: Tue, 12 May 2020 13:29:30 -0700 Subject: [PATCH 04/36] Fixed num_lu_harvest_cats variable bug Changes to be committed: modified: biogeochem/EDLoggingMortalityMod.F90 modified: main/FatesInterfaceMod.F90 --- biogeochem/EDLoggingMortalityMod.F90 | 4 ++-- main/FatesInterfaceMod.F90 | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index fc2ad0d6a9..537d3d28a2 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -45,7 +45,7 @@ module EDLoggingMortalityMod use FatesInterfaceMod , only : hlm_day_of_year use FatesInterfaceMod , only : hlm_days_per_year use FatesInterfaceMod , only : hlm_use_lu_harvest - use FatesInterfaceMod , only : hlm_use_num_lu_harvest_cats + use FatesInterfaceMod , only : hlm_num_lu_harvest_cats use FatesInterfaceMod , only : hlm_use_logging use FatesInterfaceMod , only : hlm_use_planthydro use FatesConstantsMod , only : itrue,ifalse @@ -239,7 +239,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! determine the annual hlm harvest rate for the current cohort based on patch info harvest_rate = 0._r8 - do h_index = 1,hlm_use_num_lu_harvest_cats + do h_index = 1,hlm_num_lu_harvest_cats if (use_history .eq. primaryforest) then if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. & hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 5ab2d262d0..b8229a198b 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1469,7 +1469,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_use_spitfire = unset_int hlm_use_planthydro = unset_int hlm_use_lu_harvest = unset_int - hlm_use_num_lu_harvest_cats = unset_int + hlm_num_lu_harvest_cats = unset_int hlm_use_logging = unset_int hlm_use_ed_st3 = unset_int hlm_use_ed_prescribed_phys = unset_int @@ -1527,7 +1527,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if ( (hlm_use_num_lu_harvest_cats .lt. 0) ) then + if ( (hlm_num_lu_harvest_cats .lt. 0) ) then if (fates_global_verbose()) then write(fates_log(), *) 'The FATES number of hlm harvest cats must be >= 0, exiting' end if @@ -1752,10 +1752,10 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(),*) 'Transfering hlm_use_lu_harvest= ',ival,' to FATES' end if - case('use_num_lu_harvest_cats') - hlm_use_num_lu_harvest_cats = ival + case('num_lu_harvest_cats') + hlm_num_lu_harvest_cats = ival if (fates_global_verbose()) then - write(fates_log(),*) 'Transfering hlm_use_num_lu_harvest_cats= ',ival,' to FATES' + write(fates_log(),*) 'Transfering hlm_num_lu_harvest_cats= ',ival,' to FATES' end if case('use_logging') From 40a3c0a9a4455d59a32f16293810c3421b0eeaa9 Mon Sep 17 00:00:00 2001 From: adivi Date: Wed, 13 May 2020 15:35:35 -0700 Subject: [PATCH 05/36] HLM harvest data passed to FATES The harvest data from landuse.timeseries are passed to FATES for logging. The FATES logging parameters are used to guide harvest. Only area fraction data are implemented, and these fractions are applied to the FATES logging intensity parameters. The desired effect is for the given logging intensity to be applied to a fraction of the cohort determined by the HLM harvest rate. Separate fractions are applied to primary and secondary patches. This is a working version for testing. Changes to be committed: modified: biogeochem/EDLoggingMortalityMod.F90 --- biogeochem/EDLoggingMortalityMod.F90 | 36 +++++++++++++++------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 537d3d28a2..adc20c993e 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -58,7 +58,7 @@ module EDLoggingMortalityMod use PRTGenericMod , only : fnrt_organ, store_organ, repro_organ use FatesAllometryMod , only : set_root_fraction use FatesAllometryMod , only : i_biomass_rootprof_context - use FatesConstantsMod , only : primaryforest, secondaryforest, secondary_age_threshold + use FatesConstantsMod , only : primaryforest, secondaryforest, secondary_age_threshold implicit none private @@ -84,6 +84,8 @@ module EDLoggingMortalityMod public :: logging_litter_fluxes public :: logging_time public :: IsItLoggingTime + integer, parameter, public :: hlm_harvest_area_fraction = 1 + integer, parameter, public :: hlm_harvest_carbon = 2 contains @@ -199,31 +201,33 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! are moved to newly-anthro-disturbed secondary ! forest patch) - ! Parameters + ! Local variables integer :: icode ! Integer equivalent of the event code (parameter file only allows reals) - real(r89) :: harvest_rate ! the final harvest rate to apply to this cohort today + real(r8) :: harvest_rate ! the final harvest rate to apply to this cohort today + integer :: h_index ! for looping over harvest categories + + ! Parameters real(r8), parameter :: adjustment = 1.0 ! adjustment for mortality rates real(r8), parameter :: months_per_year = 12.0 - icode = int(logging_event_code) - if (logging_time) then if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees ! todo: add a logging_dbhmax parameter, and probably lower the dbhmin one to 30 cm ! todo: change the default logging_event_code to 1 september (-244) ! todo: change the default logging_direct_frac to 0.7, which is closer to a clearcut event +! todo: check secondary forest creation ! todo: check outputs against the LUH2 carbon data ! todo: eventually set up distinct harvest practices, each with a set of input paramaeters ! Pass logging rates to cohort level if (hlm_use_lu_harvest == 0) then - ! 0=use fates logging when logging_time == .true. - ! this means harvest the whole cohort + ! 0=use fates logging parameters directly when logging_time == .true. + ! this means harvest the whole cohort area harvest_rate = 1._r8 - else if (hlm_use_lu_harvest == 1) then + else if (hlm_use_lu_harvest == hlm_harvest_area_fraction) then ! 1=use area fraction from hlm ! combine forest and non-forest fracs and then apply: ! primary and secondary area fractions to the logging rates, which are fates parameters @@ -242,29 +246,28 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & do h_index = 1,hlm_num_lu_harvest_cats if (use_history .eq. primaryforest) then if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. & - hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") + hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") then harvest_rate = harvest_rate + hlm_harvest(h_index) endif - else if (use_history .eq. secondaryforest .and. + else if (use_history .eq. secondaryforest .and. & secondary_age >= secondary_age_threshold) then - if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then harvest_rate = harvest_rate + hlm_harvest(h_index) endif - else if (use_history .eq. secondaryforest .and. + else if (use_history .eq. secondaryforest .and. & secondary_age < secondary_age_threshold) then if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. & - hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") + hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then harvest_rate = harvest_rate + hlm_harvest(h_index) endif endif end do ! calculate today's harvest rate - ! the timing has already been determined by IsItLoggingTime + ! whether to harvest today has already been determined by IsItLoggingTime ! for icode == 2, icode < 0, and icode > 10000 apply the annual rate one time (no calc) ! Bad logging event flag is caught in IsItLoggingTime, so don't check it here icode = int(logging_event_code) - !harvest_rate = harvest_rate / hlm_days_per_year if(icode .eq. 1) then ! Logging is turned off - not sure why we need another switch harvest_rate = 0._r8 @@ -275,9 +278,10 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! logging event once a month if(hlm_current_day.eq.1 ) then harvest_rate = harvest_rate / months_per_year + end if end if - else if (hlm_use_lu_harvest == 2) then + else if (hlm_use_lu_harvest == hlm_harvest_carbon) then ! 2=use carbon from hlm ! not implemented yet write(fates_log(),*) 'HLM harvest carbon data not implemented yet. Exiting.' From 30af2ac433c27a69d1d4c36832d63835551c3720 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 20 May 2020 21:36:36 -0600 Subject: [PATCH 06/36] bugfixes to get it to compile and run --- main/FatesInterfaceMod.F90 | 1 + main/FatesInterfaceTypesMod.F90 | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index f566f797e2..fbc6170b08 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -363,6 +363,7 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats) ! while the logging flag signifies only that logging is occurring (which could just be FATES logging) if (hlm_use_lu_harvest .gt. 0) then allocate(bc_in%hlm_harvest(num_lu_harvest_cats)) + allocate(bc_in%hlm_harvest_catnames(num_lu_harvest_cats)) end if return diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index c26a4163b4..a24515bc5c 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -88,7 +88,7 @@ module FatesInterfaceTypesMod integer, public :: hlm_use_spitfire ! This flag signals whether or not to use SPITFIRE ! 1 = TRUE, 0 = FALSE - integer, public, protected :: hlm_use_lu_harvest ! This flag signals whether or not to use + integer, public :: hlm_use_lu_harvest ! This flag signals whether or not to use ! harvest area data from the hlm ! 0 = do not use lu harvest from hlm ! 1 = use area fraction of vegetated land @@ -97,7 +97,7 @@ module FatesInterfaceTypesMod ! If 1 or 2, it automatically sets ! hlm_use_logging to 1 - integer, public, protected :: hlm_num_lu_harvest_cats ! number of hlm harvest categories + integer, public :: hlm_num_lu_harvest_cats ! number of hlm harvest categories ! this is the first dimension of: ! harvest_rates in dynHarvestMod ! bc_in%hlm_harvest and bc_in%hlm_harvest_catnames From c1e8b35c6f45765d7530e6b7b41898355e0ba82f Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 1 Jun 2020 20:38:53 -0600 Subject: [PATCH 07/36] added disturbance rate diagnostics to site type and history vars --- biogeochem/EDPatchDynamicsMod.F90 | 35 ++++++++++++++- main/EDTypesMod.F90 | 6 +++ main/FatesHistoryInterfaceMod.F90 | 74 +++++++++++++++++++++++++++++++ 3 files changed, 114 insertions(+), 1 deletion(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 99d10c15a6..9e751a9065 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -36,6 +36,8 @@ module EDPatchDynamicsMod use EDTypesMod , only : lg_sf use EDTypesMod , only : dl_sf use EDTypesMod , only : dump_patch + use EDTypesMod , only : N_DIST_TYPES + use EDTypesMod , only : AREA_INV use FatesConstantsMod , only : rsnbl_math_prec use FatesInterfaceTypesMod , only : hlm_use_planthydro use FatesInterfaceTypesMod , only : hlm_numSWb @@ -169,6 +171,7 @@ subroutine disturbance_rates( site_in, bc_in) ! secondary forest patch) real(r8) :: dist_rate_ldist_notharvested integer :: threshold_sizeclass + integer :: i_dist !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) @@ -218,6 +221,9 @@ subroutine disturbance_rates( site_in, bc_in) ! Calculate Disturbance Rates based on the mortality rates just calculated ! --------------------------------------------------------------------------------------------- + ! zero the diagnostic disturbance rate fields + site_in%potential_disturbance_rates(1:N_DIST_TYPES) = 0._r8 + currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -260,9 +266,14 @@ subroutine disturbance_rates( site_in, bc_in) endif ! Fire Disturbance Rate - ! Fires can't burn the whole patch, as this causes /0 errors. currentPatch%disturbance_rates(dtype_ifire) = currentPatch%frac_burnt + do i_dist = 1,N_DIST_TYPES + site_in%potential_disturbance_rates(i_dist) = site_in%potential_disturbance_rates(i_dist) + & + currentPatch%disturbance_rates(i_dist) * currentPatch%area * AREA_INV + end do + + ! Fires can't burn the whole patch, as this causes /0 errors. if (debug) then if (currentPatch%disturbance_rates(dtype_ifire) > 0.98_r8)then write(fates_log(),*) 'very high fire areas', & @@ -424,6 +435,11 @@ subroutine spawn_patches( currentSite, bc_in) site_areadis_primary = 0.0_r8 site_areadis_secondary = 0.0_r8 + ! zero the diagnostic disturbance rate fields + currentSite%disturbance_rates_primary_to_primary(1:N_DIST_TYPES) = 0._r8 + currentSite%disturbance_rates_primary_to_secondary(1:N_DIST_TYPES) = 0._r8 + currentSite%disturbance_rates_secondary_to_secondary(1:N_DIST_TYPES) = 0._r8 + do while(associated(currentPatch)) @@ -449,8 +465,25 @@ subroutine spawn_patches( currentSite, bc_in) (currentPatch%disturbance_mode .ne. dtype_ilog) ) then site_areadis_primary = site_areadis_primary + currentPatch%area * currentPatch%disturbance_rate + + ! track disturbance rates to output to history + currentSite%disturbance_rates_primary_to_primary(currentPatch%disturbance_mode) = & + currentSite%disturbance_rates_primary_to_primary(currentPatch%disturbance_mode) + & + currentPatch%area * currentPatch%disturbance_rate * AREA_INV else site_areadis_secondary = site_areadis_secondary + currentPatch%area * currentPatch%disturbance_rate + + ! track disturbance rates to output to history + if (currentPatch%anthro_disturbance_label .eq. secondaryforest) then + currentSite%disturbance_rates_secondary_to_secondary(currentPatch%disturbance_mode) = & + currentSite%disturbance_rates_secondary_to_secondary(currentPatch%disturbance_mode) + & + currentPatch%area * currentPatch%disturbance_rate * AREA_INV + else + currentSite%disturbance_rates_primary_to_secondary(currentPatch%disturbance_mode) = & + currentSite%disturbance_rates_primary_to_secondary(currentPatch%disturbance_mode) + & + currentPatch%area * currentPatch%disturbance_rate * AREA_INV + endif + endif end if diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 5cb4480e4f..b3bc0bbdd9 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -764,6 +764,12 @@ module EDTypesMod ! Canopy Spread real(r8) :: spread ! dynamic canopy allometric term [unitless] + + ! site-level variables to keep track of the disturbance rates, both actual and "potential" + real(r8) :: disturbance_rates_primary_to_primary(N_DIST_TYPES) ! actual disturbance rates from primary patches to primary patches [m2/m2/day] + real(r8) :: disturbance_rates_primary_to_secondary(N_DIST_TYPES) ! actual disturbance rates from primary patches to secondary patches [m2/m2/day] + real(r8) :: disturbance_rates_secondary_to_secondary(N_DIST_TYPES) ! actual disturbance rates from secondary patches to secondary patches [m2/m2/day] + real(r8) :: potential_disturbance_rates(N_DIST_TYPES) ! "potential" disturb rates (i.e. prior to the "which is most" logic) [m2/m2/day] end type ed_site_type diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 05a236174c..1d67a956bb 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -24,6 +24,10 @@ module FatesHistoryInterfaceMod use EDTypesMod , only : num_vegtemp_mem use EDTypesMod , only : site_massbal_type use EDTypesMod , only : element_list + use EDTypesMod , only : N_DIST_TYPES + use EDTypesMod , only : dtype_ifall + use EDTypesMod , only : dtype_ifire + use EDTypesMod , only : dtype_ilog use FatesIODimensionsMod , only : fates_io_dimension_type use FatesIOVariableKindMod , only : fates_io_variable_kind_type use FatesHistoryVariableType , only : fates_history_variable_type @@ -187,6 +191,14 @@ module FatesHistoryInterfaceMod integer :: ih_canopy_biomass_si integer :: ih_understory_biomass_si + integer :: ih_disturbance_rate_p2p_si + integer :: ih_disturbance_rate_p2s_si + integer :: ih_disturbance_rate_s2s_si + integer :: ih_fire_disturbance_rate_si + integer :: ih_logging_disturbance_rate_si + integer :: ih_fall_disturbance_rate_si + integer :: ih_potential_disturbance_rate_si + ! Indices to site by size-class by age variables integer :: ih_nplant_si_scag integer :: ih_nplant_canopy_si_scag @@ -1774,6 +1786,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_agb_si => this%hvars(ih_agb_si)%r81d, & hio_canopy_biomass_si => this%hvars(ih_canopy_biomass_si)%r81d, & hio_understory_biomass_si => this%hvars(ih_understory_biomass_si)%r81d, & + hio_disturbance_rate_p2p_si => this%hvars(ih_disturbance_rate_p2p_si)%r81d, & + hio_disturbance_rate_p2s_si => this%hvars(ih_disturbance_rate_p2s_si)%r81d, & + hio_disturbance_rate_s2s_si => this%hvars(ih_disturbance_rate_s2s_si)%r81d, & + hio_fire_disturbance_rate_si => this%hvars(ih_fire_disturbance_rate_si)%r81d, & + 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_potential_disturbance_rate_si => this%hvars(ih_potential_disturbance_rate_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, & @@ -2037,6 +2056,25 @@ subroutine update_history_dyn(this,nc,nsites,sites) this%hvars(ih_h2oveg_pheno_err_si)%r81d(io_si) = sites(s)%si_hydr%h2oveg_pheno_err end if + ! output site-level disturbance rates + hio_disturbance_rate_p2p_si = sum(sites(s)%disturbance_rates_primary_to_primary(1:N_DIST_TYPES)) + hio_disturbance_rate_p2s_si = sum(sites(s)%disturbance_rates_primary_to_secondary(1:N_DIST_TYPES)) + hio_disturbance_rate_s2s_si = sum(sites(s)%disturbance_rates_secondary_to_secondary(1:N_DIST_TYPES)) + + hio_fire_disturbance_rate_si = sites(s)%disturbance_rates_primary_to_primary(dtype_ifire) + & + sites(s)%disturbance_rates_primary_to_secondary(dtype_ifire) + & + sites(s)%disturbance_rates_secondary_to_secondary(dtype_ifire) + + hio_logging_disturbance_rate_si = sites(s)%disturbance_rates_primary_to_primary(dtype_ilog) + & + sites(s)%disturbance_rates_primary_to_secondary(dtype_ilog) + & + sites(s)%disturbance_rates_secondary_to_secondary(dtype_ilog) + + hio_fall_disturbance_rate_si = sites(s)%disturbance_rates_primary_to_primary(dtype_ifall) + & + sites(s)%disturbance_rates_primary_to_secondary(dtype_ifall) + & + sites(s)%disturbance_rates_secondary_to_secondary(dtype_ifall) + + hio_potential_disturbance_rate_si = sum(sites(s)%potential_disturbance_rates(1:N_DIST_TYPES)) + ipa = 0 cpatch => sites(s)%oldest_patch do while(associated(cpatch)) @@ -4200,6 +4238,42 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_understory_biomass_si ) + ! disturbance rates + call this%set_history_var(vname='DISTURBANCE_RATE_P2P', units='m2 m-2 d-1', & + long='Disturbance rate from primary to primary lands', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_disturbance_rate_p2p_si ) + + call this%set_history_var(vname='DISTURBANCE_RATE_P2S', units='m2 m-2 d-1', & + long='Disturbance rate from primary to secondary lands', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_disturbance_rate_p2s_si ) + + call this%set_history_var(vname='DISTURBANCE_RATE_S2S', units='m2 m-2 d-1', & + long='Disturbance rate from secondary to secondary lands', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_disturbance_rate_s2s_si ) + + call this%set_history_var(vname='DISTURBANCE_RATE_FIRE', units='m2 m-2 d-1', & + long='Disturbance rate from fire', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_fire_disturbance_rate_si ) + + call this%set_history_var(vname='DISTURBANCE_RATE_LOGGING', units='m2 m-2 d-1', & + long='Disturbance rate from logging', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_logging_disturbance_rate_si ) + + call this%set_history_var(vname='DISTURBANCE_RATE_FALL', units='m2 m-2 d-1', & + long='Disturbance rate from treefall', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_fall_disturbance_rate_si ) + + call this%set_history_var(vname='DISTURBANCE_RATE_POTENTIAL', units='m2 m-2 d-1', & + long='Potential (i.e., including unresolved) disturbance rate', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_potential_disturbance_rate_si ) + ! Canopy Resistance call this%set_history_var(vname='C_STOMATA', units='umol m-2 s-1', & From 597f404bc22947c16ddb6cd007952afd0c30f1e5 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 2 Jun 2020 10:24:58 -0600 Subject: [PATCH 08/36] added site-level diagnostic and history variable for tracking primary/secondary area error from patch fusion --- biogeochem/EDPatchDynamicsMod.F90 | 22 ++++++++++++++++++++++ main/EDTypesMod.F90 | 1 + main/FatesHistoryInterfaceMod.F90 | 10 ++++++++++ 3 files changed, 33 insertions(+) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 9e751a9065..6034d5e28f 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -2087,6 +2087,7 @@ subroutine fuse_patches( csite, bc_in ) integer :: iterate !switch of patch reduction iteration scheme. 1 to keep going, 0 to stop integer :: fuse_flag !do patches get fused (1) or not (0). integer :: i_disttype !iterator over anthropogenic disturbance categories + real(r8) :: primary_land_fraction_beforefusion,primary_land_fraction_afterfusion ! !--------------------------------------------------------------------- @@ -2094,11 +2095,20 @@ subroutine fuse_patches( csite, bc_in ) profiletol = ED_val_patch_fusion_tol + primary_land_fraction_beforefusion = 0._r8 + primary_land_fraction_afterfusion = 0._r8 + nopatches(1:n_anthro_disturbance_categories) = 0 currentPatch => currentSite%youngest_patch do while(associated(currentPatch)) nopatches(currentPatch%anthro_disturbance_label) = & nopatches(currentPatch%anthro_disturbance_label) + 1 + + if (currentPatch%anthro_disturbance_label .eq. primaryforest) then + primary_land_fraction_beforefusion = primary_land_fraction_beforefusion + & + currentPatch%area * AREA_INV + endif + currentPatch => currentPatch%older enddo @@ -2294,6 +2304,18 @@ subroutine fuse_patches( csite, bc_in ) enddo !do while nopatches>maxPatchesPerSite end do ! i_disttype loop + + do while(associated(currentPatch)) + + if (currentPatch%anthro_disturbance_label .eq. primaryforest) then + primary_land_fraction_afterfusion = primary_land_fraction_afterfusion + & + currentPatch%area * AREA_INV + endif + + currentPatch => currentPatch%older + enddo + + currentSite%primary_land_patchfusion_error = primary_land_fraction_afterfusion - primary_land_fraction_beforefusion end subroutine fuse_patches diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index b3bc0bbdd9..2b2e9fd11a 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -770,6 +770,7 @@ module EDTypesMod real(r8) :: disturbance_rates_primary_to_secondary(N_DIST_TYPES) ! actual disturbance rates from primary patches to secondary patches [m2/m2/day] real(r8) :: disturbance_rates_secondary_to_secondary(N_DIST_TYPES) ! actual disturbance rates from secondary patches to secondary patches [m2/m2/day] real(r8) :: potential_disturbance_rates(N_DIST_TYPES) ! "potential" disturb rates (i.e. prior to the "which is most" logic) [m2/m2/day] + real(r8) :: primary_land_patchfusion_error ! error term in total area of primary patches associated with patch fusion [m2/m2/day] end type ed_site_type diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 1d67a956bb..fe7687dfac 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -191,6 +191,7 @@ module FatesHistoryInterfaceMod integer :: ih_canopy_biomass_si integer :: ih_understory_biomass_si + integer :: ih_primaryland_fusion_error_si integer :: ih_disturbance_rate_p2p_si integer :: ih_disturbance_rate_p2s_si integer :: ih_disturbance_rate_s2s_si @@ -1786,6 +1787,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_agb_si => this%hvars(ih_agb_si)%r81d, & hio_canopy_biomass_si => this%hvars(ih_canopy_biomass_si)%r81d, & hio_understory_biomass_si => this%hvars(ih_understory_biomass_si)%r81d, & + hio_primaryland_fusion_error_si => this%hvars(ih_primaryland_fusion_error_si)%r81d, & hio_disturbance_rate_p2p_si => this%hvars(ih_disturbance_rate_p2p_si)%r81d, & hio_disturbance_rate_p2s_si => this%hvars(ih_disturbance_rate_p2s_si)%r81d, & hio_disturbance_rate_s2s_si => this%hvars(ih_disturbance_rate_s2s_si)%r81d, & @@ -2056,6 +2058,9 @@ subroutine update_history_dyn(this,nc,nsites,sites) this%hvars(ih_h2oveg_pheno_err_si)%r81d(io_si) = sites(s)%si_hydr%h2oveg_pheno_err end if + ! error in primary lands from patch fusion + hio_primaryland_fusion_error_si = sites(s)%primary_land_patchfusion_error + ! output site-level disturbance rates hio_disturbance_rate_p2p_si = sum(sites(s)%disturbance_rates_primary_to_primary(1:N_DIST_TYPES)) hio_disturbance_rate_p2s_si = sum(sites(s)%disturbance_rates_primary_to_secondary(1:N_DIST_TYPES)) @@ -4239,6 +4244,11 @@ subroutine define_history_vars(this, initialize_variables) ivar=ivar, initialize=initialize_variables, index = ih_understory_biomass_si ) ! disturbance rates + call this%set_history_var(vname='PRIMARYLAND_PATCHFUSION_ERROR', units='m2 m-2 d-1', & + long='Error in total primary lands associated with patch fusion', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_primaryland_fusion_error_si ) + call this%set_history_var(vname='DISTURBANCE_RATE_P2P', units='m2 m-2 d-1', & long='Disturbance rate from primary to primary lands', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & From 0b6ab98e9fff54dff5695385f3279038b2371425 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 2 Jun 2020 10:26:29 -0600 Subject: [PATCH 09/36] changed patch insertion logic to put newly-disturbed primary patch after all the secondary patches, and error-check if the code tries to merge patches with different labels --- biogeochem/EDPatchDynamicsMod.F90 | 42 ++++++++++++++++++++++++++++--- 1 file changed, 38 insertions(+), 4 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 6034d5e28f..311fe762c2 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -424,6 +424,7 @@ subroutine spawn_patches( currentSite, bc_in) real(r8) :: leaf_burn_frac ! fraction of leaves burned in fire ! for both woody and grass species real(r8) :: leaf_m ! leaf mass during partial burn calculations + logical :: foundfirstprimary ! logical for finding the first primary forest patch !--------------------------------------------------------------------- storesmallcohort => null() ! storage of the smallest cohort for insertion routine @@ -1043,12 +1044,41 @@ subroutine spawn_patches( currentSite, bc_in) if ( site_areadis_primary .gt. nearzero) then currentPatch => currentSite%youngest_patch - new_patch_primary%older => currentPatch - new_patch_primary%younger => null() - currentPatch%younger => new_patch_primary - currentSite%youngest_patch => new_patch_primary + ! insert first primary patch after all the secondary patches, if there are any + if (currentPatch%anthro_disturbance_label .eq. secondaryforest ) then + foundfirstprimary = .false. + do while(associated(currentPatch) .and. .not. foundfirstprimary) + currentPatch => currentPatch%older + if (associated(currentPatch)) then + if (currentPatch%anthro_disturbance_label .eq. primaryforest) then + foundfirstprimary = .true. + endif + endif + end do + if (associated(currentPatch)) then + ! the case where we've found a youngest primary patch + new_patch_primary%older => currentPatch + new_patch_primary%younger => currentPatch%younger + currentPatch%younger%older => new_patch_primary + currentPatch%younger => new_patch_primary + else + ! the case where we haven't, and are putting a primary patch at the oldest end of the + ! linked list (not sure how this could happen, but who knows...) + new_patch_primary%older => null() + new_patch_primary%younger => currentSite%oldest_patch + currentSite%oldest_patch%older => new_patch_primary + currentSite%oldest_patch => new_patch_primary + endif + else + ! the case where there are no secondary patches at the start of the LL (prior logic) + new_patch_primary%older => currentPatch + new_patch_primary%younger => null() + currentPatch%younger => new_patch_primary + currentSite%youngest_patch => new_patch_primary + endif endif + ! insert first secondary at the start of the list if ( site_areadis_secondary .gt. nearzero) then currentPatch => currentSite%youngest_patch new_patch_secondary%older => currentPatch @@ -2363,6 +2393,10 @@ subroutine fuse_2_patches(csite, dp, rp) call rp%litter(el)%FuseLitter(rp%area,dp%area,dp%litter(el)) end do + if ( rp%anthro_disturbance_label .ne. dp%anthro_disturbance_label) then + write(fates_log(),*) 'trying to fuse patches with different anthro_disturbance_label values' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif rp%fuel_eff_moist = (dp%fuel_eff_moist*dp%area + rp%fuel_eff_moist*rp%area) * inv_sum_area rp%livegrass = (dp%livegrass*dp%area + rp%livegrass*rp%area) * inv_sum_area From 6320c5c70b1455472565970697742da0eeba57ab Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 4 Jun 2020 23:03:33 -0600 Subject: [PATCH 10/36] added logic to normalize logging rates by primary and secondary fraction --- biogeochem/EDLoggingMortalityMod.F90 | 23 ++++++++++++++++++++++- biogeochem/EDMortalityFunctionsMod.F90 | 6 ++++-- biogeochem/EDPatchDynamicsMod.F90 | 14 +++++++++++++- main/EDMainMod.F90 | 15 +++++++++++++-- 4 files changed, 52 insertions(+), 6 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 10f71558ce..9bdbfc6c6a 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -59,6 +59,7 @@ module EDLoggingMortalityMod use FatesAllometryMod , only : set_root_fraction use FatesAllometryMod , only : i_biomass_rootprof_context use FatesConstantsMod , only : primaryforest, secondaryforest, secondary_age_threshold + use FatesConstantsMod , only : fates_tiny implicit none private @@ -183,7 +184,8 @@ end subroutine IsItLoggingTime subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & lmort_collateral,lmort_infra, l_degrad, & hlm_harvest, hlm_harvest_catnames, & - use_history, secondary_age) + use_history, secondary_age, & + frac_site_primary) ! Arguments integer, intent(in) :: pft_i ! pft index @@ -200,6 +202,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! but suffer from forest degradation (i.e. they ! are moved to newly-anthro-disturbed secondary ! forest patch) + real(r8), intent(in) :: frac_site_primary ! Local variables integer :: icode ! Integer equivalent of the event code (parameter file only allows reals) @@ -263,6 +266,24 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & endif end do + ! normalize by site-level fractio primary or secondary forest fraction + ! since harvest_rate is specified in 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 (use_history .eq. primaryforest) then + if (frac_site_primary .gt. fates_tiny) then + harvest_rate = min((harvest_rate / frac_site_primary),frac_site_primary) + 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)) + else + harvest_rate = 0._r8 + endif + endif + ! calculate today's harvest rate ! whether to harvest today has already been determined by IsItLoggingTime ! for icode == 2, icode < 0, and icode > 10000 apply the annual rate one time (no calc) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 70168a59af..074194a9a0 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -208,7 +208,7 @@ end subroutine mortality_rates ! ============================================================================ - subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) + subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary) ! ! !DESCRIPTION: @@ -224,6 +224,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) type(ed_site_type), intent(inout), target :: currentSite type(ed_cohort_type),intent(inout), target :: currentCohort type(bc_in_type), intent(in) :: bc_in + real(r8), intent(in) :: frac_site_primary ! ! !LOCAL VARIABLES: real(r8) :: cmort ! starvation mortality rate (fraction per year) @@ -249,7 +250,8 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) bc_in%hlm_harvest, & bc_in%hlm_harvest_catnames, & currentCohort%patchptr%anthro_disturbance_label, & - currentCohort%patchptr%age_since_anthro_disturbance) + currentCohort%patchptr%age_since_anthro_disturbance, & + frac_site_primary) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 311fe762c2..d23815f958 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -172,12 +172,23 @@ subroutine disturbance_rates( site_in, bc_in) real(r8) :: dist_rate_ldist_notharvested integer :: threshold_sizeclass integer :: i_dist + real(r8) :: frac_site_primary !---------------------------------------------------------------------------------------------- ! 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 fractino of the site that is primary land + frac_site_primary = 0._r8 + currentPatch => site_in%oldest_patch + do while (associated(currentPatch)) + if (currentPatch%anthro_disturbance_label .eq. primaryforest) then + frac_site_primary = frac_site_primary + currentPatch%area * AREA_INV + endif + currentPatch => currentPatch%younger + end do + currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -204,7 +215,8 @@ subroutine disturbance_rates( site_in, bc_in) bc_in%hlm_harvest, & bc_in%hlm_harvest_catnames, & currentPatch%anthro_disturbance_label, & - currentPatch%age_since_anthro_disturbance) + currentPatch%age_since_anthro_disturbance, & + frac_site_primary) currentCohort%lmort_direct = lmort_direct currentCohort%lmort_collateral = lmort_collateral diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 1858dfbf03..96d41cd294 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -71,7 +71,7 @@ module EDMainMod use FatesGlobals , only : endrun => fates_endrun use ChecksBalancesMod , only : SiteMassStock use EDMortalityFunctionsMod , only : Mortality_Derivative - + use EDTypesMod , only : AREA_INV use PRTGenericMod, only : carbon12_element use PRTGenericMod, only : all_carbon_elements use PRTGenericMod, only : leaf_organ @@ -305,6 +305,17 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) real(r8) :: current_npp ! place holder for calculating npp each year in prescribed physiology mode !----------------------------------------------------------------------- + real(r8) :: frac_site_primary + + ! first calculate the fractino of the site that is primary land + frac_site_primary = 0._r8 + currentPatch => currentSite%oldest_patch + do while (associated(currentPatch)) + if (currentPatch%anthro_disturbance_label .eq. primaryforest) then + frac_site_primary = frac_site_primary + currentPatch%area * AREA_INV + endif + currentPatch => currentPatch%younger + end do ! Set a pointer to this sites carbon12 mass balance site_cmass => currentSite%mass_balance(element_pos(carbon12_element)) @@ -338,7 +349,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) ft = currentCohort%pft ! Calculate the mortality derivatives - call Mortality_Derivative( currentSite, currentCohort, bc_in ) + call Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_primary ) ! ----------------------------------------------------------------------------- ! Apply Plant Allocation and Reactive Transport From 4517c202408316ad06fe64b6017fe99a1994e3eb Mon Sep 17 00:00:00 2001 From: ckoven Date: Fri, 5 Jun 2020 14:39:17 -0600 Subject: [PATCH 11/36] adding non-woody PFT area to secondary patches at same rate --- biogeochem/EDLoggingMortalityMod.F90 | 178 +++++++++++++-------------- 1 file changed, 89 insertions(+), 89 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 9bdbfc6c6a..4fde0ec37d 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -214,101 +214,101 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & real(r8), parameter :: months_per_year = 12.0 if (logging_time) then - if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees - -! todo: add a logging_dbhmax parameter, and probably lower the dbhmin one to 30 cm -! todo: change the default logging_event_code to 1 september (-244) -! todo: change the default logging_direct_frac to 0.7, which is closer to a clearcut event -! todo: check secondary forest creation -! todo: check outputs against the LUH2 carbon data -! todo: eventually set up distinct harvest practices, each with a set of input paramaeters - - ! Pass logging rates to cohort level - - if (hlm_use_lu_harvest == 0) then - ! 0=use fates logging parameters directly when logging_time == .true. - ! this means harvest the whole cohort area - harvest_rate = 1._r8 - - else if (hlm_use_lu_harvest == hlm_harvest_area_fraction) then - ! 1=use area fraction from hlm - ! combine forest and non-forest fracs and then apply: - ! primary and secondary area fractions to the logging rates, which are fates parameters - - ! Definitions of the underlying harvest land category variables - ! these are hardcoded to match the LUH input data via landuse.timseries file (see dynHarvestMod) - ! these are fraction of vegetated area harvested, split into five land category variables - ! HARVEST_VH1 = harvest from primary forest - ! HARVEST_VH2 = harvest from primary non-forest - ! HARVEST_SH1 = harvest from secondary mature forest - ! HARVEST_SH2 = harvest from secondary young forest - ! HARVEST_SH3 = harvest from secondary non-forest (assume this is young for biomass) - - ! determine the annual hlm harvest rate for the current cohort based on patch info - harvest_rate = 0._r8 - do h_index = 1,hlm_num_lu_harvest_cats - if (use_history .eq. primaryforest) then - if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. & - hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") then - harvest_rate = harvest_rate + hlm_harvest(h_index) - endif - else if (use_history .eq. secondaryforest .and. & - secondary_age >= secondary_age_threshold) then - if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then - harvest_rate = harvest_rate + hlm_harvest(h_index) - endif - else if (use_history .eq. secondaryforest .and. & - secondary_age < secondary_age_threshold) then - if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. & - hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then - harvest_rate = harvest_rate + hlm_harvest(h_index) - endif - endif - end do - - ! normalize by site-level fractio primary or secondary forest fraction - ! since harvest_rate is specified in 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 + ! todo: add a logging_dbhmax parameter, and probably lower the dbhmin one to 30 cm + ! todo: change the default logging_event_code to 1 september (-244) + ! todo: change the default logging_direct_frac to 0.7, which is closer to a clearcut event + ! todo: check secondary forest creation + ! todo: check outputs against the LUH2 carbon data + ! todo: eventually set up distinct harvest practices, each with a set of input paramaeters + + ! Pass logging rates to cohort level + + if (hlm_use_lu_harvest == 0) then + ! 0=use fates logging parameters directly when logging_time == .true. + ! this means harvest the whole cohort area + harvest_rate = 1._r8 + + else if (hlm_use_lu_harvest == hlm_harvest_area_fraction) then + ! 1=use area fraction from hlm + ! combine forest and non-forest fracs and then apply: + ! primary and secondary area fractions to the logging rates, which are fates parameters + + ! Definitions of the underlying harvest land category variables + ! these are hardcoded to match the LUH input data via landuse.timseries file (see dynHarvestMod) + ! these are fraction of vegetated area harvested, split into five land category variables + ! HARVEST_VH1 = harvest from primary forest + ! HARVEST_VH2 = harvest from primary non-forest + ! HARVEST_SH1 = harvest from secondary mature forest + ! HARVEST_SH2 = harvest from secondary young forest + ! HARVEST_SH3 = harvest from secondary non-forest (assume this is young for biomass) + + ! determine the annual hlm harvest rate for the current cohort based on patch info + harvest_rate = 0._r8 + do h_index = 1,hlm_num_lu_harvest_cats if (use_history .eq. primaryforest) then - if (frac_site_primary .gt. fates_tiny) then - harvest_rate = min((harvest_rate / frac_site_primary),frac_site_primary) - else - harvest_rate = 0._r8 + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. & + hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") then + harvest_rate = harvest_rate + hlm_harvest(h_index) 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)) - else - harvest_rate = 0._r8 + else if (use_history .eq. secondaryforest .and. & + secondary_age >= secondary_age_threshold) then + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then + harvest_rate = harvest_rate + hlm_harvest(h_index) + endif + else if (use_history .eq. secondaryforest .and. & + secondary_age < secondary_age_threshold) then + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. & + hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then + harvest_rate = harvest_rate + hlm_harvest(h_index) endif endif + end do - ! calculate today's harvest rate - ! whether to harvest today has already been determined by IsItLoggingTime - ! for icode == 2, icode < 0, and icode > 10000 apply the annual rate one time (no calc) - ! Bad logging event flag is caught in IsItLoggingTime, so don't check it here - icode = int(logging_event_code) - if(icode .eq. 1) then - ! Logging is turned off - not sure why we need another switch + ! normalize by site-level fractio primary or secondary forest fraction + ! since harvest_rate is specified in 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 (use_history .eq. primaryforest) then + if (frac_site_primary .gt. fates_tiny) then + harvest_rate = min((harvest_rate / frac_site_primary),frac_site_primary) + else harvest_rate = 0._r8 - else if(icode .eq. 3) then - ! Logging event every day - this may not work due to the mortality exclusivity - harvest_rate = harvest_rate / hlm_days_per_year - else if(icode .eq. 4) then - ! logging event once a month - if(hlm_current_day.eq.1 ) then - harvest_rate = harvest_rate / months_per_year - end if + 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)) + else + harvest_rate = 0._r8 + endif + endif + + ! calculate today's harvest rate + ! whether to harvest today has already been determined by IsItLoggingTime + ! for icode == 2, icode < 0, and icode > 10000 apply the annual rate one time (no calc) + ! Bad logging event flag is caught in IsItLoggingTime, so don't check it here + icode = int(logging_event_code) + if(icode .eq. 1) then + ! Logging is turned off - not sure why we need another switch + harvest_rate = 0._r8 + else if(icode .eq. 3) then + ! Logging event every day - this may not work due to the mortality exclusivity + harvest_rate = harvest_rate / hlm_days_per_year + else if(icode .eq. 4) then + ! logging event once a month + if(hlm_current_day.eq.1 ) then + harvest_rate = harvest_rate / months_per_year end if + end if - else if (hlm_use_lu_harvest == hlm_harvest_carbon) then - ! 2=use carbon from hlm - ! not implemented yet - write(fates_log(),*) 'HLM harvest carbon data not implemented yet. Exiting.' - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif + else if (hlm_use_lu_harvest == hlm_harvest_carbon) then + ! 2=use carbon from hlm + ! not implemented yet + write(fates_log(),*) 'HLM harvest carbon data not implemented yet. Exiting.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees + if (dbh >= logging_dbhmin ) then lmort_direct = harvest_rate * logging_direct_frac * adjustment l_degrad = 0._r8 @@ -337,11 +337,11 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & lmort_collateral = 0._r8 endif - else + else ! non-woody plants in the canopy still become moved to secondary patch at same rate lmort_direct = 0.0_r8 lmort_collateral = 0.0_r8 - lmort_infra = 0.0_r8 - l_degrad = 0.0_r8 + lmort_infra = harvest_rate * logging_mechanical_frac * adjustment + l_degrad = harvest_rate * logging_direct_frac * adjustment end if else lmort_direct = 0.0_r8 From 42e3c2e77f77a52c4ba42efe1b9ffb24dec01aba Mon Sep 17 00:00:00 2001 From: ckoven Date: Fri, 5 Jun 2020 17:38:22 -0600 Subject: [PATCH 12/36] added logic to disturb a fraction of bare-ground patch area, if any exists, during logging --- biogeochem/EDLoggingMortalityMod.F90 | 183 +++++++++++++++------------ biogeochem/EDPatchDynamicsMod.F90 | 57 +++++++-- main/EDMainMod.F90 | 12 +- 3 files changed, 156 insertions(+), 96 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 4fde0ec37d..28a6fef229 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -85,6 +85,7 @@ module EDLoggingMortalityMod public :: logging_litter_fluxes public :: logging_time public :: IsItLoggingTime + public :: get_harvest_rate_area integer, parameter, public :: hlm_harvest_area_fraction = 1 integer, parameter, public :: hlm_harvest_carbon = 2 @@ -205,22 +206,16 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & real(r8), intent(in) :: frac_site_primary ! Local variables - integer :: icode ! Integer equivalent of the event code (parameter file only allows reals) real(r8) :: harvest_rate ! the final harvest rate to apply to this cohort today - integer :: h_index ! for looping over harvest categories - - ! Parameters - real(r8), parameter :: adjustment = 1.0 ! adjustment for mortality rates - real(r8), parameter :: months_per_year = 12.0 + ! todo: add a logging_dbhmax parameter, and probably lower the dbhmin one to 30 cm + ! todo: change the default logging_event_code to 1 september (-244) + ! todo: change the default logging_direct_frac to 0.7, which is closer to a clearcut event + ! todo: check secondary forest creation + ! todo: check outputs against the LUH2 carbon data + ! todo: eventually set up distinct harvest practices, each with a set of input paramaeters + if (logging_time) then - ! todo: add a logging_dbhmax parameter, and probably lower the dbhmin one to 30 cm - ! todo: change the default logging_event_code to 1 september (-244) - ! todo: change the default logging_direct_frac to 0.7, which is closer to a clearcut event - ! todo: check secondary forest creation - ! todo: check outputs against the LUH2 carbon data - ! todo: eventually set up distinct harvest practices, each with a set of input paramaeters - ! Pass logging rates to cohort level if (hlm_use_lu_harvest == 0) then @@ -241,64 +236,9 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! HARVEST_SH1 = harvest from secondary mature forest ! HARVEST_SH2 = harvest from secondary young forest ! HARVEST_SH3 = harvest from secondary non-forest (assume this is young for biomass) - - ! determine the annual hlm harvest rate for the current cohort based on patch info - harvest_rate = 0._r8 - do h_index = 1,hlm_num_lu_harvest_cats - if (use_history .eq. primaryforest) then - if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. & - hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") then - harvest_rate = harvest_rate + hlm_harvest(h_index) - endif - else if (use_history .eq. secondaryforest .and. & - secondary_age >= secondary_age_threshold) then - if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then - harvest_rate = harvest_rate + hlm_harvest(h_index) - endif - else if (use_history .eq. secondaryforest .and. & - secondary_age < secondary_age_threshold) then - if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. & - hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then - harvest_rate = harvest_rate + hlm_harvest(h_index) - endif - endif - end do - - ! normalize by site-level fractio primary or secondary forest fraction - ! since harvest_rate is specified in 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 (use_history .eq. primaryforest) then - if (frac_site_primary .gt. fates_tiny) then - harvest_rate = min((harvest_rate / frac_site_primary),frac_site_primary) - 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)) - else - harvest_rate = 0._r8 - endif - endif - - ! calculate today's harvest rate - ! whether to harvest today has already been determined by IsItLoggingTime - ! for icode == 2, icode < 0, and icode > 10000 apply the annual rate one time (no calc) - ! Bad logging event flag is caught in IsItLoggingTime, so don't check it here - icode = int(logging_event_code) - if(icode .eq. 1) then - ! Logging is turned off - not sure why we need another switch - harvest_rate = 0._r8 - else if(icode .eq. 3) then - ! Logging event every day - this may not work due to the mortality exclusivity - harvest_rate = harvest_rate / hlm_days_per_year - else if(icode .eq. 4) then - ! logging event once a month - if(hlm_current_day.eq.1 ) then - harvest_rate = harvest_rate / months_per_year - end if - end if + + call get_harvest_rate_area (use_history, hlm_harvest_catnames, hlm_harvest, & + frac_site_primary, secondary_age, harvest_rate) else if (hlm_use_lu_harvest == hlm_harvest_carbon) then ! 2=use carbon from hlm @@ -310,18 +250,18 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees if (dbh >= logging_dbhmin ) then - lmort_direct = harvest_rate * logging_direct_frac * adjustment + lmort_direct = harvest_rate * logging_direct_frac l_degrad = 0._r8 else lmort_direct = 0.0_r8 - l_degrad = harvest_rate * logging_direct_frac * adjustment + l_degrad = harvest_rate * logging_direct_frac end if if (dbh >= logging_dbhmax_infra) then lmort_infra = 0.0_r8 - l_degrad = l_degrad + harvest_rate * logging_mechanical_frac * adjustment + l_degrad = l_degrad + harvest_rate * logging_mechanical_frac else - lmort_infra = harvest_rate * logging_mechanical_frac * adjustment + lmort_infra = harvest_rate * logging_mechanical_frac end if !damage rates for size class < & > threshold_size need to be specified seperately @@ -332,7 +272,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! for collateral damage, even understory collateral damage. if (canopy_layer .eq. 1) then - lmort_collateral = harvest_rate * logging_collateral_frac * adjustment + lmort_collateral = harvest_rate * logging_collateral_frac else lmort_collateral = 0._r8 endif @@ -340,8 +280,8 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & else ! non-woody plants in the canopy still become moved to secondary patch at same rate lmort_direct = 0.0_r8 lmort_collateral = 0.0_r8 - lmort_infra = harvest_rate * logging_mechanical_frac * adjustment - l_degrad = harvest_rate * logging_direct_frac * adjustment + lmort_infra = harvest_rate * logging_mechanical_frac + l_degrad = harvest_rate * logging_direct_frac end if else lmort_direct = 0.0_r8 @@ -354,6 +294,93 @@ end subroutine LoggingMortality_frac ! ============================================================================ + subroutine get_harvest_rate_area (use_history, hlm_harvest_catnames, hlm_harvest, & + frac_site_primary, secondary_age, harvest_rate) + + + ! ------------------------------------------------------------------------------------------- + ! + ! DESCRIPTION: + ! get the area-based harvest rates based on info passed to FATES from the bioundary conditions in. + ! assumes logging_time == true + + ! Arguments + real(r8), intent(in) :: hlm_harvest(:) ! annual harvest rate per hlm category + character(len=64), intent(in) :: hlm_harvest_catnames(:) ! names of hlm harvest categories + integer, intent(in) :: use_history ! patch level anthro_disturbance_label + real(r8), intent(in) :: secondary_age ! patch level age_since_anthro_disturbance + real(r8), intent(in) :: frac_site_primary + real(r8), intent(out) :: harvest_rate + + ! Local Variables + integer :: h_index ! for looping over harvest categories + integer :: icode ! Integer equivalent of the event code (parameter file only allows reals) + + ! Parameters + real(r8), parameter :: months_per_year = 12.0 + + ! determine the annual hlm harvest rate for the current cohort based on patch info + harvest_rate = 0._r8 + do h_index = 1,hlm_num_lu_harvest_cats + if (use_history .eq. primaryforest) then + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. & + hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") then + harvest_rate = harvest_rate + hlm_harvest(h_index) + endif + else if (use_history .eq. secondaryforest .and. & + secondary_age >= secondary_age_threshold) then + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then + harvest_rate = harvest_rate + hlm_harvest(h_index) + endif + else if (use_history .eq. secondaryforest .and. & + secondary_age < secondary_age_threshold) then + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. & + hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then + harvest_rate = harvest_rate + hlm_harvest(h_index) + endif + endif + end do + + ! normalize by site-level fractio primary or secondary forest fraction + ! since harvest_rate is specified in 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 (use_history .eq. primaryforest) then + if (frac_site_primary .gt. fates_tiny) then + harvest_rate = min((harvest_rate / frac_site_primary),frac_site_primary) + 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)) + else + harvest_rate = 0._r8 + endif + endif + + ! calculate today's harvest rate + ! whether to harvest today has already been determined by IsItLoggingTime + ! for icode == 2, icode < 0, and icode > 10000 apply the annual rate one time (no calc) + ! Bad logging event flag is caught in IsItLoggingTime, so don't check it here + icode = int(logging_event_code) + if(icode .eq. 1) then + ! Logging is turned off - not sure why we need another switch + harvest_rate = 0._r8 + else if(icode .eq. 3) then + ! Logging event every day - this may not work due to the mortality exclusivity + harvest_rate = harvest_rate / hlm_days_per_year + else if(icode .eq. 4) then + ! logging event once a month + if(hlm_current_day.eq.1 ) then + harvest_rate = harvest_rate / months_per_year + end if + end if + + end subroutine get_harvest_rate_area + + ! ============================================================================ + subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site_areadis) ! ------------------------------------------------------------------------------------------- diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index d23815f958..7a437d045e 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -39,6 +39,7 @@ module EDPatchDynamicsMod use EDTypesMod , only : N_DIST_TYPES use EDTypesMod , only : AREA_INV use FatesConstantsMod , only : rsnbl_math_prec + use FatesConstantsMod , only : fates_tiny use FatesInterfaceTypesMod , only : hlm_use_planthydro use FatesInterfaceTypesMod , only : hlm_numSWb use FatesInterfaceTypesMod , only : bc_in_type @@ -52,6 +53,7 @@ module EDPatchDynamicsMod use FatesPlantHydraulicsMod, only : DeallocateHydrCohort use EDLoggingMortalityMod, only : logging_litter_fluxes use EDLoggingMortalityMod, only : logging_time + use EDLoggingMortalityMod, only : get_harvest_rate_area use EDParamsMod , only : fates_mortality_disturbance_fraction use FatesAllometryMod , only : carea_allom use FatesAllometryMod , only : set_root_fraction @@ -99,6 +101,7 @@ module EDPatchDynamicsMod public :: check_patch_area public :: set_patchno private:: fuse_2_patches + public :: get_frac_site_primary character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -173,6 +176,7 @@ subroutine disturbance_rates( site_in, bc_in) integer :: threshold_sizeclass integer :: i_dist real(r8) :: frac_site_primary + real(r8) :: harvest_rate !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) @@ -180,15 +184,8 @@ subroutine disturbance_rates( site_in, bc_in) !---------------------------------------------------------------------------------------------- ! first calculate the fractino of the site that is primary land - frac_site_primary = 0._r8 - currentPatch => site_in%oldest_patch - do while (associated(currentPatch)) - if (currentPatch%anthro_disturbance_label .eq. primaryforest) then - frac_site_primary = frac_site_primary + currentPatch%area * AREA_INV - endif - currentPatch => currentPatch%younger - end do - + call get_frac_site_primary(site_in, frac_site_primary) + currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -280,6 +277,18 @@ subroutine disturbance_rates( site_in, bc_in) ! Fire Disturbance Rate currentPatch%disturbance_rates(dtype_ifire) = currentPatch%frac_burnt + ! for non-closed-canopy areas subject to logging, add an additional increment of area disturbed + ! equivalent to the fradction loged to account for transfer of interstitial ground area to new secondary lands + if ( currentPatch%disturbance_rates(dtype_ilog) .gt. fates_tiny .and. & + (currentPatch%area - currentPatch%total_canopy_area) .gt. fates_tiny ) then + + call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, bc_in%hlm_harvest, & + frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate) + + currentPatch%disturbance_rates(dtype_ilog) = currentPatch%disturbance_rates(dtype_ilog) + & + (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate * AREA_INV + endif + do i_dist = 1,N_DIST_TYPES site_in%potential_disturbance_rates(i_dist) = site_in%potential_disturbance_rates(i_dist) + & currentPatch%disturbance_rates(i_dist) * currentPatch%area * AREA_INV @@ -2805,4 +2814,34 @@ function countPatches( nsites, sites ) result ( totNumPatches ) end function countPatches + ! ===================================================================================== + + subroutine get_frac_site_primary(site_in, frac_site_primary) + + ! + ! !DESCRIPTION: + ! Calculate how much of a site is primary land + ! + ! !USES: + use EDTypesMod , only : ed_site_type + ! + ! !ARGUMENTS: + type(ed_site_type) , intent(in), target :: site_in + real(r8) , intent(out) :: frac_site_primary + + ! !LOCAL VARIABLES: + type (ed_patch_type), pointer :: currentPatch + + frac_site_primary = 0._r8 + currentPatch => site_in%oldest_patch + do while (associated(currentPatch)) + if (currentPatch%anthro_disturbance_label .eq. primaryforest) then + frac_site_primary = frac_site_primary + currentPatch%area * AREA_INV + endif + currentPatch => currentPatch%younger + end do + + end subroutine get_frac_site_primary + end module EDPatchDynamicsMod + diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 96d41cd294..adb803e639 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -68,6 +68,7 @@ module EDMainMod use FatesAllometryMod , only : h_allom,tree_sai,tree_lai use FatesPlantHydraulicsMod , only : UpdateSizeDepRhizHydStates use EDLoggingMortalityMod , only : IsItLoggingTime + use EDPatchDynamicsMod , only : get_frac_site_primary use FatesGlobals , only : endrun => fates_endrun use ChecksBalancesMod , only : SiteMassStock use EDMortalityFunctionsMod , only : Mortality_Derivative @@ -307,15 +308,8 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) !----------------------------------------------------------------------- real(r8) :: frac_site_primary - ! first calculate the fractino of the site that is primary land - frac_site_primary = 0._r8 - currentPatch => currentSite%oldest_patch - do while (associated(currentPatch)) - if (currentPatch%anthro_disturbance_label .eq. primaryforest) then - frac_site_primary = frac_site_primary + currentPatch%area * AREA_INV - endif - currentPatch => currentPatch%younger - end do + + call get_frac_site_primary(currentSite, frac_site_primary) ! Set a pointer to this sites carbon12 mass balance site_cmass => currentSite%mass_balance(element_pos(carbon12_element)) From 7847c61121cbf851335b53a92195bc58b0590cd2 Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 8 Jun 2020 14:22:22 -0600 Subject: [PATCH 13/36] added calculation of total_canopy_area prior to bare-ground harvest-associated disturbance --- biogeochem/EDPatchDynamicsMod.F90 | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 7a437d045e..58d4cccf00 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -233,6 +233,20 @@ subroutine disturbance_rates( site_in, bc_in) ! zero the diagnostic disturbance rate fields site_in%potential_disturbance_rates(1:N_DIST_TYPES) = 0._r8 + ! summarize how open the canopy is prior to resolving the disturbance + currentPatch => site_in%oldest_patch + do while (associated(currentPatch)) + currentPatch%total_canopy_area = 0._r8 + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + if(currentCohort%canopy_layer==1)then + currentPatch%total_canopy_area = currentPatch%total_canopy_area + currentCohort%c_area + end if + currentCohort => currentCohort%taller + end do + currentPatch => currentPatch%younger + end do + currentPatch => site_in%oldest_patch do while (associated(currentPatch)) From 29dabcca2ac2baff7b5ebb861fd571f51f6a69af Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 8 Jun 2020 16:42:43 -0600 Subject: [PATCH 14/36] added diagnostic history variable of primary forest rates from BCs --- main/FatesHistoryInterfaceMod.F90 | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index fe7687dfac..478d4a18de 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -199,6 +199,7 @@ module FatesHistoryInterfaceMod integer :: ih_logging_disturbance_rate_si integer :: ih_fall_disturbance_rate_si integer :: ih_potential_disturbance_rate_si + integer :: ih_harvest_primary_bcin_si ! Indices to site by size-class by age variables integer :: ih_nplant_si_scag @@ -1638,7 +1639,7 @@ end subroutine update_history_cbal ! ==================================================================================== - subroutine update_history_dyn(this,nc,nsites,sites) + subroutine update_history_dyn(this,nc,nsites,sites,bc_in) ! --------------------------------------------------------------------------------- ! This is the call to update the history IO arrays that are expected to only change @@ -1665,6 +1666,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) integer , intent(in) :: nc ! clump index integer , intent(in) :: nsites type(ed_site_type) , intent(inout), target :: sites(nsites) + type(bc_in_type) , intent(in) :: bc_in(nsites) ! Locals type(litter_type), pointer :: litt_c ! Pointer to the carbon12 litter pool @@ -1795,6 +1797,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) 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_potential_disturbance_rate_si => this%hvars(ih_potential_disturbance_rate_si)%r81d, & + hio_harvest_primary_bcin_si => this%hvars(ih_harvest_primary_bcin_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, & @@ -2080,6 +2083,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_potential_disturbance_rate_si = sum(sites(s)%potential_disturbance_rates(1:N_DIST_TYPES)) + hio_harvest_primary_bcin_si = sum(bc_in(s)%hlm_harvest(1:2)) + ipa = 0 cpatch => sites(s)%oldest_patch do while(associated(cpatch)) @@ -4284,6 +4289,11 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_potential_disturbance_rate_si ) + call this%set_history_var(vname='HARVEST_PRIMARY_RATE_BCIN', units='m2 m-2 d-1', & + long='harvest disturbance rate input from HLM', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_harvest_primary_bcin_si ) + ! Canopy Resistance call this%set_history_var(vname='C_STOMATA', units='umol m-2 s-1', & From 6b574f9415181087e65ea79a50cac17246537342 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 9 Jun 2020 12:51:26 -0600 Subject: [PATCH 15/36] fixed error on logging-related disturbance of bare-ground area --- biogeochem/EDPatchDynamicsMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 58d4cccf00..33c949327d 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -293,14 +293,14 @@ subroutine disturbance_rates( site_in, bc_in) ! for non-closed-canopy areas subject to logging, add an additional increment of area disturbed ! equivalent to the fradction loged to account for transfer of interstitial ground area to new secondary lands - if ( currentPatch%disturbance_rates(dtype_ilog) .gt. fates_tiny .and. & + if ( logging_time .and. & (currentPatch%area - currentPatch%total_canopy_area) .gt. fates_tiny ) then call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, bc_in%hlm_harvest, & frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate) currentPatch%disturbance_rates(dtype_ilog) = currentPatch%disturbance_rates(dtype_ilog) + & - (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate * AREA_INV + (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area endif do i_dist = 1,N_DIST_TYPES From 3afe8edb9c5cf7ef15c7a900e36e6b6acfbfce32 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 9 Jun 2020 13:35:17 -0600 Subject: [PATCH 16/36] adding patch area conservation error term to history output --- biogeochem/EDPatchDynamicsMod.F90 | 5 ++++- main/EDMainMod.F90 | 3 +++ main/EDTypesMod.F90 | 1 + main/FatesHistoryInterfaceMod.F90 | 10 ++++++++++ 4 files changed, 18 insertions(+), 1 deletion(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 33c949327d..3ccd344242 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1161,7 +1161,7 @@ subroutine check_patch_area( currentSite ) ! !USES: ! ! !ARGUMENTS: - type(ed_site_type), intent(in), target :: currentSite + type(ed_site_type), intent(inout), target :: currentSite ! ! !LOCAL VARIABLES: real(r8) :: areatot @@ -1216,6 +1216,9 @@ subroutine check_patch_area( currentSite ) currentSite%mass_balance(el)%patch_resize_err + mass_gain end do + + currentSite%patch_area_conservation_error = currentSite%patch_area_conservation_error + & + area_site-areatot largestPatch%area = largestPatch%area + (area_site-areatot) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index adb803e639..4c5286bb0a 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -241,6 +241,9 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in) ! Patch dynamics sub-routines: fusion, new patch creation (spwaning), termination. !********************************************************************************* + ! zero site-level patch area error tracking + currentSite%patch_area_conservation_error = 0._r8 + ! make new patches from disturbed land if ( hlm_use_ed_st3.eq.ifalse ) then call spawn_patches(currentSite, bc_in) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 2b2e9fd11a..a2b5df5892 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -771,6 +771,7 @@ module EDTypesMod real(r8) :: disturbance_rates_secondary_to_secondary(N_DIST_TYPES) ! actual disturbance rates from secondary patches to secondary patches [m2/m2/day] real(r8) :: potential_disturbance_rates(N_DIST_TYPES) ! "potential" disturb rates (i.e. prior to the "which is most" logic) [m2/m2/day] real(r8) :: primary_land_patchfusion_error ! error term in total area of primary patches associated with patch fusion [m2/m2/day] + real(r8) :: patch_area_conservation_error ! error term in total area of patches associated with patch area conservation [m2/m2/day] end type ed_site_type diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 478d4a18de..c6a916a4a1 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -191,6 +191,7 @@ module FatesHistoryInterfaceMod integer :: ih_canopy_biomass_si integer :: ih_understory_biomass_si + integer :: ih_patch_area_conservation_error_si integer :: ih_primaryland_fusion_error_si integer :: ih_disturbance_rate_p2p_si integer :: ih_disturbance_rate_p2s_si @@ -1789,6 +1790,7 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_agb_si => this%hvars(ih_agb_si)%r81d, & hio_canopy_biomass_si => this%hvars(ih_canopy_biomass_si)%r81d, & hio_understory_biomass_si => this%hvars(ih_understory_biomass_si)%r81d, & + hio_patch_area_conservation_error_si => this%hvars(ih_patch_area_conservation_error_si)%r81d, & hio_primaryland_fusion_error_si => this%hvars(ih_primaryland_fusion_error_si)%r81d, & hio_disturbance_rate_p2p_si => this%hvars(ih_disturbance_rate_p2p_si)%r81d, & hio_disturbance_rate_p2s_si => this%hvars(ih_disturbance_rate_p2s_si)%r81d, & @@ -2064,6 +2066,9 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) ! error in primary lands from patch fusion hio_primaryland_fusion_error_si = sites(s)%primary_land_patchfusion_error + ! total error in patch area fixed during conservation check + hio_patch_area_conservation_error_si = sites(s)%patch_area_conservation_error + ! output site-level disturbance rates hio_disturbance_rate_p2p_si = sum(sites(s)%disturbance_rates_primary_to_primary(1:N_DIST_TYPES)) hio_disturbance_rate_p2s_si = sum(sites(s)%disturbance_rates_primary_to_secondary(1:N_DIST_TYPES)) @@ -4249,6 +4254,11 @@ subroutine define_history_vars(this, initialize_variables) ivar=ivar, initialize=initialize_variables, index = ih_understory_biomass_si ) ! disturbance rates + call this%set_history_var(vname='PATCH_AREA_CONSERVATION_ERROR', units='m2 m-2 d-1', & + long='Error in total patch area fixed to ensure conservation', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_patch_area_conservation_error_si ) + call this%set_history_var(vname='PRIMARYLAND_PATCHFUSION_ERROR', units='m2 m-2 d-1', & long='Error in total primary lands associated with patch fusion', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & From 2bf49b8dcf1aa6cf6d3bf3f7a5d65bfa12a4179d Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 9 Jun 2020 13:35:28 -0600 Subject: [PATCH 17/36] Revert "adding patch area conservation error term to history output" This reverts commit 3afe8edb9c5cf7ef15c7a900e36e6b6acfbfce32. --- biogeochem/EDPatchDynamicsMod.F90 | 5 +---- main/EDMainMod.F90 | 3 --- main/EDTypesMod.F90 | 1 - main/FatesHistoryInterfaceMod.F90 | 10 ---------- 4 files changed, 1 insertion(+), 18 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 3ccd344242..33c949327d 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1161,7 +1161,7 @@ subroutine check_patch_area( currentSite ) ! !USES: ! ! !ARGUMENTS: - type(ed_site_type), intent(inout), target :: currentSite + type(ed_site_type), intent(in), target :: currentSite ! ! !LOCAL VARIABLES: real(r8) :: areatot @@ -1216,9 +1216,6 @@ subroutine check_patch_area( currentSite ) currentSite%mass_balance(el)%patch_resize_err + mass_gain end do - - currentSite%patch_area_conservation_error = currentSite%patch_area_conservation_error + & - area_site-areatot largestPatch%area = largestPatch%area + (area_site-areatot) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 4c5286bb0a..adb803e639 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -241,9 +241,6 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in) ! Patch dynamics sub-routines: fusion, new patch creation (spwaning), termination. !********************************************************************************* - ! zero site-level patch area error tracking - currentSite%patch_area_conservation_error = 0._r8 - ! make new patches from disturbed land if ( hlm_use_ed_st3.eq.ifalse ) then call spawn_patches(currentSite, bc_in) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index a2b5df5892..2b2e9fd11a 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -771,7 +771,6 @@ module EDTypesMod real(r8) :: disturbance_rates_secondary_to_secondary(N_DIST_TYPES) ! actual disturbance rates from secondary patches to secondary patches [m2/m2/day] real(r8) :: potential_disturbance_rates(N_DIST_TYPES) ! "potential" disturb rates (i.e. prior to the "which is most" logic) [m2/m2/day] real(r8) :: primary_land_patchfusion_error ! error term in total area of primary patches associated with patch fusion [m2/m2/day] - real(r8) :: patch_area_conservation_error ! error term in total area of patches associated with patch area conservation [m2/m2/day] end type ed_site_type diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index c6a916a4a1..478d4a18de 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -191,7 +191,6 @@ module FatesHistoryInterfaceMod integer :: ih_canopy_biomass_si integer :: ih_understory_biomass_si - integer :: ih_patch_area_conservation_error_si integer :: ih_primaryland_fusion_error_si integer :: ih_disturbance_rate_p2p_si integer :: ih_disturbance_rate_p2s_si @@ -1790,7 +1789,6 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_agb_si => this%hvars(ih_agb_si)%r81d, & hio_canopy_biomass_si => this%hvars(ih_canopy_biomass_si)%r81d, & hio_understory_biomass_si => this%hvars(ih_understory_biomass_si)%r81d, & - hio_patch_area_conservation_error_si => this%hvars(ih_patch_area_conservation_error_si)%r81d, & hio_primaryland_fusion_error_si => this%hvars(ih_primaryland_fusion_error_si)%r81d, & hio_disturbance_rate_p2p_si => this%hvars(ih_disturbance_rate_p2p_si)%r81d, & hio_disturbance_rate_p2s_si => this%hvars(ih_disturbance_rate_p2s_si)%r81d, & @@ -2066,9 +2064,6 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) ! error in primary lands from patch fusion hio_primaryland_fusion_error_si = sites(s)%primary_land_patchfusion_error - ! total error in patch area fixed during conservation check - hio_patch_area_conservation_error_si = sites(s)%patch_area_conservation_error - ! output site-level disturbance rates hio_disturbance_rate_p2p_si = sum(sites(s)%disturbance_rates_primary_to_primary(1:N_DIST_TYPES)) hio_disturbance_rate_p2s_si = sum(sites(s)%disturbance_rates_primary_to_secondary(1:N_DIST_TYPES)) @@ -4254,11 +4249,6 @@ subroutine define_history_vars(this, initialize_variables) ivar=ivar, initialize=initialize_variables, index = ih_understory_biomass_si ) ! disturbance rates - call this%set_history_var(vname='PATCH_AREA_CONSERVATION_ERROR', units='m2 m-2 d-1', & - long='Error in total patch area fixed to ensure conservation', use_default='active', & - avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & - ivar=ivar, initialize=initialize_variables, index = ih_patch_area_conservation_error_si ) - call this%set_history_var(vname='PRIMARYLAND_PATCHFUSION_ERROR', units='m2 m-2 d-1', & long='Error in total primary lands associated with patch fusion', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & From a9bb15ed140857b53e552fa2865b18a1cf9283ab Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 9 Jun 2020 13:50:36 -0600 Subject: [PATCH 18/36] Revert "added diagnostic history variable of primary forest rates from BCs" This reverts commit 29dabcca2ac2baff7b5ebb861fd571f51f6a69af. --- main/FatesHistoryInterfaceMod.F90 | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 478d4a18de..fe7687dfac 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -199,7 +199,6 @@ module FatesHistoryInterfaceMod integer :: ih_logging_disturbance_rate_si integer :: ih_fall_disturbance_rate_si integer :: ih_potential_disturbance_rate_si - integer :: ih_harvest_primary_bcin_si ! Indices to site by size-class by age variables integer :: ih_nplant_si_scag @@ -1639,7 +1638,7 @@ end subroutine update_history_cbal ! ==================================================================================== - subroutine update_history_dyn(this,nc,nsites,sites,bc_in) + subroutine update_history_dyn(this,nc,nsites,sites) ! --------------------------------------------------------------------------------- ! This is the call to update the history IO arrays that are expected to only change @@ -1666,7 +1665,6 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) integer , intent(in) :: nc ! clump index integer , intent(in) :: nsites type(ed_site_type) , intent(inout), target :: sites(nsites) - type(bc_in_type) , intent(in) :: bc_in(nsites) ! Locals type(litter_type), pointer :: litt_c ! Pointer to the carbon12 litter pool @@ -1797,7 +1795,6 @@ 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_potential_disturbance_rate_si => this%hvars(ih_potential_disturbance_rate_si)%r81d, & - hio_harvest_primary_bcin_si => this%hvars(ih_harvest_primary_bcin_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, & @@ -2083,8 +2080,6 @@ subroutine update_history_dyn(this,nc,nsites,sites,bc_in) hio_potential_disturbance_rate_si = sum(sites(s)%potential_disturbance_rates(1:N_DIST_TYPES)) - hio_harvest_primary_bcin_si = sum(bc_in(s)%hlm_harvest(1:2)) - ipa = 0 cpatch => sites(s)%oldest_patch do while(associated(cpatch)) @@ -4289,11 +4284,6 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_potential_disturbance_rate_si ) - call this%set_history_var(vname='HARVEST_PRIMARY_RATE_BCIN', units='m2 m-2 d-1', & - long='harvest disturbance rate input from HLM', use_default='active', & - avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & - ivar=ivar, initialize=initialize_variables, index = ih_harvest_primary_bcin_si ) - ! Canopy Resistance call this%set_history_var(vname='C_STOMATA', units='umol m-2 s-1', & From 3d32b1abb1ce8ff11edf6807a786639fe31364fc Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 9 Jun 2020 21:43:05 -0600 Subject: [PATCH 19/36] added site-level diagnostic variable for carbon flux of harvested wood --- biogeochem/EDPatchDynamicsMod.F90 | 14 ++++++++++++++ main/EDTypesMod.F90 | 1 + main/FatesHistoryInterfaceMod.F90 | 9 +++++++++ 3 files changed, 24 insertions(+) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 33c949327d..10eae42f75 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -82,6 +82,8 @@ module EDPatchDynamicsMod use FatesInterfaceTypesMod, only : hlm_parteh_mode use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp + use SFParamsMod, only : SF_VAL_CWD_FRAC + use EDParamsMod, only : logging_event_code ! CIME globals use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) @@ -186,6 +188,8 @@ subroutine disturbance_rates( site_in, bc_in) ! first calculate the fractino of the site that is primary land call get_frac_site_primary(site_in, frac_site_primary) + site_in%harvest_carbon_flux = 0._r8 + currentPatch => site_in%oldest_patch do while (associated(currentPatch)) @@ -220,6 +224,15 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort%lmort_infra = lmort_infra currentCohort%l_degrad = l_degrad + if (currentCohort%canopy_layer==1) then + site_in%harvest_carbon_flux = site_in%harvest_carbon_flux + & + currentCohort%lmort_direct * & + ( currentCohort%prt%GetState(sapw_organ, all_carbon_elements) + & + currentCohort%prt%GetState(struct_organ, all_carbon_elements)) * & + EDPftvarcon_inst%allom_agb_frac(currentCohort%pft) * & + SF_val_CWD_frac(ncwd) + endif + currentCohort => currentCohort%taller end do currentPatch%disturbance_mode = fates_unset_int @@ -2370,6 +2383,7 @@ subroutine fuse_patches( csite, bc_in ) end do ! i_disttype loop + currentPatch => currentSite%youngest_patch do while(associated(currentPatch)) if (currentPatch%anthro_disturbance_label .eq. primaryforest) then diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 2b2e9fd11a..74934befac 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -771,6 +771,7 @@ module EDTypesMod real(r8) :: disturbance_rates_secondary_to_secondary(N_DIST_TYPES) ! actual disturbance rates from secondary patches to secondary patches [m2/m2/day] real(r8) :: potential_disturbance_rates(N_DIST_TYPES) ! "potential" disturb rates (i.e. prior to the "which is most" logic) [m2/m2/day] real(r8) :: primary_land_patchfusion_error ! error term in total area of primary patches associated with patch fusion [m2/m2/day] + real(r8) :: harvest_carbon_flux ! diagnostic site level flux of carbon as harvested plants [kg C / m2 / day] end type ed_site_type diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index fe7687dfac..24eef8e9ac 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -199,6 +199,7 @@ module FatesHistoryInterfaceMod integer :: ih_logging_disturbance_rate_si integer :: ih_fall_disturbance_rate_si integer :: ih_potential_disturbance_rate_si + integer :: ih_harvest_carbonflux_si ! Indices to site by size-class by age variables integer :: ih_nplant_si_scag @@ -1795,6 +1796,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) 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_potential_disturbance_rate_si => this%hvars(ih_potential_disturbance_rate_si)%r81d, & + hio_harvest_carbonflux_si => this%hvars(ih_harvest_carbonflux_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, & @@ -2080,6 +2082,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_potential_disturbance_rate_si = sum(sites(s)%potential_disturbance_rates(1:N_DIST_TYPES)) + hio_harvest_carbonflux_si = sites(s)%harvest_carbon_flux + ipa = 0 cpatch => sites(s)%oldest_patch do while(associated(cpatch)) @@ -4284,6 +4288,11 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_potential_disturbance_rate_si ) + call this%set_history_var(vname='HARVEST_CARBON_FLUX', units='kg C m-2 d-1', & + long='Harvest carbon flux', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_harvest_carbonflux_si ) + ! Canopy Resistance call this%set_history_var(vname='C_STOMATA', units='umol m-2 s-1', & From 0cd417edaff979b4869d353061c3cb8508441df6 Mon Sep 17 00:00:00 2001 From: adivi Date: Fri, 12 Jun 2020 16:10:10 -0700 Subject: [PATCH 20/36] Updated fates logging calc to account for all secondary land creation Conversion to secondary land due to harvest is based on the area harvested, not just the logged crown. Also, the harvest_carbon_flux diagnostic has been updated to align with the wood_product output. Note that these are in different units and that monthly average output during harvest month is different for each. Changes to be committed: modified: biogeochem/EDLoggingMortalityMod.F90 modified: biogeochem/EDPatchDynamicsMod.F90 --- biogeochem/EDLoggingMortalityMod.F90 | 12 +++++++----- biogeochem/EDPatchDynamicsMod.F90 | 10 ++++++---- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 28a6fef229..4cf4cd0a6c 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -210,10 +210,10 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! todo: add a logging_dbhmax parameter, and probably lower the dbhmin one to 30 cm ! todo: change the default logging_event_code to 1 september (-244) - ! todo: change the default logging_direct_frac to 0.7, which is closer to a clearcut event - ! todo: check secondary forest creation + ! todo: change the default logging_direct_frac to 1.0 for cmip inputs ! todo: check outputs against the LUH2 carbon data ! todo: eventually set up distinct harvest practices, each with a set of input paramaeters + ! todo: implement harvested carbon inputs if (logging_time) then ! Pass logging rates to cohort level @@ -247,14 +247,16 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & call endrun(msg=errMsg(sourcefile, __LINE__)) endif + ! transfer of area to secondary land is based on overall area affected, not just logged crown area + ! l_degrad accounts for the affected area between logged crowns if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees if (dbh >= logging_dbhmin ) then lmort_direct = harvest_rate * logging_direct_frac - l_degrad = 0._r8 + l_degrad = harvest_rate * (1._r8 - logging_direct_frac) else lmort_direct = 0.0_r8 - l_degrad = harvest_rate * logging_direct_frac + l_degrad = harvest_rate end if if (dbh >= logging_dbhmax_infra) then @@ -281,7 +283,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & lmort_direct = 0.0_r8 lmort_collateral = 0.0_r8 lmort_infra = harvest_rate * logging_mechanical_frac - l_degrad = harvest_rate * logging_direct_frac + l_degrad = harvest_rate end if else lmort_direct = 0.0_r8 diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 10eae42f75..11aab52140 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -84,6 +84,7 @@ module EDPatchDynamicsMod use PRTGenericMod, only : prt_cnp_flex_allom_hyp use SFParamsMod, only : SF_VAL_CWD_FRAC use EDParamsMod, only : logging_event_code + use EDParamsMod, only : logging_export_frac ! CIME globals use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) @@ -223,14 +224,15 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort%lmort_collateral = lmort_collateral currentCohort%lmort_infra = lmort_infra currentCohort%l_degrad = l_degrad - - if (currentCohort%canopy_layer==1) then + + ! estimate the wood product (trunk_product_site) + if (currentCohort%canopy_layer>=1) then site_in%harvest_carbon_flux = site_in%harvest_carbon_flux + & - currentCohort%lmort_direct * & + currentCohort%lmort_direct * currentCohort%n * & ( currentCohort%prt%GetState(sapw_organ, all_carbon_elements) + & currentCohort%prt%GetState(struct_organ, all_carbon_elements)) * & EDPftvarcon_inst%allom_agb_frac(currentCohort%pft) * & - SF_val_CWD_frac(ncwd) + SF_val_CWD_frac(ncwd) * logging_export_frac endif currentCohort => currentCohort%taller From 4f9a485ab34103197f201e032c11f4f38ae6d2cd Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 17 Jun 2020 17:01:25 -0600 Subject: [PATCH 21/36] added dbhmax logic to logging code --- biogeochem/EDLoggingMortalityMod.F90 | 6 ++++-- main/EDParamsMod.F90 | 1 - 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 5953212982..1f67902069 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -33,6 +33,7 @@ module EDLoggingMortalityMod use EDParamsMod , only : logging_export_frac use EDParamsMod , only : logging_event_code use EDParamsMod , only : logging_dbhmin + use EDParamsMod , only : logging_dbhmax use EDParamsMod , only : logging_collateral_frac use EDParamsMod , only : logging_direct_frac use EDParamsMod , only : logging_mechanical_frac @@ -207,7 +208,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! Local variables real(r8) :: harvest_rate ! the final harvest rate to apply to this cohort today - ! todo: add a logging_dbhmax parameter, and probably lower the dbhmin one to 30 cm + ! todo: probably lower the dbhmin default value to 30 cm ! todo: change the default logging_event_code to 1 september (-244) ! todo: change the default logging_direct_frac to 1.0 for cmip inputs ! todo: check outputs against the LUH2 carbon data @@ -250,7 +251,8 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! l_degrad accounts for the affected area between logged crowns if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees - if (dbh >= logging_dbhmin ) then + if (dbh >= logging_dbhmin .and. .not. & + ((logging_dbhmax < fates_check_param_set) .and. (dbh < logging_dbhmax )) ) then lmort_direct = harvest_rate * logging_direct_frac l_degrad = harvest_rate * (1._r8 - logging_direct_frac) else diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 883f23a605..440c20c921 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -126,7 +126,6 @@ module EDParamsMod real(r8),protected,public :: logging_dbhmax ! Maximum dbh at which logging is applied (cm) ! Typically associated with fire suppression - ! (THIS PARAMETER IS NOT USED YET) character(len=param_string_length),parameter,public :: logging_name_dbhmax = "fates_logging_dbhmax" From a317ef5adabed57137e6d512d844ae34dc3fd1e7 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 17 Jun 2020 22:44:22 -0600 Subject: [PATCH 22/36] bugfix and changed value of use_this_pft when not on restart to 1. --- biogeochem/EDLoggingMortalityMod.F90 | 1 + main/FatesRestartInterfaceMod.F90 | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 1f67902069..7ebe2dd4aa 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -60,6 +60,7 @@ module EDLoggingMortalityMod use FatesAllometryMod , only : set_root_fraction use FatesConstantsMod , only : primaryforest, secondaryforest, secondary_age_threshold use FatesConstantsMod , only : fates_tiny + use FatesConstantsMod, only : fates_check_param_set implicit none private diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 164a1cae35..6f32df2382 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -1008,7 +1008,7 @@ subroutine define_restart_vars(this, initialize_variables) call this%set_restart_var(vname='fates_use_this_pft', vtype=cohort_int, & !should this be cohort_int as above? long_name='in fixed biogeog mode, is pft in gridcell?', & - units='0/1', flushval = flushzero, & + units='0/1', flushval = flushone, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_use_this_pft_sift) call this%set_restart_var(vname='fates_area_pft', vtype=cohort_r8, & From c2332b92a16a1e199b0711493c83fa37610c28a4 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 6 May 2020 10:16:31 -0600 Subject: [PATCH 23/36] added PFTcanopycrownarea variable --- main/FatesHistoryInterfaceMod.F90 | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 24eef8e9ac..e84afdd18e 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -433,7 +433,7 @@ module FatesHistoryInterfaceMod integer :: ih_recruitment_si_pft integer :: ih_mortality_si_pft integer :: ih_crownarea_si_pft - + integer :: ih_canopycrownarea_si_pft ! indices to (site x patch-age) variables integer :: ih_area_si_age @@ -1752,6 +1752,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_recruitment_si_pft => this%hvars(ih_recruitment_si_pft)%r82d, & hio_mortality_si_pft => this%hvars(ih_mortality_si_pft)%r82d, & hio_crownarea_si_pft => this%hvars(ih_crownarea_si_pft)%r82d, & + hio_canopycrownarea_si_pft => this%hvars(ih_canopycrownarea_si_pft)%r82d, & hio_nesterov_fire_danger_si => this%hvars(ih_nesterov_fire_danger_si)%r81d, & hio_spitfire_ros_si => this%hvars(ih_spitfire_ros_si)%r81d, & hio_fire_ros_area_product_si=> this%hvars(ih_fire_ros_area_product_si)%r81d, & @@ -2246,6 +2247,12 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_crownarea_si_pft(io_si, ft) = hio_crownarea_si_pft(io_si, ft) + & ccohort%c_area + if (ccohort%canopy_layer .eq. 1) then + ! Update PFT canopy crown area + hio_canopycrownarea_si_pft(io_si, ft) = hio_canopycrownarea_si_pft(io_si, ft) + & + ccohort%c_area + end if + ! update total biomass per age bin hio_biomass_si_age(io_si,cpatch%age_class) = hio_biomass_si_age(io_si,cpatch%age_class) & + total_c * ccohort%n * AREA_INV @@ -3936,6 +3943,11 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_crownarea_si_pft ) + call this%set_history_var(vname='PFTcanopycrownarea', units='m2/ha', & + long='total PFT-level canopy-layer crown area', use_default='inactive', & + avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_canopycrownarea_si_pft ) + call this%set_history_var(vname='PFTnindivs', units='indiv / m2', & long='total PFT level number of individuals', use_default='active', & avgflag='A', vtype=site_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & From b78f410bdcf8014046de98338fdd9d9c3c1c52f5 Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 22 Jun 2020 17:04:26 -0600 Subject: [PATCH 24/36] adding more diagnostics to dump_patch --- main/EDTypesMod.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index e5c720580a..0e843df123 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -948,6 +948,8 @@ subroutine dump_patch(cpatch) write(fates_log(),*) 'pa%c_stomata = ',cpatch%c_stomata write(fates_log(),*) 'pa%c_lblayer = ',cpatch%c_lblayer write(fates_log(),*) 'pa%disturbance_rate = ',cpatch%disturbance_rate + write(fates_log(),*) 'pa%disturbance_rates = ',cpatch%disturbance_rates(:) + write(fates_log(),*) 'pa%anthro_disturbance_label = ',cpatch%anthro_disturbance_label write(fates_log(),*) '----------------------------------------' do el = 1,num_elements write(fates_log(),*) 'element id: ',element_list(el) From 2905a9ba33e97319639dc5806ed7d212bccbb9f2 Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 22 Jun 2020 17:05:11 -0600 Subject: [PATCH 25/36] adding more logic in fuse_patches to handle case where only tiny amounts of primary or secondary patch land exists --- biogeochem/EDPatchDynamicsMod.F90 | 57 ++++++++++++++++++++++++------- 1 file changed, 44 insertions(+), 13 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 85897faf9e..85773838fc 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -495,6 +495,7 @@ subroutine spawn_patches( currentSite, bc_in) if(currentPatch%disturbance_rate>1.0_r8) then write(fates_log(),*) 'patch disturbance rate > 1 ?',currentPatch%disturbance_rate + call dump_patch(currentPatch) call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -2581,6 +2582,7 @@ subroutine terminate_patches(currentSite) integer, parameter :: max_cycles = 10 ! After 10 loops through ! You should had fused integer :: count_cycles + logical :: gotfused real(r8) areatot ! variable for checking whether the total patch area is wrong. !--------------------------------------------------------------------- @@ -2601,6 +2603,8 @@ subroutine terminate_patches(currentSite) if ( .not.associated(currentPatch,currentSite%youngest_patch) .or. & currentPatch%area <= min_patch_area_forced ) then + gotfused = .false. + if(associated(currentPatch%older) )then if(debug) & @@ -2612,25 +2616,52 @@ subroutine terminate_patches(currentSite) ! it will be returned by the subroutine as de-referenced olderPatch => currentPatch%older - call fuse_2_patches(currentSite, olderPatch, currentPatch) - - ! The fusion process has updated the "older" pointer on currentPatch - ! for us. + + if (currentPatch%anthro_disturbance_label .eq. olderPatch%anthro_disturbance_label) then + + call fuse_2_patches(currentSite, olderPatch, currentPatch) - ! This logic checks to make sure that the younger patch is not the youngest - ! patch. As mentioned earlier, we try not to fuse it. + ! The fusion process has updated the "older" pointer on currentPatch + ! for us. - elseif( associated(currentPatch%younger) ) then + ! This logic checks to make sure that the younger patch is not the youngest + ! patch. As mentioned earlier, we try not to fuse it. + + gotfused = .true. + else + if (count_cycles .gt. 0) then + ! if we're having an incredibly hard time fusing patches because of their differing anthropogenic disturbance labels, + ! since the size is so small, let's sweep the problem under the rug and change the tiny patch's label to that of its older sibling + currentPatch%anthro_disturbance_label = olderPatch%anthro_disturbance_label + call fuse_2_patches(currentSite, olderPatch, currentPatch) + gotfused = .true. + endif + endif + endif + + if( .not. gotfused .and. associated(currentPatch%younger) ) then if(debug) & - write(fates_log(),*) 'fusing to younger patch because oldest one is too small', & - currentPatch%area + write(fates_log(),*) 'fusing to younger patch because oldest one is too small', & + currentPatch%area youngerPatch => currentPatch%younger - call fuse_2_patches(currentSite, youngerPatch, currentPatch) - - ! The fusion process has updated the "younger" pointer on currentPatch - + + if (currentPatch%anthro_disturbance_label .eq. youngerPatch% anthro_disturbance_label) then + + call fuse_2_patches(currentSite, youngerPatch, currentPatch) + + ! The fusion process has updated the "younger" pointer on currentPatch + + else + if (count_cycles .gt. 0) then + ! if we're having an incredibly hard time fusing patches because of their differing anthropogenic disturbance labels, + ! since the size is so small, let's sweep the problem under the rug and change the tiny patch's label to that of its younger sibling + currentPatch%anthro_disturbance_label = youngerPatch%anthro_disturbance_label + call fuse_2_patches(currentSite, youngerPatch, currentPatch) + gotfused = .true. + endif + endif endif endif endif From 5e6fa21f6d44dfe3d1813d3bc149e3525c668e96 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 23 Jun 2020 11:06:30 -0600 Subject: [PATCH 26/36] bugfix on history-outputting of site-level diagnostic variables --- main/FatesHistoryInterfaceMod.F90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index e84afdd18e..eba5b7de37 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -2062,28 +2062,28 @@ subroutine update_history_dyn(this,nc,nsites,sites) end if ! error in primary lands from patch fusion - hio_primaryland_fusion_error_si = sites(s)%primary_land_patchfusion_error + hio_primaryland_fusion_error_si(io_si) = sites(s)%primary_land_patchfusion_error ! output site-level disturbance rates - hio_disturbance_rate_p2p_si = sum(sites(s)%disturbance_rates_primary_to_primary(1:N_DIST_TYPES)) - hio_disturbance_rate_p2s_si = sum(sites(s)%disturbance_rates_primary_to_secondary(1:N_DIST_TYPES)) - hio_disturbance_rate_s2s_si = sum(sites(s)%disturbance_rates_secondary_to_secondary(1:N_DIST_TYPES)) + hio_disturbance_rate_p2p_si(io_si) = sum(sites(s)%disturbance_rates_primary_to_primary(1:N_DIST_TYPES)) + hio_disturbance_rate_p2s_si(io_si) = sum(sites(s)%disturbance_rates_primary_to_secondary(1:N_DIST_TYPES)) + hio_disturbance_rate_s2s_si(io_si) = sum(sites(s)%disturbance_rates_secondary_to_secondary(1:N_DIST_TYPES)) - hio_fire_disturbance_rate_si = sites(s)%disturbance_rates_primary_to_primary(dtype_ifire) + & + hio_fire_disturbance_rate_si(io_si) = sites(s)%disturbance_rates_primary_to_primary(dtype_ifire) + & sites(s)%disturbance_rates_primary_to_secondary(dtype_ifire) + & sites(s)%disturbance_rates_secondary_to_secondary(dtype_ifire) - hio_logging_disturbance_rate_si = sites(s)%disturbance_rates_primary_to_primary(dtype_ilog) + & + hio_logging_disturbance_rate_si(io_si) = sites(s)%disturbance_rates_primary_to_primary(dtype_ilog) + & sites(s)%disturbance_rates_primary_to_secondary(dtype_ilog) + & sites(s)%disturbance_rates_secondary_to_secondary(dtype_ilog) - hio_fall_disturbance_rate_si = sites(s)%disturbance_rates_primary_to_primary(dtype_ifall) + & + hio_fall_disturbance_rate_si(io_si) = sites(s)%disturbance_rates_primary_to_primary(dtype_ifall) + & sites(s)%disturbance_rates_primary_to_secondary(dtype_ifall) + & sites(s)%disturbance_rates_secondary_to_secondary(dtype_ifall) - hio_potential_disturbance_rate_si = sum(sites(s)%potential_disturbance_rates(1:N_DIST_TYPES)) + hio_potential_disturbance_rate_si(io_si) = sum(sites(s)%potential_disturbance_rates(1:N_DIST_TYPES)) - hio_harvest_carbonflux_si = sites(s)%harvest_carbon_flux + hio_harvest_carbonflux_si(io_si) = sites(s)%harvest_carbon_flux ipa = 0 cpatch => sites(s)%oldest_patch From cbc4ca001f77edf68b6651082d8150fd9e1e40b2 Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 25 Jun 2020 17:38:29 -0600 Subject: [PATCH 27/36] fixing conservation of patch%age_since_anthro_disturbance during patch fusion --- biogeochem/EDPatchDynamicsMod.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 85773838fc..4f7ab336c3 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -2436,6 +2436,8 @@ subroutine fuse_2_patches(csite, dp, rp) inv_sum_area = 1.0_r8/(dp%area + rp%area) rp%age = (dp%age * dp%area + rp%age * rp%area) * inv_sum_area + rp%age_since_anthro_disturbance = (dp%age_since_anthro_disturbance * dp%area & + + rp%age_since_anthro_disturbance * rp%area) * inv_sum_area rp%age_class = get_age_class_index(rp%age) From bb04f896c8179d17fd21d969cebd072ee93f7516 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Wed, 1 Jul 2020 11:27:24 -0700 Subject: [PATCH 28/36] Apply suggestions from code review Co-authored-by: Rosie Fisher --- biogeochem/EDLoggingMortalityMod.F90 | 17 +++++++++-------- biogeochem/EDPatchDynamicsMod.F90 | 10 +++++----- main/FatesHistoryInterfaceMod.F90 | 2 +- 3 files changed, 15 insertions(+), 14 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 7ebe2dd4aa..7d16870600 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -87,8 +87,8 @@ module EDLoggingMortalityMod public :: logging_time public :: IsItLoggingTime public :: get_harvest_rate_area - integer, parameter, public :: hlm_harvest_area_fraction = 1 - integer, parameter, public :: hlm_harvest_carbon = 2 + integer, parameter, public :: hlm_harvest_area_fraction = 1 ! Code for harvesting by area + integer, parameter, public :: hlm_harvest_carbon = 2 ! Code for harvesting based on carbon extracted. contains @@ -225,19 +225,20 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & harvest_rate = 1._r8 else if (hlm_use_lu_harvest == hlm_harvest_area_fraction) then + ! We are harvesting based on areal fraction, not carbon/biomass terms. ! 1=use area fraction from hlm ! combine forest and non-forest fracs and then apply: ! primary and secondary area fractions to the logging rates, which are fates parameters ! Definitions of the underlying harvest land category variables ! these are hardcoded to match the LUH input data via landuse.timseries file (see dynHarvestMod) - ! these are fraction of vegetated area harvested, split into five land category variables + ! these are fractions of vegetated area harvested, split into five land category variables ! HARVEST_VH1 = harvest from primary forest ! HARVEST_VH2 = harvest from primary non-forest ! HARVEST_SH1 = harvest from secondary mature forest ! HARVEST_SH2 = harvest from secondary young forest ! HARVEST_SH3 = harvest from secondary non-forest (assume this is young for biomass) - + ! Get the area-based harvest rates based on info passed to FATES from the bioundary condition call get_harvest_rate_area (use_history, hlm_harvest_catnames, hlm_harvest, & frac_site_primary, secondary_age, harvest_rate) @@ -255,7 +256,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & if (dbh >= logging_dbhmin .and. .not. & ((logging_dbhmax < fates_check_param_set) .and. (dbh < logging_dbhmax )) ) then lmort_direct = harvest_rate * logging_direct_frac - l_degrad = harvest_rate * (1._r8 - logging_direct_frac) + l_degrad = harvest_rate * (1._r8 - logging_direct_frac) ! fraction passed to 'degraded' forest. else lmort_direct = 0.0_r8 l_degrad = harvest_rate @@ -323,7 +324,7 @@ subroutine get_harvest_rate_area (use_history, hlm_harvest_catnames, hlm_harvest ! Parameters real(r8), parameter :: months_per_year = 12.0 - ! determine the annual hlm harvest rate for the current cohort based on patch info + ! Loop around harvest categories to determine the annual hlm harvest rate for the current cohort based on patch history info harvest_rate = 0._r8 do h_index = 1,hlm_num_lu_harvest_cats if (use_history .eq. primaryforest) then @@ -345,8 +346,8 @@ subroutine get_harvest_rate_area (use_history, hlm_harvest_catnames, hlm_harvest endif end do - ! normalize by site-level fractio primary or secondary forest fraction - ! since harvest_rate is specified in fraction of the gridcell + ! Normalize by site-level primary or secondary forest fraction + ! since harvest_rate is specified as a fraction of the gridcell ! also need to put a cap so as not to harvest more primary or secondary area than there is in a gridcell if (use_history .eq. primaryforest) then if (frac_site_primary .gt. fates_tiny) then diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 4f7ab336c3..2c0b128a35 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -247,7 +247,7 @@ subroutine disturbance_rates( site_in, bc_in) ! zero the diagnostic disturbance rate fields site_in%potential_disturbance_rates(1:N_DIST_TYPES) = 0._r8 - ! summarize how open the canopy is prior to resolving the disturbance + ! Recalculate total canopy area prior to resolving the disturbance currentPatch => site_in%oldest_patch do while (associated(currentPatch)) currentPatch%total_canopy_area = 0._r8 @@ -309,14 +309,14 @@ subroutine disturbance_rates( site_in, bc_in) ! equivalent to the fradction loged to account for transfer of interstitial ground area to new secondary lands if ( logging_time .and. & (currentPatch%area - currentPatch%total_canopy_area) .gt. fates_tiny ) then - +! The canopy is NOT closed. call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, bc_in%hlm_harvest, & frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate) currentPatch%disturbance_rates(dtype_ilog) = currentPatch%disturbance_rates(dtype_ilog) + & (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area endif - +! Sum of disturbance rates for different classes of disturbance across all patches in this site. do i_dist = 1,N_DIST_TYPES site_in%potential_disturbance_rates(i_dist) = site_in%potential_disturbance_rates(i_dist) + & currentPatch%disturbance_rates(i_dist) * currentPatch%area * AREA_INV @@ -1120,7 +1120,7 @@ subroutine spawn_patches( currentSite, bc_in) currentSite%oldest_patch => new_patch_primary endif else - ! the case where there are no secondary patches at the start of the LL (prior logic) + ! the case where there are no secondary patches at the start of the linked list (prior logic) new_patch_primary%older => currentPatch new_patch_primary%younger => null() currentPatch%younger => new_patch_primary @@ -2634,6 +2634,7 @@ subroutine terminate_patches(currentSite) if (count_cycles .gt. 0) then ! if we're having an incredibly hard time fusing patches because of their differing anthropogenic disturbance labels, ! since the size is so small, let's sweep the problem under the rug and change the tiny patch's label to that of its older sibling + ! and then allow them to fuse together. currentPatch%anthro_disturbance_label = olderPatch%anthro_disturbance_label call fuse_2_patches(currentSite, olderPatch, currentPatch) gotfused = .true. @@ -2905,4 +2906,3 @@ subroutine get_frac_site_primary(site_in, frac_site_primary) end subroutine get_frac_site_primary end module EDPatchDynamicsMod - diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index eba5b7de37..a022e64624 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -4290,7 +4290,7 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_logging_disturbance_rate_si ) - call this%set_history_var(vname='DISTURBANCE_RATE_FALL', units='m2 m-2 d-1', & + call this%set_history_var(vname='DISTURBANCE_RATE_TREEFALL', units='m2 m-2 d-1', & long='Disturbance rate from treefall', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_fall_disturbance_rate_si ) From dafd55ff6aeb6aaf9a0535107f10a4d21d0a54c8 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Wed, 1 Jul 2020 14:21:12 -0700 Subject: [PATCH 29/36] Apply suggestions from code review --- main/FatesConstantsMod.F90 | 2 +- main/FatesInterfaceTypesMod.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 0aaa946a8d..1fd278240b 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -34,7 +34,7 @@ module FatesConstantsMod integer, parameter, public :: n_anthro_disturbance_categories = 2 integer, parameter, public :: primaryforest = 1 integer, parameter, public :: secondaryforest = 2 - integer, parameter, public :: secondary_age_threshold = 94 ! less than this value is young + real(fates_r8), parameter, public :: secondary_age_threshold = 94._fates_r8 ! less than this value is young secondary land ! based on average age of global ! secondary 1900s land in hurtt-2011 diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index 2488d7db6b..6d47909fa3 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -97,7 +97,7 @@ module FatesInterfaceTypesMod ! If 1 or 2, it automatically sets ! hlm_use_logging to 1 - integer, public :: hlm_num_lu_harvest_cats ! number of hlm harvest categories + integer, public :: hlm_num_lu_harvest_cats ! number of hlm harvest categories (e.g. primary forest harvest, secondary young forest harvest, etc.) ! this is the first dimension of: ! harvest_rates in dynHarvestMod ! bc_in%hlm_harvest and bc_in%hlm_harvest_catnames From 7afdb0475c46210691402e728866cd9a1d848ed7 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 1 Jul 2020 16:37:47 -0600 Subject: [PATCH 30/36] made a bunch of minor changes based on code review --- biogeochem/EDLoggingMortalityMod.F90 | 36 ++++++++++++-------------- biogeochem/EDMortalityFunctionsMod.F90 | 2 +- biogeochem/EDPatchDynamicsMod.F90 | 6 ++--- main/FatesConstantsMod.F90 | 5 +++- main/FatesInterfaceMod.F90 | 2 +- main/FatesInterfaceTypesMod.F90 | 6 ++--- 6 files changed, 28 insertions(+), 29 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 7d16870600..6d2a96fcb8 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -60,6 +60,7 @@ module EDLoggingMortalityMod use FatesAllometryMod , only : set_root_fraction use FatesConstantsMod , only : primaryforest, secondaryforest, secondary_age_threshold use FatesConstantsMod , only : fates_tiny + use FatesConstantsMod , only : months_per_year use FatesConstantsMod, only : fates_check_param_set implicit none @@ -185,17 +186,17 @@ end subroutine IsItLoggingTime subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & lmort_collateral,lmort_infra, l_degrad, & - hlm_harvest, hlm_harvest_catnames, & - use_history, secondary_age, & + hlm_harvest_rates, hlm_harvest_catnames, & + patch_anthro_disturbance_label, secondary_age, & frac_site_primary) ! Arguments integer, intent(in) :: pft_i ! pft index real(r8), intent(in) :: dbh ! diameter at breast height (cm) integer, intent(in) :: canopy_layer ! canopy layer of this cohort - real(r8), intent(in) :: hlm_harvest(:) ! annual harvest rate per hlm category + real(r8), intent(in) :: hlm_harvest_rates(:) ! annual harvest rate per hlm category character(len=64), intent(in) :: hlm_harvest_catnames(:) ! names of hlm harvest categories - integer, intent(in) :: use_history ! patch level anthro_disturbance_label + integer, intent(in) :: patch_anthro_disturbance_label ! patch level anthro_disturbance_label real(r8), intent(in) :: secondary_age ! patch level age_since_anthro_disturbance real(r8), intent(out) :: lmort_direct ! direct (harvestable) mortality fraction real(r8), intent(out) :: lmort_collateral ! collateral damage mortality fraction @@ -239,8 +240,8 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! HARVEST_SH2 = harvest from secondary young forest ! HARVEST_SH3 = harvest from secondary non-forest (assume this is young for biomass) ! Get the area-based harvest rates based on info passed to FATES from the bioundary condition - call get_harvest_rate_area (use_history, hlm_harvest_catnames, hlm_harvest, & - frac_site_primary, secondary_age, harvest_rate) + call get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_catnames, & + hlm_harvest_rates, frac_site_primary, secondary_age, harvest_rate) else if (hlm_use_lu_harvest == hlm_harvest_carbon) then ! 2=use carbon from hlm @@ -299,7 +300,7 @@ end subroutine LoggingMortality_frac ! ============================================================================ - subroutine get_harvest_rate_area (use_history, hlm_harvest_catnames, hlm_harvest, & + subroutine get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_catnames, hlm_harvest_rates, & frac_site_primary, secondary_age, harvest_rate) @@ -310,9 +311,9 @@ subroutine get_harvest_rate_area (use_history, hlm_harvest_catnames, hlm_harvest ! assumes logging_time == true ! Arguments - real(r8), intent(in) :: hlm_harvest(:) ! annual harvest rate per hlm category + real(r8), intent(in) :: hlm_harvest_rates(:) ! annual harvest rate per hlm category character(len=64), intent(in) :: hlm_harvest_catnames(:) ! names of hlm harvest categories - integer, intent(in) :: use_history ! patch level anthro_disturbance_label + integer, intent(in) :: patch_anthro_disturbance_label ! patch level anthro_disturbance_label real(r8), intent(in) :: secondary_age ! patch level age_since_anthro_disturbance real(r8), intent(in) :: frac_site_primary real(r8), intent(out) :: harvest_rate @@ -321,27 +322,24 @@ subroutine get_harvest_rate_area (use_history, hlm_harvest_catnames, hlm_harvest integer :: h_index ! for looping over harvest categories integer :: icode ! Integer equivalent of the event code (parameter file only allows reals) - ! Parameters - real(r8), parameter :: months_per_year = 12.0 - ! Loop around harvest categories to determine the annual hlm harvest rate for the current cohort based on patch history info harvest_rate = 0._r8 do h_index = 1,hlm_num_lu_harvest_cats - if (use_history .eq. primaryforest) then + if (patch_anthro_disturbance_label .eq. primaryforest) then if(hlm_harvest_catnames(h_index) .eq. "HARVEST_VH1" .or. & hlm_harvest_catnames(h_index) .eq. "HARVEST_VH2") then - harvest_rate = harvest_rate + hlm_harvest(h_index) + harvest_rate = harvest_rate + hlm_harvest_rates(h_index) endif - else if (use_history .eq. secondaryforest .and. & + else if (patch_anthro_disturbance_label .eq. secondaryforest .and. & secondary_age >= secondary_age_threshold) then if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then - harvest_rate = harvest_rate + hlm_harvest(h_index) + harvest_rate = harvest_rate + hlm_harvest_rates(h_index) endif - else if (use_history .eq. secondaryforest .and. & + else if (patch_anthro_disturbance_label .eq. secondaryforest .and. & secondary_age < secondary_age_threshold) then if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. & hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then - harvest_rate = harvest_rate + hlm_harvest(h_index) + harvest_rate = harvest_rate + hlm_harvest_rates(h_index) endif endif end do @@ -349,7 +347,7 @@ subroutine get_harvest_rate_area (use_history, hlm_harvest_catnames, hlm_harvest ! Normalize by site-level primary or secondary forest fraction ! since harvest_rate is specified as a fraction of the gridcell ! also need to put a cap so as not to harvest more primary or secondary area than there is in a gridcell - if (use_history .eq. primaryforest) then + if (patch_anthro_disturbance_label .eq. primaryforest) then if (frac_site_primary .gt. fates_tiny) then harvest_rate = min((harvest_rate / frac_site_primary),frac_site_primary) else diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 074194a9a0..07f8ea7024 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -247,7 +247,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr currentCohort%lmort_collateral, & currentCohort%lmort_infra, & currentCohort%l_degrad, & - bc_in%hlm_harvest, & + bc_in%hlm_harvest_rates, & bc_in%hlm_harvest_catnames, & currentCohort%patchptr%anthro_disturbance_label, & currentCohort%patchptr%age_since_anthro_disturbance, & diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 2c0b128a35..2273eee6af 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -213,7 +213,7 @@ subroutine disturbance_rates( site_in, bc_in) call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_layer, & lmort_direct,lmort_collateral,lmort_infra,l_degrad,& - bc_in%hlm_harvest, & + bc_in%hlm_harvest_rates, & bc_in%hlm_harvest_catnames, & currentPatch%anthro_disturbance_label, & currentPatch%age_since_anthro_disturbance, & @@ -310,8 +310,8 @@ subroutine disturbance_rates( site_in, bc_in) if ( logging_time .and. & (currentPatch%area - currentPatch%total_canopy_area) .gt. fates_tiny ) then ! The canopy is NOT closed. - call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, bc_in%hlm_harvest, & - frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate) + call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, & + bc_in%hlm_harvest_rates, frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate) currentPatch%disturbance_rates(dtype_ilog) = currentPatch%disturbance_rates(dtype_ilog) + & (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 1fd278240b..b0c201f4ff 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -137,7 +137,10 @@ module FatesConstantsMod ! Conversion: years per day. assume HLM uses 365 day calendar. ! If we need to link to 365.25-day-calendared HLM, rewire to pass through interface real(fates_r8), parameter, public :: years_per_day = 1.0_fates_r8/365.00_fates_r8 - + + ! Conversion: months per year + real(fates_r8), parameter, public :: months_per_year = 12.0_fates_r8 + ! Physical constants ! universal gas constant [J/K/kmol] diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index bb893c44dc..aea2f67f65 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -363,7 +363,7 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats) ! harvest flag denote data from hlm, ! while the logging flag signifies only that logging is occurring (which could just be FATES logging) if (hlm_use_lu_harvest .gt. 0) then - allocate(bc_in%hlm_harvest(num_lu_harvest_cats)) + allocate(bc_in%hlm_harvest_rates(num_lu_harvest_cats)) allocate(bc_in%hlm_harvest_catnames(num_lu_harvest_cats)) end if diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index 6d47909fa3..7909a237bb 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -100,7 +100,7 @@ module FatesInterfaceTypesMod integer, public :: hlm_num_lu_harvest_cats ! number of hlm harvest categories (e.g. primary forest harvest, secondary young forest harvest, etc.) ! this is the first dimension of: ! harvest_rates in dynHarvestMod - ! bc_in%hlm_harvest and bc_in%hlm_harvest_catnames + ! bc_in%hlm_harvest_rates and bc_in%hlm_harvest_catnames integer, public :: hlm_use_logging ! This flag signals whether or not to use ! the logging module @@ -449,9 +449,7 @@ module FatesInterfaceTypesMod ! Land use ! --------------------------------------------------------------------------------- - logical :: hlm_do_harvest_today ! denotes whether hlm harvest data - ! are available for today - real(r8),allocatable :: hlm_harvest(:) ! annual harvest rate per cat from hlm for a site + real(r8),allocatable :: hlm_harvest_rates(:) ! annual harvest rate per cat from hlm for a site character(len=64), allocatable :: hlm_harvest_catnames(:) ! names of hlm_harvest d1 ! Fixed biogeography mode From 3807a364955414d3dea1787b1884bf6b4e58b67f Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 1 Jul 2020 20:08:53 -0600 Subject: [PATCH 31/36] moved harvest unit specification to bc_in structure --- biogeochem/EDLoggingMortalityMod.F90 | 15 +++++++++------ biogeochem/EDMortalityFunctionsMod.F90 | 1 + biogeochem/EDPatchDynamicsMod.F90 | 1 + main/FatesConstantsMod.F90 | 4 ++++ main/FatesInterfaceTypesMod.F90 | 3 +++ 5 files changed, 18 insertions(+), 6 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 6d2a96fcb8..228076b197 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -61,6 +61,8 @@ module EDLoggingMortalityMod use FatesConstantsMod , only : primaryforest, secondaryforest, secondary_age_threshold use FatesConstantsMod , only : fates_tiny use FatesConstantsMod , only : months_per_year + use FatesConstantsMod , only : hlm_harvest_area_fraction + use FatesConstantsMod , only : hlm_harvest_carbon use FatesConstantsMod, only : fates_check_param_set implicit none @@ -88,8 +90,6 @@ module EDLoggingMortalityMod public :: logging_time public :: IsItLoggingTime public :: get_harvest_rate_area - integer, parameter, public :: hlm_harvest_area_fraction = 1 ! Code for harvesting by area - integer, parameter, public :: hlm_harvest_carbon = 2 ! Code for harvesting based on carbon extracted. contains @@ -187,6 +187,7 @@ end subroutine IsItLoggingTime subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & lmort_collateral,lmort_infra, l_degrad, & hlm_harvest_rates, hlm_harvest_catnames, & + hlm_harvest_units, & patch_anthro_disturbance_label, secondary_age, & frac_site_primary) @@ -196,6 +197,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & integer, intent(in) :: canopy_layer ! canopy layer of this cohort real(r8), intent(in) :: hlm_harvest_rates(:) ! annual harvest rate per hlm category character(len=64), intent(in) :: hlm_harvest_catnames(:) ! names of hlm harvest categories + integer, intent(in) :: hlm_harvest_units ! unit type of hlm harvest rates: [area vs. mass] integer, intent(in) :: patch_anthro_disturbance_label ! patch level anthro_disturbance_label real(r8), intent(in) :: secondary_age ! patch level age_since_anthro_disturbance real(r8), intent(out) :: lmort_direct ! direct (harvestable) mortality fraction @@ -220,12 +222,12 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & if (logging_time) then ! Pass logging rates to cohort level - if (hlm_use_lu_harvest == 0) then + if (hlm_use_lu_harvest == ifalse) then ! 0=use fates logging parameters directly when logging_time == .true. ! this means harvest the whole cohort area harvest_rate = 1._r8 - else if (hlm_use_lu_harvest == hlm_harvest_area_fraction) then + 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 ! combine forest and non-forest fracs and then apply: @@ -239,11 +241,12 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! HARVEST_SH1 = harvest from secondary mature forest ! HARVEST_SH2 = harvest from secondary young forest ! HARVEST_SH3 = harvest from secondary non-forest (assume this is young for biomass) - ! Get the area-based harvest rates based on info passed to FATES from the bioundary condition + + ! Get the area-based harvest rates based on info passed to FATES from the boundary condition call get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_catnames, & hlm_harvest_rates, frac_site_primary, secondary_age, harvest_rate) - else if (hlm_use_lu_harvest == hlm_harvest_carbon) then + else if (hlm_use_lu_harvest == itrue .and. hlm_harvest_units == hlm_harvest_carbon) then ! 2=use carbon from hlm ! not implemented yet write(fates_log(),*) 'HLM harvest carbon data not implemented yet. Exiting.' diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 07f8ea7024..fa0b933fc5 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -249,6 +249,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in, frac_site_pr currentCohort%l_degrad, & bc_in%hlm_harvest_rates, & bc_in%hlm_harvest_catnames, & + bc_in%hlm_harvest_units, & currentCohort%patchptr%anthro_disturbance_label, & currentCohort%patchptr%age_since_anthro_disturbance, & frac_site_primary) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 2273eee6af..9e893cc838 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -215,6 +215,7 @@ subroutine disturbance_rates( site_in, bc_in) lmort_direct,lmort_collateral,lmort_infra,l_degrad,& bc_in%hlm_harvest_rates, & bc_in%hlm_harvest_catnames, & + bc_in%hlm_harvest_units, & currentPatch%anthro_disturbance_label, & currentPatch%age_since_anthro_disturbance, & frac_site_primary) diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index b0c201f4ff..e553918028 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -37,6 +37,10 @@ module FatesConstantsMod real(fates_r8), parameter, public :: secondary_age_threshold = 94._fates_r8 ! less than this value is young secondary land ! based on average age of global ! secondary 1900s land in hurtt-2011 + + ! integer labels for specifying harvest units + integer, parameter, public :: hlm_harvest_area_fraction = 1 ! Code for harvesting by area + integer, parameter, public :: hlm_harvest_carbon = 2 ! Code for harvesting based on carbon extracted. ! Error Tolerances diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index 7909a237bb..2035e1a517 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -450,8 +450,11 @@ module FatesInterfaceTypesMod ! Land use ! --------------------------------------------------------------------------------- real(r8),allocatable :: hlm_harvest_rates(:) ! annual harvest rate per cat from hlm for a site + character(len=64), allocatable :: hlm_harvest_catnames(:) ! names of hlm_harvest d1 + integer :: hlm_harvest_units ! what units are the harvest rates specified in? [area vs carbon] + ! Fixed biogeography mode real(r8), allocatable :: pft_areafrac(:) ! Fractional area of the FATES column occupied by each PFT From d08dfd4665884b8c6b0ed5c6f683de2bbec646ea Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 2 Jul 2020 10:52:38 -0600 Subject: [PATCH 32/36] clarifying some variable names and documentation in spawn_patches --- biogeochem/EDPatchDynamicsMod.F90 | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 9e893cc838..9274d21ba0 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -474,7 +474,7 @@ subroutine spawn_patches( currentSite, bc_in) real(r8) :: leaf_burn_frac ! fraction of leaves burned in fire ! for both woody and grass species real(r8) :: leaf_m ! leaf mass during partial burn calculations - logical :: foundfirstprimary ! logical for finding the first primary forest patch + logical :: found_youngest_primary ! logical for finding the first primary forest patch !--------------------------------------------------------------------- storesmallcohort => null() ! storage of the smallest cohort for insertion routine @@ -1091,18 +1091,19 @@ subroutine spawn_patches( currentSite, bc_in) !*************************/ !** INSERT NEW PATCH(ES) INTO LINKED LIST - !**********`***************/ + !*************************/ if ( site_areadis_primary .gt. nearzero) then currentPatch => currentSite%youngest_patch - ! insert first primary patch after all the secondary patches, if there are any + ! insert new youngest primary patch after all the secondary patches, if there are any. + ! this requires first finding the current youngest primary to insert the new one ahead of if (currentPatch%anthro_disturbance_label .eq. secondaryforest ) then - foundfirstprimary = .false. - do while(associated(currentPatch) .and. .not. foundfirstprimary) + found_youngest_primary = .false. + do while(associated(currentPatch) .and. .not. found_youngest_primary) currentPatch => currentPatch%older if (associated(currentPatch)) then if (currentPatch%anthro_disturbance_label .eq. primaryforest) then - foundfirstprimary = .true. + found_youngest_primary = .true. endif endif end do @@ -1113,7 +1114,8 @@ subroutine spawn_patches( currentSite, bc_in) currentPatch%younger%older => new_patch_primary currentPatch%younger => new_patch_primary else - ! the case where we haven't, and are putting a primary patch at the oldest end of the + ! the case where we haven't, because the patches are all secondaary, + ! and are putting a primary patch at the oldest end of the ! linked list (not sure how this could happen, but who knows...) new_patch_primary%older => null() new_patch_primary%younger => currentSite%oldest_patch From 8d9f189bcf4f696f128b0391c50f930ddcf911ae Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 6 Jul 2020 17:32:25 -0600 Subject: [PATCH 33/36] replacing valeu of cpatch%age_since_anthro_disturbance for primary patches with fates_unset_r8 --- biogeochem/EDPatchDynamicsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 9274d21ba0..99964a5273 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -2003,7 +2003,7 @@ subroutine create_patch(currentSite, new_patch, age, areap, label) if (label .eq. secondaryforest) then new_patch%age_since_anthro_disturbance = age else - new_patch%age_since_anthro_disturbance = -1._r8 ! replace with fates_unset_r8 when possible + new_patch%age_since_anthro_disturbance = fates_unset_r8 endif ! This new value will be generated when the calculate disturbance From 568648ff73c8765ed9615b6e40b58639c32a010d Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 21 Jul 2020 18:22:19 -0600 Subject: [PATCH 34/36] changes to land degradation logic to fix bugs --- biogeochem/EDLoggingMortalityMod.F90 | 23 +++++++++++++---------- biogeochem/EDPatchDynamicsMod.F90 | 22 +++++++++++++--------- 2 files changed, 26 insertions(+), 19 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 228076b197..117154932f 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -257,41 +257,44 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! l_degrad accounts for the affected area between logged crowns if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees + ! direct logging rates, based on dbh min and max criteria if (dbh >= logging_dbhmin .and. .not. & ((logging_dbhmax < fates_check_param_set) .and. (dbh < logging_dbhmax )) ) then lmort_direct = harvest_rate * logging_direct_frac - l_degrad = harvest_rate * (1._r8 - logging_direct_frac) ! fraction passed to 'degraded' forest. + else lmort_direct = 0.0_r8 - l_degrad = harvest_rate end if + ! infrastructure (roads, skid trails, etc) mortality rates if (dbh >= logging_dbhmax_infra) then lmort_infra = 0.0_r8 - l_degrad = l_degrad + harvest_rate * logging_mechanical_frac else lmort_infra = harvest_rate * logging_mechanical_frac end if - !damage rates for size class < & > threshold_size need to be specified seperately ! Collateral damage to smaller plants below the direct logging size threshold ! will be applied via "understory_death" via the disturbance algorithm - ! Important: Degredation rates really only have an impact when - ! applied to the canopy layer. So we don't add to degredation - ! for collateral damage, even understory collateral damage. - if (canopy_layer .eq. 1) then lmort_collateral = harvest_rate * logging_collateral_frac else lmort_collateral = 0._r8 endif - else ! non-woody plants in the canopy still become moved to secondary patch at same rate + 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 - l_degrad = harvest_rate end if + + ! the area occupied by all plants in the canopy that aren't killed is still disturbed at the harvest rate + if (canopy_layer .eq. 1) then + l_degrad = harvest_rate * (1._r8 - & + (logging_direct_frac + logging_collateral_frac + logging_mechanical_frac)) ! fraction passed to 'degraded' forest. + else + l_degrad = 0._r8 + endif + else lmort_direct = 0.0_r8 lmort_collateral = 0.0_r8 diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 99964a5273..b78fdeecc3 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -297,15 +297,6 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort => currentCohort%taller enddo !currentCohort - ! fraction of the logging disturbance rate that is non-harvested - if (currentPatch%disturbance_rates(dtype_ilog) .gt. nearzero) then - currentPatch%fract_ldist_not_harvested = dist_rate_ldist_notharvested / & - currentPatch%disturbance_rates(dtype_ilog) - endif - - ! Fire Disturbance Rate - currentPatch%disturbance_rates(dtype_ifire) = currentPatch%frac_burnt - ! for non-closed-canopy areas subject to logging, add an additional increment of area disturbed ! equivalent to the fradction loged to account for transfer of interstitial ground area to new secondary lands if ( logging_time .and. & @@ -316,6 +307,10 @@ subroutine disturbance_rates( site_in, bc_in) currentPatch%disturbance_rates(dtype_ilog) = currentPatch%disturbance_rates(dtype_ilog) + & (currentPatch%area - currentPatch%total_canopy_area) * 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 endif ! Sum of disturbance rates for different classes of disturbance across all patches in this site. do i_dist = 1,N_DIST_TYPES @@ -323,6 +318,15 @@ subroutine disturbance_rates( site_in, bc_in) currentPatch%disturbance_rates(i_dist) * currentPatch%area * AREA_INV end do + ! fraction of the logging disturbance rate that is non-harvested + if (currentPatch%disturbance_rates(dtype_ilog) .gt. nearzero) then + currentPatch%fract_ldist_not_harvested = dist_rate_ldist_notharvested / & + currentPatch%disturbance_rates(dtype_ilog) + endif + + ! Fire Disturbance Rate + currentPatch%disturbance_rates(dtype_ifire) = currentPatch%frac_burnt + ! Fires can't burn the whole patch, as this causes /0 errors. if (debug) then if (currentPatch%disturbance_rates(dtype_ifire) > 0.98_r8)then From abd57ec513a7d270e5f3b135e3fe32a5d5ed5be1 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 22 Jul 2020 12:20:26 -0600 Subject: [PATCH 35/36] bugfixes upon bugfixes --- biogeochem/EDLoggingMortalityMod.F90 | 3 +-- biogeochem/EDPatchDynamicsMod.F90 | 16 +++++++++------- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 117154932f..bda8c8cb71 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -289,8 +289,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! 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 * (1._r8 - & - (logging_direct_frac + logging_collateral_frac + logging_mechanical_frac)) ! fraction passed to 'degraded' forest. + l_degrad = harvest_rate - (lmort_direct + lmort_infra + lmort_collateral) ! fraction passed to 'degraded' forest. else l_degrad = 0._r8 endif diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index b78fdeecc3..0c12e08185 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -301,22 +301,18 @@ subroutine disturbance_rates( site_in, bc_in) ! equivalent to the fradction loged to account for transfer of interstitial ground area to new secondary lands if ( logging_time .and. & (currentPatch%area - currentPatch%total_canopy_area) .gt. fates_tiny ) then -! The canopy is NOT closed. + ! The canopy is NOT closed. + call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, & bc_in%hlm_harvest_rates, frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate) currentPatch%disturbance_rates(dtype_ilog) = currentPatch%disturbance_rates(dtype_ilog) + & (currentPatch%area - currentPatch%total_canopy_area) * harvest_rate / currentPatch%area - ! Non-harvested part of the logging disturbance rate + ! 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 endif -! Sum of disturbance rates for different classes of disturbance across all patches in this site. - do i_dist = 1,N_DIST_TYPES - site_in%potential_disturbance_rates(i_dist) = site_in%potential_disturbance_rates(i_dist) + & - currentPatch%disturbance_rates(i_dist) * currentPatch%area * AREA_INV - end do ! fraction of the logging disturbance rate that is non-harvested if (currentPatch%disturbance_rates(dtype_ilog) .gt. nearzero) then @@ -327,6 +323,12 @@ subroutine disturbance_rates( site_in, bc_in) ! Fire Disturbance Rate currentPatch%disturbance_rates(dtype_ifire) = currentPatch%frac_burnt + ! calculate a disgnostic sum of disturbance rates for different classes of disturbance across all patches in this site. + do i_dist = 1,N_DIST_TYPES + site_in%potential_disturbance_rates(i_dist) = site_in%potential_disturbance_rates(i_dist) + & + currentPatch%disturbance_rates(i_dist) * currentPatch%area * AREA_INV + end do + ! Fires can't burn the whole patch, as this causes /0 errors. if (debug) then if (currentPatch%disturbance_rates(dtype_ifire) > 0.98_r8)then From 8872a912bec6b22b793c4cf9e2f8d3cd07e69364 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 22 Jul 2020 16:53:14 -0600 Subject: [PATCH 36/36] adding runtime check on logging parameter values --- main/EDPftvarcon.F90 | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index c116375468..d816b21824 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -2118,6 +2118,7 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! ----------------------------------------------------------------------------------- use FatesConstantsMod , only : fates_check_param_set use FatesConstantsMod , only : itrue, ifalse + use EDParamsMod , only : logging_mechanical_frac, logging_collateral_frac, logging_direct_frac ! Argument logical, intent(in) :: is_master ! Only log if this is the master proc @@ -2172,7 +2173,13 @@ subroutine FatesCheckParams(is_master, parteh_mode) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - + ! logging parameters, make sure they make sense + if ( (logging_mechanical_frac + logging_collateral_frac + logging_direct_frac) .gt. 1._r8) then + write(fates_log(),*) 'the sum of logging_mechanical_frac + logging_collateral_frac + logging_direct_frac' + write(fates_log(),*) 'must be less than or equal to 1.' + write(fates_log(),*) 'Exiting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif do ipft = 1,npft