diff --git a/biogeochem/FatesCohortMod.F90 b/biogeochem/FatesCohortMod.F90 index eb3f9ae65d..85c47175ac 100644 --- a/biogeochem/FatesCohortMod.F90 +++ b/biogeochem/FatesCohortMod.F90 @@ -289,6 +289,9 @@ module FatesCohortMod procedure :: Init procedure :: NanValues procedure :: ZeroValues + procedure :: InitRaw + procedure :: UpdateDerivedState + procedure :: CreateBare procedure :: Create procedure :: Copy procedure :: FreeMemory @@ -559,12 +562,14 @@ end subroutine ZeroValues !=========================================================================== - subroutine Create(this, prt, pft, nn, height, coage, dbh, status, & - ctrim, carea, clayer, crowndamage, spread, can_tlai, elongf_leaf, & - elongf_fnrt, elongf_stem) + !=========================================================================== + + subroutine InitRaw(this, prt, pft, nn, height, coage, dbh, status, & + ctrim, clayer, crowndamage, elongf_leaf, elongf_fnrt, elongf_stem) ! ! DESCRIPTION: - ! set up values for a newly created cohort + ! Stage 1 cohort initialization: sets up memory allocation and identity properties + ! without running allometry routines or querying parameter infrastructure. ! ARGUMENTS class(fates_cohort_type), intent(inout), target :: this ! cohort object @@ -576,24 +581,16 @@ subroutine Create(this, prt, pft, nn, height, coage, dbh, status, & real(r8), intent(in) :: nn ! number of individuals in cohort [/m2] real(r8), intent(in) :: height ! cohort height [m] real(r8), intent(in) :: coage ! cohort age [yr] - real(r8), intent(in) :: dbh ! cohort diameter at breat height [cm] - real(r8), intent(in) :: ctrim ! fraction of the maximum leaf biomass - real(r8), intent(in) :: spread ! how spread crowns are in horizontal space - real(r8), intent(in) :: carea ! area of cohort, for SP mode [m2] - real(r8), intent(in) :: can_tlai(:) ! patch-level total LAI of each canopy layer - real(r8), intent(in) :: elongf_leaf ! leaf elongation factor [fraction] - real(r8), intent(in) :: elongf_fnrt ! fine-root "elongation factor" [fraction] - real(r8), intent(in) :: elongf_stem ! stem "elongation factor" [fraction] + real(r8), intent(in) :: dbh ! cohort diameter at breast height [cm] + real(r8), intent(in) :: ctrim ! fraction of maximum leaf biomass + real(r8), intent(in), optional :: elongf_leaf ! leaf elongation factor [fraction] + real(r8), intent(in), optional :: elongf_fnrt ! fine-root "elongation factor" [fraction] + real(r8), intent(in), optional :: elongf_stem ! stem "elongation factor" [fraction] - ! LOCAL VARIABLES: - integer :: iage ! loop counter for leaf age classes - real(r8) :: leaf_c ! total leaf carbon [kgC] - real(r8) :: treesai ! stem area index within crown [m2/m2] - - ! initialize cohort + ! initialize cohort memory call this%Init(prt) - - ! set values + + ! set base values this%pft = pft this%crowndamage = crowndamage this%canopy_layer = clayer @@ -604,26 +601,53 @@ subroutine Create(this, prt, pft, nn, height, coage, dbh, status, & this%dbh = dbh this%coage = coage this%canopy_trim = ctrim - this%efleaf_coh = elongf_leaf - this%effnrt_coh = elongf_fnrt - this%efstem_coh = elongf_stem - - ! This routine may be called during restarts, and at this point in the call sequence - ! the actual cohort data is unknown, as this is really only used for allocation - ! In these cases, testing if things like biomass are reasonable is premature - ! However, in this part of the code, we will pass in nominal values for size, number and type + + ! Default elongation factors to 1.0 if optional arguments omitted + if (present(elongf_leaf)) then + this%efleaf_coh = elongf_leaf + else + this%efleaf_coh = 1.0_r8 + endif + + if (present(elongf_fnrt)) then + this%effnrt_coh = elongf_fnrt + else + this%effnrt_coh = 1.0_r8 + endif + + if (present(elongf_stem)) then + this%efstem_coh = elongf_stem + else + this%efstem_coh = 1.0_r8 + endif + if (this%dbh <= 0._r8 .or. this%n == 0._r8 .or. this%pft == 0) then - write(fates_log(),*) 'FATES: something is zero in cohort%Create', & + write(fates_log(),*) 'FATES: something is zero in cohort%InitRaw', & this%dbh, this%n, this%pft call endrun(msg=errMsg(sourcefile, __LINE__)) endif + end subroutine InitRaw + + !=========================================================================== + + subroutine UpdateDerivedState(this, spread, carea, can_tlai) + ! + ! DESCRIPTION: + ! Stage 2 cohort initialization: calculates all physical and derived allometric states. + + ! ARGUMENTS + class(fates_cohort_type), intent(inout) :: this + real(r8), intent(in) :: spread ! crown spread factor + real(r8), intent(in) :: carea ! cohort area for SP mode [m2] + real(r8), intent(in) :: can_tlai(:) ! patch-level canopy layer total LAI + + ! LOCAL VARIABLES: + real(r8) :: leaf_c ! total leaf carbon [kgC] + real(r8) :: treesai ! stem area index within crown [m2/m2] + ! Initialize the leaf to fineroot biomass ratio. - ! For C-only, this will stay constant, for nutrient-enabled this will be - ! dynamic. In both cases, new cohorts are initialized with the minimum. - ! This works in the nutrient enabled case because cohorts are also - ! initialized with full stores, which match with minimum fineroot biomass - this%l2fr = prt_params%allom_l2fr(pft) + this%l2fr = prt_params%allom_l2fr(this%pft) if (hlm_parteh_mode == carbon_nitrogen_phosphorus) then this%cx_int = 0._r8 ! Assume balanced N,P/C stores ie log(1) = 0 @@ -632,48 +656,135 @@ subroutine Create(this, prt, pft, nn, height, coage, dbh, status, & this%cnp_limiter = 0 ! Assume limitations are unknown end if - ! This sets things like vcmax25top, that depend on the leaf age fractions - ! (which are defined by PARTEH) + ! Set biophysical rates depending on leaf age fractions call this%UpdateCohortBioPhysRates() - ! calculate size classes + ! Calculate size classes call sizetype_class_index(this%dbh, this%pft, this%size_class, & this%size_by_pft_class) - ! If cohort age tracking is off we call this here once, just so everything - ! is in the first bin. This makes it easier to copy and terminate cohorts - ! later. - ! We don't need to update this ever if cohort age tracking is off + ! Calculate age classes call coagetype_class_index(this%coage, this%pft, this%coage_class, & this%coage_by_pft_class) - ! asssign or calculate canopy extent and depth + ! Assign or calculate canopy extent and depth if (hlm_use_sp .eq. ifalse) then call carea_allom(this%dbh, this%n, spread, this%pft, this%crowndamage, & this%c_area) else - ! set this from previously precision-controlled value in SP mode this%c_area = carea endif - ! Query PARTEH for the leaf carbon [kg] - leaf_c = this%prt%GetState(leaf_organ, carbon12_element) + ! Query PARTEH for leaf carbon [kg] + if (associated(this%prt)) then + leaf_c = this%prt%GetState(leaf_organ, carbon12_element) + else + leaf_c = 0._r8 + endif - call tree_lai_sai(leaf_c, this%pft, this%c_area, this%n, & - this%canopy_layer, can_tlai, this%vcmax25top, this%dbh, this%crowndamage, & + call tree_lai_sai(leaf_c, this%pft, this%c_area, this%n, & + this%canopy_layer, can_tlai, this%vcmax25top, this%dbh, this%crowndamage, & this%canopy_trim, this%efstem_coh, 2, this%treelai, treesai) if (hlm_use_sp .eq. ifalse) then this%treesai = treesai end if - call this%InitPRTBoundaryConditions() + end subroutine UpdateDerivedState + + !=========================================================================== + + subroutine Create(this, prt, pft, nn, height, coage, dbh, status, & + ctrim, carea, clayer, crowndamage, spread, can_tlai, elongf_leaf, & + elongf_fnrt, elongf_stem) + ! + ! DESCRIPTION: + ! Production cohort constructor shim executing two-stage creation. + + ! ARGUMENTS + class(fates_cohort_type), intent(inout), target :: this ! cohort object + class(prt_vartypes), intent(inout), pointer :: prt ! The allocated PARTEH object + integer, intent(in) :: pft ! cohort Plant Functional Type + integer, intent(in) :: crowndamage ! cohort damage class + integer, intent(in) :: clayer ! canopy status of cohort [canopy/understory] + integer, intent(in) :: status ! growth status of cohort [leaves on/off] + real(r8), intent(in) :: nn ! number of individuals in cohort [/m2] + real(r8), intent(in) :: height ! cohort height [m] + real(r8), intent(in) :: coage ! cohort age [yr] + real(r8), intent(in) :: dbh ! cohort diameter at breast height [cm] + real(r8), intent(in) :: ctrim ! fraction of maximum leaf biomass + real(r8), intent(in) :: spread ! how spread crowns are in horizontal space + real(r8), intent(in) :: carea ! area of cohort, for SP mode [m2] + real(r8), intent(in) :: can_tlai(:) ! patch-level total LAI of each canopy layer + real(r8), intent(in) :: elongf_leaf ! leaf elongation factor [fraction] + real(r8), intent(in) :: elongf_fnrt ! fine-root "elongation factor" [fraction] + real(r8), intent(in) :: elongf_stem ! stem "elongation factor" [fraction] + + call this%InitRaw(prt, pft, nn, height, coage, dbh, status, & + ctrim, clayer, crowndamage, elongf_leaf, elongf_fnrt, elongf_stem) + + call this%UpdateDerivedState(spread, carea, can_tlai) + end subroutine Create !=========================================================================== + subroutine CreateBare(this, prt, pft, nn, height, coage, dbh, status, & + ctrim, carea, clayer, crowndamage, treelai, treesai, vcmax25top) + ! + ! DESCRIPTION: + ! Lightweight constructor for unit tests and synthetic fixtures that sets up a valid + ! cohort without calling allometry routines or querying global parameter infrastructure. + + ! ARGUMENTS + class(fates_cohort_type), intent(inout), target :: this ! cohort object + class(prt_vartypes), intent(inout), pointer :: prt ! The allocated PARTEH object + integer, intent(in) :: pft ! cohort Plant Functional Type + real(r8), intent(in) :: nn ! number of individuals in cohort [/m2] + real(r8), intent(in) :: height ! cohort height [m] + real(r8), intent(in), optional :: coage ! cohort age [yr] + real(r8), intent(in), optional :: dbh ! cohort diameter at breast height [cm] + integer, intent(in), optional :: status ! growth status [leaves on/off] + real(r8), intent(in), optional :: ctrim ! fraction of maximum leaf biomass + real(r8), intent(in), optional :: carea ! area of cohort [m2] + integer, intent(in), optional :: clayer ! canopy layer status + integer, intent(in), optional :: crowndamage ! crown damage class + real(r8), intent(in), optional :: treelai ! leaf area index [m2/m2] + real(r8), intent(in), optional :: treesai ! stem area index [m2/m2] + real(r8), intent(in), optional :: vcmax25top ! max carboxylation rate at top + + real(r8) :: local_coage, local_dbh, local_ctrim, local_carea + integer :: local_status, local_clayer, local_crowndamage + + ! Default nominal fallback values for optional unit test arguments + local_coage = 0.0_r8; if (present(coage)) local_coage = coage + ! Default nominal dbh = 1.0 cm for synthetic testing + local_dbh = 1.0_r8; if (present(dbh)) local_dbh = dbh + ! Default status = 2 (leaves on) for synthetic testing + local_status = 2; if (present(status)) local_status = status + ! Default canopy trim = 1.0 (untrimmed) for synthetic testing + local_ctrim = 1.0_r8; if (present(ctrim)) local_ctrim = ctrim + ! Default crown area = 1.0 m2 for synthetic testing + local_carea = 1.0_r8; if (present(carea)) local_carea = carea + ! Default canopy layer = 1 (upper canopy) for synthetic testing + local_clayer = 1; if (present(clayer)) local_clayer = clayer + ! Default crown damage = 1 (undamaged) for synthetic testing + local_crowndamage = 1; if (present(crowndamage)) local_crowndamage = crowndamage + + call this%InitRaw(prt, pft, nn, height, local_coage, local_dbh, local_status, & + local_ctrim, local_clayer, local_crowndamage) + + if (present(carea)) this%c_area = carea + if (present(treelai)) this%treelai = treelai + if (present(treesai)) this%treesai = treesai + if (present(vcmax25top)) this%vcmax25top = vcmax25top + + end subroutine CreateBare + + !=========================================================================== + subroutine Copy(this, copyCohort) ! ! DESCRIPTION: diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index b4f93ba5c9..19941bea87 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -79,13 +79,8 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) ! Accumulate fluxes from hourly to daily values. ! _tstep fluxes are KgC/indiv/timestep _acc are KgC/indiv/day - ccohort%gpp_acc = ccohort%gpp_acc + ccohort%gpp_tstep - ccohort%resp_m_acc = ccohort%resp_m_acc + ccohort%resp_m_tstep - - ccohort%sym_nfix_daily = ccohort%sym_nfix_daily + ccohort%sym_nfix_tstep - ! weighted mean of D13C by gpp - if((ccohort%gpp_acc + ccohort%gpp_tstep) .eq. 0.0_r8) then + if((ccohort%gpp_acc + ccohort%gpp_tstep) <= 0.0_r8) then ccohort%c13disc_acc = 0.0_r8 else ccohort%c13disc_acc = ((ccohort%c13disc_acc * ccohort%gpp_acc) + & @@ -93,6 +88,11 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) (ccohort%gpp_acc + ccohort%gpp_tstep) endif + ccohort%gpp_acc = ccohort%gpp_acc + ccohort%gpp_tstep + ccohort%resp_m_acc = ccohort%resp_m_acc + ccohort%resp_m_tstep + + ccohort%sym_nfix_daily = ccohort%sym_nfix_daily + ccohort%sym_nfix_tstep + do iv=1,ccohort%nv if(ccohort%year_net_uptake(iv) == 999._r8)then ! note that there were leaves in this layer this year. ccohort%year_net_uptake(iv) = 0._r8 diff --git a/testing/CMakeLists.txt b/testing/CMakeLists.txt index 99814f9e75..4aceffb763 100644 --- a/testing/CMakeLists.txt +++ b/testing/CMakeLists.txt @@ -17,3 +17,4 @@ add_subdirectory(tests/unit/count_cohorts_test fates_count_cohorts_utest) add_subdirectory(tests/unit/fire_equations_test fates_fire_equations_utest) add_subdirectory(tests/unit/quadratic_roots_test fates_quadratic_roots_utest) add_subdirectory(tests/unit/great_circle_test fates_great_circle_utest) +add_subdirectory(tests/unit/ed_accumulate_fluxes_test ed_accumulate_fluxes_test) diff --git a/testing/config/unit.cfg b/testing/config/unit.cfg index ce73c6fdbb..4f2140715d 100644 --- a/testing/config/unit.cfg +++ b/testing/config/unit.cfg @@ -25,3 +25,5 @@ test_dir = fates_quadratic_roots_utest [great_circle] test_dir = fates_great_circle_utest +[ed_accumulate_fluxes] +test_dir = ed_accumulate_fluxes_test diff --git a/testing/tests/fortran_shr/FatesFactoryMod.F90 b/testing/tests/fortran_shr/FatesFactoryMod.F90 index 5a3dfb6cbd..348cf7757d 100644 --- a/testing/tests/fortran_shr/FatesFactoryMod.F90 +++ b/testing/tests/fortran_shr/FatesFactoryMod.F90 @@ -13,7 +13,8 @@ module FatesFactoryMod use FatesGlobals, only : endrun => fates_endrun use FatesCohortMod, only : fates_cohort_type use FatesPatchMod, only : fates_patch_type - use EDTypesMod, only : init_spread_inventory + use EDTypesMod, only : init_spread_inventory, ed_site_type + use FatesInterfaceTypesMod, only : bc_in_type, bc_out_type use FatesRadiationMemMod, only : num_swb use EDParamsMod, only : vai_top_bin_width use EDParamsMod, only : vai_width_increase_factor @@ -390,10 +391,10 @@ subroutine CohortFactory(cohort, pft, can_lai, dbh, number, crown_damage, status ! initialize the PRT object call PRTFactory(prt, pft, c_struct, c_leaf, c_fnrt, c_sapw, c_store) - ! create the cohort - call cohort%Create(prt, pft, number_local, height, age_local, dbh_local, & - status_local, canopy_trim_local, can_area, canopy_layer_local, crown_damage_local, & - init_spread_inventory, can_lai, elongf_leaf, elongf_fnrt, elongf_stem) + ! create the cohort using lightweight CreateBare constructor + call cohort%CreateBare(prt=prt, pft=pft, nn=number_local, height=height, coage=age_local, & + dbh=dbh_local, status=status_local, ctrim=canopy_trim_local, carea=can_area, & + clayer=canopy_layer_local, crowndamage=crown_damage_local) end subroutine CohortFactory @@ -514,6 +515,7 @@ subroutine CreateTestPatchList(patch, heights, dbhs) ! LOCALS: type(fates_cohort_type), pointer :: cohort, next_cohort ! cohort objects + class(prt_vartypes), pointer :: null_prt => null() ! null pointer fallback for CreateBare integer :: num_cohorts ! number of cohorts to add to list integer :: i ! looping index @@ -529,15 +531,22 @@ subroutine CreateTestPatchList(patch, heights, dbhs) ! initialize first cohort allocate(cohort) - cohort%height = heights(1) - if (present(dbhs)) cohort%dbh = dbhs(1) + ! Default nominal PFT = 1 and density = 0.1 /m2 for synthetic test patch list + if (present(dbhs)) then + call cohort%CreateBare(prt=null_prt, pft=1, nn=0.1_r8, height=heights(1), dbh=dbhs(1)) + else + call cohort%CreateBare(prt=null_prt, pft=1, nn=0.1_r8, height=heights(1)) + endif patch%shortest => cohort ! initialize the rest of the cohorts do i = 2, num_cohorts allocate(next_cohort) - next_cohort%height = heights(i) - if (present(dbhs)) next_cohort%dbh = dbhs(i) + if (present(dbhs)) then + call next_cohort%CreateBare(prt=null_prt, pft=1, nn=0.1_r8, height=heights(i), dbh=dbhs(i)) + else + call next_cohort%CreateBare(prt=null_prt, pft=1, nn=0.1_r8, height=heights(i)) + endif cohort%taller => next_cohort next_cohort%shorter => cohort cohort => next_cohort @@ -545,5 +554,86 @@ subroutine CreateTestPatchList(patch, heights, dbhs) patch%tallest => cohort end subroutine CreateTestPatchList - + + !--------------------------------------------------------------------------------------- + ! DESCRIPTION: + ! Helper for setting transient timestep flux rates and accumulators on a cohort for unit testing. + ! Parameter choices represent nominal physical flux fixtures (e.g. gpp_tstep, resp_m_tstep) + ! to verify rate accumulation without invoking numerical solver machinery. + subroutine SetTestFluxes(cohort, gpp_acc, gpp_tstep, resp_m_acc, resp_m_tstep, & + sym_nfix_daily, sym_nfix_tstep, c13disc_acc, c13disc_clm, & + year_net_uptake, ts_net_uptake, nv) + type(fates_cohort_type), pointer, intent(inout) :: cohort + real(r8), optional, intent(in) :: gpp_acc ! baseline accumulated GPP [kgC/m2] + real(r8), optional, intent(in) :: gpp_tstep ! current timestep GPP flux rate [kgC/m2/s] + real(r8), optional, intent(in) :: resp_m_acc ! baseline accumulated maintenance respiration [kgC/m2] + real(r8), optional, intent(in) :: resp_m_tstep ! current timestep maintenance respiration [kgC/m2/s] + real(r8), optional, intent(in) :: sym_nfix_daily ! daily accumulated symbiotic N-fixation [gN/m2] + real(r8), optional, intent(in) :: sym_nfix_tstep ! current timestep symbiotic N-fixation [gN/m2/s] + real(r8), optional, intent(in) :: c13disc_acc ! baseline C13 discrimination accumulator [per mil] + real(r8), optional, intent(in) :: c13disc_clm ! canopy-level C13 discrimination rate [per mil] + real(r8), optional, intent(in) :: year_net_uptake ! annual net uptake rate fixture [kgC/m2/yr] + real(r8), optional, intent(in) :: ts_net_uptake ! timestep net uptake rate fixture [kgC/m2/s] + integer, optional, intent(in) :: nv ! number of active canopy leaf layers + + if (present(gpp_acc)) cohort%gpp_acc = gpp_acc + if (present(gpp_tstep)) cohort%gpp_tstep = gpp_tstep + if (present(resp_m_acc)) cohort%resp_m_acc = resp_m_acc + if (present(resp_m_tstep)) cohort%resp_m_tstep = resp_m_tstep + if (present(sym_nfix_daily)) cohort%sym_nfix_daily = sym_nfix_daily + if (present(sym_nfix_tstep)) cohort%sym_nfix_tstep = sym_nfix_tstep + if (present(c13disc_acc)) cohort%c13disc_acc = c13disc_acc + if (present(c13disc_clm)) cohort%c13disc_clm = c13disc_clm + if (present(nv)) cohort%nv = nv + + if (present(year_net_uptake)) then + cohort%year_net_uptake(1) = year_net_uptake + end if + if (present(ts_net_uptake)) then + cohort%ts_net_uptake(1) = ts_net_uptake + end if + + end subroutine SetTestFluxes + + !--------------------------------------------------------------------------------------- + ! DESCRIPTION: + ! High-level factory building a site containing a patch with a cohort linked list. + ! Magic constants: + ! - local_pft = 1: default Plant Functional Type index for synthetic test fixtures + ! - filter_photo_pa(1) = 3: FATES boundary condition flag indicating active photosynthesis + subroutine CreateTestSite(sites, patch, cohort, bc_in, heights, dbhs, pft) + type(ed_site_type), target, intent(inout) :: sites(:) + type(fates_patch_type), pointer, intent(out) :: patch + type(fates_cohort_type), pointer, intent(out) :: cohort + type(bc_in_type), intent(inout) :: bc_in(:) + real(r8), intent(in) :: heights(:) + real(r8), optional, intent(in) :: dbhs(:) + integer, optional, intent(in) :: pft + + integer :: local_pft + + ! Default to PFT index 1 when optional pft argument is omitted + local_pft = 1; if (present(pft)) local_pft = pft + + allocate(patch) + ! Create patch and cohort linked list via CreateTestPatchList + if (present(dbhs)) then + call CreateTestPatchList(patch, heights, dbhs=dbhs) + else + call CreateTestPatchList(patch, heights) + end if + + ! Set default single-patch topology and label + patch%patchno = 1 + patch%nocomp_pft_label = local_pft + sites(1)%oldest_patch => patch + cohort => patch%shortest + + if (allocated(bc_in(1)%filter_photo_pa)) deallocate(bc_in(1)%filter_photo_pa) + allocate(bc_in(1)%filter_photo_pa(1)) + ! FATES photosynthesis filter flag 3 indicates active photosynthesis on vegetated patch + bc_in(1)%filter_photo_pa(1) = 3 + + end subroutine CreateTestSite + end module FatesFactoryMod diff --git a/testing/tests/unit/ed_accumulate_fluxes_test/CMakeLists.txt b/testing/tests/unit/ed_accumulate_fluxes_test/CMakeLists.txt new file mode 100644 index 0000000000..d74dbff971 --- /dev/null +++ b/testing/tests/unit/ed_accumulate_fluxes_test/CMakeLists.txt @@ -0,0 +1,5 @@ +set(pfunit_sources test_EDAccumulateFluxes.pf) + +add_pfunit_ctest(EDAccumulateFluxes + TEST_SOURCES "${pfunit_sources}" + LINK_LIBRARIES fates csm_share) diff --git a/testing/tests/unit/ed_accumulate_fluxes_test/test_EDAccumulateFluxes.pf b/testing/tests/unit/ed_accumulate_fluxes_test/test_EDAccumulateFluxes.pf new file mode 100644 index 0000000000..8998849819 --- /dev/null +++ b/testing/tests/unit/ed_accumulate_fluxes_test/test_EDAccumulateFluxes.pf @@ -0,0 +1,681 @@ +!> Summary: Unit tests for EDAccumulateFluxesMod covering normal, edge-case, multi-cohort, multi-patch, multi-site, multi-PFT, zero-layer, and subnormal floating point behaviors. +!> This module contains fortified unit tests for the EDAccumulateFluxesMod production module, +!> which aggregates and tracks plant carbon, water, and nutrient fluxes (e.g., GPP, respiration) over time. +module test_EDAccumulateFluxes + use FatesConstantsMod, only : r8 => fates_r8, nocomp_bareground + use EDTypesMod, only : ed_site_type + use FatesPatchMod, only : fates_patch_type + use FatesCohortMod, only : fates_cohort_type + use FatesInterfaceTypesMod, only : bc_in_type, bc_out_type + use EDAccumulateFluxesMod, only : AccumulateFluxes_ED + use FatesFactoryMod, only : CreateTestSite, SetTestFluxes + use funit + implicit none + + @TestCase + type, extends(TestCase) :: TestEDAccumulateFluxes + contains + end type TestEDAccumulateFluxes + + contains + + ! Test normal accumulation of GPP, respiration, N-fixation, C13 discrimination, and net uptake + @Test + subroutine test_AccumulateFluxes_ED_Normal(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for isolating single-location flux integration + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch + type(fates_cohort_type), pointer :: cohort + real(r8) :: dt_time + + call CreateTestSite(sites, patch, cohort, bc_in, heights=[10.0_r8]) + + call SetTestFluxes(cohort, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=1) + + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Comprehensive assertions for accumulated variables + @assertEqual(1.5_r8, cohort%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.3_r8, cohort%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(0.07_r8, cohort%sym_nfix_daily, tolerance=1.e-8_r8) + ! c13disc_acc = (0.1*1.0 + 0.2*0.5)/(1.0+0.5) = 0.2 / 1.5 = 0.13333333_r8 + @assertEqual(0.13333333333333333_r8, cohort%c13disc_acc, tolerance=1.e-8_r8) + @assertEqual(0.7_r8, cohort%year_net_uptake(1), tolerance=1.e-8_r8) + + ! Assert input timestep fluxes remain unmutated + @assertEqual(0.5_r8, cohort%gpp_tstep, tolerance=1.e-8_r8) + @assertEqual(0.1_r8, cohort%resp_m_tstep, tolerance=1.e-8_r8) + @assertEqual(0.02_r8, cohort%sym_nfix_tstep, tolerance=1.e-8_r8) + + deallocate(cohort) + deallocate(patch) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_Normal + + ! Test divide-by-zero protection in C13 discrimination when total GPP accumulates to zero + @Test + subroutine test_AccumulateFluxes_ED_c13ZeroDenom(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for isolating zero-GPP denominator branch + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch + type(fates_cohort_type), pointer :: cohort + real(r8) :: dt_time + + call CreateTestSite(sites, patch, cohort, bc_in, heights=[10.0_r8]) + + call SetTestFluxes(cohort, gpp_acc=1.0_r8, gpp_tstep=-1.0_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=1) + + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + @assertEqual(0.0_r8, cohort%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.0_r8, cohort%c13disc_acc, tolerance=1.e-8_r8) + ! Ensure other fluxes still accumulate correctly under zero GPP + @assertEqual(0.3_r8, cohort%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(0.07_r8, cohort%sym_nfix_daily, tolerance=1.e-8_r8) + @assertEqual(0.7_r8, cohort%year_net_uptake(1), tolerance=1.e-8_r8) + + deallocate(cohort) + deallocate(patch) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_c13ZeroDenom + + ! Test reset behavior when year_net_uptake contains the 999.0 uninitialized sentinel value + @Test + subroutine test_AccumulateFluxes_ED_YearNetUptake999(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for sentinel net uptake reset validation + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch + type(fates_cohort_type), pointer :: cohort + real(r8) :: dt_time + + call CreateTestSite(sites, patch, cohort, bc_in, heights=[10.0_r8]) + + call SetTestFluxes(cohort, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.2_r8, sym_nfix_tstep=0.2_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=999.0_r8, & + ts_net_uptake=0.5_r8, nv=1) + + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! year_net_uptake(1) should reset from sentinel 999.0 to 0, then add 0.5 + @assertEqual(0.5_r8, cohort%year_net_uptake(1), tolerance=1.e-8_r8) + @assertEqual(1.5_r8, cohort%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.3_r8, cohort%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(0.4_r8, cohort%sym_nfix_daily, tolerance=1.e-8_r8) + + deallocate(cohort) + deallocate(patch) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_YearNetUptake999 + + ! Test skip conditions for bare ground patches and inactive photosynthesis filter values + @Test + subroutine test_AccumulateFluxes_ED_SkipConditions(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for verifying patch skip conditions + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch_bare, patch_skip + type(fates_cohort_type), pointer :: cohort_bare, cohort_skip + real(r8) :: dt_time + + allocate(patch_bare) + allocate(cohort_bare) + call cohort_bare%CreateBare(pft=1, nn=0.1_r8, height=10.0_r8) + allocate(patch_skip) + allocate(cohort_skip) + call cohort_skip%CreateBare(pft=1, nn=0.1_r8, height=10.0_r8) + + ! Patch 1: Bare ground patch label + patch_bare%patchno = 1 + patch_bare%nocomp_pft_label = nocomp_bareground + call SetTestFluxes(cohort_bare, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.2_r8, sym_nfix_tstep=0.2_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=1) + cohort_bare%taller => null() + patch_bare%shortest => cohort_bare + + ! Patch 2: filter_photo_pa != 3 (inactive photosynthesis filter) + patch_skip%patchno = 2 + patch_skip%nocomp_pft_label = 1 + call SetTestFluxes(cohort_skip, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.2_r8, sym_nfix_tstep=0.2_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=1) + cohort_skip%taller => null() + patch_skip%shortest => cohort_skip + + patch_bare%younger => patch_skip + patch_skip%younger => null() + sites(1)%oldest_patch => patch_bare + + allocate(bc_in(1)%filter_photo_pa(2)) + ! FATES photosynthesis filter flag 3 indicates active photosynthesis on vegetated patch + bc_in(1)%filter_photo_pa(1) = 3 + ! FATES photosynthesis filter flag 1 indicates inactive photosynthesis skip condition + bc_in(1)%filter_photo_pa(2) = 1 + + ! Timestep duration of 1 hour (3600s) + 0.2s offset (unreferenced dummy argument in AccumulateFluxes_ED) + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Verify NO fields were mutated for bare ground + @assertEqual(1.0_r8, cohort_bare%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.2_r8, cohort_bare%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(0.2_r8, cohort_bare%sym_nfix_daily, tolerance=1.e-8_r8) + @assertEqual(0.1_r8, cohort_bare%c13disc_acc, tolerance=1.e-8_r8) + @assertEqual(0.2_r8, cohort_bare%year_net_uptake(1), tolerance=1.e-8_r8) + + ! Verify NO fields were mutated for inactive filter patch + @assertEqual(1.0_r8, cohort_skip%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.2_r8, cohort_skip%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(0.2_r8, cohort_skip%sym_nfix_daily, tolerance=1.e-8_r8) + @assertEqual(0.1_r8, cohort_skip%c13disc_acc, tolerance=1.e-8_r8) + @assertEqual(0.2_r8, cohort_skip%year_net_uptake(1), tolerance=1.e-8_r8) + + deallocate(cohort_skip) + deallocate(patch_skip) + deallocate(cohort_bare) + deallocate(patch_bare) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_SkipConditions + + ! Test cohort linked list traversal for multiple cohorts (understory to canopy) on a single patch + @Test + subroutine test_AccumulateFluxes_ED_MultiCohort(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for cohort linked list traversal + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch + type(fates_cohort_type), pointer :: cohort_under, cohort_canopy + real(r8) :: dt_time + + call CreateTestSite(sites, patch, cohort_under, bc_in, heights=[5.0_r8, 15.0_r8]) + cohort_canopy => cohort_under%taller + + call SetTestFluxes(cohort_under, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=1) + + call SetTestFluxes(cohort_canopy, gpp_acc=2.0_r8, gpp_tstep=1.0_r8, resp_m_acc=0.4_r8, & + resp_m_tstep=0.2_r8, sym_nfix_daily=0.10_r8, sym_nfix_tstep=0.04_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.4_r8, & + ts_net_uptake=1.0_r8, nv=1) + + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Verify full accumulation on understory cohort + @assertEqual(1.5_r8, cohort_under%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.3_r8, cohort_under%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(0.07_r8, cohort_under%sym_nfix_daily, tolerance=1.e-8_r8) + @assertEqual(0.7_r8, cohort_under%year_net_uptake(1), tolerance=1.e-8_r8) + + ! Verify full accumulation on canopy cohort + @assertEqual(3.0_r8, cohort_canopy%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.6_r8, cohort_canopy%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(0.14_r8, cohort_canopy%sym_nfix_daily, tolerance=1.e-8_r8) + @assertEqual(1.4_r8, cohort_canopy%year_net_uptake(1), tolerance=1.e-8_r8) + + deallocate(cohort_canopy) + deallocate(cohort_under) + deallocate(patch) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_MultiCohort + + ! Test vertical canopy layer loop accumulation across multiple leaf layers (nv > 1) + @Test + subroutine test_AccumulateFluxes_ED_MultiLayerCanopy(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for multi-layer canopy testing + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch + type(fates_cohort_type), pointer :: cohort + real(r8) :: dt_time + + call CreateTestSite(sites, patch, cohort, bc_in, heights=[10.0_r8]) + + call SetTestFluxes(cohort, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, nv=3) + cohort%year_net_uptake(1) = 0.1_r8 + cohort%year_net_uptake(2) = 999.0_r8 + cohort%year_net_uptake(3) = 0.3_r8 + + cohort%ts_net_uptake(1) = 0.5_r8 + cohort%ts_net_uptake(2) = 0.4_r8 + cohort%ts_net_uptake(3) = 0.2_r8 + + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Verify layer 1, layer 2 (reset from 999.0), and layer 3 all accumulated + @assertEqual(0.6_r8, cohort%year_net_uptake(1), tolerance=1.e-8_r8) + @assertEqual(0.4_r8, cohort%year_net_uptake(2), tolerance=1.e-8_r8) + @assertEqual(0.5_r8, cohort%year_net_uptake(3), tolerance=1.e-8_r8) + @assertEqual(1.5_r8, cohort%gpp_acc, tolerance=1.e-8_r8) + + deallocate(cohort) + deallocate(patch) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_MultiLayerCanopy + + ! Test patch age-distribution linked list traversal for multiple active vegetated patches + @Test + subroutine test_AccumulateFluxes_ED_MultiPatchActive(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for patch age sequence traversal + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch_old, patch_young + type(fates_cohort_type), pointer :: cohort_old, cohort_young + real(r8) :: dt_time + + allocate(patch_old) + allocate(cohort_old) + call cohort_old%CreateBare(pft=1, nn=0.1_r8, height=15.0_r8) + allocate(patch_young) + allocate(cohort_young) + call cohort_young%CreateBare(pft=1, nn=0.1_r8, height=5.0_r8) + + patch_old%patchno = 1 + patch_old%nocomp_pft_label = 1 + call SetTestFluxes(cohort_old, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=1) + cohort_old%taller => null() + patch_old%shortest => cohort_old + + patch_young%patchno = 2 + patch_young%nocomp_pft_label = 1 + call SetTestFluxes(cohort_young, gpp_acc=2.0_r8, gpp_tstep=1.0_r8, resp_m_acc=0.4_r8, & + resp_m_tstep=0.2_r8, sym_nfix_daily=0.10_r8, sym_nfix_tstep=0.04_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.4_r8, & + ts_net_uptake=1.0_r8, nv=1) + cohort_young%taller => null() + patch_young%shortest => cohort_young + + patch_old%younger => patch_young + patch_young%younger => null() + sites(1)%oldest_patch => patch_old + + allocate(bc_in(1)%filter_photo_pa(2)) + ! FATES photosynthesis filter flag 3 indicates active photosynthesis on vegetated patch + bc_in(1)%filter_photo_pa(1) = 3 + bc_in(1)%filter_photo_pa(2) = 3 + ! Timestep duration of 1 hour (3600s) + 0.2s offset (unreferenced dummy argument in AccumulateFluxes_ED) + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Verify accumulation occurred on both oldest and younger patches + @assertEqual(1.5_r8, cohort_old%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.3_r8, cohort_old%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(3.0_r8, cohort_young%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.6_r8, cohort_young%resp_m_acc, tolerance=1.e-8_r8) + + deallocate(cohort_young) + deallocate(patch_young) + deallocate(cohort_old) + deallocate(patch_old) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_MultiPatchActive + + ! Test near-zero floating point denominator handling in C13 discrimination calculation + @Test + subroutine test_AccumulateFluxes_ED_c13NearZeroDenom(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for near-zero subnormal float accumulation + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch + type(fates_cohort_type), pointer :: cohort + real(r8) :: dt_time + + call CreateTestSite(sites, patch, cohort, bc_in, heights=[10.0_r8]) + + call SetTestFluxes(cohort, gpp_acc=1.0e-16_r8, gpp_tstep=-1.0e-16_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=1) + + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! c13disc_acc should fall back cleanly to 0.0_r8 without NaN or overflow + @assertEqual(0.0_r8, cohort%c13disc_acc, tolerance=1.e-8_r8) + @assertEqual(0.3_r8, cohort%resp_m_acc, tolerance=1.e-8_r8) + + deallocate(cohort) + deallocate(patch) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_c13NearZeroDenom + + ! Test multi-site iteration (nsites > 1) with fortified assertions for all accumulated state fields + @Test + subroutine test_AccumulateFluxes_ED_MultiSite(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Two site configuration for testing multi-site outer loop iteration + integer, parameter :: nsites = 2 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch1, patch2 + type(fates_cohort_type), pointer :: cohort1, cohort2 + real(r8) :: dt_time + + ! Site 1 + call CreateTestSite(sites(1:1), patch1, cohort1, bc_in(1:1), heights=[10.0_r8]) + call SetTestFluxes(cohort1, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=1) + + ! Site 2 + call CreateTestSite(sites(2:2), patch2, cohort2, bc_in(2:2), heights=[10.0_r8]) + call SetTestFluxes(cohort2, gpp_acc=5.0_r8, gpp_tstep=2.0_r8, resp_m_acc=1.0_r8, & + resp_m_tstep=0.4_r8, sym_nfix_daily=0.1_r8, sym_nfix_tstep=0.05_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=1.0_r8, & + ts_net_uptake=2.0_r8, nv=1) + + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Fortified assertions: Verify ALL state variables accumulated across site 1 and site 2 + @assertEqual(1.5_r8, cohort1%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.3_r8, cohort1%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(0.07_r8, cohort1%sym_nfix_daily, tolerance=1.e-8_r8) + @assertEqual(0.13333333333333333_r8, cohort1%c13disc_acc, tolerance=1.e-8_r8) + @assertEqual(0.7_r8, cohort1%year_net_uptake(1), tolerance=1.e-8_r8) + + @assertEqual(7.0_r8, cohort2%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(1.4_r8, cohort2%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(0.15_r8, cohort2%sym_nfix_daily, tolerance=1.e-8_r8) + @assertEqual(0.12857142857142856_r8, cohort2%c13disc_acc, tolerance=1.e-8_r8) + @assertEqual(3.0_r8, cohort2%year_net_uptake(1), tolerance=1.e-8_r8) + + deallocate(cohort1) + deallocate(patch1) + deallocate(cohort2) + deallocate(patch2) + deallocate(bc_in(1)%filter_photo_pa) + deallocate(bc_in(2)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_MultiSite + + ! Test empty cohort linked list (patch with shortest => null) with fortified state assertions + @Test + subroutine test_AccumulateFluxes_ED_NullCohort(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for null cohort list traversal + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch + real(r8) :: dt_time + + allocate(patch) + patch%patchno = 1 + patch%nocomp_pft_label = 1 + patch%shortest => null() ! empty cohort list + patch%younger => null() + sites(1)%oldest_patch => patch + + allocate(bc_in(1)%filter_photo_pa(1)) + ! FATES photosynthesis filter flag 3 indicates active photosynthesis on vegetated patch + bc_in(1)%filter_photo_pa(1) = 3 + ! Timestep duration of 1 hour (3600s) + 0.2s offset (unreferenced dummy argument in AccumulateFluxes_ED) + dt_time = 3600.2_r8 + + ! Execution must complete cleanly without crash/segfault + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Fortified assertions: Verify patch structure remains intact and cohort list stays null + @assertTrue(associated(sites(1)%oldest_patch)) + @assertFalse(associated(sites(1)%oldest_patch%shortest)) + @assertEqual(1, sites(1)%oldest_patch%patchno) + deallocate(patch) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_NullCohort + + ! Test protection when total GPP accumulation sums to a negative value with fortified assertions + @Test + subroutine test_AccumulateFluxes_ED_c13NegativeDenom(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for negative GPP accumulation branch + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch + type(fates_cohort_type), pointer :: cohort + real(r8) :: dt_time + + call CreateTestSite(sites, patch, cohort, bc_in, heights=[10.0_r8]) + + call SetTestFluxes(cohort, gpp_acc=0.0_r8, gpp_tstep=-0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=1) + + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Fortified assertions: Verify c13disc_acc falls back to 0.0_r8 while other state fields accumulate + @assertEqual(0.0_r8, cohort%c13disc_acc, tolerance=1.e-8_r8) + @assertEqual(-0.5_r8, cohort%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.3_r8, cohort%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(0.07_r8, cohort%sym_nfix_daily, tolerance=1.e-8_r8) + @assertEqual(0.7_r8, cohort%year_net_uptake(1), tolerance=1.e-8_r8) + + deallocate(cohort) + deallocate(patch) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_c13NegativeDenom + + ! NEW: Test zero leaf layer count (nv = 0) to ensure loop boundary protection + @Test + subroutine test_AccumulateFluxes_ED_ZeroLayers(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for zero canopy leaf layer boundary test + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch + type(fates_cohort_type), pointer :: cohort + real(r8) :: dt_time + + call CreateTestSite(sites, patch, cohort, bc_in, heights=[10.0_r8]) + + call SetTestFluxes(cohort, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=0) + + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Verify cohort-level GPP and respiration accumulated despite nv=0 + @assertEqual(1.5_r8, cohort%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.3_r8, cohort%resp_m_acc, tolerance=1.e-8_r8) + ! Net uptake loop (1 to nv) should not execute for nv = 0 + @assertEqual(0.2_r8, cohort%year_net_uptake(1), tolerance=1.e-8_r8) + + deallocate(cohort) + deallocate(patch) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_ZeroLayers + + ! NEW: Test multi-PFT cohort mix on a single patch (PFT 1 and PFT 2) + @Test + subroutine test_AccumulateFluxes_ED_MultiPFT(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for testing multi-PFT cohort mixes + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch + type(fates_cohort_type), pointer :: cohort_pft1, cohort_pft2 + real(r8) :: dt_time + + call CreateTestSite(sites, patch, cohort_pft1, bc_in, heights=[10.0_r8, 10.0_r8], pft=1) + cohort_pft2 => cohort_pft1%taller + cohort_pft2%pft = 2 + + call SetTestFluxes(cohort_pft1, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, & + resp_m_tstep=0.1_r8, sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.2_r8, & + ts_net_uptake=0.5_r8, nv=1) + + call SetTestFluxes(cohort_pft2, gpp_acc=3.0_r8, gpp_tstep=1.5_r8, resp_m_acc=0.6_r8, & + resp_m_tstep=0.3_r8, sym_nfix_daily=0.15_r8, sym_nfix_tstep=0.06_r8, & + c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, year_net_uptake=0.6_r8, & + ts_net_uptake=1.5_r8, nv=1) + + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Verify accumulation on both PFT 1 and PFT 2 cohorts + @assertEqual(1.5_r8, cohort_pft1%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.3_r8, cohort_pft1%resp_m_acc, tolerance=1.e-8_r8) + @assertEqual(4.5_r8, cohort_pft2%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(0.9_r8, cohort_pft2%resp_m_acc, tolerance=1.e-8_r8) + + deallocate(cohort_pft2) + deallocate(cohort_pft1) + deallocate(patch) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_MultiPFT + + ! NEW: Test multi-cohort linked list combined with multi-patch age sequence + @Test + subroutine test_AccumulateFluxes_ED_MultiCohortMultiPatch(this) + class(TestEDAccumulateFluxes), intent(inout) :: this + ! Single site configuration for multi-cohort multi-patch traversal + integer, parameter :: nsites = 1 + type(ed_site_type), target :: sites(nsites) + type(bc_in_type) :: bc_in(nsites) + type(bc_out_type) :: bc_out(nsites) + type(fates_patch_type), pointer :: patch1, patch2 + type(fates_cohort_type), pointer :: c1_p1, c2_p1, c1_p2, c2_p2 + real(r8) :: dt_time + + allocate(patch1) + allocate(patch2) + allocate(c1_p1) + call c1_p1%CreateBare(pft=1, nn=0.1_r8, height=5.0_r8) + allocate(c2_p1) + call c2_p1%CreateBare(pft=1, nn=0.1_r8, height=15.0_r8) + allocate(c1_p2) + call c1_p2%CreateBare(pft=1, nn=0.1_r8, height=5.0_r8) + allocate(c2_p2) + call c2_p2%CreateBare(pft=1, nn=0.1_r8, height=15.0_r8) + + ! Patch 1 with 2 cohorts + patch1%patchno = 1 + patch1%nocomp_pft_label = 1 + call SetTestFluxes(c1_p1, gpp_acc=1.0_r8, gpp_tstep=0.5_r8, resp_m_acc=0.2_r8, resp_m_tstep=0.1_r8, & + sym_nfix_daily=0.05_r8, sym_nfix_tstep=0.02_r8, c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, & + year_net_uptake=0.2_r8, ts_net_uptake=0.5_r8, nv=1) + + call SetTestFluxes(c2_p1, gpp_acc=2.0_r8, gpp_tstep=1.0_r8, resp_m_acc=0.4_r8, resp_m_tstep=0.2_r8, & + sym_nfix_daily=0.1_r8, sym_nfix_tstep=0.04_r8, c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, & + year_net_uptake=0.4_r8, ts_net_uptake=1.0_r8, nv=1) + + c1_p1%taller => c2_p1 + c2_p1%taller => null() + patch1%shortest => c1_p1 + + ! Patch 2 with 2 cohorts + patch2%patchno = 2 + patch2%nocomp_pft_label = 1 + call SetTestFluxes(c1_p2, gpp_acc=3.0_r8, gpp_tstep=1.5_r8, resp_m_acc=0.6_r8, resp_m_tstep=0.3_r8, & + sym_nfix_daily=0.15_r8, sym_nfix_tstep=0.06_r8, c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, & + year_net_uptake=0.6_r8, ts_net_uptake=1.5_r8, nv=1) + + call SetTestFluxes(c2_p2, gpp_acc=4.0_r8, gpp_tstep=2.0_r8, resp_m_acc=0.8_r8, resp_m_tstep=0.4_r8, & + sym_nfix_daily=0.20_r8, sym_nfix_tstep=0.08_r8, c13disc_acc=0.1_r8, c13disc_clm=0.2_r8, & + year_net_uptake=0.8_r8, ts_net_uptake=2.0_r8, nv=1) + + c1_p2%taller => c2_p2 + c2_p2%taller => null() + patch2%shortest => c1_p2 + + patch1%younger => patch2 + patch2%younger => null() + sites(1)%oldest_patch => patch1 + + allocate(bc_in(1)%filter_photo_pa(2)) + ! FATES photosynthesis filter flag 3 indicates active photosynthesis on vegetated patch + bc_in(1)%filter_photo_pa(1) = 3 + bc_in(1)%filter_photo_pa(2) = 3 + ! Timestep duration of 1 hour (3600s) + 0.2s offset (unreferenced dummy argument in AccumulateFluxes_ED) + dt_time = 3600.2_r8 + + call AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) + + ! Verify accumulation across all 4 cohorts over 2 patches + @assertEqual(1.5_r8, c1_p1%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(3.0_r8, c2_p1%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(4.5_r8, c1_p2%gpp_acc, tolerance=1.e-8_r8) + @assertEqual(6.0_r8, c2_p2%gpp_acc, tolerance=1.e-8_r8) + + deallocate(c2_p2) + deallocate(c1_p2) + deallocate(c2_p1) + deallocate(c1_p1) + deallocate(patch2) + deallocate(patch1) + deallocate(bc_in(1)%filter_photo_pa) + end subroutine test_AccumulateFluxes_ED_MultiCohortMultiPatch + +end module test_EDAccumulateFluxes