diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 83294dfb28..bda8c8cb71 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 @@ -44,6 +45,8 @@ module EDLoggingMortalityMod use FatesInterfaceTypesMod , only : hlm_model_day use FatesInterfaceTypesMod , only : hlm_day_of_year use FatesInterfaceTypesMod , only : hlm_days_per_year + use FatesInterfaceTypesMod , only : hlm_use_lu_harvest + use FatesInterfaceTypesMod , only : hlm_num_lu_harvest_cats use FatesInterfaceTypesMod , only : hlm_use_logging use FatesInterfaceTypesMod , only : hlm_use_planthydro use FatesConstantsMod , only : itrue,ifalse @@ -55,6 +58,12 @@ module EDLoggingMortalityMod use PRTGenericMod , only : sapw_organ, struct_organ, leaf_organ use PRTGenericMod , only : fnrt_organ, store_organ, repro_organ use FatesAllometryMod , only : set_root_fraction + use FatesConstantsMod , only : primaryforest, secondaryforest, secondary_age_threshold + use FatesConstantsMod , only : 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 private @@ -80,6 +89,7 @@ module EDLoggingMortalityMod public :: logging_litter_fluxes public :: logging_time public :: IsItLoggingTime + public :: get_harvest_rate_area contains @@ -104,10 +114,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 + ! 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 @@ -168,12 +185,21 @@ 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_rates, hlm_harvest_catnames, & + hlm_harvest_units, & + 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_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 real(r8), intent(out) :: lmort_collateral ! collateral damage mortality fraction real(r8), intent(out) :: lmort_infra ! infrastructure mortality fraction @@ -181,49 +207,93 @@ 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 + real(r8) :: harvest_rate ! the final harvest rate to apply to this cohort today - ! Parameters - real(r8), parameter :: adjustment = 1.0 ! adjustment for mortality rates - + ! 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 + ! 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 + + if (hlm_use_lu_harvest == ifalse) then + ! 0=use fates logging parameters directly when logging_time == .true. + ! this means harvest the whole cohort area + harvest_rate = 1._r8 + + else if (hlm_use_lu_harvest == itrue .and. hlm_harvest_units == hlm_harvest_area_fraction) then + ! We are harvesting based on areal fraction, not carbon/biomass terms. + ! 1=use area fraction from hlm + ! 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 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 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 == 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.' + 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 + + ! 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 - ! Pass logging rates to cohort level - - if (dbh >= logging_dbhmin ) then - lmort_direct = 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 end if - + + ! infrastructure (roads, skid trails, etc) mortality rates if (dbh >= logging_dbhmax_infra) then lmort_infra = 0.0_r8 - l_degrad = l_degrad + logging_mechanical_frac * adjustment else - lmort_infra = 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 ! 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 = logging_collateral_frac * adjustment + lmort_collateral = harvest_rate * logging_collateral_frac else - lmort_collateral = 0._r8 + lmort_collateral = 0._r8 endif - else + else ! non-woody plants still killed by infrastructure 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 end if + + ! the area occupied by all plants in the canopy that aren't killed is still disturbed at the harvest rate + if (canopy_layer .eq. 1) then + l_degrad = harvest_rate - (lmort_direct + lmort_infra + lmort_collateral) ! fraction passed to 'degraded' forest. + else + l_degrad = 0._r8 + endif + else lmort_direct = 0.0_r8 lmort_collateral = 0.0_r8 @@ -235,6 +305,90 @@ end subroutine LoggingMortality_frac ! ============================================================================ + subroutine get_harvest_rate_area (patch_anthro_disturbance_label, hlm_harvest_catnames, hlm_harvest_rates, & + 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_rates(:) ! annual harvest rate per hlm category + character(len=64), intent(in) :: hlm_harvest_catnames(:) ! names of hlm harvest categories + integer, intent(in) :: patch_anthro_disturbance_label ! patch level anthro_disturbance_label + real(r8), intent(in) :: secondary_age ! patch level age_since_anthro_disturbance + real(r8), intent(in) :: 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) + + ! 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 (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_rates(h_index) + endif + else if (patch_anthro_disturbance_label .eq. secondaryforest .and. & + secondary_age >= secondary_age_threshold) then + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH1") then + harvest_rate = harvest_rate + hlm_harvest_rates(h_index) + endif + else if (patch_anthro_disturbance_label .eq. secondaryforest .and. & + secondary_age < secondary_age_threshold) then + if(hlm_harvest_catnames(h_index) .eq. "HARVEST_SH2" .or. & + hlm_harvest_catnames(h_index) .eq. "HARVEST_SH3") then + harvest_rate = harvest_rate + hlm_harvest_rates(h_index) + endif + endif + end do + + ! Normalize by site-level primary or secondary forest fraction + ! since harvest_rate is specified as a fraction of the gridcell + ! also need to put a cap so as not to harvest more primary or secondary area than there is in a gridcell + if (patch_anthro_disturbance_label .eq. primaryforest) then + if (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/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 2c5cbdd0d5..fa0b933fc5 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) @@ -245,7 +246,13 @@ 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_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 8475f7762a..af0208cf43 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -36,7 +36,10 @@ 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 FatesConstantsMod , only : fates_tiny use FatesInterfaceTypesMod , only : hlm_use_planthydro use FatesInterfaceTypesMod , only : hlm_numSWb use FatesInterfaceTypesMod , only : bc_in_type @@ -50,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 @@ -77,6 +81,9 @@ 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 + use EDParamsMod, only : logging_export_frac ! CIME globals use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) @@ -96,6 +103,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__ @@ -168,12 +176,20 @@ subroutine disturbance_rates( site_in, bc_in) ! secondary forest patch) real(r8) :: dist_rate_ldist_notharvested 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) ! And the same rates in understory plants have already been applied to %dndt !---------------------------------------------------------------------------------------------- + ! 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)) @@ -196,13 +212,29 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort%asmort = asmort 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_rates, & + bc_in%hlm_harvest_catnames, & + bc_in%hlm_harvest_units, & + currentPatch%anthro_disturbance_label, & + currentPatch%age_since_anthro_disturbance, & + frac_site_primary) currentCohort%lmort_direct = lmort_direct currentCohort%lmort_collateral = lmort_collateral currentCohort%lmort_infra = lmort_infra currentCohort%l_degrad = l_degrad - + + ! 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%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) * logging_export_frac + endif + currentCohort => currentCohort%taller end do currentPatch%disturbance_mode = fates_unset_int @@ -213,6 +245,23 @@ 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 + + ! Recalculate total canopy area 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)) @@ -248,6 +297,23 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort => currentCohort%taller enddo !currentCohort + ! 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. & + (currentPatch%area - currentPatch%total_canopy_area) .gt. fates_tiny ) then + ! The canopy is NOT closed. + + call get_harvest_rate_area (currentPatch%anthro_disturbance_label, bc_in%hlm_harvest_catnames, & + bc_in%hlm_harvest_rates, frac_site_primary, currentPatch%age_since_anthro_disturbance, harvest_rate) + + 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 + ! 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 / & @@ -255,9 +321,15 @@ 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 + ! 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 write(fates_log(),*) 'very high fire areas', & @@ -408,6 +480,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 :: found_youngest_primary ! logical for finding the first primary forest patch !--------------------------------------------------------------------- storesmallcohort => null() ! storage of the smallest cohort for insertion routine @@ -419,11 +492,17 @@ 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)) 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 @@ -444,8 +523,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 @@ -1001,16 +1097,47 @@ subroutine spawn_patches( currentSite, bc_in) !*************************/ !** INSERT NEW PATCH(ES) INTO LINKED LIST - !**********`***************/ + !*************************/ 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 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 + 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 + found_youngest_primary = .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, 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 + 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 linked list (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 @@ -1883,7 +2010,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 @@ -2049,6 +2176,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 ! !--------------------------------------------------------------------- @@ -2056,11 +2184,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 @@ -2256,6 +2393,19 @@ subroutine fuse_patches( csite, bc_in ) enddo !do while nopatches>maxPatchesPerSite end do ! i_disttype loop + + currentPatch => currentSite%youngest_patch + 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 @@ -2296,6 +2446,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) @@ -2303,6 +2455,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 @@ -2438,6 +2594,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. !--------------------------------------------------------------------- @@ -2458,6 +2615,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) & @@ -2469,25 +2628,53 @@ 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 + ! and then allow them to fuse together. + 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 @@ -2699,4 +2886,33 @@ 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 0fbcb88b59..63f2fa6778 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -68,10 +68,11 @@ 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 - + use EDTypesMod , only : AREA_INV use PRTGenericMod, only : carbon12_element use PRTGenericMod, only : all_carbon_elements use PRTGenericMod, only : leaf_organ @@ -305,6 +306,10 @@ 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 + + + 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)) @@ -338,7 +343,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 diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 071ffcd69f..dfaa4d00e4 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -125,7 +125,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" diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 356c1e03c0..6e61c45fd8 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -2116,6 +2116,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 @@ -2170,7 +2171,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 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index a769911db2..0e843df123 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -766,6 +766,14 @@ 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] + 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 @@ -940,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) diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 34518d8ac8..e553918028 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -34,6 +34,13 @@ module FatesConstantsMod integer, parameter, public :: n_anthro_disturbance_categories = 2 integer, parameter, public :: primaryforest = 1 integer, parameter, public :: secondaryforest = 2 + 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 @@ -134,7 +141,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/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 05a236174c..a022e64624 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,16 @@ 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 + integer :: ih_fire_disturbance_rate_si + 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 integer :: ih_nplant_canopy_si_scag @@ -419,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 @@ -1738,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, & @@ -1774,6 +1789,15 @@ 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, & + 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_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, & @@ -2037,6 +2061,30 @@ 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(io_si) = sites(s)%primary_land_patchfusion_error + + ! output site-level disturbance rates + 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(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(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(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(io_si) = sum(sites(s)%potential_disturbance_rates(1:N_DIST_TYPES)) + + hio_harvest_carbonflux_si(io_si) = sites(s)%harvest_carbon_flux + ipa = 0 cpatch => sites(s)%oldest_patch do while(associated(cpatch)) @@ -2199,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 @@ -3889,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, & @@ -4200,6 +4259,52 @@ 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='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, & + 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_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 ) + + 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 ) + + 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', & diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 5b440dc5cf..f67387e1a5 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -49,6 +49,7 @@ module FatesInterfaceMod use EDParamsMod , only : ED_val_history_coageclass_bin_edges use CLMFatesParamInterfaceMod , only : FatesReadParameters use PRTAllometricCarbonMod , only : InitPRTGlobalAllometricCarbon + use decompMod , only : bounds_type ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -233,7 +234,7 @@ end subroutine zero_bcs ! =========================================================================== - subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in) + subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats) ! --------------------------------------------------------------------------------- ! Allocate and Initialze the FATES boundary condition vectors @@ -243,9 +244,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 + integer,intent(in) :: num_lu_harvest_cats + ! Allocate input boundaries - bc_in%nlevsoil = nlevsoil_in if(nlevsoil_in > numlevsoil_max) then @@ -360,8 +362,16 @@ 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 - allocate(bc_in%pft_areafrac(maxpft)) + ! Land use + + ! 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_rates(num_lu_harvest_cats)) + allocate(bc_in%hlm_harvest_catnames(num_lu_harvest_cats)) + end if + allocate(bc_in%pft_areafrac(maxpft)) return end subroutine allocate_bcin @@ -991,6 +1001,8 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_parteh_mode = unset_int hlm_spitfire_mode = unset_int hlm_use_planthydro = unset_int + hlm_use_lu_harvest = unset_int + hlm_num_lu_harvest_cats = unset_int hlm_use_cohort_age_tracking = unset_int hlm_use_logging = unset_int hlm_use_ed_st3 = unset_int @@ -1044,6 +1056,20 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(), *) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' end if + 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 lu_harvest flag must be 0 or 1 or 2, exiting' + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + 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 + 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' @@ -1316,13 +1342,24 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(),*) 'Transfering hlm_use_planthydro= ',ival,' to FATES' end if + case('use_lu_harvest') + hlm_use_lu_harvest = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering hlm_use_lu_harvest= ',ival,' to FATES' + end if + + case('num_lu_harvest_cats') + hlm_num_lu_harvest_cats = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering hlm_num_lu_harvest_cats= ',ival,' to FATES' + end if + case('use_cohort_age_tracking') hlm_use_cohort_age_tracking = ival if (fates_global_verbose()) then write(fates_log(),*) 'Transfering hlm_use_cohort_age_tracking= ',ival,' to FATES' end if - case('use_logging') hlm_use_logging = ival if (fates_global_verbose()) then diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index fc90aed1dc..7146cc2272 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -89,9 +89,29 @@ module FatesInterfaceTypesMod ! See namelist_definition_clm4_5.xml ! ignitions: 1=constant, >1=external data sources (lightning and/or anthropogenic) + 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 + ! from from hlm + ! 2 = use carbon from hlm + ! If 1 or 2, it automatically sets + ! hlm_use_logging to 1 + + 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_rates and bc_in%hlm_harvest_catnames integer, public :: hlm_use_logging ! This flag signals whether or not to use ! the logging module + ! If hlm_use_lu_harvest is zero, + ! then logging is determined by + ! the fates parameter file + ! If hlm_use_lu_harvest is non-zero, + ! then this flag is automatically + ! set to 1 and logging is determined + ! by the lu harvest input from the hlm integer, public :: hlm_use_planthydro ! This flag signals whether or not to use ! plant hydraulics (bchristo/xu methods) @@ -433,7 +453,15 @@ module FatesInterfaceTypesMod 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 + ! --------------------------------------------------------------------------------- + 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 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, &