diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 30a41e8790..8b0ca90135 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1134,8 +1134,8 @@ subroutine spawn_patches( currentSite, bc_in ) endif - currentSite%mass_balance(el)%burn_flux_to_atm = & - currentSite%mass_balance(el)%burn_flux_to_atm + & + currentSite%mass_balance(el)%burn_flux_to_atm(i_disturbance_type) = & + currentSite%mass_balance(el)%burn_flux_to_atm(i_disturbance_type) + & leaf_burn_frac * leaf_m * nc%n ! This term increments the loss flux from surviving trees @@ -2042,7 +2042,7 @@ subroutine TransLitterNewPatch(currentSite, & if (debug) then - burn_flux0 = site_mass%burn_flux_to_atm + burn_flux0 = site_mass%burn_flux_to_atm(dist_type) litter_stock0 = curr_litt%GetTotalLitterMass()*currentPatch%area + & new_litt%GetTotalLitterMass()*newPatch%area end if @@ -2063,7 +2063,7 @@ subroutine TransLitterNewPatch(currentSite, & new_litt%ag_cwd(c) = new_litt%ag_cwd(c) + donatable_mass*donate_m2 curr_litt%ag_cwd(c) = curr_litt%ag_cwd(c) + donatable_mass*retain_m2 - site_mass%burn_flux_to_atm = site_mass%burn_flux_to_atm + burned_mass + site_mass%burn_flux_to_atm(dist_type) = site_mass%burn_flux_to_atm(dist_type) + burned_mass ! Transfer below ground CWD (none burns) do sl = 1,currentSite%nlevsoil @@ -2091,7 +2091,7 @@ subroutine TransLitterNewPatch(currentSite, & new_litt%leaf_fines(dcmpy) = new_litt%leaf_fines(dcmpy) + donatable_mass*donate_m2 curr_litt%leaf_fines(dcmpy) = curr_litt%leaf_fines(dcmpy) + donatable_mass*retain_m2 - site_mass%burn_flux_to_atm = site_mass%burn_flux_to_atm + burned_mass + site_mass%burn_flux_to_atm(dist_type) = site_mass%burn_flux_to_atm(dist_type) + burned_mass ! Transfer root fines (none burns) do sl = 1,currentSite%nlevsoil @@ -2122,7 +2122,7 @@ subroutine TransLitterNewPatch(currentSite, & ! EDMainMod start triggering. ! -------------------------------------------------------------------------- if (debug) then - burn_flux1 = site_mass%burn_flux_to_atm + burn_flux1 = site_mass%burn_flux_to_atm(dist_type) litter_stock1 = curr_litt%GetTotalLitterMass()*remainder_area + & new_litt%GetTotalLitterMass()*newPatch%area error = (litter_stock1 - litter_stock0) + (burn_flux1-burn_flux0) @@ -2301,7 +2301,7 @@ subroutine fire_litter_fluxes(currentSite, currentPatch, & donatable_mass*retain_m2*dcmpy_frac end do - site_mass%burn_flux_to_atm = site_mass%burn_flux_to_atm + burned_mass + site_mass%burn_flux_to_atm(dtype_ifire) = site_mass%burn_flux_to_atm(dtype_ifire) + burned_mass call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil, & bc_in%max_rooting_depth_index_col) @@ -2363,7 +2363,7 @@ subroutine fire_litter_fluxes(currentSite, currentPatch, & donatable_mass = donatable_mass * (1.0_r8-currentCohort%fraction_crown_burned) burned_mass = num_dead_trees * SF_val_CWD_frac_adj(c) * bstem * & currentCohort%fraction_crown_burned - site_mass%burn_flux_to_atm = site_mass%burn_flux_to_atm + burned_mass + site_mass%burn_flux_to_atm(dtype_ifire) = site_mass%burn_flux_to_atm(dtype_ifire) + burned_mass endif new_litt%ag_cwd(c) = new_litt%ag_cwd(c) + donatable_mass * donate_m2 curr_litt%ag_cwd(c) = curr_litt%ag_cwd(c) + donatable_mass * retain_m2 @@ -2768,7 +2768,8 @@ subroutine landusechange_litter_fluxes(currentSite, currentPatch, & donatable_mass*retain_m2*dcmpy_frac end do - site_mass%burn_flux_to_atm = site_mass%burn_flux_to_atm + burned_mass + site_mass%burn_flux_to_atm(dtype_ilandusechange) = & + site_mass%burn_flux_to_atm(dtype_ilandusechange) + burned_mass call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil, & bc_in%max_rooting_depth_index_col) @@ -2828,8 +2829,8 @@ subroutine landusechange_litter_fluxes(currentSite, currentPatch, & burned_mass = num_dead_trees * SF_val_CWD_frac(c) * bstem * & EDPftvarcon_inst%landusechange_frac_burned(pft) - site_mass%burn_flux_to_atm = site_mass%burn_flux_to_atm + burned_mass - + site_mass%burn_flux_to_atm(dtype_ilandusechange) = & + site_mass%burn_flux_to_atm(dtype_ilandusechange) + burned_mass else ! all other pools can end up as timber products or burn or go to litter donatable_mass = donatable_mass * (1.0_r8-EDPftvarcon_inst%landusechange_frac_exported(pft)) * & (1.0_r8-EDPftvarcon_inst%landusechange_frac_burned(pft)) @@ -2841,7 +2842,8 @@ subroutine landusechange_litter_fluxes(currentSite, currentPatch, & woodproduct_mass = num_dead_trees * SF_val_CWD_frac(c) * bstem * & EDPftvarcon_inst%landusechange_frac_exported(pft) - site_mass%burn_flux_to_atm = site_mass%burn_flux_to_atm + burned_mass + site_mass%burn_flux_to_atm(dtype_ilandusechange) = & + site_mass%burn_flux_to_atm(dtype_ilandusechange) + burned_mass ! Amount of trunk mass exported off site [kg/m2] elflux_diags%exported_harvest = elflux_diags%exported_harvest + & diff --git a/main/ChecksBalancesMod.F90 b/main/ChecksBalancesMod.F90 index f041389fd3..d93aa88a21 100644 --- a/main/ChecksBalancesMod.F90 +++ b/main/ChecksBalancesMod.F90 @@ -331,7 +331,7 @@ subroutine CheckIntegratedMassPools(site) ibal%iflux_litter = ibal%iflux_litter + & tot_litter_input - & (site_mass%frag_out*area_inv - ediag%tot_seed_turnover) - & - (site_mass%burn_flux_to_atm*area_inv - ediag%burned_liveveg) + (sum(site_mass%burn_flux_to_atm(:))*area_inv - ediag%burned_liveveg) ediag%err_liveveg = ibal%iflux_liveveg - ibal%state_liveveg diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index df21329d30..f07a3c75de 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -917,9 +917,8 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting ) bc_out%seed_c_si = bc_out%seed_c_si * g_per_kg * AREA_INV ! Set boundary condition to HLM for carbon loss to atm from fires and grazing - ! [kgC/ha/day]*[ha/m2]*[day/s] = [kg/m2/s] - - bc_out%fire_closs_to_atm_si = site_cmass%burn_flux_to_atm * area_inv * days_per_sec + ! [kgC/ha/day]*[ha/m2]*[day/s] = [kg/m2/s] + bc_out%fire_closs_to_atm_si = sum(site_cmass%burn_flux_to_atm(:)) * area_inv * days_per_sec bc_out%grazing_closs_to_atm_si = site_cmass%herbivory_flux_out * area_inv * days_per_sec end subroutine ed_update_site @@ -1009,9 +1008,10 @@ subroutine TotalBalanceCheck (currentSite, call_index, is_restarting ) site_mass%flux_generic_in + & site_mass%patch_resize_err + flux_out = sum(site_mass%wood_product_harvest(:)) + & sum(site_mass%wood_product_landusechange(:)) + & - site_mass%burn_flux_to_atm + & + sum(site_mass%burn_flux_to_atm(:)) + & site_mass%seed_out + & site_mass%flux_generic_out + & site_mass%frag_out + & @@ -1045,7 +1045,7 @@ subroutine TotalBalanceCheck (currentSite, call_index, is_restarting ) write(fates_log(),*) 'wood_product_harvest: ',site_mass%wood_product_harvest(:) write(fates_log(),*) 'wood_product_landusechange: ',site_mass%wood_product_landusechange(:) write(fates_log(),*) 'error from patch resizing: ',site_mass%patch_resize_err - write(fates_log(),*) 'burn_flux_to_atm: ',site_mass%burn_flux_to_atm + write(fates_log(),*) 'burn_flux_to_atm: ',site_mass%burn_flux_to_atm(:) write(fates_log(),*) 'seed_out: ',site_mass%seed_out write(fates_log(),*) 'flux_generic_out: ',site_mass%flux_generic_out write(fates_log(),*) 'frag_out: ',site_mass%frag_out diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 72a8614992..5aca294987 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -298,7 +298,7 @@ module EDTypesMod real(r8) :: wood_product_landusechange(maxpft) ! Total mass exported as wood product from land use change [kg/site/day] - real(r8) :: burn_flux_to_atm ! Total mass burned and exported to the atmosphere [kg/site/day] + real(r8) :: burn_flux_to_atm(n_dist_types) ! Total mass burned and exported to the atmosphere [kg/site/day] real(r8) :: flux_generic_in ! Used for prescribed or artificial input fluxes ! and initialization [kg/site/day] @@ -308,7 +308,7 @@ module EDTypesMod ! due to re-sizing patches when area math starts to lose ! precision - real(r8) :: herbivory_flux_out ! loss of element due to grazing (and/or browsing) by herbivores + real(r8) :: herbivory_flux_out ! loss of element due to grazing (and/or browsing) by herbivores [kg/site/day] contains @@ -717,7 +717,7 @@ subroutine ZeroMassBalFlux(this) this%frag_out = 0._r8 this%wood_product_harvest(:) = 0._r8 this%wood_product_landusechange(:) = 0._r8 - this%burn_flux_to_atm = 0._r8 + this%burn_flux_to_atm(:) = 0._r8 this%flux_generic_in = 0._r8 this%flux_generic_out = 0._r8 this%patch_resize_err = 0._r8 diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index ee8cc1ed3c..81b5a6f028 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -9,6 +9,7 @@ module FatesHistoryInterfaceMod use FatesConstantsMod , only : mg_per_kg use FatesConstantsMod , only : pi_const use FatesConstantsMod , only : nearzero + use FatesConstantsMod , only : rsnbl_math_prec use FatesConstantsMod , only : t_water_freeze_k_1atm use FatesConstantsMod , only : n_term_mort_types use FatesConstantsMod , only : i_term_mort_type_cstarv @@ -41,6 +42,7 @@ module FatesHistoryInterfaceMod use FatesConstantsMod , only : dtype_ifall use FatesConstantsMod , only : dtype_ifire use FatesConstantsMod , only : dtype_ilog + use FatesConstantsMod , only : dtype_ilandusechange use FatesIODimensionsMod , only : fates_io_dimension_type use FatesIOVariableKindMod , only : fates_io_variable_kind_type use FatesIOVariableKindMod , only : site_int @@ -138,6 +140,7 @@ module FatesHistoryInterfaceMod use FatesSizeAgeTypeIndicesMod, only : get_cdamagesize_class_index use FatesSizeAgeTypeIndicesMod, only : get_cdamagesizepft_class_index use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index + use FatesInterfaceTypesMod , only : hlm_use_luh implicit none private ! By default everything is private @@ -306,6 +309,7 @@ module FatesHistoryInterfaceMod integer :: ih_seedling_pool_si ! carbon only integer :: ih_ba_weighted_height_si integer :: ih_ca_weighted_height_si + integer :: ih_patch_weighted_95thpctile_height_si integer :: ih_seeds_in_local_elem integer :: ih_seeds_in_extern_elem integer :: ih_seed_decay_elem @@ -362,6 +366,12 @@ module FatesHistoryInterfaceMod integer :: ih_burnedarea_si_landuse integer :: ih_gpp_si_landuse integer :: ih_npp_si_landuse + integer :: ih_tveg_si_landuse + integer :: ih_tsa_si_landuse + integer :: ih_sw_abs_si_landuse + integer :: ih_lw_net_si_landuse + integer :: ih_shflux_si_landuse + integer :: ih_lhflux_si_landuse ! land use by land use variables integer :: ih_disturbance_rate_si_lulu @@ -406,6 +416,8 @@ module FatesHistoryInterfaceMod integer :: ih_vis_rad_err_si integer :: ih_nir_rad_err_si integer :: ih_fire_c_to_atm_si + integer :: ih_fire_c_to_atm_landusechange_si + integer :: ih_fire_livec_to_atm_si integer :: ih_interr_liveveg_elem integer :: ih_interr_litter_elem integer :: ih_cbal_err_fates_si @@ -675,14 +687,9 @@ module FatesHistoryInterfaceMod integer :: ih_biomass_si_age integer :: ih_c_stomata_si_age integer :: ih_c_lblayer_si_age - integer :: ih_agesince_anthrodist_si integer :: ih_agesince_anthrodist_si_age - integer :: ih_secondarylands_area_si_age - integer :: ih_primarylands_area_si_age - integer :: ih_area_burnt_si_age - integer :: ih_primarylands_fracarea_si - integer :: ih_secondarylands_fracarea_si integer :: ih_secondarylands_fracarea_si_age + integer :: ih_secondary_agb_si_agesinceanthro integer :: ih_primarylands_fracarea_si_age integer :: ih_fracarea_burnt_si_age integer :: ih_rx_fracarea_burnt_si_age @@ -693,6 +700,7 @@ module FatesHistoryInterfaceMod integer :: ih_rx_intensity_si_age integer :: ih_nonrx_intensity_si_age + ! indices to (site x height) variables integer :: ih_canopy_height_dist_si_height integer :: ih_leaf_height_dist_si_height @@ -858,6 +866,7 @@ module FatesHistoryInterfaceMod procedure :: update_history_hifrq_sitelevel procedure :: update_history_hifrq_subsite procedure :: update_history_hifrq_subsite_ageclass + procedure :: update_history_hifrq_landuse procedure :: update_history_hydraulics procedure :: update_history_nutrflux @@ -2429,7 +2438,8 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) real(r8) :: leaf_herbivory ! mass of leaves eaten by herbivores [kg/yr] real(r8) :: n_perm2 ! abundance per m2 real(r8) :: patch_fracarea ! Fraction of area for this patch - + real(r8) :: crown_area_covered ! accumulator variable for patch crown area + associate( hio_npatches_si => this%hvars(ih_npatches_si)%r81d, & hio_ncohorts_si => this%hvars(ih_ncohorts_si)%r81d, & hio_ncl_si => this%hvars(ih_ncl_si)%r81d, & @@ -2440,6 +2450,7 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) hio_fates_fraction_si => this%hvars(ih_fates_fraction_si)%r81d, & hio_ba_weighted_height_si => this%hvars(ih_ba_weighted_height_si)%r81d, & hio_ca_weighted_height_si => this%hvars(ih_ca_weighted_height_si)%r81d, & + hio_patch_weighted_95thpctile_height_si => this%hvars(ih_patch_weighted_95thpctile_height_si)%r81d, & hio_canopy_spread_si => this%hvars(ih_canopy_spread_si)%r81d, & hio_nesterov_fire_danger_si => this%hvars(ih_nesterov_fire_danger_si)%r81d, & hio_rx_burn_window_si => this%hvars(ih_rx_burn_window_si)%r81d, & @@ -2496,6 +2507,8 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) hio_canopy_mortality_crownarea_si => this%hvars(ih_canopy_mortality_crownarea_si)%r81d, & hio_ustory_mortality_crownarea_si => this%hvars(ih_understory_mortality_crownarea_si)%r81d, & hio_fire_c_to_atm_si => this%hvars(ih_fire_c_to_atm_si)%r81d, & + hio_fire_livec_to_atm_si => this%hvars(ih_fire_livec_to_atm_si)%r81d, & + hio_fire_c_to_atm_landusechange_si => this%hvars(ih_fire_c_to_atm_landusechange_si)%r81d, & hio_demotion_carbonflux_si => this%hvars(ih_demotion_carbonflux_si)%r81d, & hio_promotion_carbonflux_si => this%hvars(ih_promotion_carbonflux_si)%r81d, & hio_canopy_mortality_carbonflux_si => this%hvars(ih_canopy_mortality_carbonflux_si)%r81d, & @@ -2540,7 +2553,11 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) ! Total carbon lost to atmosphere from burning (kgC/site/day -> kgC/m2/s) hio_fire_c_to_atm_si(io_si) = & - sites(s)%mass_balance(element_pos(carbon12_element))%burn_flux_to_atm * & + sum(sites(s)%mass_balance(element_pos(carbon12_element))%burn_flux_to_atm(:)) * & + ha_per_m2 * days_per_sec + + hio_fire_c_to_atm_landusechange_si(io_si) = & + sites(s)%mass_balance(element_pos(carbon12_element))%burn_flux_to_atm(dtype_ilandusechange) * & ha_per_m2 * days_per_sec ! damage variables - site level - this needs to be OUT of the patch loop @@ -2677,6 +2694,9 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) sum(elflux_diags_c%root_litter_input(:))) * & AREA_INV * days_per_sec + hio_fire_livec_to_atm_si(io_si) = & + elflux_diags_c%burned_liveveg * ha_per_m2 * days_per_sec + ! Loop through patches to sum up diagonistics cpatch => sites(s)%oldest_patch patchloop: do while(associated(cpatch)) @@ -2953,7 +2973,7 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) hio_npp_stor_si(io_si) = hio_npp_stor_si(io_si) + & store_m_net_alloc * n_perm2 / days_per_year / sec_per_day - leaf_herbivory = ccohort%prt%GetHerbivory(leaf_organ, carbon12_element) * days_per_year !cdkcdk + leaf_herbivory = ccohort%prt%GetHerbivory(leaf_organ, carbon12_element) * days_per_year hio_grazing_si(io_si) = hio_grazing_si(io_si) + leaf_herbivory * n_perm2 / days_per_year / sec_per_day ! Woody State Variables (basal area growth increment) @@ -3006,6 +3026,24 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) ccohort => ccohort%taller enddo cohortloop ! cohort loop + ! mean 95th percentile height. loop through cohorts on patch again, this time from tallest to shortest + crown_area_covered = 0._r8 + ccohort => cpatch%tallest + do while(associated(ccohort)) + if (ccohort%canopy_layer .eq. 1) then ! ignore anything not in the canopy + if (((ccohort%c_area + crown_area_covered)/cpatch%area) .ge. 0.05_r8 ) then + hio_patch_weighted_95thpctile_height_si(io_si) = & + hio_patch_weighted_95thpctile_height_si(io_si) + ccohort%height * cpatch%area/AREA + ccohort => null() ! exit the cohort loop + else + crown_area_covered = crown_area_covered + ccohort%c_area + ccohort => ccohort%shorter + endif + else + ccohort => ccohort%shorter + endif + enddo ! cohort loop + cpatch => cpatch%younger end do patchloop !patch loop @@ -3300,9 +3338,6 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) hio_yesterdaycanopylevel_understory_si_scls => this%hvars(ih_yesterdaycanopylevel_understory_si_scls)%r82d, & hio_fracarea_si => this%hvars(ih_fracarea_si)%r81d, & hio_canopy_fracarea_si => this%hvars(ih_canopy_fracarea_si)%r81d, & - hio_agesince_anthrodist_si => this%hvars(ih_agesince_anthrodist_si)%r81d, & - hio_primarylands_fracarea_si => this%hvars(ih_primarylands_fracarea_si)%r81d, & - hio_secondarylands_fracarea_si => this%hvars(ih_secondarylands_fracarea_si)%r81d, & hio_fracarea_si_landuse => this%hvars(ih_fracarea_si_landuse)%r82d, & hio_npp_si_landuse => this%hvars(ih_npp_si_landuse)%r82d, & hio_biomass_si_landuse => this%hvars(ih_biomass_si_landuse)%r82d, & @@ -3337,6 +3372,7 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) hio_cstarvmortality_continuous_carbonflux_si_pft => this%hvars(ih_cstarvmortality_continuous_carbonflux_si_pft)%r82d, & hio_transition_matrix_si_lulu => this%hvars(ih_transition_matrix_si_lulu)%r82d, & hio_scorch_height_si_pft => this%hvars(ih_scorch_height_si_pft)%r82d, & + hio_secondary_agb_si_agesinceanthro => this%hvars(ih_secondary_agb_si_agesinceanthro)%r82d, & hio_sapwood_area_scpf => this%hvars(ih_sapwood_area_scpf)%r82d) model_day_int = nint(hlm_model_day) @@ -3395,7 +3431,7 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) ! Total element lost to atmosphere from burning (kg/site/day -> kg/m2/s) hio_burn_flux_elem(io_si,el) = & - sites(s)%mass_balance(el)%burn_flux_to_atm * ha_per_m2 * & + sum(sites(s)%mass_balance(el)%burn_flux_to_atm(:)) * ha_per_m2 * & days_per_sec end do @@ -3443,25 +3479,6 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) cpatch%frac_burnt * cpatch%area * AREA_INV / sec_per_day end if - ! some diagnostics on secondary forest area and its age distribution - if ( cpatch%land_use_label .eq. secondaryland ) then - - hio_agesince_anthrodist_si(io_si) = & - hio_agesince_anthrodist_si(io_si) & - + cpatch%area * AREA_INV - - hio_secondarylands_fracarea_si(io_si) = & - hio_secondarylands_fracarea_si(io_si) & - + cpatch%area * AREA_INV - - else if ( cpatch%land_use_label .eq. primaryland ) then - - hio_primarylands_fracarea_si(io_si) = & - hio_primarylands_fracarea_si(io_si) & - + cpatch%area * AREA_INV - - endif - do ft = 1,numpft hio_scorch_height_si_pft(io_si,ft) = hio_scorch_height_si_pft(io_si,ft) + & cpatch%Scorch_ht(ft) * cpatch%area * AREA_INV @@ -3890,6 +3907,16 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) hio_biomass_si_scls(io_si,scls) = hio_biomass_si_scls(io_si,scls) + & total_m * ccohort%n * AREA_INV + ! agb by patch age since anthropogenic disturbance + if ( cpatch%land_use_label .eq. secondaryland ) then + + iscag_anthrodist = get_age_class_index(cpatch%age_since_anthro_disturbance) + + hio_secondary_agb_si_agesinceanthro(io_si,iscag_anthrodist) = & + hio_secondary_agb_si_agesinceanthro(io_si,iscag_anthrodist) & + + total_m * ccohort%n * prt_params%allom_agb_frac(ccohort%pft) * AREA_INV + endif + ! update size-class quantities hio_nplant_si_scls(io_si,scls) = hio_nplant_si_scls(io_si,scls) + ccohort%n / m2_per_ha @@ -4878,6 +4905,8 @@ subroutine update_history_dyn_subsite_ageclass(this,nc,nsites,sites) hio_secondarylands_fracarea_si_age(io_si,cpatch%age_class) = & hio_secondarylands_fracarea_si_age(io_si,cpatch%age_class) & + patch_area_div_site_area + + else if ( cpatch%land_use_label .eq. primaryland) then hio_primarylands_fracarea_si_age(io_si,cpatch%age_class) = & hio_primarylands_fracarea_si_age(io_si,cpatch%age_class) & @@ -5125,6 +5154,9 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,bc_out,dt_tstep) if(hlm_hist_level_hifrq>1) then call update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) call update_history_hifrq_subsite_ageclass(this,nsites,sites,dt_tstep) + if (hlm_use_luh .eq. itrue) then + call update_history_hifrq_landuse(this,nc,nsites,sites,bc_in,dt_tstep) + end if end if end if @@ -5358,6 +5390,162 @@ end subroutine update_history_hifrq_sitelevel ! =============================================================================================== + + subroutine update_history_hifrq_landuse(this,nc,nsites,sites,bc_in,dt_tstep) + + ! + ! Arguments + class(fates_history_interface_type) :: this + 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) + real(r8) , intent(in) :: dt_tstep + + ! Locals + integer :: s ! The local site index + integer :: io_si ! The site index of the IO array + + real(r8) :: landuse_statevector(n_landuse_cats) + real(r8) :: canopy_area_bylanduse(n_landuse_cats) + integer :: i_lu + logical :: foundbaregroundpatch + + type(fates_patch_type),pointer :: cpatch + type(fates_cohort_type),pointer :: ccohort + real(r8) :: dt_tstep_inv ! Time step in frequency units (/s) + real(r8) :: vegarea_per_patcharea ! temporary weighting variable (unitless) + + associate( hio_tveg_si_landuse => this%hvars(ih_tveg_si_landuse)%r82d,& + hio_gpp_si_landuse => this%hvars(ih_gpp_si_landuse)%r82d, & + hio_tsa_si_landuse => this%hvars(ih_tsa_si_landuse)%r82d,& + hio_sw_abs_si_landuse => this%hvars(ih_sw_abs_si_landuse)%r82d,& + hio_lw_net_si_landuse => this%hvars(ih_lw_net_si_landuse)%r82d,& + hio_shflux_si_landuse => this%hvars(ih_shflux_si_landuse)%r82d,& + hio_lhflux_si_landuse => this%hvars(ih_lhflux_si_landuse)%r82d) + + do_sites: do s = 1,nsites + + io_si = sites(s)%h_gid + + dt_tstep_inv = 1.0_r8/dt_tstep + + ! biophysical properties that are indexed by land use + landuse_statevector(:) = sites(s)%get_current_landuse_statevector() * AREA + + ! get the total canopy area for each land use type + canopy_area_bylanduse(:) = 0._r8 + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + if ( cpatch%land_use_label .ne. nocomp_bareground_land ) then + canopy_area_bylanduse(cpatch%land_use_label) = canopy_area_bylanduse(cpatch%land_use_label) + & + cpatch%total_canopy_area + endif + cpatch => cpatch%younger + end do + + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + if (cpatch%total_canopy_area .gt. rsnbl_math_prec) then + ! for TVEG, since it is only defined on vegetated area of vegetated patches, normalize by the total vegetated area + hio_tveg_si_landuse(io_si,cpatch%land_use_label) = hio_tveg_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%t_veg_pa(cpatch%patchno) * cpatch%total_canopy_area/canopy_area_bylanduse(cpatch%land_use_label) + + ! for the rest of these, first weight by the vegetated area of each patch over the total patch area for each land use type + vegarea_per_patcharea = cpatch%total_canopy_area/landuse_statevector(cpatch%land_use_label) + + hio_tsa_si_landuse(io_si,cpatch%land_use_label) = hio_tsa_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%t2m_pa(cpatch%patchno) * vegarea_per_patcharea + + hio_sw_abs_si_landuse(io_si,cpatch%land_use_label) = hio_sw_abs_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%swabs_pa(cpatch%patchno) * vegarea_per_patcharea + + hio_lw_net_si_landuse(io_si,cpatch%land_use_label) = hio_lw_net_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%netlw_pa(cpatch%patchno) * vegarea_per_patcharea + + hio_shflux_si_landuse(io_si,cpatch%land_use_label) = hio_shflux_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%shflux_pa(cpatch%patchno) * vegarea_per_patcharea + + hio_lhflux_si_landuse(io_si,cpatch%land_use_label) = hio_lhflux_si_landuse(io_si,cpatch%land_use_label) + & + bc_in(s)%lhflux_pa(cpatch%patchno) * vegarea_per_patcharea + endif + cpatch => cpatch%younger + end do + + ! for all the land-use indexed variables, except for TVEG, also add in the component for the unvegetated area of each land use + foundbaregroundpatch = .false. + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + if (cpatch%land_use_label .eq. nocomp_bareground_land .and. .not. foundbaregroundpatch) then + foundbaregroundpatch = .true. + do i_lu = 1, n_landuse_cats + if ( landuse_statevector(i_lu) .gt. rsnbl_math_prec ) then + hio_tsa_si_landuse(io_si,i_lu) = hio_tsa_si_landuse(io_si,i_lu) + & + bc_in(s)%t2m_pa(cpatch%patchno) * & + (landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu) + + hio_sw_abs_si_landuse(io_si,i_lu) = hio_sw_abs_si_landuse(io_si,i_lu) + & + bc_in(s)%swabs_pa(cpatch%patchno) * & + (landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu) + + hio_lw_net_si_landuse(io_si,i_lu) = hio_lw_net_si_landuse(io_si,i_lu) + & + bc_in(s)%netlw_pa(cpatch%patchno) * & + (landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu) + + hio_shflux_si_landuse(io_si,i_lu) = hio_shflux_si_landuse(io_si,i_lu) + & + bc_in(s)%shflux_pa(cpatch%patchno) * & + (landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu) + + hio_lhflux_si_landuse(io_si,i_lu) = hio_lhflux_si_landuse(io_si,i_lu) + & + bc_in(s)%lhflux_pa(cpatch%patchno) * & + (landuse_statevector(i_lu) - canopy_area_bylanduse(i_lu)) / landuse_statevector(i_lu) + end if + end do + end if + cpatch => cpatch%younger + end do + + ! instead of leaving the values for unoccupied areas as zero, set as missing values + do i_lu = 1, n_landuse_cats + + ! if a given land use type is not present, set the value as missing + if ( landuse_statevector(i_lu) .le. rsnbl_math_prec ) then + hio_tsa_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + hio_sw_abs_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + hio_lw_net_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + hio_shflux_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + hio_lhflux_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + end if + + ! for tveg, ignore if there is no vegetation present on any patches of a given land use type + if ( canopy_area_bylanduse(i_lu) .le. rsnbl_math_prec ) then + hio_tveg_si_landuse(io_si,i_lu) = hlm_hio_ignore_val + end if + + end do + + ! for GPP by land use, we need to loop over both patches and cohorts + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + ccohort => cpatch%shortest + do while(associated(ccohort)) + if ( (.not. ccohort%isnew) .and. (cpatch%land_use_label .gt. nocomp_bareground_land) ) then + hio_gpp_si_landuse(io_si,cpatch%land_use_label) = hio_gpp_si_landuse(io_si,cpatch%land_use_label) & + + ccohort%gpp_tstep * ccohort%n * dt_tstep_inv + end if + ccohort => ccohort%taller + end do + cpatch => cpatch%younger + end do + + end do do_sites + + end associate + return + end subroutine update_history_hifrq_landuse + + ! =============================================================================================== + subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) ! --------------------------------------------------------------------------------- @@ -5416,7 +5604,6 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) hio_froot_mr_understory_si_scls => this%hvars(ih_froot_mr_understory_si_scls)%r82d, & hio_resp_g_understory_si_scls => this%hvars(ih_resp_g_understory_si_scls)%r82d, & hio_resp_m_understory_si_scls => this%hvars(ih_resp_m_understory_si_scls)%r82d, & - hio_gpp_si_landuse => this%hvars(ih_gpp_si_landuse)%r82d, & hio_parsun_z_si_cnlf => this%hvars(ih_parsun_z_si_cnlf)%r82d, & hio_parsha_z_si_cnlf => this%hvars(ih_parsha_z_si_cnlf)%r82d, & hio_ts_net_uptake_si_cnlf => this%hvars(ih_ts_net_uptake_si_cnlf)%r82d, & @@ -5507,11 +5694,6 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) hio_ar_frootm_si_scpf(io_si,scpf) = hio_ar_frootm_si_scpf(io_si,scpf) + & ccohort%froot_mr * n_perm2 - if (cpatch%land_use_label .gt. nocomp_bareground_land) then - hio_gpp_si_landuse(io_si,cpatch%land_use_label) = hio_gpp_si_landuse(io_si,cpatch%land_use_label) & - + ccohort%gpp_tstep * ccohort%n * dt_tstep_inv - end if - ! accumulate fluxes on canopy- and understory- separated fluxes if (ccohort%canopy_layer .eq. 1) then @@ -6418,6 +6600,12 @@ subroutine define_history_vars(this, initialize_variables) upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & index=ih_ca_weighted_height_si) + call this%set_history_var(vname='FATES_MEAN_95PCTILE_HEIGHT', units='m', & + long='The patch-area-weighted mean of 95th percentile height of canopy plants within a given patch', use_default='active', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_patch_weighted_95thpctile_height_si) + call this%set_history_var(vname='FATES_COLD_STATUS', units='', & long='site-level cold status, 0=not cold-dec, 1=too cold for leaves, 2=not too cold', & use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & @@ -7020,6 +7208,18 @@ subroutine define_history_vars(this, initialize_variables) upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & index = ih_fire_c_to_atm_si) + call this%set_history_var(vname='FATES_FIRE_CLOSS_LIVEFUELS', units='kg m-2 s-1', & + long='carbon loss to atmosphere from live fuels only via fire in kg carbon per m2 per second', & + use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index = ih_fire_livec_to_atm_si) + + call this%set_history_var(vname='FATES_FIRE_CLOSS_LANDUSECHANGE', units='kg m-2 s-1', & + long='carbon loss to atmosphere from fire in kg carbon per m2 per second from land use change only', & + use_default='active', avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index = ih_fire_c_to_atm_landusechange_si) + call this%set_history_var(vname='FATES_CBALANCE_ERROR', & units='kg s-1', & long='total carbon error in kg carbon per second', & @@ -7411,20 +7611,6 @@ subroutine define_history_vars(this, initialize_variables) hlms='CLM:ALM', upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & index=ih_agesince_anthrodist_si_age) - call this%set_history_var(vname='FATES_SECONDARY_AREA_ANTHRO', & - units='m2 m-2', & - long='secondary forest patch area since anthropgenic disturbance', & - use_default='inactive', avgflag='A', vtype=site_r8, & - hlms='CLM:ALM', upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & - index=ih_agesince_anthrodist_si) - - call this%set_history_var(vname='FATES_SECONDARY_AREA', & - units='m2 m-2', & - long='secondary forest patch area since any kind of disturbance', & - use_default='inactive', avgflag='A', vtype=site_r8, & - hlms='CLM:ALM', upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & - index=ih_secondarylands_fracarea_si) - call this%set_history_var(vname='FATES_SECONDARY_AREA_AP', & units='m2 m-2', & long='secondary forest patch area age distribution since any kind of disturbance', & @@ -7432,12 +7618,12 @@ subroutine define_history_vars(this, initialize_variables) hlms='CLM:ALM', upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & index=ih_secondarylands_fracarea_si_age) - call this%set_history_var(vname='FATES_PRIMARY_AREA', & - units='m2 m-2', & - long='primary forest patch area since any kind of disturbance', & - use_default='inactive', avgflag='A', vtype=site_r8, & + call this%set_history_var(vname='FATES_SECONDARY_AGB_ANTHROAGE_AP', & + units='kg m-2', & + long='secondary forest patch agb as resolved by age since anthropogenic disturbance', & + use_default='inactive', avgflag='A', vtype=site_age_r8, & hlms='CLM:ALM', upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & - index=ih_primarylands_fracarea_si) + index=ih_secondary_agb_si_agesinceanthro) call this%set_history_var(vname='FATES_PRIMARY_AREA_AP', & units='m2 m-2', & @@ -8986,6 +9172,51 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & ivar=ivar, initialize=initialize_variables, index = ih_tveg_si ) + if (hlm_use_luh .eq. itrue) then + ! biophysics variables that are indexed by land use type + call this%set_history_var(vname='FATES_TVEG_LU', units='degrees Kelvin', & + long='fates instantaneous mean vegetation temperature by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_tveg_si_landuse ) + + call this%set_history_var(vname='FATES_TSA_LU', units='degrees Kelvin', & + long='fates instantaneous mean near-surface (2m) air temperature by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_tsa_si_landuse ) + + call this%set_history_var(vname='FATES_SWABS_LU', units='W m-2', & + long='fates absorbed shortwave radiation by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_sw_abs_si_landuse ) + + call this%set_history_var(vname='FATES_NETLW_LU', units='W m-2', & + long='fates net longwave flux by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_lw_net_si_landuse ) + + call this%set_history_var(vname='FATES_SHFLUX_LU', units='W m-2', & + long='fates sensible heat flux by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_shflux_si_landuse ) + + call this%set_history_var(vname='FATES_LHFLUX_LU', units='W m-2', & + long='fates latent heat flux by land use type', & + use_default='active', & + avgflag='A', vtype=site_landuse_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & + ivar=ivar, initialize=initialize_variables, index = ih_lhflux_si_landuse ) + + call this%set_history_var(vname='FATES_GPP_LU', units='kg m-2 s-1', & + long='gross primary productivity by land use type in kg carbon per m2 per second', & + use_default='inactive', avgflag='A', vtype=site_landuse_r8, & + hlms='CLM:ALM', upfreq=group_hifr_simple, ivar=ivar, initialize=initialize_variables, & + index = ih_gpp_si_landuse) + endif + call this%set_history_var(vname='FATES_VIS_RAD_ERROR', units='-', & long='mean two-stream solver error for VIS', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=group_hifr_simple, & @@ -9145,12 +9376,6 @@ subroutine define_history_vars(this, initialize_variables) hlms='CLM:ALM', upfreq=group_hifr_complx, ivar=ivar, initialize=initialize_variables, & index = ih_gpp_si_age) - call this%set_history_var(vname='FATES_GPP_LU', units='kg m-2 s-1', & - long='gross primary productivity by land use type in kg carbon per m2 per second', & - use_default='inactive', avgflag='A', vtype=site_landuse_r8, & - hlms='CLM:ALM', upfreq=group_hifr_complx, ivar=ivar, initialize=initialize_variables, & - index = ih_gpp_si_landuse) - call this%set_history_var(vname='FATES_RDARK_USTORY_SZ', & units = 'kg m-2 s-1', & long='dark respiration for understory plants in kg carbon per m2 per second by size', & diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 3ca8220ddb..7f17ea9d5c 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -327,6 +327,12 @@ subroutine zero_bcs(fates,s) fates%bc_in(s)%hksat_sisl(:) = 0.0_r8 end if + ! Diagnostic quantities for outputting FATES patch-resolved + fates%bc_in(s)%lhflux_pa(:) = 0._r8 + fates%bc_in(s)%shflux_pa(:) = 0._r8 + fates%bc_in(s)%swabs_pa(:) = 0._r8 + fates%bc_in(s)%netlw_pa(:) = 0._r8 + fates%bc_in(s)%t2m_pa(:) = 0._r8 ! Output boundaries fates%bc_out(s)%active_suction_sl(:) = .false. @@ -569,6 +575,12 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats, end if ! Land use + ! Diagnostic quantities for outputting FATES patch-resolved + allocate(bc_in%lhflux_pa(0:maxpatch_total)) + allocate(bc_in%shflux_pa(0:maxpatch_total)) + allocate(bc_in%swabs_pa(0:maxpatch_total)) + allocate(bc_in%netlw_pa(0:maxpatch_total)) + allocate(bc_in%t2m_pa(0:maxpatch_total)) ! harvest flag denote data from hlm, ! while the logging flag signifies only that logging is occurring (which could just be FATES logging) diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index c303720792..8f97353d7f 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -568,6 +568,13 @@ module FatesInterfaceTypesMod 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) + ! Diagnostic quantities for outputting FATES patch-resolved + real(r8),allocatable :: lhflux_pa(:) ! latent heat flux + real(r8),allocatable :: shflux_pa(:) ! sensible heat flux + real(r8),allocatable :: swabs_pa(:) ! shortwave absorbed radiation + real(r8),allocatable :: netlw_pa(:) ! longwave net radiation + real(r8),allocatable :: t2m_pa(:) ! 2m air temeprature + ! Land use ! --------------------------------------------------------------------------------- real(r8),allocatable :: hlm_harvest_rates(:) ! annual harvest rate per cat from hlm for a site