diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 83f10a26b5..8b5ddf79ff 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -35,6 +35,7 @@ module EDCanopyStructureMod use FatesInterfaceTypesMod , only : hlm_use_planthydro use FatesInterfaceTypesMod , only : hlm_use_cohort_age_tracking use FatesInterfaceTypesMod , only : hlm_use_sp + use FatesInterfaceTypesMod , only : hlm_use_edge_forest use FatesInterfaceTypesMod , only : numpft use FatesInterfaceTypesMod, only : bc_in_type use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage @@ -1319,6 +1320,9 @@ subroutine canopy_summarization( nsites, sites, bc_in ) use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index use EDtypesMod , only : area use FatesConstantsMod , only : itrue + use FatesEcotypesMod , only : IsPatchForest + use EDParamsMod , only : forest_tree_fraction_threshold + use FatesEdgeForestMod , only : CalculateEdgeForestArea ! !ARGUMENTS integer , intent(in) :: nsites @@ -1447,9 +1451,15 @@ subroutine canopy_summarization( nsites, sites, bc_in ) currentPatch%total_canopy_area = currentPatch%area endif + currentPatch%is_forest = IsPatchForest(currentPatch, forest_tree_fraction_threshold) + currentPatch => currentPatch%younger end do !patch loop + if (hlm_use_edge_forest == itrue) then + call CalculateEdgeForestArea(sites(s)) + end if + call leaf_area_profile(sites(s)) if(hlm_radiation_model.eq.twostr_solver) then diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 0132f7c7ef..75bf95f014 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -32,6 +32,7 @@ module FatesPatchMod use FatesRadiationMemMod, only : num_rad_stream_types use FatesInterfaceTypesMod, only : hlm_hio_ignore_val use FatesInterfaceTypesMod, only : numpft + use FatesInterfaceTypesMod, only : nlevedgeforest use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use shr_log_mod, only : errMsg => shr_log_errMsg @@ -74,6 +75,12 @@ module FatesPatchMod !--------------------------------------------------------------------------- + ! FOREST INFO + logical :: is_forest ! whether the patch is "forest" according to FATES param file criteria + real(r8), dimension(:), allocatable :: area_in_edgeforest_bins + + !--------------------------------------------------------------------------- + ! RUNNING MEANS !class(rmean_type), pointer :: t2m ! place-holder for 2m air temperature (variable window-size) class(rmean_type), pointer :: tveg24 ! 24-hour mean vegetation temperature [K] @@ -279,6 +286,7 @@ subroutine Init(this, num_swb, num_levsoil) allocate(this%sabs_dir(num_swb)) allocate(this%sabs_dif(num_swb)) allocate(this%fragmentation_scaler(num_levsoil)) + allocate(this%area_in_edgeforest_bins(nlevedgeforest)) ! initialize all values to nan call this%NanValues() @@ -453,6 +461,10 @@ subroutine NanValues(this) this%ncl_p = fates_unset_int this%land_use_label = fates_unset_int this%age_since_anthro_disturbance = nan + + ! FOREST INFO + this%is_forest = .false. + this%area_in_edgeforest_bins(:) = nan ! LEAF ORGANIZATION this%pft_agb_profile(:,:) = nan @@ -898,6 +910,7 @@ subroutine FreeMemory(this, regeneration_model, numpft) this%sabs_dir, & this%sabs_dif, & this%fragmentation_scaler, & + this%area_in_edgeforest_bins, & stat=istat, errmsg=smsg) ! These arrays are allocated via a call from EDCanopyStructureMod @@ -1227,6 +1240,7 @@ subroutine Dump(this) write(fates_log(),*) 'pa%ncl_p = ',this%ncl_p write(fates_log(),*) 'pa%total_canopy_area = ',this%total_canopy_area write(fates_log(),*) 'pa%total_tree_area = ',this%total_tree_area + write(fates_log(),*) 'pa%is_forest = ',this%is_forest write(fates_log(),*) 'pa%total_grass_area = ',this%total_grass_area write(fates_log(),*) 'pa%zstar = ',this%zstar write(fates_log(),*) 'pa%gnd_alb_dif = ',this%gnd_alb_dif(:) diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index 7e6de3323b..829d3bdbec 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -81,7 +81,7 @@ subroutine UpdateFireWeather(currentSite, bc_in) use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use FatesConstantsMod, only : sec_per_day, sec_per_min - use EDTypesMod, only : CalculateTreeGrassAreaSite + use FatesEdgeForestMod, only : CalculateTreeGrassAreaSite use FatesInterfaceTypesMod, only : hlm_use_managed_fire use SFParamsMod, only : SF_val_rxfire_tpup, SF_val_rxfire_tplw, SF_val_rxfire_rhup, & SF_val_rxfire_rhlw, SF_val_rxfire_wdup, SF_val_rxfire_wdlw diff --git a/main/CMakeLists.txt b/main/CMakeLists.txt index 8db1d06331..dbfe2fda67 100644 --- a/main/CMakeLists.txt +++ b/main/CMakeLists.txt @@ -6,6 +6,9 @@ list(APPEND clm_sources EDTypesMod.F90 EDPftvarcon.F90 FatesConstantsMod.F90 + FatesEcotypesMod.F90 + FatesEdgeForestMod.F90 + FatesEdgeForestParamsMod.F90 FatesHydraulicsMemMod.F90 FatesParametersInterface.F90 FatesUtilsMod.F90 @@ -27,6 +30,10 @@ list(APPEND fates_sources FatesSizeAgeTypeIndicesMod.F90 FatesIntegratorsMod.F90 FatesUtilsMod.F90 - FatesSynchronizedParamsMod.F90) + FatesSynchronizedParamsMod.F90 + FatesEcotypesMod.F90 + FatesEdgeForestMod.F90 + FatesEdgeForestParamsMod.F90 + ) sourcelist_to_parent(fates_sources) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 5b1ff524fc..6bf33f1306 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -213,7 +213,9 @@ module EDParamsMod integer, protected, public :: max_cohort_per_patch character(len=param_string_length), parameter, public :: maxcohort_name = "fates_maxcohort" - + ! Ecotypes parameters + real(r8),protected,public :: forest_tree_fraction_threshold ! Tree fraction above which a patch is "forest" + character(len=param_string_length),parameter,public :: forest_tree_fraction_threshold_name = "fates_forest_tree_fraction_threshold" ! Logging Control Parameters (ONLY RELEVANT WHEN USE_FATES_LOGGING = TRUE) ! ---------------------------------------------------------------------------------------------- @@ -266,6 +268,7 @@ module EDParamsMod character(len=param_string_length),parameter,public :: eca_name_plant_escalar = "fates_cnp_eca_plant_escalar" public :: FatesParamsInit + public :: FatesParamsInitForFactory public :: FatesRegisterParams public :: FatesReceiveParams public :: FatesReportParams @@ -342,6 +345,7 @@ subroutine FatesParamsInit() dev_arbitrary = nan damage_event_code = -9 damage_canopy_layer_code = -9 + forest_tree_fraction_threshold = nan landuse_grazing_carbon_use_eff = nan landuse_grazing_nitrogen_use_eff = nan landuse_grazing_phosphorus_use_eff = nan @@ -350,6 +354,21 @@ subroutine FatesParamsInit() end subroutine FatesParamsInit + !----------------------------------------------------------------------- + + subroutine FatesParamsInitForFactory() + ! Initialize some parameters that are needed for unit-testing factories + + allocate(ED_val_history_ageclass_bin_edges(7)) + ED_val_history_ageclass_bin_edges = [0, 1, 2, 5, 10, 20, 50] + + allocate(ED_val_history_sizeclass_bin_edges(13)) + ED_val_history_sizeclass_bin_edges = [0, 5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 100] + + allocate(ED_val_history_coageclass_bin_edges(13)) + ED_val_history_coageclass_bin_edges = [0, 5] + end subroutine FatesParamsInitForFactory + !----------------------------------------------------------------------- subroutine FatesRegisterParams(fates_params) ! Register the parameters we want the host to provide, and @@ -518,6 +537,9 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=damage_name_canopy_layer_code, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) + + call fates_params%RegisterParameter(name=forest_tree_fraction_threshold_name, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) call fates_params%RegisterParameter(name=name_landuse_grazing_carbon_use_eff, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -739,6 +761,9 @@ subroutine FatesReceiveParams(fates_params) data=tmpreal) damage_canopy_layer_code = nint(tmpreal) + call fates_params%RetrieveParameter(name=forest_tree_fraction_threshold_name, & + data=forest_tree_fraction_threshold) + ! parameters that are arrays of size defined within the params file and thus need allocating as well call fates_params%RetrieveParameterAllocate(name=ED_name_history_sizeclass_bin_edges, & data=ED_val_history_sizeclass_bin_edges) @@ -864,6 +889,7 @@ subroutine FatesReportParams(is_master) write(fates_log(),'(a,L2)') 'active_crown_fire = ',active_crown_fire write(fates_log(),fmt0) 'damage_event_code = ',damage_event_code write(fates_log(),fmt0) 'damage_canopy_layer_code = ', damage_canopy_layer_code + write(fates_log(),fmt0) 'forest_tree_fraction_threshold = ', forest_tree_fraction_threshold write(fates_log(),fmt0) 'landuse_grazing_carbon_use_eff = ', landuse_grazing_carbon_use_eff write(fates_log(),fmt0) 'name_landuse_grazing_nitrogen_use_eff = ', name_landuse_grazing_nitrogen_use_eff write(fates_log(),fmt0) 'name_landuse_grazing_phosphorus_use_eff = ', name_landuse_grazing_phosphorus_use_eff diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 90be4df5ec..7ab41a7a96 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -605,12 +605,12 @@ module EDTypesMod procedure, public :: get_current_landuse_statevector procedure, public :: get_secondary_young_fraction + procedure, public :: GetNumberOfPatches end type ed_site_type ! Make public necessary subroutines and functions public :: dump_site - public :: CalculateTreeGrassAreaSite public :: set_patchno contains @@ -737,39 +737,6 @@ end subroutine ZeroMassBalFlux ! ===================================================================================== - subroutine CalculateTreeGrassAreaSite(csite, tree_fraction, grass_fraction, bare_fraction) - ! - ! DESCRIPTION: - ! Calculates total grass, tree, and bare fractions for a site - - ! ARGUMENTS: - type(ed_site_type), intent(inout) :: csite ! site object - real(r8), intent(out) :: tree_fraction ! total site tree fraction - real(r8), intent(out) :: grass_fraction ! total site grass fraction - real(r8), intent(out) :: bare_fraction ! total site bare fraction - - ! LOCALS: - type(fates_patch_type), pointer :: currentPatch ! patch object - - tree_fraction = 0.0_r8 - grass_fraction = 0.0_r8 - - currentPatch => csite%oldest_patch - do while(associated(currentPatch)) - if (currentPatch%nocomp_pft_label /= nocomp_bareground) then - call currentPatch%UpdateTreeGrassArea() - tree_fraction = tree_fraction + currentPatch%total_tree_area/AREA - grass_fraction = grass_fraction + currentPatch%total_grass_area/AREA - end if - currentPatch => currentPatch%younger - end do - - ! if cover > 1.0, grasses are under the trees - grass_fraction = min(grass_fraction, 1.0_r8 - tree_fraction) - bare_fraction = 1.0_r8 - tree_fraction - grass_fraction - - end subroutine CalculateTreeGrassAreaSite - !--------------------------------------------------------------------------------------- subroutine dump_site(csite) @@ -865,4 +832,26 @@ function get_secondary_young_fraction(this) result(secondary_young_fraction) end function get_secondary_young_fraction + function GetNumberOfPatches(this) result(n_patches) + ! DESCRIPTION + ! Returns number of patches at site + ! + ! ARGUMENTS: + class(ed_site_type) :: this + ! + ! RETURN VALUE: + integer :: n_patches + ! + ! LOCAL VARIABLES: + type(fates_patch_type), pointer :: currentPatch + + n_patches = 0 + currentPatch => this%youngest_patch + do while(associated(currentPatch)) + n_patches = n_patches + 1 + currentPatch => currentPatch%older + enddo + + end function GetNumberOfPatches + end module EDTypesMod diff --git a/main/FatesEcotypesMod.F90 b/main/FatesEcotypesMod.F90 new file mode 100644 index 0000000000..911bed7db6 --- /dev/null +++ b/main/FatesEcotypesMod.F90 @@ -0,0 +1,81 @@ +module FatesEcotypesMod + + use FatesConstantsMod, only : r8 => fates_r8 + use EDTypesMod, only : ed_site_type + use FatesPatchMod, only : fates_patch_type + + implicit none + private ! By default everything is private + + ! Make public necessary subroutines and functions + public :: IsPatchForest + ! For unit testing + public :: DoesPatchHaveForest_TreeCover + public :: DoesPatchHaveForest_GrassBiomass + +contains + + ! ===================================================================================== + + function DoesPatchHaveForest_TreeCover(patchptr, forest_tree_fraction_threshold) + ! DESCRIPTION: + ! Return boolean: Is this patch "forest"? + ! + ! ARGUMENTS: + type(fates_patch_type), intent(in), pointer :: patchptr ! pointer to patch object + real(r8), intent(in) :: forest_tree_fraction_threshold ! Tree fraction above which a patch is "forest" + ! + ! RETURN VALUE + logical :: DoesPatchHaveForest_TreeCover + ! + ! LOCAL VARIABLES + real(r8) :: tree_fraction = 0._r8 + + if (patchptr%area > 0._r8) then + tree_fraction = patchptr%total_tree_area / patchptr%area + else + tree_fraction = 0._r8 + end if + + DoesPatchHaveForest_TreeCover = tree_fraction > forest_tree_fraction_threshold + + end function DoesPatchHaveForest_TreeCover + + + function DoesPatchHaveForest_GrassBiomass(patchptr, grass_biomass_threshold) + ! DESCRIPTION: + ! Return boolean: Does this patch have grass biomass above a threshold? + ! + ! ARGUMENTS: + type(fates_patch_type), intent(in), pointer :: patchptr ! pointer to patch object + real(r8), intent(in) :: grass_biomass_threshold ! Live grass biomass (kgC/m2) above which a patch is considered to "have grass" + ! + ! RETURN VALUE + logical :: DoesPatchHaveForest_GrassBiomass + + DoesPatchHaveForest_GrassBiomass = patchptr%livegrass > grass_biomass_threshold + + end function DoesPatchHaveForest_GrassBiomass + + + function IsPatchForest(patchptr, forest_tree_fraction_threshold, grass_biomass_threshold) + ! DESCRIPTION: + ! Return boolean: Is this patch forest according to tree cover and, optionally, grass biomass? + ! + ! ARGUMENTS: + type(fates_patch_type), intent(in), pointer :: patchptr ! pointer to patch object + real(r8), intent(in) :: forest_tree_fraction_threshold ! Tree fraction above which a patch is "forest" + real(r8), intent(in), optional :: grass_biomass_threshold ! Live grass biomass (kgC/m2) above which a patch is considered to "have grass" + ! + ! RETURN VALUE + logical :: IsPatchForest + + IsPatchForest = DoesPatchHaveForest_TreeCover(patchptr, forest_tree_fraction_threshold) + if (IsPatchForest .and. present(grass_biomass_threshold)) then + IsPatchForest = .not. DoesPatchHaveForest_GrassBiomass(patchptr, grass_biomass_threshold) + end if + + end function IsPatchForest + + +end module FatesEcotypesMod diff --git a/main/FatesEdgeForestMod.F90 b/main/FatesEdgeForestMod.F90 new file mode 100644 index 0000000000..4efceab94d --- /dev/null +++ b/main/FatesEdgeForestMod.F90 @@ -0,0 +1,630 @@ +module FatesEdgeForestMod + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : nocomp_bareground + use FatesConstantsMod, only : nearzero + use FatesGlobals, only : fates_log + use FatesGlobals, only : endrun => fates_endrun + use shr_log_mod, only : errMsg => shr_log_errMsg + use EDTypesMod, only : ed_site_type + use EDTypesMod, only : AREA + use FatesPatchMod, only : fates_patch_type + use FatesEcotypesMod, only : IsPatchForest + use FatesUtilsMod, only : is_param_set + + implicit none + private ! By default everything is private + + ! Make public necessary subroutines and functions + public :: CalculateTreeGrassAreaSite + public :: CalculateEdgeForestArea + ! Public for unit testing + public :: indexx + public :: GetFracEdgeForestInEachBin + public :: GetFracEdgeForestInEachBin_norm_numerator + public :: GetFracEdgeForestInEachBin_norm_denominator + public :: GetFracEdgeForestInEachBin_norm + public :: GetFracEdgeForestInEachBin_quadratic + public :: AssignPatchToBins + +contains + + ! ===================================================================================== + + subroutine CalculateTreeGrassAreaSite(csite, tree_fraction, grass_fraction, bare_fraction) + ! + ! DESCRIPTION: + ! Calculates total grass, tree, and bare fractions for a site + + use FatesEcotypesMod, only : IsPatchForest + use EDParamsMod, only : forest_tree_fraction_threshold + + ! ARGUMENTS: + type(ed_site_type), intent(inout), target :: csite ! site object + real(r8), intent(out) :: tree_fraction ! total site tree fraction + real(r8), intent(out) :: grass_fraction ! total site grass fraction + real(r8), intent(out) :: bare_fraction ! total site bare fraction + + ! LOCALS: + type(fates_patch_type), pointer :: currentPatch ! patch object + + tree_fraction = 0.0_r8 + grass_fraction = 0.0_r8 + + currentPatch => csite%oldest_patch + do while(associated(currentPatch)) + if (currentPatch%nocomp_pft_label /= nocomp_bareground) then + call currentPatch%UpdateTreeGrassArea() + tree_fraction = tree_fraction + currentPatch%total_tree_area/AREA + grass_fraction = grass_fraction + currentPatch%total_grass_area/AREA + currentPatch%is_forest = IsPatchForest(currentPatch, forest_tree_fraction_threshold) + end if + currentPatch => currentPatch%younger + end do + + ! if cover > 1.0, grasses are under the trees + grass_fraction = min(grass_fraction, 1.0_r8 - tree_fraction) + bare_fraction = 1.0_r8 - tree_fraction - grass_fraction + + ! Must come after patch loop with IsPatchForest() call + call CalculateEdgeForestArea(csite) + + end subroutine CalculateTreeGrassAreaSite + + + subroutine GetNumberOfForestPatches(site, n_forest_patches, area_forest_patches, area_site) + ! DESCRIPTION + ! Returns number and area of forest patches at site + ! + ! ARGUMENTS: + type(ed_site_type), pointer, intent(in) :: site + integer, intent(out) :: n_forest_patches + real(r8), intent(out) :: area_forest_patches + real(r8), intent(out) :: area_site + ! + ! LOCAL VARIABLES: + type(fates_patch_type), pointer :: currentPatch + + n_forest_patches = 0 + area_forest_patches = 0._r8 + area_site = 0._r8 + currentPatch => site%youngest_patch + do while(associated(currentPatch)) + area_site = area_site + currentPatch%area + if (currentPatch%is_forest) then + n_forest_patches = n_forest_patches + 1 + area_forest_patches = area_forest_patches + currentPatch%area + end if + + currentPatch => currentPatch%older + enddo + + end subroutine GetNumberOfForestPatches + + + subroutine RankForestEdgeProximity(site, indices, index_forestpatches_to_allpatches) + ! DESCRIPTION: + ! Rank forest patches by their proximity to edge, using age as a proxy. + ! + ! ARGUMENTS: + type(ed_site_type), pointer, intent(in) :: site + integer, dimension(:), intent(inout) :: indices ! Indices to use if you want to sort patches + integer, dimension(:), intent(inout) :: index_forestpatches_to_allpatches ! Array with length (number of patches in gridcell), values 0 if not forest and otherwise an index corresponding to which number forest patch this is + ! + ! LOCAL VARIABLES: + real(r8), dimension(:), allocatable :: array ! Array to be index-sorted. + integer :: n_forest_patches ! Number of patches in above arrays + type(fates_patch_type), pointer :: currentPatch + integer :: f ! index of current forest patch + integer :: p ! index of patch + + ! Skip sites with no forest patches + n_forest_patches = size(indices) + if (n_forest_patches == 0) then + return + end if + + ! Allocate arrays + allocate(array(1:n_forest_patches)) + + ! Fill arrays + f = 0 + p = 0 + index_forestpatches_to_allpatches(:) = 0 + currentPatch => site%oldest_patch + patchloop: do while(associated(currentPatch)) + p = p + 1 + if (.not. currentPatch%is_forest) then + currentPatch => currentPatch%younger + cycle + end if + + f = f + 1 + index_forestpatches_to_allpatches(p) = f + + ! Fill with patch age. + ! TODO: Add other options. Biomass? Woody biomass? + array(f) = currentPatch%age + + currentPatch => currentPatch%younger + end do patchloop + + ! Get indices of sorted forest patches + call indexx(array, indices) + + ! Clean up + deallocate(array) + end subroutine RankForestEdgeProximity + + function GetFracEdgeForestInEachBin_norm_numerator(x_in, A, mu, sigma, lognorm) + ! DESCRIPTION + ! Gets numerator at xof normal-like function (Gaussian if lognorm==.true., lognormal otherwise) + ! + ! ARGUMENTS: + real(r8), intent(in) :: x_in + real(r8), intent(in) :: A ! Amplitude + real(r8), intent(in) :: mu ! Center + real(r8), intent(in) :: sigma ! Width + logical, intent(in) :: lognorm ! Gaussian function if true, lognormal otherwise + ! + ! RETURN VALUE: + real(r8) :: GetFracEdgeForestInEachBin_norm_numerator + ! + ! LOCAL VARIABLES: + real(r8) :: x ! either x_in or its log + + if (lognorm) then + x = log(x_in) + else + x = x_in + end if + + GetFracEdgeForestInEachBin_norm_numerator = A * exp(-(x - mu)**2 / (2*sigma**2)) + end function GetFracEdgeForestInEachBin_norm_numerator + + function GetFracEdgeForestInEachBin_norm_denominator(x, sigma, lognorm) + ! DESCRIPTION + ! Gets denominator at x of normal-like function (Gaussian if lognorm==.true., lognormal otherwise) + ! + ! ARGUMENTS: + use FatesConstantsMod, only : pi => pi_const + real(r8), intent(in) :: x + real(r8), intent(in) :: sigma ! Width + logical, intent(in) :: lognorm ! Gaussian function if true, lognormal otherwise + ! + ! RETURN VALUE: + real(r8) :: GetFracEdgeForestInEachBin_norm_denominator + + GetFracEdgeForestInEachBin_norm_denominator = sigma * sqrt(2*pi) + if (lognorm) then + GetFracEdgeForestInEachBin_norm_denominator = GetFracEdgeForestInEachBin_norm_denominator * x + end if + end function GetFracEdgeForestInEachBin_norm_denominator + + function GetFracEdgeForestInEachBin_norm(x, A, mu, sigma, lognorm) + ! DESCRIPTION + ! Gets value at x of normal-like function (Gaussian if lognorm==.true., lognormal otherwise) + ! + ! ARGUMENTS: + real(r8), intent(in) :: x + real(r8), intent(in) :: A ! Amplitude + real(r8), intent(in) :: mu ! Center + real(r8), intent(in) :: sigma ! Width + logical, intent(in) :: lognorm ! Gaussian function if true, lognormal otherwise + ! + ! RETURN VALUE: + real(r8) :: GetFracEdgeForestInEachBin_norm + + GetFracEdgeForestInEachBin_norm = GetFracEdgeForestInEachBin_norm_numerator(x, A, mu, sigma, lognorm) / GetFracEdgeForestInEachBin_norm_denominator(x, sigma, lognorm) + end function GetFracEdgeForestInEachBin_norm + + function GetFracEdgeForestInEachBin_quadratic(x, a, b, c) + ! DESCRIPTION + ! Gets value at x of quadratic function + ! + ! ARGUMENTS: + real(r8), intent(in) :: x + real(r8), intent(in) :: a, b, c ! Parameters + ! + ! RETURN VALUE: + real(r8) :: GetFracEdgeForestInEachBin_quadratic + + GetFracEdgeForestInEachBin_quadratic = a*(x**2) + b*x + c + end function GetFracEdgeForestInEachBin_quadratic + + subroutine GetFracEdgeForestInEachBin(x, nlevedgeforest, efb_gaussian_amplitudes, efb_gaussian_sigmas, efb_gaussian_centers, efb_lognormal_amplitudes, efb_lognormal_sigmas, efb_lognormal_centers, efb_quadratic_a, efb_quadratic_b, efb_quadratic_c, fraction_forest_in_bin, norm) + ! DESCRIPTION: + ! Get the fraction of forest in each bin. + ! + ! USES + use FatesConstantsMod, only : pi => pi_const + ! + ! ARGUMENTS + real(r8), intent(in) :: x ! Independent variable in the fit + integer, intent(in) :: nlevedgeforest + real(r8), dimension(:), intent(in) :: efb_gaussian_amplitudes + real(r8), dimension(:), intent(in) :: efb_gaussian_sigmas + real(r8), dimension(:), intent(in) :: efb_gaussian_centers + real(r8), dimension(:), intent(in) :: efb_lognormal_amplitudes + real(r8), dimension(:), intent(in) :: efb_lognormal_sigmas + real(r8), dimension(:), intent(in) :: efb_lognormal_centers + real(r8), dimension(:), intent(in) :: efb_quadratic_a + real(r8), dimension(:), intent(in) :: efb_quadratic_b + real(r8), dimension(:), intent(in) :: efb_quadratic_c + real(r8), dimension(:), intent(out) :: fraction_forest_in_bin + logical, optional, intent(in) :: norm + ! + ! LOCAL VARIABLES + integer :: b ! Bin index + real(r8) :: A ! Amplitude + real(r8) :: mu ! Center + real(r8) :: sigma ! Sigma + logical :: lognorm + logical :: do_norm + ! Error checking + real(r8), parameter :: tol = 1.e-9_r8 ! fraction of total forest area + real(r8) :: err_chk + + ! Initialize + fraction_forest_in_bin(:) = 0._r8 + + ! If the cell is nearly 0% forest, put any forest in the first edge bin (closest to edge) + if (x < nearzero) then + fraction_forest_in_bin(1) = 1._r8 + return + end if + + ! If the cell is (nearly) 100% forest, it's all "deep forest" + if (1._r8 - x < nearzero) then + fraction_forest_in_bin(nlevedgeforest) = 1._r8 + return + end if + + if (present(norm)) then + do_norm = norm + else + do_norm = .true. + end if + + binloop: do b = 1, nlevedgeforest + + if (is_param_set(efb_gaussian_amplitudes(b)) .or. is_param_set(efb_lognormal_amplitudes(b))) then + ! Gaussian or Lognormal + lognorm = is_param_set(efb_lognormal_amplitudes(b)) + if (lognorm) then + A = efb_lognormal_amplitudes(b) + mu = efb_lognormal_centers(b) + sigma = efb_lognormal_sigmas(b) + else + A = efb_gaussian_amplitudes(b) + mu = efb_gaussian_centers(b) + sigma = efb_gaussian_sigmas(b) + end if + fraction_forest_in_bin(b) = GetFracEdgeForestInEachBin_norm(x, A, mu, sigma, lognorm) + + else if (is_param_set(efb_quadratic_a(b))) then + ! Quadratic + fraction_forest_in_bin(b) = GetFracEdgeForestInEachBin_quadratic(x, efb_quadratic_a(b), efb_quadratic_b(b), efb_quadratic_c(b)) + + else + call endrun("Unrecognized bin fit type") + end if + end do binloop + + ! Account for fit errors by normalizing to 1 + if (do_norm) then + fraction_forest_in_bin(:) = fraction_forest_in_bin(:) / sum(fraction_forest_in_bin) + + err_chk = sum(fraction_forest_in_bin) - 1._r8 + if (abs(err_chk) > tol) then + write(fates_log(),*) "ERROR: bin fractions don't sum to 1; actual minus expected = ",err_chk + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + end if + + + end subroutine GetFracEdgeForestInEachBin + + + subroutine AssignPatchToBins(fraction_forest_in_each_bin, area_forest_patches, patch_area, nlevedgeforest, tol, sum_forest_bins_so_far_m2, area_in_edgeforest_bins_m2) + ! DESCRIPTION + ! Given one patch in a site, assign its area to edge bin(s). + ! + ! ARGUMENTS + real(r8), dimension(:), pointer, intent(in) :: fraction_forest_in_each_bin ! Fraction of site's forest area in each edge bin + real(r8), intent(in) :: area_forest_patches ! Total forest area in the site + real(r8), intent(in) :: patch_area ! Area of this patch + integer, intent(in) :: nlevedgeforest ! Number of edge forest bins, including "deep forest" + real(r8), intent(in) :: tol ! Tolerance for checking total area assigned to bins + real(r8), intent(inout) :: sum_forest_bins_so_far_m2 ! How much of the site's forest area has been assigned? + real(r8), dimension(:), intent(out) :: area_in_edgeforest_bins_m2 ! Area of this patch in each edge bin (m2) + ! + ! LOCAL VARIABLES + real(r8) :: remaining_to_assign_from_patch_m2 ! How much of this patch's area still needs to be assigned + real(r8) :: remaining_to_assign_to_bin_m2 ! How much of a given bin's area still needs to be assigned + integer :: b + ! For checks + real(r8) :: err_chk + + area_in_edgeforest_bins_m2(:) = 0._r8 + remaining_to_assign_from_patch_m2 = patch_area + binloop: do b = 1, nlevedgeforest + + ! How much area is left for this bin? + remaining_to_assign_to_bin_m2 = sum(fraction_forest_in_each_bin(1:b))*area_forest_patches - sum_forest_bins_so_far_m2 + if (remaining_to_assign_to_bin_m2 <= 0) then + cycle + end if + + ! Assign area + area_in_edgeforest_bins_m2(b) = min(remaining_to_assign_from_patch_m2, remaining_to_assign_to_bin_m2) + remaining_to_assign_from_patch_m2 = remaining_to_assign_from_patch_m2 - area_in_edgeforest_bins_m2(b) + + ! Update accounting + sum_forest_bins_so_far_m2 = sum_forest_bins_so_far_m2 + area_in_edgeforest_bins_m2(b) + + if (remaining_to_assign_from_patch_m2 == 0._r8) then + exit + end if + end do binloop + + ! Check that this patch's complete area was assigned (and no more) + err_chk = remaining_to_assign_from_patch_m2 + if (abs(err_chk) > tol) then + write(fates_log(),*) "ERROR: not enough or too much patch area was assigned to bins (check 1); remainder = ",err_chk + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + err_chk = patch_area - sum(area_in_edgeforest_bins_m2) + if (abs(err_chk) > tol) then + write(fates_log(),*) "ERROR: not enough or too much patch area was assigned to bins (check 2); remainder = ",err_chk + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + end subroutine AssignPatchToBins + + + subroutine AssignPatchesToBins(site, indices, index_forestpatches_to_allpatches, fraction_forest_in_each_bin, n_forest_patches, n_patches, area_forest_patches) + ! DESCRIPTION + ! Loops through forest patches from nearest to farthest from edge, assigning their + ! area to edge bin(s). + ! + ! USES: + use FatesInterfaceTypesMod, only : nlevedgeforest + ! ARGUMENTS + type(ed_site_type), pointer, intent(in) :: site + integer, dimension(:), intent(in) :: indices ! Indices to use if you want to sort patches + integer, dimension(:), intent(in) :: index_forestpatches_to_allpatches ! Array with length (number of patches in gridcell), values 0 if not forest and otherwise an index corresponding to which number forest patch this is + real(r8), dimension(:), pointer, intent(in) :: fraction_forest_in_each_bin + integer, intent(in) :: n_forest_patches ! Number of forest patches + integer, intent(in) :: n_patches ! Number of patches in site + real(r8), intent(in) :: area_forest_patches + ! + ! LOCAL VARIABLES + integer :: f, i, p, b + type(fates_patch_type), pointer :: currentPatch + real(r8) :: sum_forest_bins_so_far_m2 + ! For checks + real(r8), dimension(nlevedgeforest) :: bin_area_sums + real(r8), parameter :: tol = 1.e-9_r8 ! m2 + real(r8) :: err_chk + + sum_forest_bins_so_far_m2 = 0._r8 + forestpatchloop: do f = 1, n_forest_patches + + ! Get the i'th patch (which is the f'th forest patch) + i = indices(f) + currentPatch => site%oldest_patch + allpatchloop: do p = 1, n_patches + if (index_forestpatches_to_allpatches(p) == i) then + exit + end if + currentPatch => currentPatch%younger + end do allpatchloop + if ((.not. associated(currentPatch)) .and. (.not. p == n_patches)) then + write(fates_log(),*) "ERROR: i'th patch (f'th forest patch) not found." + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + + ! Make sure this is a forest patch + if (.not. currentPatch%is_forest) then + write(fates_log(),*) "ERROR: unexpected non-forest patch in forestpatchloop" + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + + ! Assign this patch's area + call AssignPatchToBins(fraction_forest_in_each_bin, area_forest_patches, currentPatch%area, nlevedgeforest, tol, sum_forest_bins_so_far_m2, currentPatch%area_in_edgeforest_bins) + + end do forestpatchloop + + ! More checks + bin_area_sums(:) = 0._r8 + currentPatch => site%oldest_patch + allpatchloop_check: do while (associated(currentPatch)) + if (currentPatch%is_forest) then + + ! Check that all area of each forest patch is assigned + err_chk = sum(currentPatch%area_in_edgeforest_bins) - currentPatch%area + if (abs(err_chk) > tol) then + write(fates_log(),*) "ERROR: unexpected patch forest bin sum (check 3); actual minus expected = ",err_chk + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + + ! Accumulate site-wide area in each bin + binloop_check4a: do b = 1, nlevedgeforest + bin_area_sums(b) = bin_area_sums(b) + currentPatch%area_in_edgeforest_bins(b) / area_forest_patches + end do binloop_check4a + + end if + currentPatch => currentPatch%younger + end do allpatchloop_check + ! Check that fraction in each bin is what was expected + binloop_check4b: do b = 1, nlevedgeforest + err_chk = bin_area_sums(b) - fraction_forest_in_each_bin(b) + if (abs(err_chk) > tol) then + write(fates_log(),*) "ERROR: unexpected bin sum (check 4); actual minus expected = ",err_chk + call endrun(msg=errMsg(__FILE__, __LINE__)) + end if + end do binloop_check4b + end subroutine AssignPatchesToBins + + + subroutine CalculateEdgeForestArea(site) + ! DESCRIPTION: + ! Loop through forest patches in decreasing order of proximity, calculating the + ! area of each patch that is in each edge bin. + ! + ! USES: + use FatesInterfaceTypesMod, only : nlevedgeforest + use FatesEdgeForestParamsMod, only : ED_val_edgeforest_gaussian_amplitude, ED_val_edgeforest_gaussian_sigma, ED_val_edgeforest_gaussian_center + use FatesEdgeForestParamsMod, only : ED_val_edgeforest_lognormal_amplitude, ED_val_edgeforest_lognormal_sigma, ED_val_edgeforest_lognormal_center + use FatesEdgeForestParamsMod, only : ED_val_edgeforest_quadratic_a, ED_val_edgeforest_quadratic_b, ED_val_edgeforest_quadratic_c + use FatesEdgeForestParamsMod, only : ED_val_edgeforest_bin_edges + ! + ! ARGUMENTS: + type(ed_site_type), pointer, intent(in) :: site + ! + ! LOCAL VARIABLES: + integer, dimension(:), allocatable :: indices ! Indices to use if you want to sort patches + integer, dimension(:), allocatable :: index_forestpatches_to_allpatches ! Array with length (number of patches in gridcell), values 0 if not forest and otherwise an index corresponding to which number forest patch this is + integer :: n_forest_patches ! Number of forest patches + integer :: n_patches ! Number of patches in site + real(r8) :: area_forest_patches + real(r8) :: area_site + real(r8) :: frac_forest + real(r8), dimension(nlevedgeforest), target :: fraction_forest_in_each_bin + + ! Skip sites with no forest patches + call GetNumberOfForestPatches(site, n_forest_patches, area_forest_patches, area_site) + if (n_forest_patches == 0) then + return + end if + + ! Allocate arrays + allocate(indices(1:n_forest_patches)) + n_patches = site%GetNumberOfPatches() + allocate(index_forestpatches_to_allpatches(1:n_patches)) + + ! Rank forest patches by their proximity to edge + call RankForestEdgeProximity(site, indices, index_forestpatches_to_allpatches) + + ! Get fraction of forest area in each bin + frac_forest = area_forest_patches / area_site + call GetFracEdgeForestInEachBin(frac_forest, nlevedgeforest, ED_val_edgeforest_gaussian_amplitude, ED_val_edgeforest_gaussian_sigma, ED_val_edgeforest_gaussian_center, ED_val_edgeforest_lognormal_amplitude, ED_val_edgeforest_lognormal_sigma, ED_val_edgeforest_lognormal_center, ED_val_edgeforest_quadratic_a, ED_val_edgeforest_quadratic_b, ED_val_edgeforest_quadratic_c, fraction_forest_in_each_bin) + + ! Assign patches to bins + call AssignPatchesToBins(site, indices, index_forestpatches_to_allpatches, fraction_forest_in_each_bin, n_forest_patches, n_patches, area_forest_patches) + + ! Clean up + deallocate(indices) + deallocate(index_forestpatches_to_allpatches) + end subroutine CalculateEdgeForestArea + + + + + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + ! + ! The following two subroutines perform an index-sort of an array. + ! They are a GPL-licenced replacement for the Numerical Recipes routine indexx. + ! They are not derived from any NR code, but are based on a quicksort routine by + ! Michael Lamont (http://linux.wku.edu/~lamonml/kb.html), originally written + ! in C, and issued under the GNU General Public License. The conversion to + ! Fortran 90, and modification to do an index sort was done by Ian Rutt. + ! + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + subroutine indexx(array, index) + + ! Performs an index sort of \texttt{array} and returns the result in + ! \texttt{index}. The order of elements in \texttt{array} is unchanged. + ! + ! This is a GPL-licenced replacement for the Numerical Recipes routine indexx. + ! It is not derived from any NR code, but are based on a quicksort routine by + ! Michael Lamont (http://linux.wku.edu/~lamonml/kb.html), originally written + ! in C, and issued under the GNU General Public License. The conversion to + ! Fortran 90, and modification to do an index sort was done by Ian Rutt. + + real(r8), dimension(:) :: array ! Array to be indexed. + integer, dimension(:) :: index ! Index of elements of patch_array + integer :: i + + if (size(array) /= size(index)) then + write(fates_log(),*) 'ERROR: INDEXX size mismatch.' + call endrun(msg=errMsg(__FILE__, __LINE__)) + else if (size(array) == 0) then + write(fates_log(),*) 'ERROR: INDEXX array size 0.' + call endrun(msg=errMsg(__FILE__, __LINE__)) + endif + + do i=1,size(index) + index(i)=i + enddo + + call q_sort_index(array,index,1,size(array)) + + end subroutine indexx + +!============================================================== + + recursive subroutine q_sort_index(numbers,index,left,right) + + !> This is the recursive subroutine actually used by sort_patches. + !> + !> This is a GPL-licenced replacement for the Numerical Recipes routine indexx. + !> It is not derived from any NR code, but are based on a quicksort routine by + !> Michael Lamont (http://linux.wku.edu/~lamonml/kb.html), originally written + !> in C, and issued under the GNU General Public License. The conversion to + !> Fortran 90, and modification to do an index sort was done by Ian Rutt. + + implicit none + + real(r8), dimension(:) :: numbers !> Numbers being sorted + integer, dimension(:) :: index !> Returned index + integer :: left, right !> Limit of sort region + + integer :: ll,rr + integer :: pv_int,l_hold, r_hold,pivpos + real(r8) :: pivot + + ll=left + rr=right + + l_hold = ll + r_hold = rr + pivot = numbers(index(ll)) + pivpos=index(ll) + + do + if (.not.(ll < rr)) exit + + do + if (.not.((numbers(index(rr)) >= pivot) .and. (ll < rr))) exit + rr=rr-1 + enddo + + if (ll /= rr) then + index(ll) = index(rr) + ll=ll+1 + endif + + do + if (.not.((numbers(index(ll)) <= pivot) .and. (ll < rr))) exit + ll=ll+1 + enddo + + if (ll /= rr) then + index(rr) = index(ll) + rr=rr-1 + endif + enddo + + index(ll) = pivpos + pv_int = ll + ll = l_hold + rr = r_hold + if (ll < pv_int) call q_sort_index(numbers, index,ll, pv_int-1) + if (rr > pv_int) call q_sort_index(numbers, index,pv_int+1, rr) + + end subroutine q_sort_index + + + end module FatesEdgeForestMod diff --git a/main/FatesEdgeForestParamsMod.F90 b/main/FatesEdgeForestParamsMod.F90 new file mode 100644 index 0000000000..6c9d8778c6 --- /dev/null +++ b/main/FatesEdgeForestParamsMod.F90 @@ -0,0 +1,224 @@ +module FatesEdgeForestParamsMod + ! + ! module that deals with reading the edge forest parameter file + ! + use FatesConstantsMod, only : r8 => fates_r8 + use FatesInterfaceTypesMod, only : nlevedgeforest + use FatesParametersInterface, only : param_string_length + use FatesGlobals, only : fates_log + use FatesGlobals, only : endrun => fates_endrun + use FatesUtilsMod, only : is_param_set + use shr_log_mod, only : errMsg => shr_log_errMsg + + implicit none + private + save + + ! + ! this is what the user can use for the actual values + ! + real(r8),protected,allocatable,public :: ED_val_edgeforest_gaussian_amplitude(:) + real(r8),protected,allocatable,public :: ED_val_edgeforest_gaussian_sigma(:) + real(r8),protected,allocatable,public :: ED_val_edgeforest_gaussian_center(:) + real(r8),protected,allocatable,public :: ED_val_edgeforest_lognormal_amplitude(:) + real(r8),protected,allocatable,public :: ED_val_edgeforest_lognormal_sigma(:) + real(r8),protected,allocatable,public :: ED_val_edgeforest_lognormal_center(:) + real(r8),protected,allocatable,public :: ED_val_edgeforest_quadratic_a(:) + real(r8),protected,allocatable,public :: ED_val_edgeforest_quadratic_b(:) + real(r8),protected,allocatable,public :: ED_val_edgeforest_quadratic_c(:) + real(r8),protected,allocatable,public :: ED_val_edgeforest_bin_edges(:) + + character(len=param_string_length),parameter,public :: ED_name_edgeforest_gaussian_amplitude = "fates_edgeforest_gaussian_amplitude" + character(len=param_string_length),parameter,public :: ED_name_edgeforest_gaussian_sigma = "fates_edgeforest_gaussian_sigma" + character(len=param_string_length),parameter,public :: ED_name_edgeforest_gaussian_center = "fates_edgeforest_gaussian_center" + character(len=param_string_length),parameter,public :: ED_name_edgeforest_lognormal_amplitude = "fates_edgeforest_lognormal_amplitude" + character(len=param_string_length),parameter,public :: ED_name_edgeforest_lognormal_sigma = "fates_edgeforest_lognormal_sigma" + character(len=param_string_length),parameter,public :: ED_name_edgeforest_lognormal_center = "fates_edgeforest_lognormal_center" + character(len=param_string_length),parameter,public :: ED_name_edgeforest_quadratic_a = "fates_edgeforest_quadratic_a" + character(len=param_string_length),parameter,public :: ED_name_edgeforest_quadratic_b = "fates_edgeforest_quadratic_b" + character(len=param_string_length),parameter,public :: ED_name_edgeforest_quadratic_c = "fates_edgeforest_quadratic_c" + character(len=param_string_length),parameter,public :: ED_name_edgeforest_bin_edges = "fates_edgeforest_bin_edges" + + character(len=*), parameter, private :: sourcefile = __FILE__ + + public :: EdgeForestRegisterParams + public :: EdgeForestReceiveParams + public :: EdgeForestCheckParams + +contains + + ! ===================================================================================== + + subroutine check_all_unset(value_array, b) + real(r8), dimension(:), intent(in) :: value_array + integer, intent(in) :: b + integer :: i + do i = 1, size(value_array) + if (is_param_set(value_array(i))) then + write(fates_log(),*) 'Multiple fit types found for bin ',b + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do + end subroutine check_all_unset + + subroutine EdgeForestCheckParams(is_master) + + ! ---------------------------------------------------------------------------------- + ! + ! This subroutine performs logical checks on user supplied parameters. It cross + ! compares various parameters and will fail if they don't make sense. + ! E.g.: Each bin should have parameters for exactly one fit type. + ! ----------------------------------------------------------------------------------- + + logical, intent(in) :: is_master ! Only check if this is the master proc + + real(r8) :: gaussian_amplitude + real(r8) :: gaussian_sigma + real(r8) :: gaussian_center + real(r8) :: lognormal_amplitude + real(r8) :: lognormal_sigma + real(r8) :: lognormal_center + real(r8) :: quadratic_a + real(r8) :: quadratic_b + real(r8) :: quadratic_c + integer :: b + + if (.not. is_master) return + + ! Check each bin + do b = 1, nlevedgeforest + + gaussian_amplitude = ED_val_edgeforest_gaussian_amplitude(b) + gaussian_sigma = ED_val_edgeforest_gaussian_sigma(b) + gaussian_center = ED_val_edgeforest_gaussian_center(b) + lognormal_amplitude = ED_val_edgeforest_lognormal_amplitude(b) + lognormal_sigma = ED_val_edgeforest_lognormal_sigma(b) + lognormal_center = ED_val_edgeforest_lognormal_center(b) + quadratic_a = ED_val_edgeforest_quadratic_a(b) + quadratic_b = ED_val_edgeforest_quadratic_b(b) + quadratic_c = ED_val_edgeforest_quadratic_c(b) + + if (is_param_set(gaussian_amplitude)) then + ! Has all gaussian parameters + if (.not. (is_param_set(gaussian_center) .and. is_param_set(gaussian_sigma))) then + write(fates_log(),*) 'Not all gaussian forest edge parameters found for bin ',b + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + ! Has no other parameters + call check_all_unset( (/ lognormal_amplitude, lognormal_sigma, lognormal_center, quadratic_a, quadratic_b, quadratic_c /), b ) + + else if (is_param_set(lognormal_amplitude)) then + ! Has all lognormal parameters + if (.not. (is_param_set(lognormal_center) .and. is_param_set(lognormal_sigma))) then + write(fates_log(),*) 'Not all lognormal forest edge parameters found for bin ',b + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + ! Has no other parameters + call check_all_unset( (/ gaussian_amplitude, gaussian_sigma, gaussian_center, quadratic_a, quadratic_b, quadratic_c /), b ) + + else if (is_param_set(quadratic_a)) then + ! Has all quadratic parameters + if (.not. (is_param_set(quadratic_b) .and. is_param_set(quadratic_c))) then + write(fates_log(),*) 'Not all quadratic forest edge parameters found for bin ',b + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + ! Has no other parameters + call check_all_unset( (/ gaussian_amplitude, gaussian_sigma, gaussian_center, lognormal_amplitude, lognormal_sigma, lognormal_center /), b ) + + else + write(fates_log(),*) 'Unrecognized bin fit type for bin ',b + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do + + return + end subroutine EdgeForestCheckParams + + !----------------------------------------------------------------------- + subroutine EdgeForestRegisterParams(fates_params) + + use FatesParametersInterface, only : fates_parameters_type, dimension_name_scalar, dimension_shape_scalar + use FatesParametersInterface, only : dimension_name_edgeforest_bins + use FatesParametersInterface, only : dimension_shape_1d + + implicit none + + class(fates_parameters_type), intent(inout) :: fates_params + + character(len=param_string_length), parameter :: dim_names_edgeforest(1)= (/dimension_name_edgeforest_bins/) + + call fates_params%RegisterParameter(name=ED_name_edgeforest_gaussian_amplitude, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_edgeforest) + + call fates_params%RegisterParameter(name=ED_name_edgeforest_gaussian_sigma, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_edgeforest) + + call fates_params%RegisterParameter(name=ED_name_edgeforest_gaussian_center, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_edgeforest) + + call fates_params%RegisterParameter(name=ED_name_edgeforest_lognormal_amplitude, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_edgeforest) + + call fates_params%RegisterParameter(name=ED_name_edgeforest_lognormal_sigma, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_edgeforest) + + call fates_params%RegisterParameter(name=ED_name_edgeforest_lognormal_center, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_edgeforest) + + call fates_params%RegisterParameter(name=ED_name_edgeforest_quadratic_a, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_edgeforest) + + call fates_params%RegisterParameter(name=ED_name_edgeforest_quadratic_b, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_edgeforest) + + call fates_params%RegisterParameter(name=ED_name_edgeforest_quadratic_c, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_edgeforest) + + call fates_params%RegisterParameter(name=ED_name_edgeforest_bin_edges, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_edgeforest) + + end subroutine EdgeForestRegisterParams + + !----------------------------------------------------------------------- + subroutine EdgeForestReceiveParams(fates_params) + + use FatesParametersInterface, only : fates_parameters_type, dimension_name_scalar + + implicit none + + class(fates_parameters_type), intent(inout) :: fates_params + + call fates_params%RetrieveParameterAllocate(name=ED_name_edgeforest_gaussian_amplitude, & + data=ED_val_edgeforest_gaussian_amplitude) + + call fates_params%RetrieveParameterAllocate(name=ED_name_edgeforest_gaussian_sigma, & + data=ED_val_edgeforest_gaussian_sigma) + + call fates_params%RetrieveParameterAllocate(name=ED_name_edgeforest_gaussian_center, & + data=ED_val_edgeforest_gaussian_center) + + call fates_params%RetrieveParameterAllocate(name=ED_name_edgeforest_lognormal_amplitude, & + data=ED_val_edgeforest_lognormal_amplitude) + + call fates_params%RetrieveParameterAllocate(name=ED_name_edgeforest_lognormal_sigma, & + data=ED_val_edgeforest_lognormal_sigma) + + call fates_params%RetrieveParameterAllocate(name=ED_name_edgeforest_lognormal_center, & + data=ED_val_edgeforest_lognormal_center) + + call fates_params%RetrieveParameterAllocate(name=ED_name_edgeforest_quadratic_a, & + data=ED_val_edgeforest_quadratic_a) + + call fates_params%RetrieveParameterAllocate(name=ED_name_edgeforest_quadratic_b, & + data=ED_val_edgeforest_quadratic_b) + + call fates_params%RetrieveParameterAllocate(name=ED_name_edgeforest_quadratic_c, & + data=ED_val_edgeforest_quadratic_c) + + call fates_params%RetrieveParameterAllocate(name=ED_name_edgeforest_bin_edges, & + data=ED_val_edgeforest_bin_edges) + + end subroutine EdgeForestReceiveParams + + +end module FatesEdgeForestParamsMod diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index b9dc19bf49..fc18253705 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -50,6 +50,7 @@ module FatesHistoryInterfaceMod use FatesInterfaceTypesMod , only : hlm_use_ed_st3 use FatesInterfaceTypesMod , only : hlm_use_cohort_age_tracking use FatesInterfaceTypesMod , only : hlm_use_tree_damage + use FatesInterfaceTypesMod , only : hlm_use_edge_forest use FatesInterfaceTypesMod , only : nlevdamage use FatesInterfaceTypesMod , only : numpft use FatesInterfaceTypesMod , only : hlm_freq_day @@ -60,6 +61,7 @@ module FatesHistoryInterfaceMod use EDParamsMod , only : nlevleaf use EDParamsMod , only : ED_val_history_height_bin_edges use EDParamsMod , only : ED_val_history_ageclass_bin_edges + use EDParamsMod , only : ED_val_history_height_bin_edges use FatesInterfaceTypesMod , only : nlevsclass, nlevage use FatesInterfaceTypesMod , only : nlevheight use FatesInterfaceTypesMod , only : bc_in_type @@ -68,6 +70,7 @@ module FatesHistoryInterfaceMod use FatesInterfaceTypesMod , only : nlevcoage use FatesInterfaceTypesMod , only : hlm_use_nocomp use FatesInterfaceTypesMod , only : hlm_use_fixed_biogeog + use FatesInterfaceTypesMod , only : nlevedgeforest use FatesRadiationMemMod , only : ivis,inir use FatesInterfaceTypesMod , only : hlm_hist_level_hifrq,hlm_hist_level_dynam use FatesIOVariableKindMod, only : site_r8, site_soil_r8, site_size_pft_r8 @@ -81,6 +84,7 @@ module FatesHistoryInterfaceMod use FatesIOVariableKindMod, only : site_elcwd_r8, site_elage_r8, site_clscpf_r8 use FatesIOVariableKindMod, only : site_cdpf_r8, site_cdsc_r8, site_cdam_r8 use FatesIOVariableKindMod, only : site_landuse_r8, site_lulu_r8, site_lupft_r8 + use FatesIOVariableKindMod, only : site_edgebin_r8 use FatesConstantsMod , only : n_landuse_cats use FatesAllometryMod , only : CrownDepth use FatesAllometryMod , only : bstore_allom, bsap_allom @@ -293,6 +297,19 @@ module FatesHistoryInterfaceMod integer :: ih_trimming_si integer :: ih_fracarea_plant_si integer :: ih_fracarea_trees_si + integer :: ih_is_forest_si + integer :: ih_is_forest_si_age + integer :: ih_is_forest_pct10_si + integer :: ih_is_forest_pct25_si + integer :: ih_is_forest_pct50_si + integer :: ih_is_forest_pct75_si + integer :: ih_is_forest_pct90_si + integer :: ih_is_forest_pct10_0grass_si + integer :: ih_is_forest_pct25_0grass_si + integer :: ih_is_forest_pct50_0grass_si + integer :: ih_is_forest_pct75_0grass_si + integer :: ih_is_forest_pct90_0grass_si + integer :: ih_forest_edge_bin_area_si integer :: ih_litter_in_elem integer :: ih_litter_out_elem integer :: ih_seed_bank_elem @@ -664,6 +681,8 @@ module FatesHistoryInterfaceMod ! indices to (site x patch-age) variables integer :: ih_fracarea_si_age + integer :: ih_fracarea_plant_si_age + integer :: ih_fracarea_trees_si_age integer :: ih_lai_si_age integer :: ih_canopy_fracarea_si_age integer :: ih_gpp_si_age @@ -831,6 +850,7 @@ module FatesHistoryInterfaceMod !! THESE WERE EXPLICITLY PRIVATE WHEN TYPE WAS PUBLIC integer, private :: column_index_, levsoil_index_, levscpf_index_ integer, private :: levscls_index_, levpft_index_, levage_index_ + integer, private :: levedgeforest_index_ integer, private :: levfuel_index_, levcwdsc_index_, levscag_index_ integer, private :: levcan_index_, levcnlf_index_, levcnlfpft_index_ integer, private :: levcdpf_index_, levcdsc_index_, levcdam_index_ @@ -872,6 +892,7 @@ module FatesHistoryInterfaceMod procedure :: levcacls_index procedure :: levpft_index procedure :: levage_index + procedure :: levedgeforest_index procedure :: levfuel_index procedure :: levcwdsc_index procedure :: levcan_index @@ -908,6 +929,7 @@ module FatesHistoryInterfaceMod procedure, private :: set_levscls_index procedure, private :: set_levpft_index procedure, private :: set_levage_index + procedure, private :: set_levedgeforest_index procedure, private :: set_levfuel_index procedure, private :: set_levcwdsc_index procedure, private :: set_levcan_index @@ -953,6 +975,7 @@ subroutine Init(this, num_threads, fates_bounds) use FatesIODimensionsMod, only : column, levsoil, levscpf use FatesIODimensionsMod, only : levscls, levpft, levage + use FatesIODimensionsMod, only : levedgeforest use FatesIODimensionsMod, only : levcacls, levcapf use FatesIODimensionsMod, only : levfuel, levcwdsc, levscag use FatesIODimensionsMod, only : levscagpft, levagepft @@ -1012,6 +1035,11 @@ subroutine Init(this, num_threads, fates_bounds) call this%dim_bounds(dim_count)%Init(levage, num_threads, & fates_bounds%age_class_begin, fates_bounds%age_class_end) + dim_count = dim_count + 1 + call this%set_levedgeforest_index(dim_count) + call this%dim_bounds(dim_count)%Init(levedgeforest, num_threads, & + fates_bounds%edgeforest_class_begin, fates_bounds%edgeforest_class_end) + dim_count = dim_count + 1 call this%set_levfuel_index(dim_count) call this%dim_bounds(dim_count)%Init(levfuel, num_threads, & @@ -1165,6 +1193,10 @@ subroutine SetThreadBoundsEach(this, thread_index, thread_bounds) call this%dim_bounds(index)%SetThreadBounds(thread_index, & thread_bounds%age_class_begin, thread_bounds%age_class_end) + index = this%levedgeforest_index() + call this%dim_bounds(index)%SetThreadBounds(thread_index, & + thread_bounds%edgeforest_class_begin, thread_bounds%edgeforest_class_end) + index = this%levfuel_index() call this%dim_bounds(index)%SetThreadBounds(thread_index, & thread_bounds%fuel_begin, thread_bounds%fuel_end) @@ -1255,8 +1287,6 @@ end subroutine SetThreadBoundsEach ! =================================================================================== subroutine assemble_history_output_types(this) - - implicit none class(fates_history_interface_type), intent(inout) :: this @@ -1286,6 +1316,9 @@ subroutine assemble_history_output_types(this) call this%set_dim_indices(site_age_r8, 1, this%column_index()) call this%set_dim_indices(site_age_r8, 2, this%levage_index()) + call this%set_dim_indices(site_edgebin_r8, 1, this%column_index()) + call this%set_dim_indices(site_edgebin_r8, 2, this%levedgeforest_index()) + call this%set_dim_indices(site_fuel_r8, 1, this%column_index()) call this%set_dim_indices(site_fuel_r8, 2, this%levfuel_index()) @@ -1503,6 +1536,20 @@ integer function levage_index(this) levage_index = this%levage_index_ end function levage_index + ! ======================================================================= + subroutine set_levedgeforest_index(this, index) + implicit none + class(fates_history_interface_type), intent(inout) :: this + integer, intent(in) :: index + this%levedgeforest_index_ = index + end subroutine set_levedgeforest_index + + integer function levedgeforest_index(this) + implicit none + class(fates_history_interface_type), intent(in) :: this + levedgeforest_index = this%levedgeforest_index_ + end function levedgeforest_index + ! ======================================================================= subroutine set_levfuel_index(this, index) implicit none @@ -1975,16 +2022,6 @@ subroutine init_dim_kinds_maps(this) ! number of entries listed here. ! ! ---------------------------------------------------------------------------------- - use FatesIOVariableKindMod, only : site_r8, site_soil_r8, site_size_pft_r8 - use FatesIOVariableKindMod, only : site_size_r8, site_pft_r8, site_age_r8 - use FatesIOVariableKindMod, only : site_coage_r8, site_coage_pft_r8 - use FatesIOVariableKindMod, only : site_fuel_r8, site_cwdsc_r8, site_scag_r8 - use FatesIOVariableKindMod, only : site_scagpft_r8, site_agepft_r8 - use FatesIOVariableKindMod, only : site_can_r8, site_cnlf_r8, site_cnlfpft_r8 - use FatesIOVariableKindMod, only : site_height_r8, site_agefuel_r8 - use FatesIOVariableKindMod, only : site_elem_r8, site_elpft_r8 - use FatesIOVariableKindMod, only : site_elcwd_r8, site_elage_r8, site_clscpf_r8 - use FatesIOVariableKindMod, only : site_cdpf_r8, site_cdsc_r8, site_cdam_r8 implicit none @@ -2026,6 +2063,10 @@ subroutine init_dim_kinds_maps(this) index = index + 1 call this%dim_kinds(index)%Init(site_age_r8, 2) + ! site x forest-edge-bin class + index = index + 1 + call this%dim_kinds(index)%Init(site_edgebin_r8, 2) + ! site x fuel size class index = index + 1 call this%dim_kinds(index)%Init(site_fuel_r8, 2) @@ -2387,6 +2428,8 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) ! updated here, but not FATES_VEGC_PF. ! --------------------------------------------------------------------------------- + use FatesEcotypesMod, only : IsPatchForest + use FatesUtilsMod, only : logical_to_real ! Arguments class(fates_history_interface_type) :: this @@ -2436,7 +2479,18 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) hio_zstar_si => this%hvars(ih_zstar_si)%r81d, & hio_trimming_si => this%hvars(ih_trimming_si)%r81d, & hio_fracarea_plant_si => this%hvars(ih_fracarea_plant_si)%r81d, & - hio_fracarea_trees_si => this%hvars(ih_fracarea_trees_si)%r81d, & + hio_fracarea_trees_si => this%hvars(ih_fracarea_trees_si)%r81d, & + hio_is_forest_si => this%hvars(ih_is_forest_si)%r81d, & + hio_is_forest_pct10_si => this%hvars(ih_is_forest_pct10_si)%r81d, & + hio_is_forest_pct25_si => this%hvars(ih_is_forest_pct25_si)%r81d, & + hio_is_forest_pct50_si => this%hvars(ih_is_forest_pct50_si)%r81d, & + hio_is_forest_pct75_si => this%hvars(ih_is_forest_pct75_si)%r81d, & + hio_is_forest_pct90_si => this%hvars(ih_is_forest_pct90_si)%r81d, & + hio_is_forest_pct10_0grass_si => this%hvars(ih_is_forest_pct10_0grass_si)%r81d, & + hio_is_forest_pct25_0grass_si => this%hvars(ih_is_forest_pct25_0grass_si)%r81d, & + hio_is_forest_pct50_0grass_si => this%hvars(ih_is_forest_pct50_0grass_si)%r81d, & + hio_is_forest_pct75_0grass_si => this%hvars(ih_is_forest_pct75_0grass_si)%r81d, & + hio_is_forest_pct90_0grass_si => this%hvars(ih_is_forest_pct90_0grass_si)%r81d, & 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, & @@ -2719,6 +2773,31 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) hio_fracarea_plant_si(io_si) = hio_fracarea_plant_si(io_si) + min(cpatch%total_canopy_area,cpatch%area) * AREA_INV hio_fracarea_trees_si(io_si) = hio_fracarea_trees_si(io_si) + min(cpatch%total_tree_area,cpatch%area) * AREA_INV + ! whether patch is forest according to FATES parameter file threshold + hio_is_forest_si(io_si) = hio_is_forest_si(io_si) + & + merge(1._r8, 0._r8, cpatch%is_forest) * cpatch%area * AREA_INV + ! according to experimental definitions + hio_is_forest_pct10_si(io_si) = hio_is_forest_pct10_si(io_si) + cpatch%area * AREA_INV * & + logical_to_real(IsPatchForest(cpatch, 0.1_r8)) + hio_is_forest_pct25_si(io_si) = hio_is_forest_pct25_si(io_si) + cpatch%area * AREA_INV * & + logical_to_real(IsPatchForest(cpatch, 0.25_r8)) + hio_is_forest_pct50_si(io_si) = hio_is_forest_pct50_si(io_si) + cpatch%area * AREA_INV * & + logical_to_real(IsPatchForest(cpatch, 0.5_r8)) + hio_is_forest_pct75_si(io_si) = hio_is_forest_pct75_si(io_si) + cpatch%area * AREA_INV * & + logical_to_real(IsPatchForest(cpatch, 0.75_r8)) + hio_is_forest_pct90_si(io_si) = hio_is_forest_pct90_si(io_si) + cpatch%area * AREA_INV * & + logical_to_real(IsPatchForest(cpatch, 0.9_r8)) + hio_is_forest_pct10_0grass_si(io_si) = hio_is_forest_pct10_0grass_si(io_si) + cpatch%area * AREA_INV * & + logical_to_real(IsPatchForest(cpatch, 0.1_r8, grass_biomass_threshold=0._r8)) + hio_is_forest_pct25_0grass_si(io_si) = hio_is_forest_pct25_0grass_si(io_si) + cpatch%area * AREA_INV * & + logical_to_real(IsPatchForest(cpatch, 0.25_r8, grass_biomass_threshold=0._r8)) + hio_is_forest_pct50_0grass_si(io_si) = hio_is_forest_pct50_0grass_si(io_si) + cpatch%area * AREA_INV * & + logical_to_real(IsPatchForest(cpatch, 0.5_r8, grass_biomass_threshold=0._r8)) + hio_is_forest_pct75_0grass_si(io_si) = hio_is_forest_pct75_0grass_si(io_si) + cpatch%area * AREA_INV * & + logical_to_real(IsPatchForest(cpatch, 0.75_r8, grass_biomass_threshold=0._r8)) + hio_is_forest_pct90_0grass_si(io_si) = hio_is_forest_pct90_0grass_si(io_si) + cpatch%area * AREA_INV * & + logical_to_real(IsPatchForest(cpatch, 0.9_r8, grass_biomass_threshold=0._r8)) + ! Patch specific variables that are already calculated ! These things are all duplicated. Should they all be converted to LL or array structures RF? ! define scalar to counteract the patch albedo scaling logic for conserved quantities @@ -3084,7 +3163,7 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) type(litter_type), pointer :: litt_c ! Pointer to the carbon12 litter pool type(litter_type), pointer :: litt ! Generic pointer to any litter pool integer :: s ! site counter - integer :: ipa2 ! patch index matching host model array space + integer :: b ! edge bin counter integer :: io_si ! site's index in the history output array space integer :: el ! element index integer :: ft ! pft index @@ -3319,6 +3398,7 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) hio_cwd_ag_out_si_cwdsc => this%hvars(ih_cwd_ag_out_si_cwdsc)%r82d, & hio_cwd_bg_out_si_cwdsc => this%hvars(ih_cwd_bg_out_si_cwdsc)%r82d, & hio_crownarea_si_cnlf => this%hvars(ih_crownarea_si_cnlf)%r82d, & + hio_forest_edge_bin_area_si => this%hvars(ih_forest_edge_bin_area_si)%r82d, & hio_crownarea_cl => this%hvars(ih_crownarea_cl)%r82d) ! Break up associates for NAG compilers @@ -3431,6 +3511,13 @@ subroutine update_history_dyn_subsite(this,nc,nsites,sites,bc_in) hio_fracarea_si(io_si) = hio_fracarea_si(io_si) & + cpatch%area * AREA_INV + ! area of forest in each edge bin + if (hlm_use_edge_forest == itrue .and. cpatch%is_forest) then + binloop: do b = 1, nlevedgeforest + hio_forest_edge_bin_area_si(io_si,b) = hio_forest_edge_bin_area_si(io_si,b) + & + cpatch%area_in_edgeforest_bins(b) + end do binloop + end if ! ignore land use info on nocomp bareground (where landuse label = 0) if (cpatch%land_use_label .gt. nocomp_bareground_land) then @@ -4782,6 +4869,7 @@ subroutine update_history_dyn_subsite_ageclass(this,nc,nsites,sites) associate( & hio_lai_si_age => this%hvars(ih_lai_si_age)%r82d, & + hio_is_forest_si_age => this%hvars(ih_is_forest_si_age)%r82d, & hio_ncl_si_age => this%hvars(ih_ncl_si_age)%r82d, & hio_scorch_height_si_agepft => this%hvars(ih_scorch_height_si_agepft)%r82d, & hio_zstar_si_age => this%hvars(ih_zstar_si_age)%r82d, & @@ -4808,6 +4896,8 @@ subroutine update_history_dyn_subsite_ageclass(this,nc,nsites,sites) hio_nplant_canopy_si_scag => this%hvars(ih_nplant_canopy_si_scag)%r82d, & hio_nplant_understory_si_scag => this%hvars(ih_nplant_understory_si_scag)%r82d, & hio_fracarea_si_age => this%hvars(ih_fracarea_si_age)%r82d, & + hio_fracarea_plant_si_age => this%hvars(ih_fracarea_plant_si_age)%r82d, & + hio_fracarea_trees_si_age => this%hvars(ih_fracarea_trees_si_age)%r82d, & hio_agesince_anthrodist_si_age => this%hvars(ih_agesince_anthrodist_si_age)%r82d, & hio_primarylands_fracarea_si_age => this%hvars(ih_primarylands_fracarea_si_age)%r82d, & hio_secondarylands_fracarea_si_age => this%hvars(ih_secondarylands_fracarea_si_age)%r82d, & @@ -4866,6 +4956,10 @@ subroutine update_history_dyn_subsite_ageclass(this,nc,nsites,sites) + cpatch%zstar * patch_area_div_site_area end if + ! whether patch is forest according to FATES parameter file threshold + hio_is_forest_si_age(io_si,cpatch%age_class) = hio_is_forest_si_age(io_si,cpatch%age_class) + & + merge(1._r8, 0._r8, cpatch%is_forest) * patch_area_div_site_area + ! some diagnostics on secondary forest area and its age distribution if ( cpatch%land_use_label .eq. secondaryland ) then @@ -4900,16 +4994,24 @@ subroutine update_history_dyn_subsite_ageclass(this,nc,nsites,sites) hio_fire_intensity_si_age(io_si, cpatch%age_class) = hio_fire_intensity_si_age(io_si,cpatch%age_class) + & cpatch%FI * J_per_kJ & ! [kJ/m/s] -> [J/m/s] * cpatch%frac_burnt * patch_area_div_site_area - ! hio_fire_rate_of_spread_front_si_age(io_si, cpatch%age_class) = hio_fire_rate_of_spread_si_age(io_si, cpatch%age_class) + & - ! cpatch%ros_front * cpatch*frac_burnt * patch_area_div_site_area - hio_rx_intensity_si_age(io_si, cpatch%age_class) = hio_rx_intensity_si_age(io_si, cpatch%age_class) + & cpatch%rx_FI * J_per_kJ & ! [kJ/m/s] -> [J/m/s] * cpatch%rx_frac_burnt * patch_area_div_site_area - hio_nonrx_intensity_si_age(io_si, cpatch%age_class) = hio_nonrx_intensity_si_age(io_si, cpatch%age_class) + & cpatch%nonrx_FI * J_per_kJ & ! [kJ/m/s] -> [J/m/s] * cpatch%nonrx_frac_burnt * patch_area_div_site_area + ! hio_fire_rate_of_spread_front_si_age(io_si, cpatch%age_class) = hio_fire_rate_of_spread_si_age(io_si, cpatch%age_class) + & + ! cpatch%ros_front * cpatch*frac_burnt * patch_area_div_site_area + + ! Weighted by site-wide plant or tree canopy area + hio_fracarea_plant_si_age(io_si,cpatch%age_class) = & + hio_fracarea_plant_si_age(io_si,cpatch%age_class) + & + min(cpatch%total_canopy_area,cpatch%area) * & + AREA_INV + hio_fracarea_trees_si_age(io_si,cpatch%age_class) = & + hio_fracarea_trees_si_age(io_si,cpatch%age_class) + & + min(cpatch%total_tree_area,cpatch%area) * & + AREA_INV ! Weighted by cohort canopy area relative to site area ccohort => cpatch%shortest @@ -5428,7 +5530,7 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,dt_tstep) hio_parsun_si_can => this%hvars(ih_parsun_si_can)%r82d, & hio_parsha_si_can => this%hvars(ih_parsha_si_can)%r82d, & hio_laisun_si_can => this%hvars(ih_laisun_si_can)%r82d, & - hio_laisha_si_can => this%hvars(ih_laisha_si_can)%r82d ) + hio_laisha_si_can => this%hvars(ih_laisha_si_can)%r82d) ! THIS CAN BE REMOVED WHEN BOTH CTSM AND E3SM CALL FLUSH_ALL_HVARS @@ -6296,19 +6398,8 @@ subroutine define_history_vars(this, initialize_variables) ! a real. The applied flush value will use the NINT() intrinsic function ! --------------------------------------------------------------------------------- - use FatesIOVariableKindMod, only : site_r8, site_soil_r8, site_size_pft_r8 - use FatesIOVariableKindMod, only : site_size_r8, site_pft_r8, site_age_r8 - use FatesIOVariableKindMod, only : site_coage_pft_r8, site_coage_r8 - use FatesIOVariableKindMod, only : site_height_r8, site_agefuel_r8 use FatesInterfaceTypesMod, only : hlm_use_planthydro - use FatesIOVariableKindMod, only : site_fuel_r8, site_cwdsc_r8, site_scag_r8 - use FatesIOVariableKindMod, only : site_can_r8, site_cnlf_r8, site_cnlfpft_r8 - use FatesIOVariableKindMod, only : site_cdsc_r8, site_cdpf_r8, site_cdam_r8 - use FatesIOVariableKindMod, only : site_scagpft_r8, site_agepft_r8 - use FatesIOVariableKindMod, only : site_elem_r8, site_elpft_r8, site_clscpf_r8 - use FatesIOVariableKindMod, only : site_elcwd_r8, site_elage_r8 - implicit none @@ -6327,6 +6418,7 @@ subroutine define_history_vars(this, initialize_variables) ! patch age (site_age_r8) : AP ! canopy layer (site_can_r8) : CL ! coarse woody debris size (site_cwdsc_r8) : DC + ! forest edge bin (site_edgebin_r8): EB ! element (site_elem_r8) : EL ! leaf layer : LL ! fuel class (site_fuel_r8) : FC @@ -6378,12 +6470,101 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_r8, hlms='CLM:ALM', upfreq=group_dyna_simple, ivar=ivar, & initialize=initialize_variables, index=ih_fracarea_plant_si) + call this%set_history_var(vname='FATES_AREA_PLANTS_AP', units='m2 m-2', & + long='area occupied by all plants per m2 land area (by patch age)', use_default='active', & + avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', upfreq=group_dyna_complx, ivar=ivar, & + initialize=initialize_variables, index=ih_fracarea_plant_si_age) + call this%set_history_var(vname='FATES_AREA_TREES', units='m2 m-2', & long='area occupied by woody plants per m2 land area', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', & upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & index=ih_fracarea_trees_si) + call this%set_history_var(vname='FATES_AREA_TREES_AP', units='m2 m-2', & + long='area occupied by woody plants per m2 land area (by patch age)', use_default='active', & + avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', & + upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & + index=ih_fracarea_trees_si_age) + + call this%set_history_var(vname='FATES_IS_FOREST', units='', & + long='whether patch is forest', use_default='inactive', & + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_si) + + call this%set_history_var(vname='FATES_IS_FOREST_AP', units='', & + long='whether patch is forest (by patch age)', use_default='inactive',& + avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', & + upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_si_age) + + call this%set_history_var(vname='FATES_IS_FOREST_PCT10', units='', & + long='whether patch is forest (10% threshold)', use_default='inactive',& + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_pct10_si) + + call this%set_history_var(vname='FATES_IS_FOREST_PCT25', units='', & + long='whether patch is forest (25% threshold)', use_default='inactive',& + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_pct25_si) + + call this%set_history_var(vname='FATES_IS_FOREST_PCT50', units='', & + long='whether patch is forest (50% threshold)', use_default='inactive',& + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_pct50_si) + + call this%set_history_var(vname='FATES_IS_FOREST_PCT75', units='', & + long='whether patch is forest (75% threshold)', use_default='inactive',& + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_pct75_si) + + call this%set_history_var(vname='FATES_IS_FOREST_PCT90', units='', & + long='whether patch is forest (90% threshold)', use_default='inactive',& + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_pct90_si) + + call this%set_history_var(vname='FATES_IS_FOREST_PCT10_0GRASS', units='', & + long='whether patch is forest (10% threshold, no grass)', use_default='inactive',& + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_pct10_0grass_si) + + call this%set_history_var(vname='FATES_IS_FOREST_PCT25_0GRASS', units='', & + long='whether patch is forest (25% threshold, no grass)', use_default='inactive',& + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_pct25_0grass_si) + + call this%set_history_var(vname='FATES_IS_FOREST_PCT50_0GRASS', units='', & + long='whether patch is forest (50% threshold, no grass)', use_default='inactive',& + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_pct50_0grass_si) + + call this%set_history_var(vname='FATES_IS_FOREST_PCT75_0GRASS', units='', & + long='whether patch is forest (75% threshold, no grass)', use_default='inactive',& + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_pct75_0grass_si) + + call this%set_history_var(vname='FATES_IS_FOREST_PCT90_0GRASS', units='', & + long='whether patch is forest (90% threshold, no grass)', use_default='inactive',& + avgflag='A', vtype=site_r8, hlms='CLM:ALM', & + upfreq=group_dyna_simple, ivar=ivar, initialize=initialize_variables, & + index=ih_is_forest_pct90_0grass_si) + + call this%set_history_var(vname='FATES_FOREST_AREA_EB', units='m2', & + long='area of forest in each edge bin', use_default='inactive', & + avgflag='A', vtype=site_edgebin_r8, hlms='CLM:ALM', & + upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables, & + index=ih_forest_edge_bin_area_si) + call this%set_history_var(vname='FATES_FRACTION', units='m2 m-2', & long='total gridcell fraction which FATES is running over', use_default='active', & avgflag='A', vtype=site_r8, hlms='CLM:ALM', & diff --git a/main/FatesHistoryVariableType.F90 b/main/FatesHistoryVariableType.F90 index 3b0bcc1b29..fc1772665f 100644 --- a/main/FatesHistoryVariableType.F90 +++ b/main/FatesHistoryVariableType.F90 @@ -8,6 +8,7 @@ module FatesHistoryVariableType use FatesIOVariableKindMod, only : fates_io_variable_kind_type use FatesIOVariableKindMod, only : site_r8, site_soil_r8, site_size_pft_r8 use FatesIOVariableKindMod, only : site_size_r8, site_pft_r8, site_age_r8 + use FatesIOVariableKindMod, only : site_edgebin_r8 use FatesIOVariableKindMod, only : site_coage_r8, site_coage_pft_r8 use FatesIOVariableKindMod, only : site_height_r8 use FatesIOVariableKindMod, only : site_fuel_r8, site_cwdsc_r8, site_scag_r8 @@ -153,6 +154,10 @@ subroutine Init(this, vname, units, long, use_default, & allocate(this%r82d(lb1:ub1, lb2:ub2)) this%r82d(:,:) = flushval + case(site_edgebin_r8) + allocate(this%r82d(lb1:ub1, lb2:ub2)) + this%r82d(:,:) = flushval + case(site_height_r8) allocate(this%r82d(lb1:ub1, lb2:ub2)) this%r82d(:,:) = flushval @@ -321,6 +326,8 @@ subroutine HFlush(this, thread, dim_bounds, dim_kinds) this%r82d(lb1:ub1, lb2:ub2) = this%flushval case(site_age_r8) this%r82d(lb1:ub1, lb2:ub2) = this%flushval + case(site_edgebin_r8) + this%r82d(lb1:ub1, lb2:ub2) = this%flushval case(site_height_r8) this%r82d(lb1:ub1, lb2:ub2) = this%flushval case(site_fuel_r8) diff --git a/main/FatesIODimensionsMod.F90 b/main/FatesIODimensionsMod.F90 index ed487d7eed..b618af1902 100644 --- a/main/FatesIODimensionsMod.F90 +++ b/main/FatesIODimensionsMod.F90 @@ -20,6 +20,7 @@ module FatesIODimensionsMod character(*), parameter, public :: levscls = 'fates_levscls' ! matches histFileMod character(*), parameter, public :: levpft = 'fates_levpft' ! matches histFileMod character(*), parameter, public :: levage = 'fates_levage' ! matches histFileMod + character(*), parameter, public :: levedgeforest = 'fates_levedge' ! matches histFileMod character(*), parameter, public :: levheight = 'fates_levheight' ! matches histFileMod character(*), parameter, public :: levfuel = 'fates_levfuel' ! matches histFileMod character(*), parameter, public :: levcwdsc = 'fates_levcwdsc' ! matches histFileMod @@ -64,6 +65,9 @@ module FatesIODimensionsMod ! levage = This is a structure that records the boundaries for the ! number of patch-age-class dimension + ! levedgeforest = This is a structure that records the boundaries for the + ! number of forest-edge-bin dimension + ! levheight = This is a structure that records the boundaries for the ! number of height dimension @@ -153,6 +157,8 @@ module FatesIODimensionsMod integer :: pft_class_end integer :: age_class_begin integer :: age_class_end + integer :: edgeforest_class_begin + integer :: edgeforest_class_end integer :: height_begin integer :: height_end integer :: fuel_begin diff --git a/main/FatesIOVariableKindMod.F90 b/main/FatesIOVariableKindMod.F90 index 75ea7dbe57..efd49e84f6 100644 --- a/main/FatesIOVariableKindMod.F90 +++ b/main/FatesIOVariableKindMod.F90 @@ -27,6 +27,7 @@ module FatesIOVariableKindMod character(*), parameter, public :: cohort_int = 'CO_INT' character(*), parameter, public :: site_pft_r8 = 'SI_PFT_R8' character(*), parameter, public :: site_age_r8 = 'SI_AGE_R8' + character(*), parameter, public :: site_edgebin_r8 = 'SI_EDGEBIN_R8' character(*), parameter, public :: site_height_r8 = 'SI_HEIGHT_R8' character(*), parameter, public :: site_fuel_r8 = 'SI_FUEL_R8' character(*), parameter, public :: site_cwdsc_r8 = 'SI_CWDSC_R8' diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index c5ce7f56aa..a791a0a0b1 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -43,6 +43,7 @@ module FatesInterfaceMod use FatesConstantsMod , only : n_crop_lu_types use FatesConstantsMod , only : n_term_mort_types use FatesConstantsMod , only : nocomp_bareground + use FatesInterfaceTypesMod , only : nlevedgeforest, hlm_use_tree_damage use FatesGlobals , only : fates_global_verbose use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun @@ -53,6 +54,7 @@ module FatesInterfaceMod use EDPftvarcon , only : FatesCheckParams use EDPftvarcon , only : EDPftvarcon_inst use SFParamsMod , only : SpitFireCheckParams + use FatesEdgeForestParamsMod , only : EdgeForestCheckParams use EDParamsMod , only : FatesReportParams use EDParamsMod , only : bgc_soil_salinity use FatesPlantHydraulicsMod , only : InitHydroGlobals @@ -64,12 +66,14 @@ module FatesInterfaceMod use EDParamsMod , only : sdlng_mdd_timescale use EDParamsMod , only : ED_val_history_sizeclass_bin_edges use EDParamsMod , only : ED_val_history_ageclass_bin_edges + use FatesEdgeForestParamsMod , only : ED_val_edgeforest_bin_edges use EDParamsMod , only : ED_val_history_height_bin_edges use EDParamsMod , only : ED_val_history_coageclass_bin_edges use FatesParametersInterface , only : fates_param_reader_type use FatesParametersInterface , only : fates_parameters_type use EDParamsMod , only : FatesRegisterParams, FatesReceiveParams use SFParamsMod , only : SpitFireRegisterParams, SpitFireReceiveParams + use FatesEdgeForestParamsMod , only : EdgeForestRegisterParams, EdgeForestReceiveParams use PRTInitParamsFATESMod , only : PRTRegisterParams, PRTReceiveParams use FatesLeafBiophysParamsMod , only : LeafBiophysRegisterParams, LeafBiophysReceiveParams,LeafBiophysReportParams use FatesSynchronizedParamsMod, only : FatesSynchronizedParamsInst @@ -972,6 +976,7 @@ subroutine SetFatesGlobalElements2(use_fates) ! Identify number of size and age class bins for history output ! assume these arrays are 1-indexed nlevage = size(ED_val_history_ageclass_bin_edges,dim=1) + nlevedgeforest = size(ED_val_edgeforest_bin_edges,dim=1) nlevheight = size(ED_val_history_height_bin_edges,dim=1) nlevcoage = size(ED_val_history_coageclass_bin_edges,dim=1) nlevdamage = size(ED_val_history_damage_bin_edges, dim=1) @@ -986,6 +991,10 @@ subroutine SetFatesGlobalElements2(use_fates) write(fates_log(), *) 'age class bins specified in parameter file must start at zero' call endrun(msg=errMsg(sourcefile, __LINE__)) endif + if ( ED_val_edgeforest_bin_edges(1) .ne. 0._r8 ) then + write(fates_log(), *) 'edge forest class bins specified in parameter file must start at zero' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif if ( ED_val_history_height_bin_edges(1) .ne. 0._r8 ) then write(fates_log(), *) 'height class bins specified in parameter file must start at zero' call endrun(msg=errMsg(sourcefile, __LINE__)) @@ -1002,6 +1011,12 @@ subroutine SetFatesGlobalElements2(use_fates) call endrun(msg=errMsg(sourcefile, __LINE__)) end if end do + do i = 2,nlevedgeforest + if ( (ED_val_edgeforest_bin_edges(i) - ED_val_edgeforest_bin_edges(i-1)) .le. 0._r8) then + write(fates_log(), *) 'edge forest class bins specified in parameter file must be monotonically increasing' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do do i = 2,nlevheight if ( (ED_val_history_height_bin_edges(i) - ED_val_history_height_bin_edges(i-1)) .le. 0._r8) then write(fates_log(), *) 'height class bins specified in parameter file must be monotonically increasing' @@ -1014,7 +1029,7 @@ subroutine SetFatesGlobalElements2(use_fates) call endrun(msg=errMsg(sourcefile, __LINE__)) end if end do - + ! Set the fates dispersal kernel mode if there are any seed dispersal parameters set. ! The validation of the parameter values is check in FatesCheckParams prior to this check. ! This is currently hard coded, but could be added as a fates parameter file option, @@ -1160,6 +1175,7 @@ subroutine fates_history_maps use EDParamsMod, only : nlevleaf use EDParamsMod, only : ED_val_history_sizeclass_bin_edges use EDParamsMod, only : ED_val_history_ageclass_bin_edges + use FatesEdgeForestParamsMod, only : ED_val_edgeforest_bin_edges use EDParamsMod, only : ED_val_history_height_bin_edges use EDParamsMod, only : ED_val_history_coageclass_bin_edges @@ -1193,6 +1209,7 @@ subroutine fates_history_maps allocate( fates_hdim_levfuel(1:num_fuel_classes )) allocate( fates_hdim_levcwdsc(1:NCWD )) allocate( fates_hdim_levage(1:nlevage )) + allocate( fates_hdim_levedge(1:nlevedgeforest )) allocate( fates_hdim_levheight(1:nlevheight )) allocate( fates_hdim_levcoage(1:nlevcoage )) allocate( fates_hdim_pfmap_levcapf(1:nlevcoage*numpft)) @@ -1233,6 +1250,7 @@ subroutine fates_history_maps ! Fill the IO array of plant size classes fates_hdim_levsclass(:) = ED_val_history_sizeclass_bin_edges(:) fates_hdim_levage(:) = ED_val_history_ageclass_bin_edges(:) + fates_hdim_levedge(:) = ED_val_edgeforest_bin_edges(:) fates_hdim_levheight(:) = ED_val_history_height_bin_edges(:) fates_hdim_levcoage(:) = ED_val_history_coageclass_bin_edges(:) fates_hdim_levleaf(:) = dlower_vai(:) @@ -1491,6 +1509,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_decomp = 'unset' hlm_nitrogen_spec = unset_int hlm_use_tree_damage = unset_int + hlm_use_edge_forest = unset_int hlm_phosphorus_spec = unset_int hlm_use_ch4 = unset_int hlm_use_vertsoilc = unset_int @@ -1688,15 +1707,16 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) if(hlm_use_tree_damage .eq. unset_int) then write(fates_log(),*) 'FATES dimension/parameter unset: hlm_use_tree_damage, exiting' call endrun(msg=errMsg(sourcefile, __LINE__)) - else - if((hlm_use_tree_damage .eq. itrue) .and. & + else if ((hlm_use_tree_damage .eq. itrue) .and. & (hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp))then - write(fates_log(),*) 'FATES tree damage (use_fates_tree_damage = .true.) is not' - write(fates_log(),*) '(yet) compatible with CNP allocation (fates_parteh_mode = 2)' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) 'FATES tree damage (use_fates_tree_damage = .true.) is not' + write(fates_log(),*) '(yet) compatible with CNP allocation (fates_parteh_mode = 2)' + call endrun(msg=errMsg(sourcefile, __LINE__)) end if - + if(hlm_use_edge_forest .eq. unset_int) then + write(fates_log(),*) 'FATES dimension/parameter unset: hlm_use_edge_forest, exiting' + call endrun(msg=errMsg(sourcefile, __LINE__)) end if if(hlm_nitrogen_spec .eq. unset_int) then @@ -1923,6 +1943,12 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) if (fates_global_verbose()) then write(fates_log(),*) 'Transfering hlm_use_tree_damage = ',ival,' to FATES' end if + + case('use_edge_forest') + hlm_use_edge_forest = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering hlm_use_edge_forest = ',ival,' to FATES' + end if case('nitrogen_spec') hlm_nitrogen_spec = ival @@ -2250,6 +2276,7 @@ subroutine FatesReportParameters(masterproc) call FatesCheckParams(masterproc) ! Check general fates parameters call PRTCheckParams(masterproc) ! Check PARTEH parameters call SpitFireCheckParams(masterproc) + call EdgeForestCheckParams(masterproc) call TransferRadParams() @@ -2686,6 +2713,7 @@ subroutine FatesReadParameters(param_reader) call fates_params%Init() ! fates_params class, in FatesParameterInterfaceMod call FatesRegisterParams(fates_params) !EDParamsMod, only operates on fates_params class call SpitFireRegisterParams(fates_params) !SpitFire Mod, only operates of fates_params class + call EdgeForestRegisterParams(fates_params) !EdgeForest Mod, only operates on fates_params class call PRTRegisterParams(fates_params) ! PRT mod, only operates on fates_params class call LeafBiophysRegisterParams(fates_params) call FatesSynchronizedParamsInst%RegisterParams(fates_params) !Synchronized params class in Synchronized params mod, only operates on fates_params class @@ -2694,6 +2722,7 @@ subroutine FatesReadParameters(param_reader) call FatesReceiveParams(fates_params) call SpitFireReceiveParams(fates_params) + call EdgeForestReceiveParams(fates_params) call PRTReceiveParams(fates_params) call LeafBiophysReceiveParams(fates_params) call FatesSynchronizedParamsInst%ReceiveParams(fates_params) diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index fabddbec1c..8fb998a656 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -160,6 +160,9 @@ module FatesInterfaceTypesMod integer, public :: hlm_use_tree_damage ! This flag signals whether or not to turn on the ! tree damage module + integer, public :: hlm_use_edge_forest ! This flag signals whether or not to turn on the + ! edge forest module + integer, public :: hlm_hydr_solver ! Switch that defines which hydraulic solver to use ! 1 = Taylor solution that solves plant fluxes with 1 layer ! sequentially placing solution on top of previous layer solves @@ -299,6 +302,7 @@ module FatesInterfaceTypesMod integer , public, allocatable :: fates_hdim_pfmap_levscpf(:) ! map of pfts into size-class x pft dimension integer , public, allocatable :: fates_hdim_scmap_levscpf(:) ! map of size-class into size-class x pft dimension real(r8), public, allocatable :: fates_hdim_levage(:) ! patch age lower bound dimension + real(r8), public, allocatable :: fates_hdim_levedge(:) ! edge forest bin lower bound dimension real(r8), public, allocatable :: fates_hdim_levheight(:) ! height lower bound dimension integer , public, allocatable :: fates_hdim_levpft(:) ! plant pft dimension integer , public, allocatable :: fates_hdim_levlanduse(:) ! land use label dimension @@ -373,6 +377,7 @@ module FatesInterfaceTypesMod integer, public :: nlevcoage ! The total number of cohort age bins output to history integer, public :: nleafage ! The total number of leaf age classes integer, public :: nlevdamage ! The total number of damage classes + integer, public :: nlevedgeforest ! The total number of forest edge bins (incl. deep forest) ! ------------------------------------------------------------------------------------- ! Structured Boundary Conditions (SITE/PATCH SCALE) diff --git a/main/FatesParametersInterface.F90 b/main/FatesParametersInterface.F90 index feebf503b7..be33cb812c 100644 --- a/main/FatesParametersInterface.F90 +++ b/main/FatesParametersInterface.F90 @@ -33,6 +33,7 @@ module FatesParametersInterface character(len=*), parameter, public :: dimension_name_leaf_age = 'fates_leafage_class' character(len=*), parameter, public :: dimension_name_history_size_bins = 'fates_history_size_bins' character(len=*), parameter, public :: dimension_name_history_age_bins = 'fates_history_age_bins' + character(len=*), parameter, public :: dimension_name_edgeforest_bins = 'fates_edgeforest_bins' character(len=*), parameter, public :: dimension_name_history_height_bins = 'fates_history_height_bins' character(len=*), parameter, public :: dimension_name_history_coage_bins = 'fates_history_coage_bins' character(len=*), parameter, public :: dimension_name_hlm_pftno = 'fates_hlm_pftno' diff --git a/main/FatesUtilsMod.F90 b/main/FatesUtilsMod.F90 index 03537bd226..7067c86bbf 100644 --- a/main/FatesUtilsMod.F90 +++ b/main/FatesUtilsMod.F90 @@ -6,6 +6,7 @@ module FatesUtilsMod use FatesConstantsMod, only : r8 => fates_r8 use FatesGlobals, only : fates_log use FatesConstantsMod, only : nearzero + use FatesConstantsMod, only : fates_check_param_set use FatesGlobals, only : endrun => fates_endrun use shr_log_mod , only : errMsg => shr_log_errMsg @@ -21,6 +22,8 @@ module FatesUtilsMod public :: QuadraticRootsNSWC public :: QuadraticRootsSridharachary public :: ArrayNint + public :: is_param_set + public :: logical_to_real character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -299,5 +302,24 @@ subroutine QuadraticRootsSridharachary(a,b,c,root1,root2,err) end subroutine QuadraticRootsSridharachary + function is_param_set(param) + use shr_infnan_mod, only : isnan => shr_infnan_isnan + + real(r8), intent(in) :: param + + logical :: is_param_set + + is_param_set = .not. isnan(param) + if (is_param_set) then + is_param_set = param < fates_check_param_set + end if + end function is_param_set + + function logical_to_real(logical_in) result(real_out) + logical, intent(in) :: logical_in + real(r8) :: real_out + real_out = merge(1._r8, 0._r8, logical_in) + end function logical_to_real + ! ====================================================================================== end module FatesUtilsMod diff --git a/main/test/CMakeLists.txt b/main/test/CMakeLists.txt new file mode 100644 index 0000000000..72d2b9fc61 --- /dev/null +++ b/main/test/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(edge_forest_test) diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index e899d3ff7b..11e371addb 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -1,6 +1,7 @@ netcdf fates_params_default { dimensions: fates_NCWD = 4 ; + fates_edgeforest_bins = 6 ; fates_history_age_bins = 7 ; fates_history_coage_bins = 2 ; fates_history_damage_bins = 2 ; @@ -258,6 +259,36 @@ variables: double fates_dev_arbitrary_pft(fates_pft) ; fates_dev_arbitrary_pft:units = "unknown" ; fates_dev_arbitrary_pft:long_name = "Unassociated pft dimensioned free parameter that developers can use for testing arbitrary new hypotheses" ; + double fates_edgeforest_bin_edges(fates_edgeforest_bins) ; + fates_edgeforest_bin_edges:units = "m" ; + fates_edgeforest_bin_edges:long_name = "Boundaries of forest edge bins (for each bin, include value closest to zero)" ; + double fates_edgeforest_gaussian_amplitude(fates_edgeforest_bins) ; + fates_edgeforest_gaussian_amplitude:units = "unitless" ; + fates_edgeforest_gaussian_amplitude:long_name = "Amplitudes for calculating forest area in each edge bin (gaussian fit)" ; + double fates_edgeforest_gaussian_sigma(fates_edgeforest_bins) ; + fates_edgeforest_gaussian_sigma:units = "unitless" ; + fates_edgeforest_gaussian_sigma:long_name = "Sigmas for calculating forest area in each edge bin (gaussian fit)" ; + double fates_edgeforest_gaussian_center(fates_edgeforest_bins) ; + fates_edgeforest_gaussian_center:units = "unitless" ; + fates_edgeforest_gaussian_center:long_name = "Centers for calculating forest area in each edge bin (gaussian fit)" ; + double fates_edgeforest_lognormal_amplitude(fates_edgeforest_bins) ; + fates_edgeforest_lognormal_amplitude:units = "unitless" ; + fates_edgeforest_lognormal_amplitude:long_name = "Amplitudes for calculating forest area in each edge bin (lognormal fit)" ; + double fates_edgeforest_lognormal_sigma(fates_edgeforest_bins) ; + fates_edgeforest_lognormal_sigma:units = "unitless" ; + fates_edgeforest_lognormal_sigma:long_name = "Sigmas for calculating forest area in each edge bin (lognormal fit)" ; + double fates_edgeforest_lognormal_center(fates_edgeforest_bins) ; + fates_edgeforest_lognormal_center:units = "unitless" ; + fates_edgeforest_lognormal_center:long_name = "Centers for calculating forest area in each edge bin (lognormal fit)" ; + double fates_edgeforest_quadratic_a(fates_edgeforest_bins) ; + fates_edgeforest_quadratic_a:units = "unitless" ; + fates_edgeforest_quadratic_a:long_name = "x^2 coefficient for calculating forest area in each edge bin (quadratic fit)" ; + double fates_edgeforest_quadratic_b(fates_edgeforest_bins) ; + fates_edgeforest_quadratic_b:units = "unitless" ; + fates_edgeforest_quadratic_b:long_name = "x^1 coefficient for calculating forest area in each edge bin (quadratic fit)" ; + double fates_edgeforest_quadratic_c(fates_edgeforest_bins) ; + fates_edgeforest_quadratic_c:units = "unitless" ; + fates_edgeforest_quadratic_c:long_name = "x^0 coefficient for calculating forest area in each edge bin (quadratic fit)" ; double fates_fire_alpha_SH(fates_pft) ; fates_fire_alpha_SH:units = "m / (kw/m)**(2/3)" ; fates_fire_alpha_SH:long_name = "spitfire parameter, alpha scorch height, Equation 16 Thonicke et al 2010" ; @@ -783,6 +814,9 @@ variables: double fates_fire_threshold ; fates_fire_threshold:units = "kW/m" ; fates_fire_threshold:long_name = "spitfire parameter, fire intensity threshold for tracking fires that spread" ; + double fates_forest_tree_fraction_threshold ; + fates_forest_tree_fraction_threshold:units = "m2/m2" ; + fates_forest_tree_fraction_threshold:long_name = "Tree fraction above which patch is considered 'forest'" ; double fates_frag_cwd_fcel ; fates_frag_cwd_fcel:units = "unitless" ; fates_frag_cwd_fcel:long_name = "Cellulose fraction for CWD" ; @@ -1205,6 +1239,26 @@ data: fates_dev_arbitrary_pft = _, _, _, _, _, _, _, _, _, _, _, _, _, _ ; + fates_edgeforest_gaussian_amplitude = 0.37033665947126704, 0.25548576693720165, _, _, _, _ ; + + fates_edgeforest_gaussian_sigma = 0.4329573815987602, 0.42239011110175917, _, _, _, _ ; + + fates_edgeforest_gaussian_center = -0.35596450526414164, -0.28183332150023604, _, _, _, _ ; + + fates_edgeforest_lognormal_amplitude = _, _, _, _, _, 14.283227336273354 ; + + fates_edgeforest_lognormal_sigma = _, _, _, _, _, 1.26210031507715 ; + + fates_edgeforest_lognormal_center = _, _, _, _, _, 2.1958393621871597 ; + + fates_edgeforest_quadratic_a = _, _, 0.2549551300197741, 0.07685044819893726, 0.035189070666016925, _ ; + + fates_edgeforest_quadratic_b = _, _, -0.5457222474679617, -0.19438641157435982, -0.300528731650077, _ ; + + fates_edgeforest_quadratic_c = _, _, 0.29299184857665717, 0.11825507562859365, 0.2669694066063096, _ ; + + fates_edgeforest_bin_edges = 0, 30, 60, 120, 150, 300 ; + fates_fire_alpha_SH = 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2 ; @@ -1784,6 +1838,8 @@ data: fates_fire_threshold = 50 ; + fates_forest_tree_fraction_threshold = 0.5 ; + fates_frag_cwd_fcel = 0.76 ; fates_frag_cwd_flig = 0.24 ; diff --git a/testing/CMakeLists.txt b/testing/CMakeLists.txt index 3b788d51fa..c08e021e05 100644 --- a/testing/CMakeLists.txt +++ b/testing/CMakeLists.txt @@ -7,9 +7,12 @@ add_subdirectory(functional_testing/fire/fuel fates_fuel_ftest) add_subdirectory(functional_testing/fire/ros fates_ros_ftest) add_subdirectory(functional_testing/patch fates_patch_ftest) add_subdirectory(functional_testing/fire/mortality fates_firemort_ftest) +add_subdirectory(functional_testing/edge_forest fates_edge_forest_ftest) ## Unit tests add_subdirectory(unit_testing/fire_weather_test fates_fire_weather_utest) +add_subdirectory(unit_testing/ecotypes_test fates_ecotypes_utest) +add_subdirectory(unit_testing/edge_forest_test fates_edge_forest_utest) add_subdirectory(unit_testing/fire_fuel_test fates_fire_fuel_utest) add_subdirectory(unit_testing/sort_cohorts_test fates_sort_cohorts_utest) add_subdirectory(unit_testing/insert_cohort_test fates_insert_cohort_utest) diff --git a/testing/functional_testing/edge_forest/CMakeLists.txt b/testing/functional_testing/edge_forest/CMakeLists.txt new file mode 100644 index 0000000000..371f4d911c --- /dev/null +++ b/testing/functional_testing/edge_forest/CMakeLists.txt @@ -0,0 +1,24 @@ +set(edgeforest_sources FatesTestEdgeForest.F90) + +set(NETCDF_C_DIR ${NETCDF_C_PATH}) +set(NETCDF_FORTRAN_DIR ${NETCDF_F_PATH}) + +FIND_PATH(NETCDFC_FOUND libnetcdf.a ${NETCDF_C_DIR}/lib) +FIND_PATH(NETCDFF_FOUND libnetcdff.a ${NETCDF_FORTRAN_DIR}/lib) + + +include_directories(${NETCDF_C_DIR}/include + ${NETCDF_FORTRAN_DIR}/include) + +link_directories(${NETCDF_C_DIR}/lib + ${NETCDF_FORTRAN_DIR}/lib + ${PFUNIT_TOP_DIR}/lib) + +add_executable(FATES_edge_forest_exe ${edgeforest_sources}) + +target_link_libraries(FATES_edge_forest_exe + netcdf + netcdff + fates + csm_share + funit) diff --git a/testing/functional_testing/edge_forest/FatesTestEdgeForest.F90 b/testing/functional_testing/edge_forest/FatesTestEdgeForest.F90 new file mode 100644 index 0000000000..e3a903e234 --- /dev/null +++ b/testing/functional_testing/edge_forest/FatesTestEdgeForest.F90 @@ -0,0 +1,253 @@ +program FatesTestEdgeForest + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesUtilsMod, only : is_param_set + use FatesEdgeForestMod, only : GetFracEdgeForestInEachBin_norm, GetFracEdgeForestInEachBin_quadratic + use FatesEdgeForestMod, only : GetFracEdgeForestInEachBin + use FatesEdgeForestParamsMod, only : ED_val_edgeforest_gaussian_amplitude, ED_val_edgeforest_gaussian_sigma,ED_val_edgeforest_gaussian_center + use FatesEdgeForestParamsMod, only : ED_val_edgeforest_lognormal_amplitude, ED_val_edgeforest_lognormal_sigma,ED_val_edgeforest_lognormal_center + use FatesEdgeForestParamsMod, only : ED_val_edgeforest_quadratic_a, ED_val_edgeforest_quadratic_b,ED_val_edgeforest_quadratic_c + use FatesEdgeForestParamsMod, only : ED_val_edgeforest_bin_edges + use FatesUnitTestParamReaderMod, only :fates_unit_test_param_reader + use FatesArgumentUtils, only : command_line_arg + use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=), isnan => shr_infnan_isnan + + implicit none + + ! LOCALS: + type(fates_unit_test_param_reader) :: param_reader ! param reader instance + character(len=:), allocatable :: param_file ! input parameter file + integer :: e, i ! looping indices + integer :: n_bins ! number of edge forest bins in the parameter file + integer :: n_frac_forest ! size of frac_forest array + real(r8), allocatable :: frac_forest(:) ! fraction forest in site + real(r8), allocatable :: frac_in_bin_gaussian(:) ! output: fraction of forest in a bin with Gaussian fit + real(r8), allocatable :: frac_in_bin_lognormal(:) ! output: fraction of forest in a bin with lognormal fit + real(r8), allocatable :: frac_in_bin_quadratic(:) ! output: fraction of forest in a bin with quadratic fit + real(r8), allocatable :: frac_in_every_bin(:,:) ! output: fraction of forest in each bin + real(r8), allocatable :: frac_in_every_bin_norm(:,:) ! output: fraction of forest in each bin (normalized) + ! Edge bin parameters + real(r8) :: amplitude, mu, sigma, a, b, c + logical :: bin_found + + ! CONSTANTS: + logical, parameter :: debug = .false. + character(len=*), parameter :: out_file = 'edge_forest_out.nc' ! output file + real(r8), parameter :: min_frac_forest = 0._r8 ! minimum fraction forest to calculate + real(r8), parameter :: max_frac_forest = 1.0_r8 ! maximum fraction forest to calculate + real(r8), parameter :: frac_forest_inc = 0.001_r8 ! fraction forest increment to use + + interface + + subroutine WriteEdgeForestData(out_file, n_frac_forest, n_bins, frac_forest, frac_in_bin_gaussian, & + frac_in_bin_lognormal, frac_in_bin_quadratic, frac_in_every_bin, frac_in_every_bin_norm) + + use FatesUnitTestIOMod, only : OpenNCFile, RegisterNCDims, CloseNCFile + use FatesUnitTestIOMod, only : WriteVar, RegisterVar + use FatesUnitTestIOMod, only : type_double, type_int + use FatesConstantsMod, only : r8 => fates_r8 + implicit none + + character(len=*), intent(in) :: out_file + integer, intent(in) :: n_frac_forest + integer, intent(in) :: n_bins + real(r8), intent(in) :: frac_forest(:) + real(r8), intent(in) :: frac_in_bin_gaussian(:) + real(r8), intent(in) :: frac_in_bin_lognormal(:) + real(r8), intent(in) :: frac_in_bin_quadratic(:) + real(r8), intent(in) :: frac_in_every_bin(:,:) + real(r8), intent(in) :: frac_in_every_bin_norm(:,:) + end subroutine WriteEdgeForestData + + end interface + + ! read in parameter file name from command line + param_file = command_line_arg(1) + + ! read in parameter file + call param_reader%Init(param_file) + call param_reader%RetrieveParameters() + + ! determine sizes of arrays + n_bins = size(ED_val_edgeforest_bin_edges, dim=1) + n_frac_forest = int((max_frac_forest - min_frac_forest)/frac_forest_inc + 1) + + ! allocate arrays + allocate(frac_forest(n_frac_forest)) + allocate(frac_in_bin_gaussian(n_frac_forest)) + allocate(frac_in_bin_lognormal(n_frac_forest)) + allocate(frac_in_bin_quadratic(n_frac_forest)) + allocate(frac_in_every_bin(n_frac_forest, n_bins)) + allocate(frac_in_every_bin_norm(n_frac_forest, n_bins)) + + ! initialize frac_forest array + do i = 1, n_frac_forest + frac_forest(i) = min_frac_forest + frac_forest_inc*(i-1) + end do + + ! if debugging, print read-in parameters for all bins + if (debug) then + do e = 1, n_bins + + ! Gaussian + if (is_param_set(ED_val_edgeforest_gaussian_amplitude(e))) then + amplitude = ED_val_edgeforest_gaussian_amplitude(e) + mu = ED_val_edgeforest_gaussian_center(e) + sigma = ED_val_edgeforest_gaussian_sigma(e) + write(*, '(a, i2, a)') "Gaussian (bin ",e,"):" + write(*, '(a, E15.6)') " amplitude: ",amplitude + write(*, '(a, E15.6)') " mu: ",mu + write(*, '(a, E15.6)') " sigma: ",sigma + + ! Lognormal + else if (is_param_set(ED_val_edgeforest_lognormal_amplitude(e))) then + bin_found = .true. + amplitude = ED_val_edgeforest_lognormal_amplitude(e) + mu = ED_val_edgeforest_lognormal_center(e) + sigma = ED_val_edgeforest_lognormal_sigma(e) + write(*, '(a, i2, a)') "Lognormal (bin ",e,"):" + write(*, '(a, E15.6)') " amplitude: ",amplitude + write(*, '(a, E15.6)') " mu: ",mu + write(*, '(a, E15.6)') " sigma: ",sigma + + ! Quadratic + else if (is_param_set(ED_val_edgeforest_quadratic_a(e))) then + bin_found = .true. + a = ED_val_edgeforest_quadratic_a(e) + b = ED_val_edgeforest_quadratic_b(e) + c = ED_val_edgeforest_quadratic_c(e) + write(*, '(a, i2, a)') "Quadratic (bin ",e,"):" + write(*, '(a, E15.6)') " a: ",a + write(*, '(a, E15.6)') " b: ",b + write(*, '(a, E15.6)') " c: ",c + end if + end do + end if + + ! calculate fraction in every bin + do i = 1, n_frac_forest + + ! Pre-normalization + call GetFracEdgeForestInEachBin(frac_forest(i), n_bins, & + ED_val_edgeforest_gaussian_amplitude, ED_val_edgeforest_gaussian_sigma, ED_val_edgeforest_gaussian_center, & + ED_val_edgeforest_lognormal_amplitude, ED_val_edgeforest_lognormal_sigma, ED_val_edgeforest_lognormal_center, & + ED_val_edgeforest_quadratic_a, ED_val_edgeforest_quadratic_b, ED_val_edgeforest_quadratic_c, & + frac_in_every_bin(i,:), .false.) + + ! Normalized + call GetFracEdgeForestInEachBin(frac_forest(i), n_bins, & + ED_val_edgeforest_gaussian_amplitude, ED_val_edgeforest_gaussian_sigma, ED_val_edgeforest_gaussian_center, & + ED_val_edgeforest_lognormal_amplitude, ED_val_edgeforest_lognormal_sigma, ED_val_edgeforest_lognormal_center, & + ED_val_edgeforest_quadratic_a, ED_val_edgeforest_quadratic_b, ED_val_edgeforest_quadratic_c, & + frac_in_every_bin_norm(i,:), .true.) + end do + + ! write out data to netcdf file + call WriteEdgeForestData(out_file, n_frac_forest, n_bins, frac_forest, frac_in_bin_gaussian, & + frac_in_bin_lognormal, frac_in_bin_quadratic, frac_in_every_bin, frac_in_every_bin_norm) + +end program FatesTestEdgeForest + +! ---------------------------------------------------------------------------------------- + +subroutine WriteEdgeForestData(out_file, n_frac_forest, n_bins, frac_forest, frac_in_bin_gaussian, & + frac_in_bin_lognormal, frac_in_bin_quadratic, frac_in_every_bin, frac_in_every_bin_norm) + ! + ! DESCRIPTION: + ! Writes out data from the edge forest test + ! + use FatesUnitTestIOMod, only : OpenNCFile, RegisterNCDims, CloseNCFile + use FatesUnitTestIOMod, only : WriteVar + use FatesUnitTestIOMod, only : RegisterVar + use FatesUnitTestIOMod, only : EndNCDef + use FatesUnitTestIOMod, only : type_double, type_int + use FatesConstantsMod, only : r8 => fates_r8 + use FatesEdgeForestParamsMod, only : ED_val_edgeforest_bin_edges + implicit none + + ! ARGUMENTS + character(len=*), intent(in) :: out_file + integer, intent(in) :: n_frac_forest + integer, intent(in) :: n_bins + real(r8), intent(in) :: frac_forest(:) + real(r8), intent(in) :: frac_in_bin_gaussian(:) + real(r8), intent(in) :: frac_in_bin_lognormal(:) + real(r8), intent(in) :: frac_in_bin_quadratic(:) + real(r8), intent(in) :: frac_in_every_bin(:,:) + real(r8), intent(in) :: frac_in_every_bin_norm(:,:) + + ! LOCALS: + ! TODO: Set a local parameter for number of dimensions (2) + integer :: i ! looping index + integer :: ncid ! netcdf file id + character(len=24) :: dim_names(2) ! dimension name(s) + integer :: dimIDs(2) ! dimension ID(s) + integer :: fracforestID, binID ! variable ID(s) for dimensions + integer :: gaussianID, lognormalID, quadraticID + integer :: everybinID, everybinnormID + + ! dimension name(s) + dim_names = [character(len=24) :: 'frac_forest', 'bin'] + + ! open file + call OpenNCFile(trim(out_file), ncid, 'readwrite') + + ! register dimensions + call RegisterNCDims(ncid, dim_names, (/n_frac_forest, n_bins/), 2, dimIDs) + + ! register fraction forest + call RegisterVar(ncid, dim_names(1), dimIDs(1:1), type_double, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: 'unitless', 'Fraction of site that is forest'], 2, fracforestID) + + ! register bin + call RegisterVar(ncid, dim_names(2), dimIDs(2:2), type_double, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: 'unitless', 'FATES edge bin min. distance to nonforest'], 2, binID) + + ! register frac_in_bin_gaussian + call RegisterVar(ncid, 'frac_in_bin_gaussian', dimIDs(1:1), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'frac_forest', 'unitless', 'Fraction of forest in first bin with Gaussian fit'], & + 3, gaussianID) + + ! register frac_in_bin_lognormal + call RegisterVar(ncid, 'frac_in_bin_lognormal', dimIDs(1:1), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'frac_forest', 'unitless', 'Fraction of forest in first bin with lognormal fit'], & + 3, lognormalID) + + ! register frac_in_bin_quadratic + call RegisterVar(ncid, 'frac_in_bin_quadratic', dimIDs(1:1), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'frac_forest', 'unitless', 'Fraction of forest in first bin with quadratic fit'], & + 3, quadraticID) + + ! register frac_in_every_bin + call RegisterVar(ncid, 'frac_in_every_bin', dimIDs(1:2), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'frac_forest bin', 'unitless', 'Fraction of forest in every bin'], & + 3, everybinID) + + ! register frac_in_every_bin_norm + call RegisterVar(ncid, 'frac_in_every_bin_norm', dimIDs(1:2), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'frac_forest bin', 'unitless', 'Fraction of forest in every bin (normalized)'], & + 3, everybinnormID) + + ! finish defining variables + call EndNCDef(ncid) + + ! write out data + call WriteVar(ncid, fracforestID, frac_forest(:)) + call WriteVar(ncid, binID, ED_val_edgeforest_bin_edges(:)) + call WriteVar(ncid, gaussianID, frac_in_bin_gaussian(:)) + call WriteVar(ncid, lognormalID, frac_in_bin_lognormal(:)) + call WriteVar(ncid, quadraticID, frac_in_bin_quadratic(:)) + call WriteVar(ncid, everybinID, frac_in_every_bin(:,:)) + call WriteVar(ncid, everybinnormID, frac_in_every_bin_norm(:,:)) + + ! close the file + call CloseNCFile(ncid) + +end subroutine WriteEdgeForestData diff --git a/testing/functional_testing/edge_forest/edge_forest_test.py b/testing/functional_testing/edge_forest/edge_forest_test.py new file mode 100644 index 0000000000..68e7485e06 --- /dev/null +++ b/testing/functional_testing/edge_forest/edge_forest_test.py @@ -0,0 +1,173 @@ +""" +Concrete class for running the edge forest functional tests for FATES. +""" + +import os +import xarray as xr +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +from utils import round_up +from utils_plotting import sample_colormap, blank_plot +from functional_class import FunctionalTest + + +class EdgeForestTest(FunctionalTest): + """Quadratic test class""" + + name = "edge_forest" + + def __init__(self, test_dict): + super().__init__( + EdgeForestTest.name, + test_dict["test_dir"], + test_dict["test_exe"], + test_dict["out_file"], + test_dict["use_param_file"], + test_dict["other_args"], + ) + self.plot = True + + def plot_output(self, run_dir: str, save_figs: bool, plot_dir: str): + """Plots all edge forest plots + + Args: + run_dir (str): run directory + out_file (str): output file name + save_figs (bool): whether or not to save the figures + plot_dir (str): plot directory to save the figures to + """ + + # read in edge forest data + edge_forest_dat = xr.open_dataset(os.path.join(run_dir, self.out_file)) + + # Plot all bins + da = edge_forest_dat["frac_in_every_bin"] + da_norm = edge_forest_dat["frac_in_every_bin_norm"] + self.plot_edge_forest_frac_allbins( + da, + da_norm, + "Fraction of forest in each edge bin", + da.attrs["units"], + save_figs, + plot_dir, + ) + + # Plot individual bins + plot_dict = { + "frac_in_bin_gaussian": { + "varname": "Fraction of forest in first bin with Gaussian fit", + "units": "unitless", + }, + "frac_in_bin_lognormal": { + "varname": "Fraction of forest in first bin with lognormal fit", + "units": "unitless", + }, + "frac_in_bin_quadratic": { + "varname": "Fraction of forest in first bin with quadratic fit", + "units": "unitless", + }, + } + for plot, attributes in plot_dict.items(): + self.plot_edge_forest_frac_onebin( + edge_forest_dat[plot], + attributes["varname"], + attributes["units"], + save_figs, + plot_dir, + ) + + @staticmethod + def plot_edge_forest_frac_allbins( + data: xr.Dataset, data_norm: xr.Dataset, varname: str, units: str, save_fig: bool, plot_dir: str = None + ): + """Plot the fraction of forest in all bins + + Args: + data (xarray DataArray): the data array of the variable to plot + data_norm (xarray DataArray): the data array of the normalized variable to plot + var (str): variable name (for data structure) + varname (str): variable name for plot labels + units (str): variable units for plot labels + save_fig (bool): whether or not to write out plot + plot_dir (str): if saving figure, where to write to + """ + + x = data["frac_forest"] + max_x = x.max() + + y = data + y_norm = data_norm + max_y = round_up(y.max()) + + blank_plot(max_x, 0.0, max_y, 0.0, draw_horizontal_lines=True) + + bins = data["bin"].values + + for b, this_bin in enumerate(bins): + plt.plot( + x.values, + y.sel(bin=this_bin).values, + color=sample_colormap("jet_r", len(bins), b), + label=str(this_bin), + ) + + for b, this_bin in enumerate(bins): + plt.plot( + x.values, + y_norm.sel(bin=this_bin).values, + color=sample_colormap("jet_r", len(bins), b), + label="__nolegend__", + linestyle="--", + ) + + plt.xlabel("Frac. forest in site", fontsize=11) + plt.ylabel(f"{varname} ({units})", fontsize=11) + plt.title(f"Simulated {varname} for input parameter file\n(dashed=adjusted)", fontsize=11) + plt.legend(loc="best", title="Bin distance to\nnonforest (m)", alignment="left") + + if save_fig: + fig_name = os.path.join(plot_dir, f"edge_forest_plot_{data.name}.png") + plt.savefig(fig_name) + + @staticmethod + def plot_edge_forest_frac_onebin( + data: xr.Dataset, varname: str, units: str, save_fig: bool, plot_dir: str = None + ): + """Plot the fraction of forest in a given bin + + Args: + data (xarray DataArray): the data array of the variable to plot + var (str): variable name (for data structure) + varname (str): variable name for plot labels + units (str): variable units for plot labels + save_fig (bool): whether or not to write out plot + plot_dir (str): if saving figure, where to write to + """ + + # This is left over from AllometryTest, which had two dimensions + data_frame = pd.DataFrame( + { + "frac_forest": np.tile(data.frac_forest, 1), + data.name: data.values.flatten(), + } + ) + + max_forest = data_frame["frac_forest"].max() + max_var = round_up(data_frame[data.name].max()) + + blank_plot(max_forest, 0.0, max_var, 0.0, draw_horizontal_lines=True) + + plt.plot( + data_frame.frac_forest.values, + data_frame[data.name].values, + lw=2, + ) + + plt.xlabel("Frac. forest in site", fontsize=11) + plt.ylabel(f"{varname} ({units})", fontsize=11) + plt.title(f"Simulated {varname} for input parameter file", fontsize=11) + + if save_fig: + fig_name = os.path.join(plot_dir, f"edge_forest_plot_{data.name}.png") + plt.savefig(fig_name) diff --git a/testing/functional_testing/fire/mortality/fire_mortality_test.py b/testing/functional_testing/fire/mortality/fire_mortality_test.py index 58e6055d53..aa9a7a833d 100644 --- a/testing/functional_testing/fire/mortality/fire_mortality_test.py +++ b/testing/functional_testing/fire/mortality/fire_mortality_test.py @@ -7,7 +7,7 @@ import pandas as pd import matplotlib.pyplot as plt from functional_class import FunctionalTest -from utils import blank_plot, get_color_palette +from utils_plotting import blank_plot, get_color_palette class FireMortTest(FunctionalTest): """Fire mortality test class""" diff --git a/testing/functional_tests.cfg b/testing/functional_tests.cfg index fec8e075d6..a4aface964 100644 --- a/testing/functional_tests.cfg +++ b/testing/functional_tests.cfg @@ -40,3 +40,10 @@ test_exe = FATES_firemort_exe out_file = fire_mortality_out.nc use_param_file = True other_args = [] + +[edge_forest] +test_dir = fates_edge_forest_ftest +test_exe = FATES_edge_forest_exe +out_file = edge_forest_out.nc +use_param_file = True +other_args = [] diff --git a/testing/load_functional_tests.py b/testing/load_functional_tests.py index 3b72dcb2b7..e565466d0c 100644 --- a/testing/load_functional_tests.py +++ b/testing/load_functional_tests.py @@ -8,3 +8,4 @@ from functional_testing.fire.ros.ros_test import ROSTest from functional_testing.patch.patch_test import PatchTest from functional_testing.fire.mortality.fire_mortality_test import FireMortTest +from functional_testing.edge_forest.edge_forest_test import EdgeForestTest diff --git a/testing/testing_shr/FatesFactoryMod.F90 b/testing/testing_shr/FatesFactoryMod.F90 index 5f2cdee5fb..9f868e02ba 100644 --- a/testing/testing_shr/FatesFactoryMod.F90 +++ b/testing/testing_shr/FatesFactoryMod.F90 @@ -23,10 +23,13 @@ module FatesFactoryMod use EDParamsMod, only : nclmax use EDParamsMod, only : photo_temp_acclim_timescale use EDParamsMod, only : photo_temp_acclim_thome_time + use EDParamsMod, only : FatesParamsInitForFactory use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm use FatesRunningMeanMod, only : moving_ema_window, fixed_window use EDCohortDynamicsMod, only : InitPRTObject use PRTParametersMod, only : prt_params + use EDPftvarcon, only : EDPftvarcon_inst + use FatesParameterDerivedMod, only : param_derived use PRTGenericMod, only : element_pos use PRTGenericMod, only : num_elements use PRTGenericMod, only : element_list @@ -64,12 +67,74 @@ module FatesFactoryMod implicit none public :: GetSyntheticPatch + public :: InitializeParams public :: InitializeGlobals contains !--------------------------------------------------------------------------------------- + subroutine InitializeParams() + ! + ! DESCRIPTION: + ! Initialize parameters needed for running factory that usually come from the parameter file + + ! LOCALS: + integer, parameter :: n_pfts = 14 + integer, parameter :: n_leafage_class = 1 + + ! Things from parameter file for prt_params + call FatesParamsInitForFactory() + allocate(prt_params%allom_agb_frac(n_pfts)); prt_params%allom_agb_frac = [0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 1., 1., 1.] + allocate(prt_params%allom_agb1(n_pfts)); prt_params%allom_agb1 = [0.0673, 0.1364012, 0.0393057, 0.2653695, 0.0673, 0.0728698, 0.06896, 0.06896, 0.06896, 0.06896, 0.06896, 0.001, 0.001, 0.003] + allocate(prt_params%allom_agb2(n_pfts)); prt_params%allom_agb2 = [0.976, 0.9449041, 1.087335, 0.8321321, 0.976, 1.0373211, 0.572, 0.572, 0.572, 0.5289883, 0.6853945, 1.6592, 1.6592, 1.3456] + allocate(prt_params%allom_agb3(n_pfts)); prt_params%allom_agb3 = [1.94, 1.94, 1.94, 1.94, 1.94, 1.94, 1.94, 1.94, 1.94, 2.1010352, 1.7628613, 1.248, 1.248, 1.869] + allocate(prt_params%allom_agb4(n_pfts)); prt_params%allom_agb4 = [0.931, 0.931, 0.931, 0.931, 0.931, 0.931, 0.931, 0.931, 0.931, 0.931, 0.931, -999.9, -999.9, -999.9] + allocate(prt_params%allom_amode(n_pfts)); prt_params%allom_amode = [3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 5, 5, 5] + allocate(prt_params%allom_blca_expnt_diff(n_pfts)); prt_params%allom_blca_expnt_diff = [-0.12, -0.34, -0.32, -0.22, -0.12, -0.35, 0., 0., 0., 0., 0., -0.487, -0.487, -0.259] + allocate(prt_params%allom_cmode(n_pfts)); prt_params%allom_cmode(:) = 1 + allocate(prt_params%allom_d2bl1(n_pfts)); prt_params%allom_d2bl1 = [0.04, 0.07, 0.07, 0.01, 0.04, 0.07, 0.07, 0.07, 0.07, 0.0481934, 0.0481934, 0.0004, 0.0004, 0.0012] + allocate(prt_params%allom_d2bl2(n_pfts)); prt_params%allom_d2bl2 = [1.6019679, 1.5234373, 1.3051237, 1.9621397, 1.6019679, 1.3998939, 1.3, 1.3, 1.3, 1.0600586, 1.7176758, 1.7092, 1.7092, 1.5879] + allocate(prt_params%allom_d2bl3(n_pfts)); prt_params%allom_d2bl3 = [0.55, 0.55, 0.55, 0.55, 0.55, 0.55, 0.55, 0.55, 0.55, 0.55, 0.55, 0.3417, 0.3417, 0.9948] + allocate(prt_params%allom_d2ca_coefficient_max(n_pfts)); prt_params%allom_d2ca_coefficient_max = [0.2715891, 0.3693718, 1.0787259, 0.0579297, 0.2715891, 1.1553612, 0.6568464, 0.6568464, 0.6568464, 0.4363427, 0.3166497, 0.0408, 0.0408, 0.0862] + allocate(prt_params%allom_d2ca_coefficient_min(n_pfts)); prt_params%allom_d2ca_coefficient_min = prt_params%allom_d2ca_coefficient_max + allocate(prt_params%allom_d2h1(n_pfts)); prt_params%allom_d2h1 = [78.4087704, 306.842667, 106.8745821, 104.3586841, 78.4087704, 31.4557047, 0.64, 0.64, 0.64, 0.8165625, 0.778125, 0.1812, 0.1812, 0.3353] + allocate(prt_params%allom_d2h2(n_pfts)); prt_params%allom_d2h2 = [0.8124383, 0.752377, 0.9471302, 1.1146973, 0.8124383, 0.9734088, 0.37, 0.37, 0.37, 0.2316113, 0.4027002, 0.6384, 0.6384, 0.4235] + allocate(prt_params%allom_d2h3(n_pfts)); prt_params%allom_d2h3 = [47.6666164, 196.6865691, 93.9790461, 160.6835089, 47.6666164, 16.5928174, -999.9, -999.9, -999.9, -999.9, -999.9, -999.9, -999.9, -999.9] + allocate(prt_params%allom_dbh_maxheight(n_pfts)); prt_params%allom_dbh_maxheight = [1000., 1000., 1000., 1000., 1000., 1000., 3., 3., 2., 2.4, 1.9, 20., 20., 30.] + allocate(prt_params%allom_fmode(n_pfts)); prt_params%allom_fmode(:) = 1 + allocate(prt_params%allom_hmode(n_pfts)); prt_params%allom_hmode = [5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 1, 3, 3, 3] + allocate(prt_params%allom_l2fr(n_pfts)); prt_params%allom_l2fr = [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 0.67, 0.67, 1.41] + allocate(prt_params%allom_la_per_sa_int(n_pfts)); prt_params%allom_la_per_sa_int(:) = 0.8 + allocate(prt_params%allom_la_per_sa_slp(n_pfts)); prt_params%allom_la_per_sa_slp(:) = 0. + allocate(prt_params%allom_lmode(n_pfts)); prt_params%allom_lmode = [2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 5, 5, 5] + allocate(prt_params%allom_sai_scaler(n_pfts)); prt_params%allom_sai_scaler(:) = 0.1 + allocate(prt_params%allom_smode(n_pfts)); prt_params%allom_smode = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2] + allocate(prt_params%allom_stmode(n_pfts)); prt_params%allom_stmode(:) = 1 + allocate(prt_params%c2b(n_pfts)); prt_params%c2b(:) = 2 + allocate(prt_params%cushion(n_pfts)); prt_params%cushion = [1.2, 1.2, 1.2, 1.2, 2.4, 1.2, 1.2, 2.4, 1.2, 1.5, 1.4, 1.2, 1.2, 1.2] + allocate(prt_params%leaf_long(n_pfts,n_leafage_class)); prt_params%leaf_long = reshape([1.5, 4., 1., 1.5, 1., 1., 1.5, 1., 1., 1.5, 1., 1., 1., 1.], shape(prt_params%leaf_long)) + allocate(prt_params%leafn_vert_scaler_coeff1(n_pfts)); prt_params%leafn_vert_scaler_coeff1(:) = 0.00963 + allocate(prt_params%leafn_vert_scaler_coeff2(n_pfts)); prt_params%leafn_vert_scaler_coeff2(:) = 2.43 + allocate(prt_params%phen_fnrt_drop_fraction(n_pfts)); prt_params%phen_fnrt_drop_fraction(:) = 0. + allocate(prt_params%phen_leaf_habit(n_pfts)); prt_params%phen_leaf_habit = [1, 1, 2, 1, 3, 2, 1, 3, 2, 1, 2, 2, 3, 3] + allocate(prt_params%phen_stem_drop_fraction(n_pfts)); prt_params%phen_stem_drop_fraction(:) = 0. + allocate(prt_params%slamax(n_pfts)); prt_params%slamax = [0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.0954, 0.012, 0.03, 0.03, 0.012, 0.032, 0.05, 0.05, 0.05] + allocate(prt_params%slatop(n_pfts)); prt_params%slatop = [0.012, 0.005, 0.024, 0.009, 0.03, 0.03, 0.012, 0.03, 0.03, 0.01, 0.032, 0.027, 0.05, 0.05] + allocate(prt_params%wood_density(n_pfts)); prt_params%wood_density = [0.548327, 0.44235, 0.454845, 0.754336, 0.548327, 0.566452, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7] + + ! Things from parameter file for EDPftvarcon_inst + allocate(EDPftvarcon_inst%damage_frac(n_pfts)) + EDPftvarcon_inst%damage_frac = [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, & + 0.01, 0.01, 0.01, 0.01] + + ! Derived parameters + call param_derived%Init(n_pfts) + + end subroutine InitializeParams + + !--------------------------------------------------------------------------------------- + subroutine InitializeGlobals(step_size) ! ! DESCRIPTION: @@ -483,6 +548,10 @@ subroutine GetSyntheticPatch(patch_data, num_levsoil, patch) ! create the patch call PatchFactory(patch, patch_age, patch_data%area, num_swb, numpft, num_levsoil) + + ! Apply patch variables + patch%total_tree_area = patch_data%total_tree_area + patch%livegrass = patch_data%livegrass ! add cohorts do i = 1, patch_data%num_cohorts diff --git a/testing/testing_shr/FatesUnitTestParamReaderMod.F90 b/testing/testing_shr/FatesUnitTestParamReaderMod.F90 index d8b655136b..b0598bfa99 100644 --- a/testing/testing_shr/FatesUnitTestParamReaderMod.F90 +++ b/testing/testing_shr/FatesUnitTestParamReaderMod.F90 @@ -8,6 +8,7 @@ module FatesUnitTestParamReaderMod use FatesParametersInterface, only : dimension_shape_scalar, dimension_shape_1d, dimension_shape_2d use EDParamsMod, only : FatesRegisterParams, FatesReceiveParams use SFParamsMod, only : SpitFireRegisterParams, SpitFireReceiveParams + use FatesEdgeForestParamsMod, only : EdgeForestRegisterParams, EdgeForestReceiveParams use PRTInitParamsFatesMod, only : PRTRegisterParams, PRTReceiveParams use PRTParametersMod, only : prt_params use FatesInterfaceTypesMod, only : nleafage @@ -124,6 +125,7 @@ subroutine RetrieveParameters(this) call FatesRegisterParams(fates_params) call SpitFireRegisterParams(fates_params) + call EdgeForestRegisterParams(fates_params) call PRTRegisterParams(fates_params) call FatesSynchronizedParamsInst%RegisterParams(fates_params) call EDPftvarcon_inst%Register(fates_pft_params) @@ -133,6 +135,7 @@ subroutine RetrieveParameters(this) call FatesReceiveParams(fates_params) call SpitFireReceiveParams(fates_params) + call EdgeForestReceiveParams(fates_params) call PRTReceiveParams(fates_params) call FatesSynchronizedParamsInst%ReceiveParams(fates_params) call EDPftvarcon_inst%Receive(fates_pft_params) diff --git a/testing/testing_shr/SyntheticPatchTypes.F90 b/testing/testing_shr/SyntheticPatchTypes.F90 index 64747e5b6d..0401ddfaf2 100644 --- a/testing/testing_shr/SyntheticPatchTypes.F90 +++ b/testing/testing_shr/SyntheticPatchTypes.F90 @@ -9,6 +9,8 @@ module SyntheticPatchTypes private integer, parameter :: chunk_size = 10 + real(r8), public :: forest_tree_fraction_threshold = 0.5 + real(r8), public :: grass_biomass_threshold = 0._r8 ! patch data type to hold data about the synthetic patches type, public :: synthetic_patch_type @@ -20,6 +22,8 @@ module SyntheticPatchTypes real(r8), allocatable :: dbhs(:) ! cohort dbhs [cm] real(r8), allocatable :: densities(:) ! cohort densities [/m2] integer, allocatable :: canopy_layers(:) ! canopy layers + real(r8) :: total_tree_area ! tree cover area [m2] + real(r8) :: livegrass ! total aboveground grass biomass in patch [kgC/m2] integer, allocatable :: pft_ids(:) ! pft ids contains @@ -49,7 +53,7 @@ module SyntheticPatchTypes contains subroutine InitSyntheticPatchData(this, patch_id, patch_name, area, ages, dbhs, & - densities, pft_ids, canopy_layers) + densities, pft_ids, canopy_layers, total_tree_area, livegrass) ! ! DESCRIPTION: ! Initializes a synthetic patch with input characteristics @@ -65,6 +69,8 @@ subroutine InitSyntheticPatchData(this, patch_id, patch_name, area, ages, dbhs, real(r8), intent(in) :: densities(:) ! cohort densities [/m2] integer, intent(in) :: pft_ids(:) ! pft ids integer, intent(in) :: canopy_layers(:) ! canopy layers of cohorts + real(r8), intent(in) :: total_tree_area ! tree cover area [m2] + real(r8), intent(in) :: livegrass ! total aboveground grass biomass in patch [kgC/m2] ! LOCALS: integer :: num_cohorts ! number of cohorts on patch @@ -84,6 +90,8 @@ subroutine InitSyntheticPatchData(this, patch_id, patch_name, area, ages, dbhs, this%num_cohorts = num_cohorts this%patch_id = patch_id this%area = area + this%total_tree_area = total_tree_area + this%livegrass = livegrass do i = 1, num_cohorts this%ages(i) = ages(i) @@ -97,8 +105,8 @@ end subroutine InitSyntheticPatchData ! -------------------------------------------------------------------------------------- - subroutine AddPatch(this, patch_id, patch_name, area, ages, dbhs, densities, pft_ids, & - canopy_layers) + subroutine AddPatch(this, patch_name, area, ages, dbhs, densities, pft_ids, & + canopy_layers, total_tree_area, livegrass) ! ! DESCRIPTION: ! Adds a synthetic patch data to a dynamic array @@ -106,7 +114,6 @@ subroutine AddPatch(this, patch_id, patch_name, area, ages, dbhs, densities, pft ! ARGUMENTS: class(synthetic_patch_array_type), intent(inout) :: this ! array of synthetic patches - integer, intent(in) :: patch_id ! patch id character(len=*), intent(in) :: patch_name ! name of patch real(r8), intent(in) :: area ! patch area real(r8), intent(in) :: ages(:) ! cohort ages [yr] @@ -114,6 +121,8 @@ subroutine AddPatch(this, patch_id, patch_name, area, ages, dbhs, densities, pft real(r8), intent(in) :: densities(:) ! cohort densities [/m2] integer, intent(in) :: pft_ids(:) ! pft ids integer, intent(in) :: canopy_layers(:) ! canopy layers + real(r8), intent(in) :: total_tree_area ! tree cover area [m2] + real(r8), intent(in) :: livegrass ! total aboveground grass biomass in patch [kgC/m2] ! LOCALS: type(synthetic_patch_type) :: patch_data ! synthetic patch data @@ -137,8 +146,8 @@ subroutine AddPatch(this, patch_id, patch_name, area, ages, dbhs, densities, pft this%num_patches = 1 end if - call patch_data%InitSyntheticPatchData(patch_id, patch_name, area, ages, dbhs, & - densities, pft_ids, canopy_layers) + call patch_data%InitSyntheticPatchData(this%num_patches, patch_name, area, ages, dbhs, & + densities, pft_ids, canopy_layers, total_tree_area, livegrass) this%patches(this%num_patches) = patch_data @@ -189,7 +198,7 @@ end function PatchDataPosition ! -------------------------------------------------------------------------------------- - subroutine GetSyntheticPatchData(this) + subroutine GetSyntheticPatchData(this, patch_name_in) ! ! DESCRIPTION: ! Returns an array of hard-coded synthetic patch data @@ -198,41 +207,90 @@ subroutine GetSyntheticPatchData(this) ! ARGUMENTS: class(synthetic_patch_array_type), intent(inout) :: this ! array of synthetic patches + character(len=*), intent(in), optional :: patch_name_in + ! + ! LOCALS: + logical :: add_all + real(r8) :: patch_area = 500 + + add_all = .not. present(patch_name_in) + + if (add_all .or. patch_name_in == 'tropical') then + call this%AddPatch(patch_name='tropical', area=patch_area, & + ages=(/100.0_r8, 80.0_r8, 40.0_r8, 20.0_r8/), & + dbhs=(/60.0_r8, 50.0_r8, 25.0_r8, 10.0_r8/), & + densities=(/0.005_r8, 0.008_r8, 0.02_r8, 0.017_r8/), & + pft_ids=(/1, 1, 1, 1/), & + canopy_layers=(/1, 1, 2, 2/), & + total_tree_area=patch_area, & + livegrass=0._r8) + end if - call this%AddPatch(patch_id=1, patch_name='tropical', area=500.0_r8, & - ages=(/100.0_r8, 80.0_r8, 40.0_r8, 20.0_r8/), & - dbhs=(/60.0_r8, 50.0_r8, 25.0_r8, 10.0_r8/), & - densities=(/0.005_r8, 0.008_r8, 0.02_r8, 0.017_r8/), & - pft_ids=(/1, 1, 1, 1/), & - canopy_layers=(/1, 1, 2, 2/)) - - call this%AddPatch(patch_id=2, patch_name='evergreen', area=500.0_r8, & - ages=(/50.0_r8, 50.0_r8/), & - dbhs=(/30.0_r8, 25.0_r8/), & - densities=(/0.015_r8, 0.015_r8/), & - pft_ids=(/2, 2/), & - canopy_layers=(/1, 1/)) + if (add_all .or. patch_name_in == 'evergreen') then + call this%AddPatch(patch_name='evergreen', area=patch_area, & + ages=(/50.0_r8, 50.0_r8/), & + dbhs=(/30.0_r8, 25.0_r8/), & + densities=(/0.015_r8, 0.015_r8/), & + pft_ids=(/2, 2/), & + canopy_layers=(/1, 1/), & + total_tree_area=patch_area, & + livegrass=0._r8) + end if + + if (add_all .or. patch_name_in == 'savannah') then + call this%AddPatch(patch_name='savannah', area=patch_area, & + ages=(/20.0_r8, 1.0_r8/), & + dbhs=(/15.0_r8, 1.0_r8/), & + densities=(/0.015_r8, 0.015_r8/), & + pft_ids=(/5, 14/), & + canopy_layers=(/1, 2/), & + total_tree_area=patch_area * forest_tree_fraction_threshold, & + livegrass=grass_biomass_threshold) + end if + + if (add_all .or. patch_name_in == 'savannah_woody') then + call this%AddPatch(patch_name='savannah_woody', area=patch_area, & + ages=(/20.0_r8, 1.0_r8/), & + dbhs=(/15.0_r8, 1.0_r8/), & + densities=(/0.030_r8, 0.0075_r8/), & + pft_ids=(/5, 14/), & + canopy_layers=(/1, 2/), & + total_tree_area=patch_area * (forest_tree_fraction_threshold + 0.1), & + livegrass=grass_biomass_threshold - 1.0_r8) + end if - call this%AddPatch(patch_id=3, patch_name='savannah', area=500.0_r8, & - ages=(/20.0_r8, 1.0_r8/), & - dbhs=(/15.0_r8, 1.0_r8/), & - densities=(/0.015_r8, 0.015_r8/), & - pft_ids=(/5, 14/), & - canopy_layers=(/1, 2/)) + if (add_all .or. patch_name_in == 'savannah_grassy') then + call this%AddPatch(patch_name='savannah_grassy', area=patch_area, & + ages=(/20.0_r8, 1.0_r8/), & + dbhs=(/15.0_r8, 1.0_r8/), & + densities=(/0.0075_r8, 0.030_r8/), & + pft_ids=(/5, 14/), & + canopy_layers=(/1, 2/), & + total_tree_area=patch_area * forest_tree_fraction_threshold - 0.1, & + livegrass=grass_biomass_threshold + 1.0_r8) + end if - call this%AddPatch(patch_id=4, patch_name='grassland', area=500.0_r8, & - ages=(/1.0_r8, 2.0_r8/), & - dbhs=(/1.0_r8, 1.0_r8/), & - densities=(/0.015_r8, 0.015_r8/), & - pft_ids=(/13, 13/), & - canopy_layers=(/1, 1/)) + if (add_all .or. patch_name_in == 'grassland') then + call this%AddPatch(patch_name='grassland', area=patch_area, & + ages=(/1.0_r8, 2.0_r8/), & + dbhs=(/1.0_r8, 1.0_r8/), & + densities=(/0.015_r8, 0.015_r8/), & + pft_ids=(/13, 13/), & + canopy_layers=(/1, 1/), & + total_tree_area=0._r8, & + livegrass=grass_biomass_threshold + 2.0_r8) + end if - call this%AddPatch(patch_id=5, patch_name='temperate', area=500.0_r8, & - ages=(/80.0_r8, 50.0_r8, 20.0_r8, 5.0_r8/), & - dbhs=(/50.0_r8, 30.0_r8, 15.0_r8, 3.0_r8/), & - densities=(/0.005_r8, 0.01_r8, 0.015_r8, 0.005_r8/), & - pft_ids=(/6, 2, 2, 9/), & - canopy_layers=(/1, 1, 2, 2/)) + if (add_all .or. patch_name_in == 'temperate') then + call this%AddPatch(patch_name='temperate', area=patch_area, & + ages=(/80.0_r8, 50.0_r8, 20.0_r8, 5.0_r8/), & + dbhs=(/50.0_r8, 30.0_r8, 15.0_r8, 3.0_r8/), & + densities=(/0.005_r8, 0.01_r8, 0.015_r8, 0.005_r8/), & + pft_ids=(/6, 2, 2, 9/), & + canopy_layers=(/1, 1, 2, 2/), & + total_tree_area=patch_area, & + livegrass=0._r8) + end if end subroutine GetSyntheticPatchData diff --git a/testing/unit_testing/ecotypes_test/CMakeLists.txt b/testing/unit_testing/ecotypes_test/CMakeLists.txt new file mode 100644 index 0000000000..f27c5789ef --- /dev/null +++ b/testing/unit_testing/ecotypes_test/CMakeLists.txt @@ -0,0 +1,5 @@ +set(pfunit_sources test_Ecotypes.pf) + +add_pfunit_ctest(Ecotypes + TEST_SOURCES "${pfunit_sources}" + LINK_LIBRARIES fates csm_share) diff --git a/testing/unit_testing/ecotypes_test/test_Ecotypes.pf b/testing/unit_testing/ecotypes_test/test_Ecotypes.pf new file mode 100644 index 0000000000..097f4d37e2 --- /dev/null +++ b/testing/unit_testing/ecotypes_test/test_Ecotypes.pf @@ -0,0 +1,167 @@ +module test_Ecotypes + ! + ! DESCRIPTION: + ! Test the FATES ecotypes code + ! + use FatesConstantsMod, only : r8 => fates_r8 + use FatesEcotypesMod, only : DoesPatchHaveForest_TreeCover + use FatesEcotypesMod, only : DoesPatchHaveForest_GrassBiomass + use FatesEcotypesMod, only : IsPatchForest + use FatesFactoryMod, only : InitializeParams, InitializeGlobals, GetSyntheticPatch + use FatesPatchMod, only : fates_patch_type + use SyntheticPatchTypes, only : synthetic_patch_array_type + use SyntheticPatchTypes, only : forest_tree_fraction_threshold + use SyntheticPatchTypes, only : grass_biomass_threshold + use FatesArgumentUtils, only : command_line_arg + use funit + + implicit none + + ! LOCALS: + integer :: i + logical :: already_initialized = .false. + + ! CONSTANTS: + integer, parameter :: num_levsoil = 10 ! number of soil layers + real(r8), parameter :: step_size = 1800.0_r8 ! step-size [s] + + @TestCase + type, extends(TestCase) :: TestEcotypes + + type(synthetic_patch_array_type) :: patch_data + type(fates_patch_type), pointer :: patch + + contains + procedure :: setUp + + end type TestEcotypes + + real(r8), parameter :: tol = 1.e-7_r8 + real(r8), parameter :: nan = 0._r8 / 0._r8 + + integer, parameter :: n_to_sort = 5 + + + contains + + subroutine setUp(this) + class(TestEcotypes), intent(inout) :: this + + ! Only need to do this for the first test + if (.not. already_initialized) then + + ! initialize some global data we need + call InitializeParams() + call InitializeGlobals(step_size) + + already_initialized = .true. + end if + end subroutine setUp + + @Test + subroutine test_isforest_tropical(this) + ! Should have high tree cover and no grass, so true for all except grass check. + class(TestEcotypes), intent(inout) :: this + + call this%patch_data%GetSyntheticPatchData('tropical') + call GetSyntheticPatch(this%patch_data%patches(1), num_levsoil, this%patch) + + @assertTrue(DoesPatchHaveForest_TreeCover(this%patch, forest_tree_fraction_threshold)) + @assertTrue(IsPatchForest(this%patch, forest_tree_fraction_threshold)) + @assertFalse(DoesPatchHaveForest_GrassBiomass(this%patch, grass_biomass_threshold)) + @assertTrue(IsPatchForest(this%patch, forest_tree_fraction_threshold, grass_biomass_threshold)) + + end subroutine test_isforest_tropical + + @Test + subroutine test_isforest_evergreen(this) + ! Should have high tree cover and no grass, so true for all except grass check. + class(TestEcotypes), intent(inout) :: this + + call this%patch_data%GetSyntheticPatchData('evergreen') + call GetSyntheticPatch(this%patch_data%patches(1), num_levsoil, this%patch) + + @assertTrue(DoesPatchHaveForest_TreeCover(this%patch, forest_tree_fraction_threshold)) + @assertTrue(IsPatchForest(this%patch, forest_tree_fraction_threshold)) + @assertFalse(DoesPatchHaveForest_GrassBiomass(this%patch, grass_biomass_threshold)) + @assertTrue(IsPatchForest(this%patch, forest_tree_fraction_threshold, grass_biomass_threshold)) + + end subroutine test_isforest_evergreen + + @Test + subroutine test_isforest_savannah(this) + ! Exactly at the thresholds? All should be false. + class(TestEcotypes), intent(inout) :: this + + call this%patch_data%GetSyntheticPatchData('savannah') + call GetSyntheticPatch(this%patch_data%patches(1), num_levsoil, this%patch) + + @assertFalse(DoesPatchHaveForest_TreeCover(this%patch, forest_tree_fraction_threshold)) + @assertFalse(IsPatchForest(this%patch, forest_tree_fraction_threshold)) + @assertFalse(DoesPatchHaveForest_GrassBiomass(this%patch, grass_biomass_threshold)) + @assertFalse(IsPatchForest(this%patch, forest_tree_fraction_threshold, grass_biomass_threshold)) + + end subroutine test_isforest_savannah + + @Test + subroutine test_isforest_savannah_woody(this) + ! Should have high tree cover and little grass, so true for all except grass check. + class(TestEcotypes), intent(inout) :: this + + call this%patch_data%GetSyntheticPatchData('savannah_woody') + call GetSyntheticPatch(this%patch_data%patches(1), num_levsoil, this%patch) + + @assertTrue(DoesPatchHaveForest_TreeCover(this%patch, forest_tree_fraction_threshold)) + @assertTrue(IsPatchForest(this%patch, forest_tree_fraction_threshold)) + @assertFalse(DoesPatchHaveForest_GrassBiomass(this%patch, grass_biomass_threshold)) + @assertTrue(IsPatchForest(this%patch, forest_tree_fraction_threshold, grass_biomass_threshold)) + + end subroutine test_isforest_savannah_woody + + @Test + subroutine test_isforest_savannah_grassy(this) + ! Should have low tree cover and plenty of grass, so false for all except grass check. + class(TestEcotypes), intent(inout) :: this + + call this%patch_data%GetSyntheticPatchData('savannah_grassy') + call GetSyntheticPatch(this%patch_data%patches(1), num_levsoil, this%patch) + + @assertFalse(DoesPatchHaveForest_TreeCover(this%patch, forest_tree_fraction_threshold)) + @assertFalse(IsPatchForest(this%patch, forest_tree_fraction_threshold)) + @assertTrue(DoesPatchHaveForest_GrassBiomass(this%patch, grass_biomass_threshold)) + @assertFalse(IsPatchForest(this%patch, forest_tree_fraction_threshold, grass_biomass_threshold)) + + end subroutine test_isforest_savannah_grassy + + @Test + subroutine test_isforest_grassland(this) + ! Should have no trees and plenty of grass, so false for all except grass check. + class(TestEcotypes), intent(inout) :: this + + call this%patch_data%GetSyntheticPatchData('grassland') + call GetSyntheticPatch(this%patch_data%patches(1), num_levsoil, this%patch) + + @assertFalse(DoesPatchHaveForest_TreeCover(this%patch, forest_tree_fraction_threshold)) + @assertFalse(IsPatchForest(this%patch, forest_tree_fraction_threshold)) + @assertTrue(DoesPatchHaveForest_GrassBiomass(this%patch, grass_biomass_threshold)) + @assertFalse(IsPatchForest(this%patch, forest_tree_fraction_threshold, grass_biomass_threshold)) + + end subroutine test_isforest_grassland + + @Test + subroutine test_isforest_temperate(this) + ! Should have high tree cover and no grass, so true for all except grass check. + class(TestEcotypes), intent(inout) :: this + + call this%patch_data%GetSyntheticPatchData('temperate') + call GetSyntheticPatch(this%patch_data%patches(1), num_levsoil, this%patch) + + @assertTrue(DoesPatchHaveForest_TreeCover(this%patch, forest_tree_fraction_threshold)) + @assertTrue(IsPatchForest(this%patch, forest_tree_fraction_threshold)) + @assertFalse(DoesPatchHaveForest_GrassBiomass(this%patch, grass_biomass_threshold)) + @assertTrue(IsPatchForest(this%patch, forest_tree_fraction_threshold, grass_biomass_threshold)) + + end subroutine test_isforest_temperate + + + end module test_Ecotypes diff --git a/testing/unit_testing/edge_forest_test/CMakeLists.txt b/testing/unit_testing/edge_forest_test/CMakeLists.txt new file mode 100644 index 0000000000..63888cccbd --- /dev/null +++ b/testing/unit_testing/edge_forest_test/CMakeLists.txt @@ -0,0 +1,5 @@ +set(pfunit_sources test_EdgeForest.pf) + +add_pfunit_ctest(EdgeForest + TEST_SOURCES "${pfunit_sources}" + LINK_LIBRARIES fates csm_share) diff --git a/testing/unit_testing/edge_forest_test/test_EdgeForest.pf b/testing/unit_testing/edge_forest_test/test_EdgeForest.pf new file mode 100644 index 0000000000..0bf0b12880 --- /dev/null +++ b/testing/unit_testing/edge_forest_test/test_EdgeForest.pf @@ -0,0 +1,637 @@ +module test_EdgeForest + ! + ! DESCRIPTION: + ! Test the FATES edge forest code + ! + use FatesConstantsMod, only : r8 => fates_r8 + use FatesEdgeForestMod, only : indexx + use FatesEdgeForestMod, only : GetFracEdgeForestInEachBin + use FatesEdgeForestMod, only : GetFracEdgeForestInEachBin_norm_numerator, GetFracEdgeForestInEachBin_norm_denominator, GetFracEdgeForestInEachBin_quadratic + use FatesEdgeForestMod, only : GetFracEdgeForestInEachBin_norm + use FatesEdgeForestMod, only : AssignPatchToBins + use funit + + implicit none + + @TestCase + type, extends(TestCase) :: TestEdgeForest + + real(r8), dimension(:), allocatable :: array_to_sort + integer, dimension(:), allocatable :: sorted_indices + + contains + procedure :: setUp + procedure :: tearDown + end type TestEdgeForest + + real(r8), parameter :: tol = 1.e-7_r8 + real(r8), parameter :: nan = 0._r8 / 0._r8 + + integer, parameter :: n_to_sort = 5 + + + contains + + subroutine setUp(this) + class(TestEdgeForest), intent(inout) :: this + allocate(this%array_to_sort(n_to_sort)) + allocate(this%sorted_indices(n_to_sort)) + end subroutine setUp + + subroutine tearDown(this) + class(TestEdgeForest), intent(inout) :: this + if (allocated(this%array_to_sort)) deallocate(this%array_to_sort) + if (allocated(this%sorted_indices)) deallocate(this%sorted_indices) + end subroutine tearDown + + @Test + subroutine indexx_alreadySorted(this) + class(TestEdgeForest), intent(inout) :: this + + this%array_to_sort = (/ 1._r8, 1.62_r8, 2.72_r8, 3.14_r8, 6.28_r8 /) + + call indexx(this%array_to_sort, this%sorted_indices) + + @assertEqual((/ 1, 2, 3, 4, 5 /), this%sorted_indices) + + end subroutine indexx_alreadySorted + + @Test + subroutine indexx_reverseSorted(this) + class(TestEdgeForest), intent(inout) :: this + + this%array_to_sort = (/ 6.28_r8, 3.14_r8, 2.72_r8, 1.62_r8, 1._r8 /) + + call indexx(this%array_to_sort, this%sorted_indices) + + @assertEqual((/ 5, 4, 3, 2, 1 /), this%sorted_indices) + + end subroutine indexx_reverseSorted + + @Test + subroutine indexx_lowTie(this) + class(TestEdgeForest), intent(inout) :: this + + this%array_to_sort = (/ 1._r8, 1._r8, 2.72_r8, 3.14_r8, 6.28_r8 /) + + call indexx(this%array_to_sort, this%sorted_indices) + + @assertEqual((/ 1, 2, 3, 4, 5 /), this%sorted_indices) + + end subroutine indexx_lowTie + + @Test + subroutine indexx_highTie(this) + class(TestEdgeForest), intent(inout) :: this + + this%array_to_sort = (/ 1._r8, 1.62_r8, 2.72_r8, 3.14_r8, 3.14_r8 /) + + call indexx(this%array_to_sort, this%sorted_indices) + + @assertEqual((/ 1, 2, 3, 4, 5 /), this%sorted_indices) + + end subroutine indexx_highTie + + @Test + subroutine indexx_random(this) + class(TestEdgeForest), intent(inout) :: this + + this%array_to_sort = (/ 3._r8, 8._r8, 10._r8, 2._r8, 7._r8 /) + + call indexx(this%array_to_sort, this%sorted_indices) + + @assertEqual((/ 4, 1, 5, 2, 3 /), this%sorted_indices) + + end subroutine indexx_random + + + @Test + subroutine indexx_all_equal(this) + class(TestEdgeForest), intent(inout) :: this + + this%array_to_sort = (/ 1._r8, 1._r8, 1._r8, 1._r8, 1._r8 /) + + call indexx(this%array_to_sort, this%sorted_indices) + + @assertEqual((/ 1, 2, 3, 4, 5 /), this%sorted_indices) + + end subroutine indexx_all_equal + + + @Test + subroutine test_GetFracEdgeForestInEachBin_with_gaussian(this) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: frac_forest = 0.5_r8 + integer, parameter :: nlevedgeforest_tmp = 1 + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_amplitude = (/ 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_sigma = (/ 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_center = (/ 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_nan = (/ nan /) + logical :: norm = .false. ! DON'T normalize; we want the raw output of Gaussian + ! Output + real(r8), dimension(nlevedgeforest_tmp) :: fraction_forest_in_bin + + call GetFracEdgeForestInEachBin(frac_forest, nlevedgeforest_tmp, & + efb_params_amplitude, efb_params_sigma, efb_params_center, & + efb_params_nan, efb_params_nan, efb_params_nan, & + efb_params_nan, efb_params_nan, efb_params_nan, & + fraction_forest_in_bin, norm) + + @assertEqual(0.3520653267642879_r8, fraction_forest_in_bin(1), tol) + + end subroutine test_GetFracEdgeForestInEachBin_with_gaussian + + + @Test + subroutine test_GetFracEdgeForestInEachBin_with_lognorm(this) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: frac_forest = 0.5_r8 + integer, parameter :: nlevedgeforest_tmp = 1 + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_amplitude = (/ 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_sigma = (/ 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_center = (/ 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_nan = (/ nan /) + logical :: norm = .false. ! DON'T normalize; we want the raw output of Lognormal + ! Output + real(r8), dimension(nlevedgeforest_tmp) :: fraction_forest_in_bin + + call GetFracEdgeForestInEachBin(frac_forest, nlevedgeforest_tmp, & + efb_params_nan, efb_params_nan, efb_params_nan, & + efb_params_amplitude, efb_params_sigma, efb_params_center, & + efb_params_nan, efb_params_nan, efb_params_nan, & + fraction_forest_in_bin, norm) + + @assertEqual(0.19029780481010555, fraction_forest_in_bin(1), tol) + + end subroutine test_GetFracEdgeForestInEachBin_with_lognorm + + + @Test + subroutine test_GetFracEdgeForestInEachBin_with_quadratic(this) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: frac_forest = 0.5_r8 + integer, parameter :: nlevedgeforest_tmp = 1 + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_a = (/ 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_b = (/ 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_c = (/ 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_nan = (/ nan /) + logical :: norm = .false. ! DON'T normalize; we want the raw output of Quadratic + ! Output + real(r8), dimension(nlevedgeforest_tmp) :: fraction_forest_in_bin + + call GetFracEdgeForestInEachBin(frac_forest, nlevedgeforest_tmp, & + efb_params_nan, efb_params_nan, efb_params_nan, & + efb_params_nan, efb_params_nan, efb_params_nan, & + efb_params_a, efb_params_b, efb_params_c, & + fraction_forest_in_bin, norm) + + @assertEqual(1.75, fraction_forest_in_bin(1), tol) + + end subroutine test_GetFracEdgeForestInEachBin_with_quadratic + + + @Test + subroutine test_GetFracEdgeForestInEachBin_with_norm(this) + ! Test that normalization works correctly: If all bins have the same parameters, they should + ! get normalized to 1/nbins. + class(TestEdgeForest), intent(inout) :: this + real(r8) :: frac_forest = 0.5_r8 + integer, parameter :: nlevedgeforest_tmp = 3 + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_a = (/ 1._r8, 1._r8, 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_b = (/ 1._r8, 1._r8, 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_c = (/ 1._r8, 1._r8, 1._r8 /) + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_nan = (/ nan, nan, nan /) + logical :: norm = .true. + ! Output + real(r8), dimension(nlevedgeforest_tmp) :: fraction_forest_in_bin + real(r8) :: expected = 1._r8 / real(nlevedgeforest_tmp, r8) + + call GetFracEdgeForestInEachBin(frac_forest, nlevedgeforest_tmp, & + efb_params_nan, efb_params_nan, efb_params_nan, & + efb_params_nan, efb_params_nan, efb_params_nan, & + efb_params_a, efb_params_b, efb_params_c, & + fraction_forest_in_bin, norm) + + @assertEqual(expected, fraction_forest_in_bin(1), tol) + @assertEqual(expected, fraction_forest_in_bin(2), tol) + @assertEqual(expected, fraction_forest_in_bin(3), tol) + + end subroutine test_GetFracEdgeForestInEachBin_with_norm + + + @Test + subroutine test_GetFracEdgeForestInEachBin_x0(this) + ! If frac forest is zero, all forest should be in first edge bin (closest to edge) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: frac_forest = 0._r8 + integer, parameter :: nlevedgeforest_tmp = 3 + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_dummy = (/ 1._r8, 2._r8, 3._r8 /) + logical :: norm = .false. ! Shouldn't matter for this test + ! Output + real(r8), dimension(nlevedgeforest_tmp) :: fraction_forest_in_bin + + call GetFracEdgeForestInEachBin(frac_forest, nlevedgeforest_tmp, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + fraction_forest_in_bin, norm) + + @assertEqual( (/ 1._r8, 0._r8, 0._r8 /) , fraction_forest_in_bin) + + end subroutine test_GetFracEdgeForestInEachBin_x0 + + @Test + subroutine test_GetFracEdgeForestInEachBin_x0_norm(this) + ! If frac forest is zero, all forest should be in first edge bin (closest to edge) + ! Same test as test_GetFracEdgeForestInEachBin_x0, except norm=.true. + ! Shouldn't matter! + class(TestEdgeForest), intent(inout) :: this + real(r8) :: frac_forest = 0._r8 + integer, parameter :: nlevedgeforest_tmp = 3 + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_dummy = (/ 1._r8, 2._r8, 3._r8 /) + logical :: norm = .true. ! Shouldn't matter for this test + ! Output + real(r8), dimension(nlevedgeforest_tmp) :: fraction_forest_in_bin + + call GetFracEdgeForestInEachBin(frac_forest, nlevedgeforest_tmp, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + fraction_forest_in_bin, norm) + + @assertEqual( (/ 1._r8, 0._r8, 0._r8 /) , fraction_forest_in_bin) + + end subroutine test_GetFracEdgeForestInEachBin_x0_norm + + @Test + subroutine test_GetFracEdgeForestInEachBin_x1(this) + ! If frac forest is one, all forest should be in last edge bin ("deep forest") + class(TestEdgeForest), intent(inout) :: this + real(r8) :: frac_forest = 1._r8 + integer, parameter :: nlevedgeforest_tmp = 3 + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_dummy = (/ 1._r8, 2._r8, 3._r8 /) + logical :: norm = .false. ! Shouldn't matter for this test + ! Output + real(r8), dimension(nlevedgeforest_tmp) :: fraction_forest_in_bin + + call GetFracEdgeForestInEachBin(frac_forest, nlevedgeforest_tmp, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + fraction_forest_in_bin, norm) + + @assertEqual( (/ 0._r8, 0._r8, 1._r8 /) , fraction_forest_in_bin) + + end subroutine test_GetFracEdgeForestInEachBin_x1 + + @Test + subroutine test_GetFracEdgeForestInEachBin_x1_norm(this) + ! If frac forest is one, all forest should be in last edge bin ("deep forest"). + ! Same test as test_GetFracEdgeForestInEachBin_x1, except norm=.true. + ! Shouldn't matter! + class(TestEdgeForest), intent(inout) :: this + real(r8) :: frac_forest = 1._r8 + integer, parameter :: nlevedgeforest_tmp = 3 + real(r8), dimension(nlevedgeforest_tmp) :: efb_params_dummy = (/ 1._r8, 2._r8, 3._r8 /) + logical :: norm = .true. ! Shouldn't matter for this test + ! Output + real(r8), dimension(nlevedgeforest_tmp) :: fraction_forest_in_bin + + call GetFracEdgeForestInEachBin(frac_forest, nlevedgeforest_tmp, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + efb_params_dummy, efb_params_dummy, efb_params_dummy, & + fraction_forest_in_bin, norm) + + @assertEqual( (/ 0._r8, 0._r8, 1._r8 /) , fraction_forest_in_bin) + + end subroutine test_GetFracEdgeForestInEachBin_x1_norm + + @Test + subroutine test_gffeb_lognorm_numerator(this) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: x + real(r8) :: A ! Amplitude + real(r8) :: mu ! Center + real(r8) :: sigma ! Sigma + real(r8) :: expected + real(r8) :: actual + logical :: lognorm + + x = 10._r8 + ! Bin 4 + A = 1.601064269911203 + sigma = 0.8465094354405984 + mu = 1.973996133523811 + expected = 1.4848754270133655 + lognorm = .true. + + actual = GetFracEdgeForestInEachBin_norm_numerator(x, A, mu, sigma, lognorm) + + @assertEqual(expected, actual, tolerance=tol) + + end subroutine test_gffeb_lognorm_numerator + + @Test + subroutine test_gffeb_gaussian_numerator(this) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: x + real(r8) :: A ! Amplitude + real(r8) :: mu ! Center + real(r8) :: sigma ! Sigma + real(r8) :: expected + real(r8) :: actual + logical :: lognorm + + x = 0.87 + A = 2 + sigma = 3 + mu = 4 + expected = 1.160527865802580 + lognorm = .false. + + actual = GetFracEdgeForestInEachBin_norm_numerator(x, A, mu, sigma, lognorm) + + @assertEqual(expected, actual, tolerance=tol) + + end subroutine test_gffeb_gaussian_numerator + + @Test + subroutine test_GetFracEdgeForestInEachBin_norm_lognorm(this) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: x + real(r8) :: A ! Amplitude + real(r8) :: mu ! Center + real(r8) :: sigma ! Sigma + real(r8) :: expected + real(r8) :: actual + logical :: lognorm + + x = 0.87 + A = 2 + sigma = 3 + mu = 4 + expected = 0.11800808319557105 + lognorm = .true. + + actual = GetFracEdgeForestInEachBin_norm(x, A, mu, sigma, lognorm) + + @assertEqual(expected, actual, tolerance=tol) + + end subroutine test_GetFracEdgeForestInEachBin_norm_lognorm + + @Test + subroutine test_GetFracEdgeForestInEachBin_norm_gaussian(this) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: x + real(r8) :: A ! Amplitude + real(r8) :: mu ! Center + real(r8) :: sigma ! Sigma + real(r8) :: expected + real(r8) :: actual + logical :: lognorm + + x = 0.87 + A = 2 + sigma = 3 + mu = 4 + expected = 0.1543278780068184 + lognorm = .false. + + actual = GetFracEdgeForestInEachBin_norm(x, A, mu, sigma, lognorm) + + @assertEqual(expected, actual, tolerance=tol) + + end subroutine test_GetFracEdgeForestInEachBin_norm_gaussian + + @Test + subroutine test_gffeb_lognorm_denominator(this) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: x + real(r8) :: A ! Amplitude + real(r8) :: mu ! Center + real(r8) :: sigma ! Sigma + real(r8) :: expected + real(r8) :: actual + logical :: lognorm + + x = 10._r8 + ! Bin 4 + sigma = 0.8465094354405984 + expected = 21.21884485617329 + lognorm = .true. + + actual = GetFracEdgeForestInEachBin_norm_denominator(x, sigma, lognorm) + + @assertEqual(expected, actual, tolerance=tol) + + end subroutine test_gffeb_lognorm_denominator + + @Test + subroutine test_gffeb_gaussian_denominator(this) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: x + real(r8) :: A ! Amplitude + real(r8) :: mu ! Center + real(r8) :: sigma ! Sigma + real(r8) :: expected + real(r8) :: actual + logical :: lognorm + + x = 0.5_r8 + ! Bin 4 + sigma = 0.8465094354405984 + expected = 2.121884538742686 + lognorm = .false. + + actual = GetFracEdgeForestInEachBin_norm_denominator(x, sigma, lognorm) + + @assertEqual(expected, actual, tolerance=tol) + + end subroutine test_gffeb_gaussian_denominator + + @Test + subroutine test_GetFracEdgeForestInEachBin_quadratic(this) + class(TestEdgeForest), intent(inout) :: this + real(r8) :: x = 0.4_r8 + real(r8) :: a = 1987._r8 + real(r8) :: b = 2025._r8 + real(r8) :: c = 1984._r8 + real(r8) :: expected = 3111.92_r8 + real(r8) :: actual + + actual = GetFracEdgeForestInEachBin_quadratic(x, a, b, c) + + @assertEqual(expected, actual, tolerance=tol) + + end subroutine test_GetFracEdgeForestInEachBin_quadratic + + @Test + subroutine test_AssignPatchToBins_01(this) + ! The site has just one forest patch. All of the site's forest area is in the bin closest to edge. + class(TestEdgeForest), intent(inout) :: this + integer, parameter :: nlevedgeforest_tmp = 3 + + ! The site has 100 area of forest + real(r8) :: area_forest_patches = 100._r8 + + ! The size of this patch is 100 area + real(r8) :: patch_area = 100._r8 + + ! So far none of the site's forest area has been allocated to any bin + real(r8) :: sum_forest_bins_so_far_m2 = 0._r8 + + ! Out + real(r8), dimension(nlevedgeforest_tmp) :: area_in_edgeforest_bins + + ! All of the site's forest is in the bin closest to edge + real(r8), dimension(nlevedgeforest_tmp), target :: fraction_forest_in_each_bin + fraction_forest_in_each_bin = (/ 1._r8, 0._r8, 0._r8 /) + + call AssignPatchToBins(fraction_forest_in_each_bin, area_forest_patches, patch_area, & + nlevedgeforest_tmp, tol, sum_forest_bins_so_far_m2, area_in_edgeforest_bins) + + ! 0+100=100 area of the site's forest has been assigned after this call + @assertEqual(100._r8, sum_forest_bins_so_far_m2, tolerance=tol) + + ! All of the patch's area was assigned to the first edge bin + @assertEqual(100._r8, area_in_edgeforest_bins(1), tolerance=tol) + end subroutine test_AssignPatchToBins_01 + + @Test + subroutine test_AssignPatchToBins_02(this) + ! The site has multiple forest patches. All of the site's forest area is in the bin closest to edge. + class(TestEdgeForest), intent(inout) :: this + integer, parameter :: nlevedgeforest_tmp = 3 + + ! The site has 300 area of forest + real(r8) :: area_forest_patches = 300._r8 + + ! The size of this patch is 100 area + real(r8) :: patch_area = 100._r8 + + ! So far none of the site's forest area has been allocated to any bin + real(r8) :: sum_forest_bins_so_far_m2 = 0._r8 + + ! Out + real(r8), dimension(nlevedgeforest_tmp) :: area_in_edgeforest_bins + + ! All of the site's forest is in the bin closest to edge + real(r8), dimension(nlevedgeforest_tmp), target :: fraction_forest_in_each_bin + fraction_forest_in_each_bin = (/ 1._r8, 0._r8, 0._r8 /) + + call AssignPatchToBins(fraction_forest_in_each_bin, area_forest_patches, patch_area, & + nlevedgeforest_tmp, tol, sum_forest_bins_so_far_m2, area_in_edgeforest_bins) + + ! 0+100=100 area of the site's forest has been assigned after this call + @assertEqual(100._r8, sum_forest_bins_so_far_m2, tolerance=tol) + + ! All of the patch's area was assigned to the first edge bin + @assertEqual(100._r8, area_in_edgeforest_bins(1), tolerance=tol) + @assertEqual(0._r8, area_in_edgeforest_bins(2), tolerance=tol) + @assertEqual(0._r8, area_in_edgeforest_bins(3), tolerance=tol) + end subroutine test_AssignPatchToBins_02 + + @Test + subroutine test_AssignPatchToBins_03(this) + ! The site has just one forest patch. The site's forest split evenly across the two bins closest to edge. + class(TestEdgeForest), intent(inout) :: this + integer, parameter :: nlevedgeforest_tmp = 3 + + ! The site has 100 area of forest + real(r8) :: area_forest_patches = 100._r8 + + ! The size of this patch is 100 area + real(r8) :: patch_area = 100._r8 + + ! So far none of the site's forest area has been allocated to any bin + real(r8) :: sum_forest_bins_so_far_m2 = 0._r8 + + ! Out + real(r8), dimension(nlevedgeforest_tmp) :: area_in_edgeforest_bins + + ! The site's forest split evenly across the two bins closest to edge + real(r8), dimension(nlevedgeforest_tmp), target :: fraction_forest_in_each_bin + fraction_forest_in_each_bin = (/ 0.5_r8, 0.5_r8, 0._r8 /) + + call AssignPatchToBins(fraction_forest_in_each_bin, area_forest_patches, patch_area, & + nlevedgeforest_tmp, tol, sum_forest_bins_so_far_m2, area_in_edgeforest_bins) + + ! 0+100=100 area of the site's forest has been assigned after this call + @assertEqual(100._r8, sum_forest_bins_so_far_m2, tolerance=tol) + + ! Half of the patch's area was assigned to the first edge bin, half to the second + @assertEqual(50._r8, area_in_edgeforest_bins(1), tolerance=tol) + @assertEqual(50._r8, area_in_edgeforest_bins(2), tolerance=tol) + @assertEqual(0._r8, area_in_edgeforest_bins(3), tolerance=tol) + end subroutine test_AssignPatchToBins_03 + + @Test + subroutine test_AssignPatchToBins_04(this) + ! The site has just one forest patch. The site's forest is split evenly across the two bins farthest from edge. + class(TestEdgeForest), intent(inout) :: this + integer, parameter :: nlevedgeforest_tmp = 3 + + ! The site has 100 area of forest + real(r8) :: area_forest_patches = 100._r8 + + ! The size of this patch is 100 area + real(r8) :: patch_area = 100._r8 + + ! So far none of the site's forest area has been allocated to any bin + real(r8) :: sum_forest_bins_so_far_m2 = 0._r8 + + ! Out + real(r8), dimension(nlevedgeforest_tmp) :: area_in_edgeforest_bins + + ! The site's forest is split evenly across the two bins farthest from edge + real(r8), dimension(nlevedgeforest_tmp), target :: fraction_forest_in_each_bin + fraction_forest_in_each_bin = (/ 0._r8, 0.5_r8, 0.5_r8 /) + + call AssignPatchToBins(fraction_forest_in_each_bin, area_forest_patches, patch_area, & + nlevedgeforest_tmp, tol, sum_forest_bins_so_far_m2, area_in_edgeforest_bins) + + ! 0+100=100 area of the site's forest has been assigned after this call + @assertEqual(100._r8, sum_forest_bins_so_far_m2, tolerance=tol) + + ! Half of the patch's area was assigned to the second edge bin, half to the third + @assertEqual(0._r8, area_in_edgeforest_bins(1), tolerance=tol) + @assertEqual(50._r8, area_in_edgeforest_bins(2), tolerance=tol) + @assertEqual(50._r8, area_in_edgeforest_bins(3), tolerance=tol) + end subroutine test_AssignPatchToBins_04 + + @Test + subroutine test_AssignPatchToBins_05(this) + ! The site has multiple forest patches. Some of the site's forest area has already been assigned. The patch isn't big enough to take the rest of the site's forest area. + class(TestEdgeForest), intent(inout) :: this + integer, parameter :: nlevedgeforest_tmp = 3 + + ! The site has 100 area of forest + real(r8) :: area_forest_patches = 100._r8 + + ! The size of this patch is 50 area + real(r8) :: patch_area = 50._r8 + + ! So far 30 of the site's forest area has been allocated to any bin + real(r8) :: sum_forest_bins_so_far_m2 = 30._r8 + + ! Out + real(r8), dimension(nlevedgeforest_tmp) :: area_in_edgeforest_bins + + ! The site's forest is split evenly across the two bins farthest from edge + real(r8), dimension(nlevedgeforest_tmp), target :: fraction_forest_in_each_bin + fraction_forest_in_each_bin = (/ 0._r8, 0.5_r8, 0.5_r8 /) + + call AssignPatchToBins(fraction_forest_in_each_bin, area_forest_patches, patch_area, & + nlevedgeforest_tmp, tol, sum_forest_bins_so_far_m2, area_in_edgeforest_bins) + + ! 30+50=80 area of the site's forest has been assigned after this call + @assertEqual(80._r8, sum_forest_bins_so_far_m2, tolerance=tol) + + ! Some of the patch's area was assigned to the second edge bin, rest to the third + @assertEqual(0._r8, area_in_edgeforest_bins(1), tolerance=tol) + @assertEqual(20._r8, area_in_edgeforest_bins(2), tolerance=tol) + @assertEqual(30._r8, area_in_edgeforest_bins(3), tolerance=tol) + end subroutine test_AssignPatchToBins_05 + + end module test_EdgeForest diff --git a/testing/unit_tests.cfg b/testing/unit_tests.cfg index 426ff34874..bd4b1635ce 100644 --- a/testing/unit_tests.cfg +++ b/testing/unit_tests.cfg @@ -1,6 +1,12 @@ [fire_weather] test_dir = fates_fire_weather_utest +[ecotypes] +test_dir = fates_ecotypes_utest + +[edge_forest] +test_dir = fates_edge_forest_utest + [fire_fuel] test_dir = fates_fire_fuel_utest @@ -17,4 +23,4 @@ test_dir = fates_validate_cohorts_utest test_dir = fates_count_cohorts_utest [fire_equations] -test_dir = fates_fire_equations_utest \ No newline at end of file +test_dir = fates_fire_equations_utest diff --git a/testing/utils.py b/testing/utils.py index 55978cc7d0..186e3ab64a 100644 --- a/testing/utils.py +++ b/testing/utils.py @@ -53,6 +53,9 @@ def create_nc_from_cdl(cdl_path: str, run_dir: str) -> str: file_basename = os.path.basename(cdl_path).split(".")[-2] file_nc_name = f"{file_basename}.nc" + if not os.path.exists(run_dir): + os.makedirs(run_dir) + file_gen_command = ["ncgen -o", os.path.join(run_dir, file_nc_name), cdl_path] out = run_cmd_no_fail(" ".join(file_gen_command), combine_output=True) print(out) @@ -109,6 +112,7 @@ def config_to_dict(config_file: str) -> dict: # Define list of config file options that we expect to be paths options_that_are_paths = ["datm_file"] + options_that_are_bools = ["use_param_file"] config = configparser.ConfigParser() config.read(config_file) @@ -123,6 +127,10 @@ def config_to_dict(config_file: str) -> dict: if option in options_that_are_paths: value = get_abspath_from_config_file(value, config_file) + # If the option is one that we expect to be a boolean, convert it from string. + if option in options_that_are_bools: + value = str_to_bool(value) + # Save value to dictionary dictionary[section][option] = value @@ -172,6 +180,8 @@ def str_to_bool(val: str) -> bool: Returns: bool: True or False """ + if isinstance(val, bool): + return val if val.lower() in ("y", "yes", "t", "true", "on", "1"): return True if val.lower() in ("n", "no", "f", "false", "off", "0"): diff --git a/testing/utils_plotting.py b/testing/utils_plotting.py index df38b80617..5436cac353 100644 --- a/testing/utils_plotting.py +++ b/testing/utils_plotting.py @@ -3,6 +3,7 @@ import math import matplotlib.pyplot as plt +from matplotlib import colormaps def blank_plot( @@ -99,3 +100,11 @@ def get_color_palette(number: int) -> list: ] return colors[:number] + + +def sample_colormap(this_colormap: str, n_colors: int, i: int) -> tuple: + """ + Given a colormap and a number of desired colors evenly spaced along it, get the i'th color + """ + color = colormaps[this_colormap](i / (n_colors - 1)) + return color