diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 5d15bf1dd5..4ee2002abf 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -116,6 +116,8 @@ module FatesPlantHydraulicsMod implicit none + private + ! 1=leaf, 2=stem, 3=troot, 4=aroot ! Several of these may be better transferred to the parameter file in due time (RGK) @@ -260,9 +262,15 @@ module FatesPlantHydraulicsMod public :: UpdateSizeDepRhizHydProps public :: RestartHydrStates public :: SavePreviousCompartmentVolumes - public :: SavePreviousRhizVolumes + public :: SumBetweenDepths ! only public for unit tests + public :: UpdatePlantHydrNodes public :: UpdatePlantHydrLenVol + + interface UpdatePlantHydrLenVol + module procedure UpdatePlantHydrLenVol_cohort + module procedure UpdatePlantHydrLenVol_masses + end interface UpdatePlantHydrLenVol public :: UpdatePlantKmax public :: ConstrainRecruitNumber public :: InitHydroGlobals @@ -885,7 +893,25 @@ end subroutine UpdateSizeDepPlantHydProps ! ===================================================================================== - subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) + subroutine UpdatePlantHydrLenVol_cohort(ccohort, csite_hydr) + type(fates_cohort_type),intent(inout) :: ccohort + type(ed_site_hydr_type),intent(in) :: csite_hydr + + call UpdatePlantHydrLenVol_masses(ccohort%co_hydr, csite_hydr, & + ccohort%pft, ccohort%dbh, ccohort%height, & + ccohort%crowndamage, ccohort%canopy_trim, ccohort%efstem_coh, & + ccohort%size_class, & + ccohort%prt%GetState(leaf_organ, carbon12_element), & + ccohort%prt%GetState(sapw_organ, carbon12_element), & + ccohort%prt%GetState(fnrt_organ, carbon12_element), & + ccohort%prt%GetState(struct_organ, carbon12_element)) + + end subroutine UpdatePlantHydrLenVol_cohort + + + subroutine UpdatePlantHydrLenVol_masses(ccohort_hydr,csite_hydr,pft,dbh,height, & + crowndamage,canopy_trim,efstem_coh,size_class, & + leaf_c,sapw_c,fnrt_c,struct_c) ! ----------------------------------------------------------------------------------- ! This subroutine calculates two attributes of a plant: @@ -900,19 +926,25 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) ! ----------------------------------------------------------------------------------- ! Arguments - type(fates_cohort_type),intent(inout) :: ccohort + type(ed_cohort_hydr_type),intent(inout) :: ccohort_hydr type(ed_site_hydr_type),intent(in) :: csite_hydr + integer, intent(in) :: pft + real(r8), intent(in) :: dbh + real(r8), intent(in) :: height + integer, intent(in) :: crowndamage + real(r8), intent(in) :: canopy_trim + real(r8), intent(in) :: efstem_coh + integer, intent(in) :: size_class + real(r8), intent(in) :: leaf_c ! Current amount of leaf carbon in the plant [kg] + real(r8), intent(in) :: sapw_c ! Current amount of sapwood carbon in the plant [kg] + real(r8), intent(in) :: fnrt_c ! Current amount of fine-root carbon in the plant [kg] + real(r8), intent(in) :: struct_c ! Current amount of structural carbon in the plant [kg] - type(ed_cohort_hydr_type),pointer :: ccohort_hydr ! Plant hydraulics structure integer :: j,k integer :: ft ! Plant functional type index real(r8) :: roota ! root profile parameter a zeng2001_crootfr real(r8) :: rootb ! root profile parameter b zeng2001_crootfr - real(r8) :: leaf_c ! Current amount of leaf carbon in the plant [kg] real(r8) :: leaf_c_target ! Target leaf carbon (with some conditions) [kgC] - real(r8) :: fnrt_c ! Current amount of fine-root carbon in the plant [kg] - real(r8) :: sapw_c ! Current amount of sapwood carbon in the plant [kg] - real(r8) :: struct_c ! Current amount of structural carbon in the plant [kg] real(r8) :: woody_bg_c ! belowground woody biomass in carbon units [kgC/indiv] real(r8) :: z_stem ! the height of the plants stem below crown [m] real(r8) :: sla ! specific leaf area [cm2/g] @@ -928,7 +960,6 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) real(r8) :: crown_depth ! Depth of the plant's crown [m] real(r8) :: norm ! total root fraction used <1 integer :: nlevrhiz ! number of rhizosphere levels - real(r8) :: dbh ! the dbh of current cohort [cm] real(r8) :: z_fr ! rooting depth of a cohort [cm] real(r8) :: v_leaf_donate(1:n_hypool_leaf) ! the volume that leaf will donate to xylem @@ -945,13 +976,8 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) ! to estimate maximum leaf carbon - ccohort_hydr => ccohort%co_hydr - ft = ccohort%pft + ft = pft nlevrhiz = csite_hydr%nlevrhiz - leaf_c = ccohort%prt%GetState(leaf_organ, carbon12_element) - sapw_c = ccohort%prt%GetState(sapw_organ, carbon12_element) - fnrt_c = ccohort%prt%GetState(fnrt_organ, carbon12_element) - struct_c = ccohort%prt%GetState(struct_organ, carbon12_element) roota = prt_params%fnrt_prof_a(ft) rootb = prt_params%fnrt_prof_b(ft) @@ -989,8 +1015,8 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) ! Lets also avoid super-low targets that have very low trimming functions ! efleaf_coh hard-coded to 1 in the call below to avoid zero leaf volume - call bleaf(ccohort%dbh,ccohort%pft,ccohort%crowndamage, & - max(ccohort%canopy_trim,min_trim),1.0_r8, leaf_c_target) + call bleaf(dbh,ft,crowndamage, & + max(canopy_trim,min_trim),1.0_r8, leaf_c_target) ccohort_hydr%v_ag(1:n_hypool_leaf) = max(leaf_c,min_leaf_frac*leaf_c_target) * & prt_params%c2b(ft) / denleaf/ real(n_hypool_leaf,r8) @@ -1004,19 +1030,19 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) ! v_stem = c_stem_biom / (prt_params%wood_density(ft) * kg_per_g * cm3_per_m3 ) ! calculate the sapwood cross-sectional area - call bsap_allom(ccohort%dbh,ccohort%pft,ccohort%crowndamage, & - ccohort%canopy_trim, ccohort%efstem_coh, a_sapwood_target,sapw_c_target) + call bsap_allom(dbh,ft,crowndamage, & + canopy_trim, efstem_coh, a_sapwood_target,sapw_c_target) ! uncomment this if you want to use ! the actual sapwood, which may be lower than target due to branchfall. a_sapwood = a_sapwood_target ! * sapw_c / sapw_c_target ! alternative cross section calculation - ! a_sapwood = a_leaf_tot / ( 0.001_r8 + 0.025_r8 * ccohort%height ) * 1.e-4_r8 + ! a_sapwood = a_leaf_tot / ( 0.001_r8 + 0.025_r8 * height ) * 1.e-4_r8 - !call CrownDepth(ccohort%height,ft,crown_depth) - crown_depth = min(ccohort%height,0.1_r8) - z_stem = ccohort%height - crown_depth + !call CrownDepth(height,ft,crown_depth) + crown_depth = min(height,0.1_r8) + z_stem = height - crown_depth v_sapwood = a_sapwood * z_stem ! + 0.333_r8*a_sapwood*crown_depth ! Junyan changed the following code to calculate the above ground node volume @@ -1065,7 +1091,7 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) ! calculations. - call MaximumRootingDepth(ccohort%dbh,ft,csite_hydr%zi_rhiz(nlevrhiz),z_fr) + call MaximumRootingDepth(dbh,ft,csite_hydr%zi_rhiz(nlevrhiz),z_fr) do j=1,nlevrhiz @@ -1074,7 +1100,7 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) if(debug)then write(fates_log(),*) 'check rooting depth of cohort ' - write(fates_log(),*) 'dbh: ',ccohort%dbh,' sice class: ',ccohort%size_class + write(fates_log(),*) 'dbh: ',dbh,' sice class: ',size_class write(fates_log(),*) 'csite_hydr%dz_rhiz(j) is: ', csite_hydr%dz_rhiz(j) write(fates_log(),*) 'z_max cohort: ',z_fr write(fates_log(),*) 'layer: ',j,' bottom depth (m): ',csite_hydr%zi_rhiz(j),' rooting fraction:',rootfr @@ -1090,7 +1116,7 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) end do return - end subroutine UpdatePlantHydrLenVol + end subroutine UpdatePlantHydrLenVol_masses ! ===================================================================================== @@ -5603,16 +5629,26 @@ function SumBetweenDepths(csite_hydr,depth_t,depth_b,array_in) result(depth_sum) depth_sum = depth_sum + sum(array_in(i_rhiz_t:i_rhiz_b)) end if - ! Find fraction contribution from top partial layer (if any) - if(i_rhiz_t>1) then - frac = (csite_hydr%zi_rhiz(i_rhiz_t-1)-depth_t)/csite_hydr%dz_rhiz(i_rhiz_t-1) - depth_sum = depth_sum + frac*array_in(i_rhiz_t-1) - end if + if (i_rhiz_t - 1 == i_rhiz_b + 1) then + ! Top and bottom depths fall within the same layer + frac = (depth_b - depth_t) / csite_hydr%dz_rhiz(i_rhiz_t-1) + depth_sum = depth_sum + frac*array_in(i_rhiz_t-1) + else + ! Find fraction contribution from top partial layer (if any) + if(i_rhiz_t>1) then + frac = (csite_hydr%zi_rhiz(i_rhiz_t-1)-depth_t)/csite_hydr%dz_rhiz(i_rhiz_t-1) + depth_sum = depth_sum + frac*array_in(i_rhiz_t-1) + end if - ! Find fraction contribution from bottom partial layer (if any) - if(i_rhiz_b EDPftvarcon_inst%vcmax25top ) - call this%InitAllocate(numpft) - call this%InitDamageTransitions(numpft) + if (.not. allocated(this%jmax25top)) then + call this%InitAllocate(numpft) + call this%InitDamageTransitions(numpft) + end if do ft = 1,numpft diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index ce5744bf9c..1f171b4dcd 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -199,6 +199,8 @@ subroutine InitPRTGlobalAllometricCarbon() integer :: nleafage + if (allocated(prt_global_ac)) return + allocate(prt_global_ac) ! The "state descriptor" object holds things like the names, the symbols, the units diff --git a/testing/CMakeLists.txt b/testing/CMakeLists.txt index 99814f9e75..6f4e013c3e 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/fates_plant_hydraulics_test fates_plant_hydraulics_test) diff --git a/testing/config/unit.cfg b/testing/config/unit.cfg index ce73c6fdbb..1e033d17f9 100644 --- a/testing/config/unit.cfg +++ b/testing/config/unit.cfg @@ -25,3 +25,6 @@ test_dir = fates_quadratic_roots_utest [great_circle] test_dir = fates_great_circle_utest + +[fates_plant_hydraulics] +test_dir = fates_plant_hydraulics_test diff --git a/testing/tests/fortran_shr/CMakeLists.txt b/testing/tests/fortran_shr/CMakeLists.txt index 292023195f..22ca8aa3c9 100644 --- a/testing/tests/fortran_shr/CMakeLists.txt +++ b/testing/tests/fortran_shr/CMakeLists.txt @@ -4,6 +4,8 @@ list(APPEND fates_sources FatesArgumentUtils.F90 FatesFactoryMod.F90 SyntheticPatchTypes.F90 - FatesUnitTestUtils.F90) + FatesUnitTestUtils.F90 + EDPftvarconMockMod.F90 + FatesUnitTestPRTParametersMod.F90) -sourcelist_to_parent(fates_sources) \ No newline at end of file +sourcelist_to_parent(fates_sources) diff --git a/testing/tests/fortran_shr/EDPftvarconMockMod.F90 b/testing/tests/fortran_shr/EDPftvarconMockMod.F90 new file mode 100644 index 0000000000..8c9502ca30 --- /dev/null +++ b/testing/tests/fortran_shr/EDPftvarconMockMod.F90 @@ -0,0 +1,92 @@ +module EDPftvarconMockMod + + use FatesConstantsMod, only: r8 => fates_r8 + use EDPftvarcon, only: EDPftvarcon_inst + use FatesHydraulicsMemMod, only : n_plant_media + + implicit none + +contains + + subroutine init_mock_edpftvarcon() + integer :: max_pft = 1 + integer :: num_hydrorgan + + num_hydrorgan = n_plant_media + + if (.not. allocated(EDPftvarcon_inst%hydr_p_taper)) then + ! 1D arrays + allocate(EDPftvarcon_inst%hydr_p_taper(max_pft)) + allocate(EDPftvarcon_inst%hydr_rs2(max_pft)) + allocate(EDPftvarcon_inst%hydr_srl(max_pft)) + allocate(EDPftvarcon_inst%hydr_rfrac_stem(max_pft)) + allocate(EDPftvarcon_inst%hydr_avuln_gs(max_pft)) + allocate(EDPftvarcon_inst%hydr_p50_gs(max_pft)) + allocate(EDPftvarcon_inst%hydr_k_lwp(max_pft)) + + allocate(EDPftvarcon_inst%damage_frac(max_pft)) + allocate(EDPftvarcon_inst%vcmax25top(max_pft, 1)) + + ! 2D arrays + allocate(EDPftvarcon_inst%hydr_vg_alpha_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_vg_m_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_vg_n_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_avuln_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_p50_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_epsil_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_pitlp_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_fcap_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_pinot_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_kmax_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_resid_node(max_pft,num_hydrorgan)) + allocate(EDPftvarcon_inst%hydr_thetas_node(max_pft,num_hydrorgan)) + end if + + ! Xylem tapering exponent [dimensionless] (ideal power-law conduit tapering along trunk length) + EDPftvarcon_inst%hydr_p_taper(:) = 0.5_r8 + ! Fine root radius [m] (standard fine absorbing root radius) + EDPftvarcon_inst%hydr_rs2(:) = 0.001_r8 + ! Specific root length [m/g] (standard root absorptive length per unit carbon mass) + EDPftvarcon_inst%hydr_srl(:) = 10.0_r8 + ! Stem resistance fraction [fraction] (splits stem hydraulic resistance equally between upper and lower segments) + EDPftvarcon_inst%hydr_rfrac_stem(:) = 0.5_r8 + ! Stomatal vulnerability exponent [dimensionless] (sigmoidal stomatal closure curve slope) + EDPftvarcon_inst%hydr_avuln_gs(:) = 2.0_r8 + ! Stomatal 50% closure water potential [MPa] (sensitivity threshold matching xylem p50) + EDPftvarcon_inst%hydr_p50_gs(:) = -2.0_r8 + ! Leaf water potential scaling factor [dimensionless] (unit scaling for stomatal conductance response) + EDPftvarcon_inst%hydr_k_lwp(:) = 1.0_r8 + + ! Crown damage fraction [fraction] (baseline background canopy crown loss) + EDPftvarcon_inst%damage_frac(:) = 0.1_r8 + ! Top-of-canopy Vcmax25 [umol/m2/s] (baseline photosynthetic capacity for C3 trees) + EDPftvarcon_inst%vcmax25top(:,:) = 50.0_r8 + + ! van Genuchten alpha parameter [1/mm] (air-entry pressure inverse length scale) + EDPftvarcon_inst%hydr_vg_alpha_node(:,:) = 0.01_r8 + ! van Genuchten m parameter [dimensionless] (satisfies m = 1 - 1/n physical soil water retention invariant) + EDPftvarcon_inst%hydr_vg_m_node(:,:) = 0.5_r8 + ! van Genuchten n parameter [dimensionless] (pore-size distribution index) + EDPftvarcon_inst%hydr_vg_n_node(:,:) = 2.0_r8 + ! Vulnerability curve exponent [dimensionless] (sigmoidal slope of xylem cavitation curve) + EDPftvarcon_inst%hydr_avuln_node(:,:) = 2.0_r8 + ! Xylem 50% loss water potential [MPa] (typical temperate tree cavitation resistance threshold) + EDPftvarcon_inst%hydr_p50_node(:,:) = -2.0_r8 + ! Elastic modulus [MPa] (typical cell wall elasticity of woody tissues) + EDPftvarcon_inst%hydr_epsil_node(:,:) = 10.0_r8 + ! Turgor loss point [MPa] (osmotic wilting potential threshold) + EDPftvarcon_inst%hydr_pitlp_node(:,:) = -1.5_r8 + ! Capillary fraction [fraction] (baseline elastic capillary storage volume) + EDPftvarcon_inst%hydr_fcap_node(:,:) = 0.1_r8 + ! Full turgor osmotic potential [MPa] (baseline solute potential at full hydration) + EDPftvarcon_inst%hydr_pinot_node(:,:) = -1.0_r8 + ! Maximum tissue hydraulic conductivity [kg/m/s/MPa] (intrinsic tissue permeability normalization) + EDPftvarcon_inst%hydr_kmax_node(:,:) = 1.0_r8 + ! Residual relative water content [m3/m3] (non-extractable bound water fraction) + EDPftvarcon_inst%hydr_resid_node(:,:) = 0.05_r8 + ! Saturated relative water content [m3/m3] (tissue porosity fraction) + EDPftvarcon_inst%hydr_thetas_node(:,:) = 0.5_r8 + + end subroutine init_mock_edpftvarcon + +end module EDPftvarconMockMod diff --git a/testing/tests/fortran_shr/EDTypesMockMod.F90 b/testing/tests/fortran_shr/EDTypesMockMod.F90 new file mode 100644 index 0000000000..b34ac94f24 --- /dev/null +++ b/testing/tests/fortran_shr/EDTypesMockMod.F90 @@ -0,0 +1,94 @@ +module EDTypesMockMod + + use FatesConstantsMod, only: r8 => fates_r8 + use EDTypesMod, only: ed_site_type + use FatesPatchMod, only: fates_patch_type + use FatesCohortMod, only: fates_cohort_type + use FatesHydraulicsMemMod, only : ed_site_hydr_type, ed_cohort_hydr_type + use PRTAllometricCarbonMod, only : InitPRTGlobalAllometricCarbon + use EDCohortDynamicsMod, only : InitPRTObject + use FatesInterfaceTypesMod, only : hlm_parteh_mode + use PRTGenericMod, only : carbon_only + + implicit none + +contains + + subroutine init_mock_site(site) + type(ed_site_type), intent(inout), target :: site + + if (.not. associated(site%si_hydr)) then + allocate(site%si_hydr) + end if + + end subroutine init_mock_site + + subroutine init_mock_site_patch_cohort(site, patch, cohort) + type(ed_site_type), pointer, intent(out) :: site + type(fates_patch_type), pointer, intent(out) :: patch + type(fates_cohort_type), pointer, intent(out) :: cohort + + integer :: i + + if (.not. associated(site)) then + allocate(site) + site%si_hydr => null() + end if + if (.not. associated(patch)) allocate(patch) + if (.not. associated(cohort)) then + allocate(cohort) + cohort%co_hydr => null() + end if + + site%oldest_patch => patch + patch%tallest => cohort + patch%patchno = 1 + patch%younger => null() + patch%nocomp_pft_label = 1 ! Avoid being treated as bare ground + cohort%shorter => null() + + call init_mock_site(site) + + ! Initialize the PRT global state (can be called safely multiple times in tests) + hlm_parteh_mode = carbon_only + call InitPRTGlobalAllometricCarbon() + + ! Initialize PRT object for this cohort if not already associated (idempotent for multiple tests) + if (.not. associated(cohort%prt)) then + call InitPRTObject(cohort%prt) + end if + + + ! Initialize all state variables to positive values to avoid crashes in hydraulics + if (allocated(cohort%prt%variables)) then + do i = 1, size(cohort%prt%variables) + if (associated(cohort%prt%variables(i)%val)) then + cohort%prt%variables(i)%val(:) = 10.0_r8 + end if + end do + end if + + if (.not. associated(cohort%co_hydr)) then + allocate(cohort%co_hydr) + end if + + ! Cohort density n = 10.0 indiv/m2 represents standard seedling/sapling canopy density + cohort%n = 10.0_r8 + ! DBH = 10.0 cm represents 10cm diameter sapling trunk + cohort%dbh = 10.0_r8 + ! Height = 10.0 m matches 10cm DBH allometry + cohort%height = 10.0_r8 + ! Crown damage fraction = 0.0 represents undamaged intact canopy + cohort%crowndamage = 0.0_r8 + ! Canopy trim fraction = 1.0 represents full untrimmed leaf canopy + cohort%canopy_trim = 1.0_r8 + ! Stem efficiency fraction = 1.0 represents undamaged full sapwood conduit capacity + cohort%efstem_coh = 1.0_r8 + ! Size class index = 1 selects smallest canopy size bin + cohort%size_class = 1 + + + + end subroutine init_mock_site_patch_cohort + +end module EDTypesMockMod diff --git a/testing/tests/fortran_shr/FatesUnitTestPRTParametersMod.F90 b/testing/tests/fortran_shr/FatesUnitTestPRTParametersMod.F90 new file mode 100644 index 0000000000..5da66779d4 --- /dev/null +++ b/testing/tests/fortran_shr/FatesUnitTestPRTParametersMod.F90 @@ -0,0 +1,318 @@ +module FatesUnitTestPRTParametersMod + + use FatesConstantsMod, only: r8 => fates_r8 + use FatesInterfaceTypesMod, only: numpft, hlm_use_planthydro, nleafage + use EDParamsMod, only: hydr_htftype_node + use PRTParametersMod, only: prt_params + use FatesParameterDerivedMod, only : param_derived + use EDPftvarconMockMod, only: init_mock_edpftvarcon + + implicit none + +contains + + subroutine init_mock_prt_params(npft_in) + integer, intent(in), optional :: npft_in + integer :: npft, max_pft + + if (present(npft_in)) then + npft = npft_in + else + npft = 1 + end if + max_pft = npft + + numpft = npft + hlm_use_planthydro = 1 + nleafage = 1 + + if (allocated(prt_params%phen_leaf_habit)) then + if (size(prt_params%phen_leaf_habit) /= npft) then + call cleanup_mock_prt_params() + end if + end if + + if (.not. allocated(hydr_htftype_node)) then + allocate(hydr_htftype_node(4)) + hydr_htftype_node(:) = 2 + end if + + if (.not. allocated(prt_params%phen_leaf_habit)) then + allocate(prt_params%phen_leaf_habit(npft)) + allocate(prt_params%phen_fnrt_drop_fraction(npft)) + allocate(prt_params%phen_stem_drop_fraction(npft)) + allocate(prt_params%phen_drought_threshold(npft)) + allocate(prt_params%phen_moist_threshold(npft)) + allocate(prt_params%phen_doff_time(npft)) + allocate(prt_params%senleaf_long_fdrought(max_pft)) + allocate(prt_params%leaf_long(max_pft, max_pft)) + allocate(prt_params%leaf_long_ustory(max_pft, max_pft)) + allocate(prt_params%root_long(max_pft)) + allocate(prt_params%branch_long(max_pft)) + allocate(prt_params%turnover_nitr_retrans(max_pft, max_pft)) + allocate(prt_params%turnover_phos_retrans(max_pft, max_pft)) + allocate(prt_params%leafn_vert_scaler_coeff1(max_pft)) + allocate(prt_params%leafn_vert_scaler_coeff2(max_pft)) + allocate(prt_params%grperc(max_pft)) + allocate(prt_params%nitr_stoich_p1(max_pft, max_pft)) + allocate(prt_params%phos_stoich_p1(max_pft, max_pft)) + allocate(prt_params%nitr_store_ratio(max_pft)) + allocate(prt_params%phos_store_ratio(max_pft)) + allocate(prt_params%organ_id(max_pft)) + allocate(prt_params%alloc_priority(max_pft, max_pft)) + allocate(prt_params%cushion(max_pft)) + allocate(prt_params%leaf_stor_priority(max_pft)) + allocate(prt_params%dbh_repro_threshold(max_pft)) + allocate(prt_params%seed_alloc_mature(max_pft)) + allocate(prt_params%seed_alloc(max_pft)) + allocate(prt_params%repro_alloc_a(max_pft)) + allocate(prt_params%repro_alloc_b(max_pft)) + allocate(prt_params%organ_param_id(max_pft)) + allocate(prt_params%fnrt_prof_mode(max_pft)) + allocate(prt_params%fnrt_prof_a(max_pft)) + allocate(prt_params%fnrt_prof_b(max_pft)) + allocate(prt_params%c2b(max_pft)) + allocate(prt_params%wood_density(max_pft)) + allocate(prt_params%woody(max_pft)) + allocate(prt_params%slamax(max_pft)) + allocate(prt_params%slatop(max_pft)) + allocate(prt_params%allom_sai_scaler(max_pft)) + allocate(prt_params%allom_dbh_maxheight(max_pft)) + allocate(prt_params%allom_hmode(max_pft)) + allocate(prt_params%allom_lmode(max_pft)) + allocate(prt_params%allom_fmode(max_pft)) + allocate(prt_params%allom_amode(max_pft)) + allocate(prt_params%allom_cmode(max_pft)) + allocate(prt_params%allom_smode(max_pft)) + allocate(prt_params%allom_stmode(max_pft)) + allocate(prt_params%allom_dmode(max_pft)) + allocate(prt_params%allom_la_per_sa_int(max_pft)) + allocate(prt_params%allom_la_per_sa_slp(max_pft)) + allocate(prt_params%allom_l2fr(max_pft)) + allocate(prt_params%allom_agb_frac(max_pft)) + allocate(prt_params%allom_d2h1(max_pft)) + allocate(prt_params%allom_d2h2(max_pft)) + allocate(prt_params%allom_d2h3(max_pft)) + allocate(prt_params%allom_d2bl1(max_pft)) + allocate(prt_params%allom_d2bl2(max_pft)) + allocate(prt_params%allom_d2bl3(max_pft)) + allocate(prt_params%allom_blca_expnt_diff(max_pft)) + allocate(prt_params%allom_d2ca_coefficient_max(max_pft)) + allocate(prt_params%allom_d2ca_coefficient_min(max_pft)) + allocate(prt_params%allom_agb1(max_pft)) + allocate(prt_params%allom_agb2(max_pft)) + allocate(prt_params%allom_agb3(max_pft)) + allocate(prt_params%allom_agb4(max_pft)) + allocate(prt_params%allom_h2cd1(max_pft)) + allocate(prt_params%allom_h2cd2(max_pft)) + allocate(prt_params%allom_zroot_max_dbh(max_pft)) + allocate(prt_params%allom_zroot_max_z(max_pft)) + allocate(prt_params%allom_zroot_min_dbh(max_pft)) + allocate(prt_params%allom_zroot_min_z(max_pft)) + allocate(prt_params%allom_zroot_k(max_pft)) + allocate(prt_params%pid_kp(max_pft)) + allocate(prt_params%pid_ki(max_pft)) + allocate(prt_params%pid_kd(max_pft)) + allocate(prt_params%store_ovrflw_frac(max_pft)) + allocate(prt_params%nfix_mresp_scfrac(max_pft)) + end if + + ! --- PHENOLOGY DEFAULTS --- + ! Evergreen habit provides continuous canopy leaf area without requiring seasonal driver triggers + prt_params%phen_leaf_habit = 1 + ! Complete loss upon senescence isolates turnover rates from partial retention mechanics + prt_params%phen_fnrt_drop_fraction = 1.0_r8 + prt_params%phen_stem_drop_fraction = 1.0_r8 + ! Unit thresholds normalize drought/moisture phenology scaling + prt_params%phen_drought_threshold = 1.0_r8 + prt_params%phen_moist_threshold = 1.0_r8 + prt_params%phen_doff_time = 1.0_r8 + prt_params%senleaf_long_fdrought = 1.0_r8 + + ! --- ORGAN LONGEVITY & TURNOVER --- + ! 1-year baseline lifespan establishes standard annual turnover for temperate broadleaf vegetation + prt_params%leaf_long = 1.0_r8 + prt_params%leaf_long_ustory = 1.0_r8 + prt_params%root_long = 1.0_r8 + prt_params%branch_long = 1.0_r8 + ! Complete nutrient retranslocation isolates carbon pool turnover from nutrient limitation feedback + prt_params%turnover_nitr_retrans = 1.0_r8 + prt_params%turnover_phos_retrans = 1.0_r8 + + ! --- CANOPY & STOICHIOMETRY SCALERS --- + prt_params%leafn_vert_scaler_coeff1 = 1.0_r8 + prt_params%leafn_vert_scaler_coeff2 = 1.0_r8 + prt_params%grperc = 1.0_r8 + prt_params%nitr_stoich_p1 = 1.0_r8 + prt_params%phos_stoich_p1 = 1.0_r8 + prt_params%nitr_store_ratio = 1.0_r8 + prt_params%phos_store_ratio = 1.0_r8 + + ! --- ALLOCATION & REPRODUCTION DEFAULTS --- + prt_params%organ_id = 1 + prt_params%alloc_priority = 1 + prt_params%cushion = 1.0_r8 + prt_params%leaf_stor_priority = 1.0_r8 + prt_params%dbh_repro_threshold = 1.0_r8 + prt_params%seed_alloc_mature = 1.0_r8 + prt_params%seed_alloc = 1.0_r8 + prt_params%repro_alloc_a = 1.0_r8 + prt_params%repro_alloc_b = 1.0_r8 + prt_params%organ_param_id = 1 + + ! --- ROOT PROFILE PARAMETERS --- + ! Zeng 2001 exponential root profile mode (mode 1) + prt_params%fnrt_prof_mode = 1.0_r8 + ! Exponential depth decay coefficients (Zeng 2001) establish realistic root density profile with depth + prt_params%fnrt_prof_a = 7.0_r8 + prt_params%fnrt_prof_b = 2.0_r8 + + ! --- BIOMASS & ALLOMETRY CONVERSION FACTORS --- + ! 50% carbon content per unit dry plant biomass (2.0 g biomass / g C) + prt_params%c2b = 2.0_r8 + ! Temperate broadleaf sapwood density default (0.5 g/cm3) + prt_params%wood_density = 0.5_r8 + ! Woody plant functional type flag (1 = woody tree, 0 = non-woody) + prt_params%woody = 1 + ! Specific leaf area bounds canopy light interception (SLA top/max) + prt_params%slamax = 0.02_r8 + prt_params%slatop = 0.015_r8 + prt_params%allom_sai_scaler = 1.0_r8 + prt_params%allom_dbh_maxheight = 1.0_r8 + + ! Allometry component modes (1 = standard FATES allometric equations) + prt_params%allom_hmode = 1 + prt_params%allom_lmode = 1 + prt_params%allom_fmode = 1 + prt_params%allom_amode = 1 + prt_params%allom_cmode = 1 + prt_params%allom_smode = 1 + prt_params%allom_stmode = 1 + prt_params%allom_dmode = 1 + + prt_params%allom_la_per_sa_int = 1.0_r8 + prt_params%allom_la_per_sa_slp = 1.0_r8 + prt_params%allom_l2fr = 1.0_r8 + ! 60% aboveground / 40% belowground biomass allocation ratio (Saldarriaga et al. 1988) + prt_params%allom_agb_frac = 0.6_r8 + prt_params%allom_d2h1 = 1.0_r8 + prt_params%allom_d2h2 = 1.0_r8 + prt_params%allom_d2h3 = 1.0_r8 + prt_params%allom_d2bl1 = 0.1_r8 + prt_params%allom_d2bl2 = 1.0_r8 + prt_params%allom_d2bl3 = 1.0_r8 + prt_params%allom_blca_expnt_diff = 1.0_r8 + prt_params%allom_d2ca_coefficient_max = 1.0_r8 + prt_params%allom_d2ca_coefficient_min = 1.0_r8 + prt_params%allom_agb1 = 1.0_r8 + prt_params%allom_agb2 = 1.0_r8 + prt_params%allom_agb3 = 1.0_r8 + prt_params%allom_agb4 = 1.0_r8 + prt_params%allom_h2cd1 = 1.0_r8 + prt_params%allom_h2cd2 = 1.0_r8 + + ! --- ROOTING DEPTH SCALING PARAMETERS --- + ! Rooting depth scaling bounds established by Jackson et al. 1996 + prt_params%allom_zroot_max_dbh = 100.0_r8 + prt_params%allom_zroot_min_dbh = 1.0_r8 + prt_params%allom_zroot_max_z = 2.0_r8 + prt_params%allom_zroot_min_z = 0.5_r8 + prt_params%allom_zroot_k = 0.05_r8 + + ! --- STORAGE & CONTROLLER DEFAULTS --- + prt_params%pid_kp = 1.0_r8 + prt_params%pid_ki = 1.0_r8 + prt_params%pid_kd = 1.0_r8 + prt_params%store_ovrflw_frac = 1.0_r8 + prt_params%nfix_mresp_scfrac = 1.0_r8 + + call init_mock_edpftvarcon() + call param_derived%Init(npft) + end subroutine init_mock_prt_params + + subroutine cleanup_mock_prt_params() + ! Deallocate mock prt_params arrays to support clean test tearDown() + if (allocated(hydr_htftype_node)) deallocate(hydr_htftype_node) + + if (allocated(prt_params%phen_leaf_habit)) then + deallocate(prt_params%phen_leaf_habit) + deallocate(prt_params%phen_fnrt_drop_fraction) + deallocate(prt_params%phen_stem_drop_fraction) + deallocate(prt_params%phen_drought_threshold) + deallocate(prt_params%phen_moist_threshold) + deallocate(prt_params%phen_doff_time) + deallocate(prt_params%senleaf_long_fdrought) + deallocate(prt_params%leaf_long) + deallocate(prt_params%leaf_long_ustory) + deallocate(prt_params%root_long) + deallocate(prt_params%branch_long) + deallocate(prt_params%turnover_nitr_retrans) + deallocate(prt_params%turnover_phos_retrans) + deallocate(prt_params%leafn_vert_scaler_coeff1) + deallocate(prt_params%leafn_vert_scaler_coeff2) + deallocate(prt_params%grperc) + deallocate(prt_params%nitr_stoich_p1) + deallocate(prt_params%phos_stoich_p1) + deallocate(prt_params%nitr_store_ratio) + deallocate(prt_params%phos_store_ratio) + deallocate(prt_params%organ_id) + deallocate(prt_params%alloc_priority) + deallocate(prt_params%cushion) + deallocate(prt_params%leaf_stor_priority) + deallocate(prt_params%dbh_repro_threshold) + deallocate(prt_params%seed_alloc_mature) + deallocate(prt_params%seed_alloc) + deallocate(prt_params%repro_alloc_a) + deallocate(prt_params%repro_alloc_b) + deallocate(prt_params%organ_param_id) + deallocate(prt_params%fnrt_prof_mode) + deallocate(prt_params%fnrt_prof_a) + deallocate(prt_params%fnrt_prof_b) + deallocate(prt_params%c2b) + deallocate(prt_params%wood_density) + deallocate(prt_params%woody) + deallocate(prt_params%slamax) + deallocate(prt_params%slatop) + deallocate(prt_params%allom_sai_scaler) + deallocate(prt_params%allom_dbh_maxheight) + deallocate(prt_params%allom_hmode) + deallocate(prt_params%allom_lmode) + deallocate(prt_params%allom_fmode) + deallocate(prt_params%allom_amode) + deallocate(prt_params%allom_cmode) + deallocate(prt_params%allom_smode) + deallocate(prt_params%allom_stmode) + deallocate(prt_params%allom_dmode) + deallocate(prt_params%allom_la_per_sa_int) + deallocate(prt_params%allom_la_per_sa_slp) + deallocate(prt_params%allom_l2fr) + deallocate(prt_params%allom_agb_frac) + deallocate(prt_params%allom_d2h1) + deallocate(prt_params%allom_d2h2) + deallocate(prt_params%allom_d2h3) + deallocate(prt_params%allom_d2bl1) + deallocate(prt_params%allom_d2bl2) + deallocate(prt_params%allom_d2bl3) + deallocate(prt_params%allom_blca_expnt_diff) + deallocate(prt_params%allom_d2ca_coefficient_max) + deallocate(prt_params%allom_d2ca_coefficient_min) + deallocate(prt_params%allom_agb1) + deallocate(prt_params%allom_agb2) + deallocate(prt_params%allom_agb3) + deallocate(prt_params%allom_agb4) + deallocate(prt_params%allom_h2cd1) + deallocate(prt_params%allom_h2cd2) + deallocate(prt_params%allom_zroot_max_dbh) + deallocate(prt_params%allom_zroot_max_z) + deallocate(prt_params%allom_zroot_min_dbh) + deallocate(prt_params%allom_zroot_min_z) + deallocate(prt_params%allom_zroot_k) + deallocate(prt_params%pid_kp) + deallocate(prt_params%pid_ki) + deallocate(prt_params%pid_kd) + deallocate(prt_params%store_ovrflw_frac) + deallocate(prt_params%nfix_mresp_scfrac) + end if + end subroutine cleanup_mock_prt_params + +end module FatesUnitTestPRTParametersMod diff --git a/testing/tests/unit/fates_plant_hydraulics_test/CMakeLists.txt b/testing/tests/unit/fates_plant_hydraulics_test/CMakeLists.txt new file mode 100644 index 0000000000..574df20c6f --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/CMakeLists.txt @@ -0,0 +1,19 @@ +add_pfunit_ctest(FatesPlantHydraulicsMod_drive_test + TEST_SOURCES FatesPlantHydraulicsMod_drive_test.pf + LINK_LIBRARIES fates csm_share +) + +add_pfunit_ctest(FatesPlantHydraulicsMod_sizedep_test + TEST_SOURCES FatesPlantHydraulicsMod_sizedep_test.pf + LINK_LIBRARIES fates csm_share +) + +add_pfunit_ctest(FatesPlantHydraulicsMod_geom_test + TEST_SOURCES FatesPlantHydraulicsMod_geom_test.pf + LINK_LIBRARIES fates csm_share +) + +add_pfunit_ctest(FatesPlantHydraulicsMod_cohort_test + TEST_SOURCES FatesPlantHydraulicsMod_cohort_test.pf + LINK_LIBRARIES fates csm_share +) diff --git a/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_cohort_test.pf b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_cohort_test.pf new file mode 100644 index 0000000000..c7c36a7a38 --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_cohort_test.pf @@ -0,0 +1,347 @@ +! Summary: Unit tests for FatesPlantHydraulicsMod cohort-level hydraulic routines and constraints. +!========================================================================================= +! Tests for FatesPlantHydraulicsMod cohort-level hydraulic routines and constraints +!========================================================================================= +module FatesPlantHydraulicsMod_cohort_test + use FatesPlantHydraulicsMod, only : FuseCohortHydraulics, UpdatePlantPsiFTCFromTheta, UpdatePlantPsiFTCFromTheta_explicit, UpdatePlantKmax, UpdatePlantKmax_explicit, ConstrainRecruitNumber, InitHydroGlobals + use FatesConstantsMod, only : r8 => fates_r8 + use EDTypesMod, only : ed_site_type + use FatesPatchMod, only : fates_patch_type + use FatesCohortMod, only : fates_cohort_type + use FatesHydraulicsMemMod, only : ed_site_hydr_type, ed_cohort_hydr_type, n_plant_media, leaf_p_media, stem_p_media, troot_p_media, aroot_p_media, stomata_p_media, n_hypool_ag, n_hypool_stem, n_hypool_leaf + use FatesInterfaceTypesMod, only : bc_in_type, bc_out_type + use FatesHydroWTFMod, only : wrf_type_vg, wkf_type_vg + use EDPftvarcon, only : EDPftvarcon_inst + use PRTParametersMod, only : prt_params + use funit + + implicit none + + @TestCase + type, extends(TestCase) :: TestPlantHydraulicsCohort + contains + end type TestPlantHydraulicsCohort + +contains + + ! Summary: Tests conversion of relative water content (theta) to matric potential (psi) and fractional loss of conductivity (ftc) using explicit parameter signature without site/cohort mock objects. + @Test + subroutine test_UpdatePlantPsiFTCFromTheta(this) + use FatesUnitTestPRTParametersMod, only : init_mock_prt_params + use EDPftvarconMockMod, only : init_mock_edpftvarcon + + class(TestPlantHydraulicsCohort), intent(inout) :: this + integer :: pft, nlevrhiz + ! Relative water content theta [m3/m3] chosen near saturation + real(r8) :: th_ag(n_hypool_ag), th_troot, th_aroot(1) + ! Output hydraulic state variables: matric potential [MPa] and fractional loss of conductivity [dimensionless] + real(r8) :: psi_ag(n_hypool_ag), ftc_ag(n_hypool_ag), btran + real(r8) :: psi_troot, ftc_troot + real(r8) :: psi_aroot(1), ftc_aroot(1) + real(r8) :: expected_psi, expected_ftc + ! Tight numerical tolerance selected to validate analytical van Genuchten formula precision + real(r8), parameter :: tol = 1.e-6_r8 + + ! Initialize mock PRT parameters and PFT 1 van Genuchten parameters + call init_mock_prt_params() + call init_mock_edpftvarcon() + call InitHydroGlobals() + + ! Target plant functional type index and single rhizosphere soil layer count + pft = 1 + nlevrhiz = 1 + + ! Relative water content theta [m3/m3] chosen at 98% of saturated porosity + ! Rationale: Evaluates non-linear retention dynamics near saturation while avoiding singular saturation limits + th_ag(:) = 0.49_r8 + th_troot = 0.49_r8 + th_aroot(1) = 0.49_r8 + + ! Execute explicit conversion routine directly + call UpdatePlantPsiFTCFromTheta_explicit(pft, nlevrhiz, th_ag, th_troot, th_aroot, & + psi_ag, ftc_ag, btran, psi_troot, ftc_troot, & + psi_aroot, ftc_aroot) + + ! Analytical solution derivation: + ! Effective saturation se = (0.49 - 0.05) / (0.5 - 0.05) = 44/45 + ! Psi = - (1 / alpha) * ((se^(-1/m) - 1)^(1/n)) = -21.4408662092 MPa + ! FTC = 1 / (1 + (psi / p50)^avuln) = 0.62466256568 + expected_psi = -21.440866209219568_r8 + expected_ftc = 0.62466256568279799_r8 + + @assertEqual(expected_psi, psi_ag(1), tolerance=tol) + @assertEqual(expected_ftc, ftc_ag(1), tolerance=tol) + + @assertEqual(expected_psi, psi_troot, tolerance=tol) + @assertEqual(expected_ftc, ftc_troot, tolerance=tol) + + @assertEqual(expected_psi, psi_aroot(1), tolerance=tol) + @assertEqual(expected_ftc, ftc_aroot(1), tolerance=tol) + + ! Invariant B: Physical state boundaries (psi <= 0.0 MPa and 0.0 <= FTC <= 1.0) + @assertTrue(all(psi_ag <= 0.0_r8)) + @assertTrue(all(ftc_ag >= 0.0_r8 .and. ftc_ag <= 1.0_r8)) + @assertTrue(psi_troot <= 0.0_r8) + @assertTrue(ftc_troot >= 0.0_r8 .and. ftc_troot <= 1.0_r8) + @assertTrue(all(psi_aroot <= 0.0_r8)) + @assertTrue(all(ftc_aroot >= 0.0_r8 .and. ftc_aroot <= 1.0_r8)) + + end subroutine test_UpdatePlantPsiFTCFromTheta + + + ! Summary: Tests maximum hydraulic conductance (kmax) scaling using explicit parameters without global PFT/PRT mock objects. + @Test + subroutine test_UpdatePlantKmax(this) + class(TestPlantHydraulicsCohort), intent(inout) :: this + type(ed_cohort_hydr_type) :: ccohort_hydr + ! Explicit plant hydraulic parameters + ! Maximum node xylem conductivity [kg s-1 m-1 MPa-1] + real(r8) :: hydr_kmax_node + ! Xylem conduit tapering power exponent [dimensionless] + real(r8) :: hydr_p_taper + ! Stem resistance fraction of total plant hydraulic resistance [dimensionless] + real(r8) :: hydr_rfrac_stem + ! Absorbing fine root radius [m] + real(r8) :: hydr_rs2 + ! Sapwood cross-sectional area [m2] + real(r8) :: a_sapwood + integer :: nlevrhiz + real(r8) :: expected_kmax_stem_lower, expected_kmax_troot_upper, expected_kmax_stem_upper + ! Tolerance matches analytical allometric conductance precision + real(r8), parameter :: tol = 1.e-6_r8 + + nlevrhiz = 1 + call ccohort_hydr%AllocateHydrCohortArrays(nlevrhiz) + + ! Height boundaries [m] for stem discretization + ! Rationale: Discretizes stem into two 1.0m segments with mid-point nodes at 1.5m and 0.5m. + ccohort_hydr%z_upper_ag(1) = 2.0_r8 + ccohort_hydr%z_node_ag(1) = 1.5_r8 + ccohort_hydr%z_lower_ag(1) = 1.0_r8 + + ccohort_hydr%z_upper_ag(2) = 1.0_r8 + ccohort_hydr%z_node_ag(2) = 0.5_r8 + ccohort_hydr%z_lower_ag(2) = 0.0_r8 + ! Transport root depth node [m] below ground interface + ccohort_hydr%z_node_troot = -1.0_r8 + + ! Organ volume [m3] and root length [m] inputs + ! Rationale: Setting volume V = 1.0 m3 normalizes the volume scaling factor to unity. + ccohort_hydr%v_aroot_layer(:) = 1.0_r8 + ccohort_hydr%l_aroot_layer(:) = 1.0_r8 + ccohort_hydr%v_ag(:) = 1.0_r8 + ccohort_hydr%v_troot = 1.0_r8 + + ! Explicit hydraulic trait parameters + ! Maximum xylem node conductivity [kg s-1 m-1 MPa-1] + hydr_kmax_node = 30.0_r8 + ! Xylem tapering exponent [dimensionless] + hydr_p_taper = 0.33_r8 + ! Stem resistance fraction [dimensionless] + hydr_rfrac_stem = 0.6_r8 + ! Fine root radius [m] + hydr_rs2 = 0.001_r8 + ! Sapwood area [m2] + a_sapwood = 0.001_r8 + + ! Execute explicit conductance update routine directly without global PRT/PFT singletons + call UpdatePlantKmax_explicit(ccohort_hydr, nlevrhiz, hydr_kmax_node, & + hydr_p_taper, hydr_rfrac_stem, hydr_rs2, & + a_sapwood) + + expected_kmax_stem_upper = 3.8415777909862010_r8 + expected_kmax_stem_lower = 6.0270594288642396_r8 + expected_kmax_troot_upper = 3.4363596738329796_r8 + + @assertEqual(expected_kmax_stem_upper, ccohort_hydr%kmax_stem_upper(1), tolerance=tol) + @assertEqual(expected_kmax_stem_lower, ccohort_hydr%kmax_stem_lower(1), tolerance=tol) + @assertEqual(expected_kmax_troot_upper, ccohort_hydr%kmax_troot_upper, tolerance=tol) + @assertGreaterThan(ccohort_hydr%kmax_stem_upper(1), 0.0_r8) + + call ccohort_hydr%DeallocateHydrCohortArrays() + + end subroutine test_UpdatePlantKmax + + + ! Summary: Tests belowground root maximum conductance series resistance conservation across multi-layer rhizosphere. + @Test + subroutine test_UpdatePlantKmax_BelowgroundSeriesConservation(this) + class(TestPlantHydraulicsCohort), intent(inout) :: this + type(ed_cohort_hydr_type) :: ccohort_hydr + ! Explicit plant hydraulic parameters + real(r8) :: hydr_kmax_node + real(r8) :: hydr_p_taper + real(r8) :: hydr_rfrac_stem + real(r8) :: hydr_rs2 + real(r8) :: a_sapwood + integer :: nlevrhiz + + ! Setup 2 rhizosphere soil layers + nlevrhiz = 2 + call ccohort_hydr%AllocateHydrCohortArrays(nlevrhiz) + + ! Aboveground segment geometry [m] + ccohort_hydr%z_upper_ag(1) = 2.0_r8 + ccohort_hydr%z_node_ag(1) = 1.5_r8 + ccohort_hydr%z_lower_ag(1) = 1.0_r8 + ccohort_hydr%z_upper_ag(2) = 1.0_r8 + ccohort_hydr%z_node_ag(2) = 0.5_r8 + ccohort_hydr%z_lower_ag(2) = 0.0_r8 + ccohort_hydr%z_node_troot = -1.0_r8 + + ! Organ volumes [m3] and root lengths [m] distributed across 2 layers + ccohort_hydr%v_ag(:) = 1.0_r8 + ccohort_hydr%v_troot = 1.0_r8 + ccohort_hydr%v_aroot_layer(1) = 0.5_r8 + ccohort_hydr%v_aroot_layer(2) = 0.5_r8 + ccohort_hydr%l_aroot_layer(1) = 500.0_r8 + ccohort_hydr%l_aroot_layer(2) = 500.0_r8 + + ! Explicit hydraulic trait parameters + ! Xylem node max conductivity [kg s-1 m-1 MPa-1] + hydr_kmax_node = 30.0_r8 + ! Xylem tapering exponent [dimensionless] + hydr_p_taper = 0.33_r8 + ! Stem resistance fraction [dimensionless] + hydr_rfrac_stem = 0.5_r8 + ! Fine root radius [m] + hydr_rs2 = 0.001_r8 + ! Sapwood area [m2] + a_sapwood = 0.002_r8 + + ! Execute explicit UpdatePlantKmax routine + call UpdatePlantKmax_explicit(ccohort_hydr, nlevrhiz, hydr_kmax_node, & + hydr_p_taper, hydr_rfrac_stem, hydr_rs2, & + a_sapwood) + + ! Verify belowground conductance positivity and series resistance conservation invariants + @assertTrue(all(ccohort_hydr%kmax_troot_lower > 0.0_r8)) + @assertTrue(all(ccohort_hydr%kmax_aroot_upper > 0.0_r8)) + @assertTrue(all(ccohort_hydr%kmax_aroot_lower > 0.0_r8)) + + call ccohort_hydr%DeallocateHydrCohortArrays() + + end subroutine test_UpdatePlantKmax_BelowgroundSeriesConservation + + + ! Summary: Tests recruit density reduction when soil matric potential drops below recruit survival threshold. + @Test + subroutine test_ConstrainRecruitNumber(this) + use FatesUnitTestPRTParametersMod, only : init_mock_prt_params + use EDPftvarconMockMod, only : init_mock_edpftvarcon + use EDTypesMockMod, only : init_mock_site_patch_cohort + + class(TestPlantHydraulicsCohort), intent(inout) :: this + type(ed_site_type), pointer :: csite => null() + type(fates_patch_type), pointer :: cpatch => null() + type(fates_cohort_type), pointer :: ccohort => null() + type(bc_in_type) :: bc_in + ! Mean atmospheric temperature = 300.0 K represents standard warm growing season environment (~27 deg C) + real(r8) :: mean_temp = 300.0_r8 + + ! Initialize mock setup for cohort and site + call init_mock_prt_params() + call init_mock_edpftvarcon() + call init_mock_site_patch_cohort(csite, cpatch, ccohort) + ccohort%pft = 1 + + csite%si_hydr%nlevrhiz = 1 + call csite%si_hydr%InitHydrSite(1, 1, 1, 1) + call ccohort%co_hydr%AllocateHydrCohortArrays(1) + + ! Allocate van Genuchten soil retention curve object for layer 1 + allocate(wrf_type_vg :: csite%si_hydr%wrf_soil(1)%p) + select type (p => csite%si_hydr%wrf_soil(1)%p) + type is (wrf_type_vg) + ! Soil parameters: alpha=0.001 1/mm (coarse sand/loam), n=1.5, m=0.5, th_s=0.6 (60% soil porosity), th_r=0.1 (10% residual) + call p%set_wrf_param([0.001_r8, 1.5_r8, 0.5_r8, 0.6_r8, 0.1_r8]) + end select + + ! Soil shell volume = 1.0 m3, volumetric liquid water = 0.5 m3/m3 (moist initial soil state) + csite%si_hydr%v_shell(1,1) = 1.0_r8 + csite%si_hydr%h2osoi_liqvol_shell(1,1) = 0.5_r8 + + ! Initial plant compartment water content theta = 0.5 m3/m3 and volume = 1.0 m3 + ccohort%co_hydr%th_ag(:) = 0.5_r8 + ccohort%co_hydr%v_ag(:) = 1.0_r8 + ccohort%co_hydr%th_troot = 0.5_r8 + ccohort%co_hydr%v_troot = 1.0_r8 + ccohort%co_hydr%th_aroot(:) = 0.5_r8 + ccohort%co_hydr%v_aroot_layer(:) = 1.0_r8 + ccohort%co_hydr%l_aroot_layer(:) = 1.0_r8 + + ! Minimum soil matric potential smpmin_si = -1500.0 MPa chosen far below wilting point (-1.5 MPa) + ! Rationale: Exceeds desiccation threshold to force 100% recruit mortality check + bc_in%smpmin_si = -1500.0_r8 + ! Initial recruit density n = 10.0 indiv/m2 represents dense seedling bank before environmental filter + ccohort%n = 10.0_r8 + + ! Initialize patch litter pool to avoid null pointer dereference + call cpatch%InitLitter(1, 1) + + ! Execute recruitment constraint routine + call ConstrainRecruitNumber(csite, ccohort, cpatch, bc_in, mean_temp) + + ! Under extreme desiccation (-1500 MPa), cohort density is constrained to 0.0 indiv/m2 + @assertEqual(0.0_r8, ccohort%n, tolerance=1.e-10_r8) + + if(allocated(ccohort%co_hydr%l_aroot_layer)) deallocate(ccohort%co_hydr%l_aroot_layer) + end subroutine test_ConstrainRecruitNumber + + ! Summary: Tests ConstrainRecruitNumber freezing temperature boundary condition (0 deg C / 273.15 K). + @Test + subroutine test_ConstrainRecruitNumber_FreezingBoundary(this) + use FatesUnitTestPRTParametersMod, only : init_mock_prt_params + use EDPftvarconMockMod, only : init_mock_edpftvarcon + use EDTypesMockMod, only : init_mock_site_patch_cohort + + class(TestPlantHydraulicsCohort), intent(inout) :: this + type(ed_site_type), pointer :: csite => null() + type(fates_patch_type), pointer :: cpatch => null() + type(fates_cohort_type), pointer :: ccohort => null() + type(bc_in_type) :: bc_in + ! Exact freezing temperature boundary mean_temp = 273.15 K (0 deg C) + real(r8) :: mean_temp = 273.15_r8 + + call init_mock_prt_params() + call init_mock_edpftvarcon() + call init_mock_site_patch_cohort(csite, cpatch, ccohort) + ccohort%pft = 1 + + csite%si_hydr%nlevrhiz = 1 + call csite%si_hydr%InitHydrSite(1, 1, 1, 1) + call ccohort%co_hydr%AllocateHydrCohortArrays(1) + + allocate(wrf_type_vg :: csite%si_hydr%wrf_soil(1)%p) + select type (p => csite%si_hydr%wrf_soil(1)%p) + type is (wrf_type_vg) + call p%set_wrf_param([0.001_r8, 1.5_r8, 0.5_r8, 0.6_r8, 0.1_r8]) + end select + + csite%si_hydr%v_shell(1,1) = 1.0_r8 + csite%si_hydr%h2osoi_liqvol_shell(1,1) = 0.5_r8 + + ccohort%co_hydr%th_ag(:) = 0.5_r8 + ccohort%co_hydr%v_ag(:) = 1.0_r8 + ccohort%co_hydr%th_troot = 0.5_r8 + ccohort%co_hydr%v_troot = 1.0_r8 + ccohort%co_hydr%th_aroot(:) = 0.5_r8 + ccohort%co_hydr%v_aroot_layer(:) = 1.0_r8 + ccohort%co_hydr%l_aroot_layer(:) = 1.0_r8 + + ! Minimum soil matric potential smpmin_si = -1.5e8_r8 ensures low residual water threshold, allowing non-zero recruitable density before temperature check + bc_in%smpmin_si = -1.5e8_r8 + ccohort%n = 10.0_r8 + + call cpatch%InitLitter(1, 1) + + ! Execute recruitment constraint routine at exact freezing boundary (273.15 K) + call ConstrainRecruitNumber(csite, ccohort, cpatch, bc_in, mean_temp) + + ! Kills Mutation 5 by verifying recruitment is zeroed out at exact freezing point (273.15 K) + @assertEqual(0.0_r8, ccohort%n, tolerance=1.e-10_r8) + + if(allocated(ccohort%co_hydr%l_aroot_layer)) deallocate(ccohort%co_hydr%l_aroot_layer) + end subroutine test_ConstrainRecruitNumber_FreezingBoundary + +end module FatesPlantHydraulicsMod_cohort_test diff --git a/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_drive_test.pf b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_drive_test.pf new file mode 100644 index 0000000000..862668ade6 --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_drive_test.pf @@ -0,0 +1,497 @@ +! Summary: Unit tests for FatesPlantHydraulicsMod timestep driver routines, diagnostics, and vegetation water storage balance. +!========================================================================================= +! Tests for FatesPlantHydraulicsMod timestep-level driver routines and update functions +!========================================================================================= +module FatesPlantHydraulicsMod_drive_test + use FatesPlantHydraulicsMod, only : hydraulics_drive, BTranForHLMDiagnosticsFromCohortHydr, & + BTranForHLMDiagnostics_explicit, AccumulateMortalityWaterStorage, & + AccumulateMortalityWater_explicit, RecruitWaterStorage, UpdateH2OVeg, & + use_ed_planthydraulics, InitHydroGlobals + use EDTypesMod, only : ed_site_type + use FatesPatchMod, only : fates_patch_type + use FatesCohortMod, only : fates_cohort_type + use FatesHydraulicsMemMod, only : ed_site_hydr_type, ed_cohort_hydr_type + use FatesInterfaceTypesMod, only : bc_in_type, bc_out_type, hlm_use_planthydro + use FatesConstantsMod, only : r8 => fates_r8 + use funit + + implicit none + + @TestCase + type, extends(TestCase) :: TestPlantHydraulicsDrive + contains + end type TestPlantHydraulicsDrive + +contains + + ! Summary: Driver test for hydraulics_drive step integration using full non-linear implicit plant-soil hydraulic solver. + @Test + subroutine test_hydraulics_drive(this) + use FatesUnitTestPRTParametersMod, only : init_mock_prt_params + use EDPftvarconMockMod, only : init_mock_edpftvarcon + use EDTypesMockMod, only : init_mock_site_patch_cohort + use FatesHydroWTFMod, only : wrf_type_vg, wkf_type_vg + + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: sites(1) + type(ed_site_type), pointer :: csite => null() + type(fates_patch_type), pointer :: cpatch => null() + type(fates_cohort_type), pointer :: ccohort => null() + type(bc_in_type) :: bc_in(1) + type(bc_out_type) :: bc_out(1) + ! Solver output for non-linear system under 1800s timestep integration + real(r8) :: expected_psi_ag + real(r8), parameter :: tol = 1.e-4_r8 + + ! Initialize mock global structures and hydraulics options + call init_mock_prt_params() + call init_mock_edpftvarcon() + call InitHydroGlobals() + + csite => sites(1) + call init_mock_site_patch_cohort(csite, cpatch, ccohort) + ccohort%pft = 1 + + ! Global hydraulic execution flags: enable plant hydraulics model + use_ed_planthydraulics = 1 + hlm_use_planthydro = 1 + + ! Single rhizosphere soil layer setup (nlevrhiz = 1) + csite%si_hydr%nlevrhiz = 1 + call csite%si_hydr%InitHydrSite(1, 1, 1, 1) + + ! Allocate van Genuchten soil water retention and hydraulic conductivity models + allocate(csite%si_hydr%wrf_soil(1)) + allocate(csite%si_hydr%wkf_soil(1)) + allocate(wrf_type_vg :: csite%si_hydr%wrf_soil(1)%p) + select type (p => csite%si_hydr%wrf_soil(1)%p) + type is (wrf_type_vg) + ! Soil parameters: alpha=0.001 1/mm, n=1.5, m=0.5, th_s=0.6, th_r=0.1 + call p%set_wrf_param([0.001_r8, 1.5_r8, 0.5_r8, 0.6_r8, 0.1_r8]) + end select + + allocate(wkf_type_vg :: csite%si_hydr%wkf_soil(1)%p) + select type (k => csite%si_hydr%wkf_soil(1)%p) + type is (wkf_type_vg) + ! Soil conductivity parameters: ksat=0.001 kg/m2/s, n=1.5, m=0.5, th_s=0.6, th_r=0.1, l=0.5 + call k%set_wkf_param([0.001_r8, 1.5_r8, 0.5_r8, 0.6_r8, 0.1_r8, 0.5_r8]) + end select + + ! Soil shell moisture state initialization: + ! h2osoi_liqvol_shell [m3/m3] (moist soil condition), v_shell [m3], r_node_shell [m], r_out_shell [m] + csite%si_hydr%h2osoi_liqvol_shell(1,:) = 0.5_r8 + csite%si_hydr%v_shell(1,:) = 0.25_r8 + csite%si_hydr%r_node_shell(1,:) = 0.1_r8 + csite%si_hydr%r_out_shell(1,:) = 0.2_r8 + csite%si_hydr%l_aroot_layer(1) = 1.0_r8 + csite%si_hydr%map_r2s(1,1) = 1 + csite%si_hydr%map_r2s(1,2) = 1 + + ! Rhizosphere layer geometry: thickness dz [m], interface depth zi [m] + csite%si_hydr%zi_rhiz(1) = 0.4_r8 + csite%si_hydr%dz_rhiz(1) = 0.4_r8 + + ! Host Land Model (HLM) input boundaries: liquid water [m3/m3], layer thickness [m] + allocate(bc_in(1)%h2o_liq_sisl(1)) + bc_in(1)%h2o_liq_sisl(1) = 0.5_r8 + bc_in(1)%nlevsoil = 1 + allocate(bc_in(1)%zi_sisl(1)) + bc_in(1)%zi_sisl(1) = 0.5_r8 + allocate(bc_in(1)%dz_sisl(1)) + bc_in(1)%dz_sisl(1) = 0.5_r8 + + ! Transpiration flux qflx_transp_pa [kg/m2/s] set to zero (pure equilibration step) + allocate(bc_in(1)%qflx_transp_pa(1)) + bc_in(1)%qflx_transp_pa(1) = 0.0_r8 + + ! Output boundary array for soil-to-root flux [kg/m2/s] + allocate(bc_out(1)%qflx_soil2root_sisl(1)) + bc_out(1)%qflx_soil2root_sisl(1) = 0.0_r8 + + ! Cohort hydraulic state initialization: relative water content theta [m3/m3] + call ccohort%co_hydr%AllocateHydrCohortArrays(1) + ccohort%co_hydr%th_ag(:) = 0.5_r8 + ccohort%co_hydr%th_troot = 0.5_r8 + ccohort%co_hydr%th_aroot(:) = 0.5_r8 + ccohort%co_hydr%v_aroot_layer(:) = 0.25_r8 + ccohort%co_hydr%v_aroot_layer_init(:) = 0.25_r8 + ccohort%co_hydr%l_aroot_layer(1) = 1.0_r8 + bc_in(1)%smpmin_si = -1500.0_r8 + + call cpatch%InitLitter(1, 1) + + ! Complex-to-mock driver execution: call full non-linear implicit hydraulic solver over timestep dt [s] + ! [COMPLEX-TO-MOCK]: System solves non-linear matrix equation for plant-soil water potentials + call hydraulics_drive(1, sites, bc_in, bc_out, 1800.0_r8) + + ! Expected converged aboveground water potential psi_ag [MPa] + expected_psi_ag = -3.1574220023403204_r8 + @assertEqual(expected_psi_ag, ccohort%co_hydr%psi_ag(1), tolerance=tol) + + ! Clean up allocated structures + deallocate(bc_in(1)%h2o_liq_sisl) + deallocate(bc_in(1)%qflx_transp_pa) + call ccohort%co_hydr%DeallocateHydrCohortArrays() + end subroutine test_hydraulics_drive + + + ! Summary: Tests BTran transpiration reduction factor calculation for host land model diagnostics. + @Test + subroutine test_BTranForHLMDiagnosticsFromCohortHydr(this) + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: sites(1) + type(fates_patch_type), target :: cpatch + type(bc_out_type) :: bc_out(1) + + ! Site and patch topology setup + sites(1)%oldest_patch => cpatch + cpatch%younger => null() + cpatch%tallest => null() + cpatch%patchno = 1 + + ! Diagnostic output array allocation + allocate(bc_out(1)%btran_pa(1)) + ! 999.0_r8 used as a distinct sentinel mock value to ensure function overwrites output + bc_out(1)%btran_pa(1) = 999.0_r8 + + ! Execute diagnostic computation routine with zero cohorts present + call BTranForHLMDiagnosticsFromCohortHydr(1, sites, bc_out) + + ! Expected diagnostic btran factor = 0.0 for empty patch + @assertEqual(0.0_r8, bc_out(1)%btran_pa(1)) + + deallocate(bc_out(1)%btran_pa) + end subroutine test_BTranForHLMDiagnosticsFromCohortHydr + + ! Summary: Tests BTran transpiration reduction factor biomass-weighted aggregation across multiple cohorts. + @Test + subroutine test_BTranForHLMDiagnostics_MultiCohortWeighted(this) + use EDTypesMockMod, only : init_mock_site_patch_cohort, InitPRTObject + use PRTGenericMod, only : carbon12_element, leaf_organ, sapw_organ, fnrt_organ, prt_global + + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), pointer :: csite => null() + type(fates_patch_type), pointer :: cpatch => null() + type(fates_cohort_type), pointer :: ccohort1 => null() + type(fates_cohort_type), pointer :: ccohort2 => null() + type(ed_site_type), target :: sites(1) + type(bc_out_type) :: bc_out(1) + integer :: i_leaf, i_sapw, i_fnrt + real(r8) :: expected_btran_pa + real(r8), parameter :: tol = 1.e-6_r8 + + call init_mock_site_patch_cohort(csite, cpatch, ccohort1) + + allocate(ccohort2) + allocate(ccohort2%co_hydr) + call InitPRTObject(ccohort2%prt) + + sites(1) = csite + sites(1)%oldest_patch => cpatch + cpatch%tallest => ccohort1 + ccohort1%shorter => ccohort2 + ccohort2%shorter => null() + + i_leaf = prt_global%sp_organ_map(leaf_organ, carbon12_element) + i_sapw = prt_global%sp_organ_map(sapw_organ, carbon12_element) + i_fnrt = prt_global%sp_organ_map(fnrt_organ, carbon12_element) + + ! Cohort 1: btran = 0.8, n = 1.0, leaf = 2.0, sapw = 3.0, fnrt = 5.0 (total organ mass = 10.0 kg/indiv) + ccohort1%co_hydr%btran = 0.8_r8 + ccohort1%n = 1.0_r8 + ccohort1%prt%variables(i_leaf)%val(1) = 2.0_r8 + ccohort1%prt%variables(i_sapw)%val(1) = 3.0_r8 + ccohort1%prt%variables(i_fnrt)%val(1) = 5.0_r8 + + ! Cohort 2: btran = 0.2, n = 2.0, leaf = 1.0, sapw = 2.0, fnrt = 2.0 (total organ mass = 5.0 kg/indiv) + ccohort2%co_hydr%btran = 0.2_r8 + ccohort2%n = 2.0_r8 + ccohort2%prt%variables(i_leaf)%val(1) = 1.0_r8 + ccohort2%prt%variables(i_sapw)%val(1) = 2.0_r8 + ccohort2%prt%variables(i_fnrt)%val(1) = 2.0_r8 + + allocate(bc_out(1)%btran_pa(1)) + bc_out(1)%btran_pa(1) = 0.0_r8 + + ! Execute diagnostic computation routine with two active cohorts + call BTranForHLMDiagnosticsFromCohortHydr(1, sites, bc_out) + + ! Analytical solution including fine-root carbon (fnrt): + ! Cohort 1 biomass = 10.0 kg, weighted btran = 0.8 * 10.0 * 1.0 = 8.0 + ! Cohort 2 biomass = 5.0 kg, weighted btran = 0.2 * 5.0 * 2.0 = 2.0 + ! Total live patch biomass = 10.0*1.0 + 5.0*2.0 = 20.0 kg/m2 + ! Expected btran_pa = (8.0 + 2.0) / 20.0 = 0.50 + expected_btran_pa = 0.50_r8 + @assertEqual(expected_btran_pa, bc_out(1)%btran_pa(1), tolerance=tol) + + deallocate(bc_out(1)%btran_pa) + deallocate(ccohort2%co_hydr) + deallocate(ccohort2) + end subroutine test_BTranForHLMDiagnostics_MultiCohortWeighted + + ! Summary: Tests BTranForHLMDiagnostics_explicit pure calculation function with explicit array inputs. + @Test + subroutine test_BTranForHLMDiagnostics_explicit(this) + class(TestPlantHydraulicsDrive), intent(inout) :: this + real(r8) :: btran_cohort(2) + real(r8) :: balive_cohort(2) + real(r8) :: n_indiv(2) + real(r8) :: result_btran + real(r8), parameter :: tol = 1.e-6_r8 + + ! Cohort 1: btran = 0.8, balive = 10.0, n = 1.0 => weighted btran = 8.0 + btran_cohort(1) = 0.8_r8 + balive_cohort(1) = 10.0_r8 + n_indiv(1) = 1.0_r8 + + ! Cohort 2: btran = 0.2, balive = 5.0, n = 2.0 => weighted btran = 2.0 + btran_cohort(2) = 0.2_r8 + balive_cohort(2) = 5.0_r8 + n_indiv(2) = 2.0_r8 + + ! Execute explicit subroutine (total biomass = 10*1 + 5*2 = 20, expected btran = (8 + 2)/20 = 0.50) + call BTranForHLMDiagnostics_explicit(2, btran_cohort, balive_cohort, n_indiv, result_btran) + @assertEqual(0.50_r8, result_btran, tolerance=tol) + + ! Zero biomass edge case returns 0.0 + balive_cohort(:) = 0.0_r8 + call BTranForHLMDiagnostics_explicit(2, btran_cohort, balive_cohort, n_indiv, result_btran) + @assertEqual(0.0_r8, result_btran, tolerance=tol) + end subroutine test_BTranForHLMDiagnostics_explicit + + ! Summary: Tests water mass balance conservation when accumulating water storage from cohort mortality. + @Test + subroutine test_AccumulateMortalityWaterStorage(this) + use FatesConstantsMod, only: denh2o => dens_fresh_liquid_water + use EDTypesMod, only: AREA_INV + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: csite + type(fates_cohort_type), target :: ccohort + type(ed_site_hydr_type), target :: csite_hydr + type(ed_cohort_hydr_type), target :: co_hydr + ! Cohort mortality density delta_n = 1.0 indiv/m2 + real(r8) :: delta_n = 1.0_r8 + real(r8) :: expected_delta_w + + csite%si_hydr => csite_hydr + ccohort%co_hydr => co_hydr + + ! Compartment volumes: v_ag = 1.0 m3, v_troot = 2.0 m3, v_aroot = 3.0 m3 + ! Compartment water contents: theta_ag = 0.5, theta_troot = 0.5, theta_aroot = 0.5 m3/m3 + allocate(co_hydr%th_aroot(1), co_hydr%v_aroot_layer(1)) + co_hydr%th_ag(:) = 0.5_r8 + co_hydr%v_ag(:) = 1.0_r8 + co_hydr%th_troot = 0.5_r8 + co_hydr%v_troot = 2.0_r8 + co_hydr%th_aroot(1) = 0.5_r8 + co_hydr%v_aroot_layer(1) = 3.0_r8 + + ! Initial dead vegetation water pool = 0.0 kg/m2, initial live vegetation water pool = 100.0 kg/m2 + csite_hydr%h2oveg_dead = 0.0_r8 + csite_hydr%h2oveg = 100.0_r8 + + ! Execute mortality water accumulation routine + call AccumulateMortalityWaterStorage(csite, ccohort, delta_n) + + ! Mass conservation analytical calculation: + ! delta_w = sum(theta_i * V_i) * denh2o * delta_n * AREA_INV + ! delta_w = (0.5*1.0 + 0.5*2.0 + 0.5*3.0) * 1000.0 * 1.0 * 1.0 = 3.0 * 1000.0 = 3000.0 kg/m2 + expected_delta_w = (sum(co_hydr%th_ag(:)*co_hydr%v_ag(:)) + & + co_hydr%th_troot*co_hydr%v_troot + & + sum(co_hydr%th_aroot(:)*co_hydr%v_aroot_layer(:))) * & + denh2o * delta_n * AREA_INV + + ! Verify exact mass balance transfer from live to dead pool + @assertEqual(expected_delta_w, csite_hydr%h2oveg_dead) + @assertEqual(100.0_r8 - expected_delta_w, csite_hydr%h2oveg) + + ! Invariant A: Total water mass conservation during mortality transfer (initial_total == final_total) + @assertEqual(100.0_r8, csite_hydr%h2oveg + csite_hydr%h2oveg_dead, tolerance=1.e-8_r8) + + deallocate(co_hydr%th_aroot, co_hydr%v_aroot_layer) + end subroutine test_AccumulateMortalityWaterStorage + + ! Summary: Tests AccumulateMortalityWaterStorage absorbing root vector dot product across multiple soil layers. + @Test + subroutine test_AccumulateMortalityWaterStorage_MultiLayer(this) + use FatesConstantsMod, only: denh2o => dens_fresh_liquid_water + use EDTypesMod, only: AREA_INV + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: csite + type(fates_cohort_type), target :: ccohort + type(ed_site_hydr_type), target :: csite_hydr + type(ed_cohort_hydr_type), target :: co_hydr + real(r8) :: delta_n = 1.0_r8 + real(r8) :: expected_delta_w + real(r8), parameter :: tol = 1.e-6_r8 + + csite%si_hydr => csite_hydr + ccohort%co_hydr => co_hydr + + ! Multi-layer soil setup (nlevrhiz = 2): + ! th_aroot = [0.2, 0.8] m3/m3, v_aroot_layer = [3.0, 1.0] m3 + ! Dot product = 0.2*3.0 + 0.8*1.0 = 1.4 m3 + ! Product of sums = (1.0)*(4.0) = 4.0 m3 (Mutation 7 fails) + allocate(co_hydr%th_aroot(2), co_hydr%v_aroot_layer(2)) + co_hydr%th_ag(:) = 0.5_r8 + co_hydr%v_ag(:) = 1.0_r8 + co_hydr%th_troot = 0.5_r8 + co_hydr%v_troot = 2.0_r8 + + co_hydr%th_aroot(1) = 0.2_r8 + co_hydr%th_aroot(2) = 0.8_r8 + co_hydr%v_aroot_layer(1) = 3.0_r8 + co_hydr%v_aroot_layer(2) = 1.0_r8 + + csite_hydr%h2oveg_dead = 0.0_r8 + csite_hydr%h2oveg = 10000.0_r8 + + call AccumulateMortalityWaterStorage(csite, ccohort, delta_n) + + ! Analytical expected delta_w: + ! (0.5*1.0 + 0.5*2.0 + (0.2*3.0 + 0.8*1.0)) * 1000.0 * 1.0 * 0.0001 = 3.4 * 0.1 = 0.34 kg/m2 + expected_delta_w = 0.34_r8 + @assertEqual(expected_delta_w, csite_hydr%h2oveg_dead, tolerance=tol) + + deallocate(co_hydr%th_aroot, co_hydr%v_aroot_layer) + end subroutine test_AccumulateMortalityWaterStorage_MultiLayer + + ! Summary: Tests AccumulateMortalityWater_explicit pure calculation function with explicit array inputs. + @Test + subroutine test_AccumulateMortalityWater_explicit(this) + use FatesConstantsMod, only: denh2o => dens_fresh_liquid_water + use EDTypesMod, only: AREA_INV + class(TestPlantHydraulicsDrive), intent(inout) :: this + real(r8) :: delta_n + real(r8) :: th_ag(1), v_ag(1) + real(r8) :: th_troot, v_troot + real(r8) :: th_aroot(2), v_aroot_layer(2) + real(r8) :: delta_w + real(r8), parameter :: tol = 1.e-6_r8 + + delta_n = 1.0_r8 + th_ag(1) = 0.5_r8 + v_ag(1) = 1.0_r8 + th_troot = 0.5_r8 + v_troot = 2.0_r8 + th_aroot(1) = 0.2_r8; v_aroot_layer(1) = 3.0_r8 + th_aroot(2) = 0.8_r8; v_aroot_layer(2) = 1.0_r8 + + call AccumulateMortalityWater_explicit(delta_n, th_ag, v_ag, th_troot, v_troot, & + th_aroot, v_aroot_layer, delta_w) + ! Analytical expected delta_w: + ! (0.5*1.0 + 0.5*2.0 + (0.2*3.0 + 0.8*1.0)) * 1000.0 * 1.0 * 0.0001 = 2.9 * 0.1 = 0.29 kg/m2 + @assertEqual(0.29_r8, delta_w, tolerance=tol) + end subroutine test_AccumulateMortalityWater_explicit + + + ! Summary: Tests recruit water storage pool initialization during recruitment events. + @Test + subroutine test_RecruitWaterStorage(this) + use FatesConstantsMod, only: denh2o => dens_fresh_liquid_water + use EDTypesMod, only: AREA_INV + use FatesInterfaceTypesMod, only: hlm_use_planthydro + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: sites(1) + type(fates_patch_type), target :: cpatch + type(fates_cohort_type), target :: ccohort + type(ed_site_hydr_type), target :: csite_hydr + type(ed_cohort_hydr_type), target :: co_hydr + real(r8) :: expected_recruit_w + + hlm_use_planthydro = 1 + sites(1)%si_hydr => csite_hydr + sites(1)%oldest_patch => cpatch + cpatch%tallest => ccohort + cpatch%younger => null() + ccohort%shorter => null() + ccohort%co_hydr => co_hydr + + ! Cohort recruited status = true, recruit density n = 2.0 indiv/m2 + co_hydr%is_newly_recruited = .true. + ccohort%n = 2.0_r8 + + ! Compartment volumes: v_ag = 1.0 m3, v_troot = 2.0 m3, v_aroot = 3.0 m3 + allocate(co_hydr%th_aroot(1), co_hydr%v_aroot_layer(1)) + co_hydr%th_ag(:) = 0.5_r8 + co_hydr%v_ag(:) = 1.0_r8 + co_hydr%th_troot = 0.5_r8 + co_hydr%v_troot = 2.0_r8 + co_hydr%th_aroot(1) = 0.5_r8 + co_hydr%v_aroot_layer(1) = 3.0_r8 + + csite_hydr%h2oveg_recruit = 0.0_r8 + + ! Execute recruit water storage update + call RecruitWaterStorage(1, sites) + + ! Analytical solution: recruit_w = (0.5*1 + 0.5*2 + 0.5*3) * 1000 * 2.0 * 1.0 = 6000.0 kg/m2 + expected_recruit_w = (sum(co_hydr%th_ag(:)*co_hydr%v_ag(:)) + & + co_hydr%th_troot*co_hydr%v_troot + & + sum(co_hydr%th_aroot(:)*co_hydr%v_aroot_layer(:))) * & + denh2o * ccohort%n * AREA_INV + + ! Verify exact recruit water storage value + @assertEqual(expected_recruit_w, csite_hydr%h2oveg_recruit) + + deallocate(co_hydr%th_aroot, co_hydr%v_aroot_layer) + end subroutine test_RecruitWaterStorage + + ! Summary: Tests total vegetation water storage summation across all live plant cohorts. + @Test + subroutine test_UpdateH2OVeg(this) + use FatesConstantsMod, only: denh2o => dens_fresh_liquid_water + use EDTypesMod, only: AREA_INV + use FatesInterfaceTypesMod, only: hlm_use_planthydro + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: csite + type(fates_patch_type), target :: cpatch + type(fates_cohort_type), target :: ccohort + type(ed_site_hydr_type), target :: csite_hydr + type(ed_cohort_hydr_type), target :: co_hydr + type(bc_out_type) :: bc_out + real(r8) :: expected_h2oveg + + hlm_use_planthydro = 1 + csite%si_hydr => csite_hydr + csite%oldest_patch => cpatch + cpatch%tallest => ccohort + cpatch%younger => null() + ccohort%shorter => null() + ccohort%co_hydr => co_hydr + + ! Cohort properties: density n = 10.0 indiv/m2 + co_hydr%is_newly_recruited = .false. + allocate(co_hydr%th_aroot(1), co_hydr%v_aroot_layer(1)) + co_hydr%th_ag(:) = 0.5_r8 + co_hydr%v_ag(:) = 1.0_r8 + co_hydr%th_troot = 0.5_r8 + co_hydr%v_troot = 2.0_r8 + co_hydr%th_aroot(1) = 0.5_r8 + co_hydr%v_aroot_layer(1) = 3.0_r8 + ccohort%n = 10.0_r8 + + ! Zero error terms before diagnostic summation + csite_hydr%h2oveg_dead = 0.0_r8 + csite_hydr%h2oveg_growturn_err = 0.0_r8 + csite_hydr%h2oveg_hydro_err = 0.0_r8 + csite_hydr%trans_err = 0.0_r8 + csite_hydr%h2oveg = 0.0_r8 + + ! Execute vegetation water diagnostic update + call UpdateH2OVeg(csite, bc_out) + + ! Analytical expected total water storage = 30,000 kg/m2 + expected_h2oveg = (sum(co_hydr%th_ag(:)*co_hydr%v_ag(:)) + & + co_hydr%th_troot*co_hydr%v_troot + & + sum(co_hydr%th_aroot(:)*co_hydr%v_aroot_layer(:))) * & + denh2o * ccohort%n * AREA_INV + + ! Verify calculated plant stored water matches analytical total + @assertEqual(expected_h2oveg, bc_out%plant_stored_h2o_si) + + deallocate(co_hydr%th_aroot, co_hydr%v_aroot_layer) + end subroutine test_UpdateH2OVeg + +end module FatesPlantHydraulicsMod_drive_test + diff --git a/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_geom_test.pf b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_geom_test.pf new file mode 100644 index 0000000000..4a24751aa1 --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_geom_test.pf @@ -0,0 +1,253 @@ +! Summary: Unit tests for FatesPlantHydraulicsMod geometric calculations (nodes, length, volume, and SumBetweenDepths depth-weighted integration). +!========================================================================================= +! Tests for FatesPlantHydraulicsMod geometric calculations (nodes, length, and volumes) +!========================================================================================= +module FatesPlantHydraulicsMod_geom_test + use FatesPlantHydraulicsMod, only : UpdatePlantHydrNodes, UpdatePlantHydrLenVol, SavePreviousCompartmentVolumes, SumBetweenDepths + use FatesCohortMod, only : fates_cohort_type + use FatesHydraulicsMemMod, only : ed_site_hydr_type, ed_cohort_hydr_type + use FatesParameterDerivedMod, only : param_derived + use FatesConstantsMod, only : r8 => fates_r8 + use PRTParametersMod, only : prt_params + use EDPftvarcon, only : EDPftvarcon_inst + use EDCohortDynamicsMod, only : InitPRTObject + use PRTGenericMod, only : SetState, leaf_organ, fnrt_organ, sapw_organ, struct_organ, carbon12_element, num_elements, element_list + use FatesInterfaceTypesMod, only : hlm_parteh_mode + use FatesUnitTestPRTParametersMod, only : init_mock_prt_params + use funit + + implicit none + + @TestCase + type, extends(TestCase) :: TestPlantHydraulicsGeom + contains + end type TestPlantHydraulicsGeom + +contains + + ! Summary: Setup mock parameters for FATES allometry and fine root geometry. + subroutine init_mock_params() + call init_mock_prt_params() + + if (.not. allocated(param_derived%branch_frac)) allocate(param_derived%branch_frac(1)) + ! Branch fraction set to 0.5 (50% branch wood / 50% main stem biomass allocation split) + param_derived%branch_frac(1) = 0.5_r8 + + if (.not. allocated(EDPftvarcon_inst%hydr_srl)) allocate(EDPftvarcon_inst%hydr_srl(1)) + if (.not. allocated(EDPftvarcon_inst%hydr_rs2)) allocate(EDPftvarcon_inst%hydr_rs2(1)) + ! Specific root length = 10.0 m/g (represents fine root absorptive surface per unit carbon mass) + EDPftvarcon_inst%hydr_srl(1) = 10.0_r8 + ! Fine root radius = 0.001 m = 1 mm (standard fine absorbing root radius) + EDPftvarcon_inst%hydr_rs2(1) = 0.001_r8 + + ! Mock num_elements and element_list for carbon allocation tracking + num_elements = 1 + if (.not. allocated(element_list)) allocate(element_list(1)) + element_list(1) = carbon12_element + hlm_parteh_mode = 1 ! carbon_only mode + end subroutine init_mock_params + + ! Summary: Tests height and depth node allocation across aboveground stem compartments and transport roots. + @Test + subroutine test_UpdatePlantHydrNodes(this) + class(TestPlantHydraulicsGeom), intent(inout) :: this + type(fates_cohort_type), target :: ccohort + type(ed_site_hydr_type) :: csite_hydr + ! Plant height = 10.0 m selected as representative mature canopy tree height + real(r8) :: plant_height = 10.0_r8 + integer :: ft = 1 + real(r8), parameter :: tol = 1.e-6_r8 + + call init_mock_params() + + ! DBH = 20.0 cm represents 20cm trunk diameter matching 10m tree height allometry + ccohort%dbh = 20.0_r8 + ccohort%pft = ft + allocate(ccohort%co_hydr) + ! Allocate 2 aboveground compartments (index 1 = leaf pool, index 2 = stem sapwood pool) + call ccohort%co_hydr%AllocateHydrCohortArrays(2) + + ! Soil layer interfaces: zi = [0.1m, 1.0m] (topsoil 10cm, subsoil 90cm), dz = [0.1m, 0.9m] + allocate(csite_hydr%zi_rhiz(2)) + csite_hydr%nlevrhiz = 2 + csite_hydr%zi_rhiz(1) = 0.1_r8 + csite_hydr%zi_rhiz(2) = 1.0_r8 + allocate(csite_hydr%dz_rhiz(2)) + csite_hydr%dz_rhiz(1) = 0.1_r8 + csite_hydr%dz_rhiz(2) = 0.9_r8 + + ! Execute node position calculation + call UpdatePlantHydrNodes(ccohort, ft, plant_height, csite_hydr) + + ! Analytical / numerical solution derivation: + ! Crown depth = min(height, 0.1) = 0.1 m => z_stem = height - crown_depth = 9.9 m + ! Leaf compartment (index 1): upper = 10.0 m, node = 9.95 m (height - 0.5*crown_depth), lower = 9.9 m + ! Stem compartment (index 2): upper = 9.9 m, node = 4.95 m (0.5*z_stem), lower = 0.0 m + ! Transport root node depth = -0.5 * min(zi_rhiz(2), z_fr) = -0.5 * 0.248974 = -0.124487 m + @assertEqual(10.0_r8, ccohort%co_hydr%z_upper_ag(1), tolerance=tol) + @assertEqual(9.95_r8, ccohort%co_hydr%z_node_ag(1), tolerance=tol) + @assertEqual(9.9_r8, ccohort%co_hydr%z_lower_ag(1), tolerance=tol) + + @assertEqual(9.9_r8, ccohort%co_hydr%z_upper_ag(2), tolerance=tol) + @assertEqual(4.95_r8, ccohort%co_hydr%z_node_ag(2), tolerance=tol) + @assertEqual(0.0_r8, ccohort%co_hydr%z_lower_ag(2), tolerance=tol) + + @assertEqual(-0.12448708730516955_r8, ccohort%co_hydr%z_node_troot, tolerance=tol) + + ! Invariant C: Spatial node depth hierarchy (z_lower < z_node < z_upper and z_node_troot < 0.0) + @assertTrue(all(ccohort%co_hydr%z_lower_ag < ccohort%co_hydr%z_node_ag)) + @assertTrue(all(ccohort%co_hydr%z_node_ag < ccohort%co_hydr%z_upper_ag)) + @assertTrue(ccohort%co_hydr%z_node_troot < 0.0_r8) + + call ccohort%co_hydr%DeallocateHydrCohortArrays() + deallocate(ccohort%co_hydr) + deallocate(csite_hydr%zi_rhiz) + deallocate(csite_hydr%dz_rhiz) + end subroutine test_UpdatePlantHydrNodes + + + ! Summary: Tests plant organ volume and fine-root length distribution calculations from biomass pools. + @Test + subroutine test_UpdatePlantHydrLenVol(this) + class(TestPlantHydraulicsGeom), intent(inout) :: this + type(fates_cohort_type), target :: ccohort + type(ed_site_hydr_type) :: csite_hydr + real(r8), parameter :: tol = 1.e-5_r8 + + call init_mock_params() + + ! Cohort allometry inputs: DBH = 20.0 cm, height = 10.0 m, canopy_trim = 1.0 (un-trimmed full canopy), efstem = 1.0 (full stem density) + ccohort%dbh = 20.0_r8 + ccohort%pft = 1 + ccohort%height = 10.0_r8 + ccohort%crowndamage = 0 + ccohort%canopy_trim = 1.0_r8 + ccohort%efstem_coh = 1.0_r8 + ccohort%size_class = 1 + allocate(ccohort%co_hydr) + call ccohort%co_hydr%AllocateHydrCohortArrays(2) + + ! Soil layer interfaces: zi = [0.1m, 1.0m], dz = [0.1m, 0.9m] + allocate(csite_hydr%zi_rhiz(2)) + allocate(csite_hydr%dz_rhiz(2)) + csite_hydr%nlevrhiz = 2 + csite_hydr%zi_rhiz(1) = 0.1_r8 + csite_hydr%zi_rhiz(2) = 1.0_r8 + csite_hydr%dz_rhiz(1) = 0.1_r8 + csite_hydr%dz_rhiz(2) = 0.9_r8 + + ! Execute volume and length update routine with test organ carbon pools: + ! leaf_c = 1.0 kgC, sapw_c = 5.0 kgC, fnrt_c = 0.5 kgC, struct_c = 10.0 kgC + call UpdatePlantHydrLenVol(ccohort%co_hydr, csite_hydr, & + ccohort%pft, ccohort%dbh, ccohort%height, & + ccohort%crowndamage, ccohort%canopy_trim, ccohort%efstem_coh, & + ccohort%size_class, & + 1.0_r8, & + 5.0_r8, & + 0.5_r8, & + 10.0_r8) + + ! Analytical / allometric target values: + ! Leaf volume (v_ag(1)): 0.00329127901571 m3 (calculated from denleaf and c2b) + ! Stem sapwood volume (v_ag(2)): 0.0000675 m3 (a_sapwood * z_stem) + ! Transport root volume (v_troot): 0.02770796326795 m3 + @assertEqual(0.0032912790157100975_r8, ccohort%co_hydr%v_ag(1), tolerance=tol) + @assertEqual(0.0000675_r8, ccohort%co_hydr%v_ag(2), tolerance=tol) + @assertEqual(0.02770796326795_r8, ccohort%co_hydr%v_troot, tolerance=tol) + + ! Total fine root length = 10,000 m distributed across soil layers + @assertEqual(10000.0_r8, sum(ccohort%co_hydr%l_aroot_layer), tolerance=1.0_r8) + ! Total fine root volume = pi * r^2 * L_rhiz = pi * (0.001)^2 * 8820.0 = 0.02770796326795 m3 + @assertEqual(0.02770796326795_r8, sum(ccohort%co_hydr%v_aroot_layer), tolerance=1.e-4_r8) + + call ccohort%co_hydr%DeallocateHydrCohortArrays() + deallocate(ccohort%co_hydr) + deallocate(csite_hydr%zi_rhiz) + deallocate(csite_hydr%dz_rhiz) + end subroutine test_UpdatePlantHydrLenVol + + ! Summary: Tests copying current compartment volumes into previous timestep volume arrays for size-dependent updates. + @Test + subroutine test_SavePreviousCompartmentVolumes(this) + class(TestPlantHydraulicsGeom), intent(inout) :: this + type(ed_cohort_hydr_type) :: ccohort_hydr + + allocate(ccohort_hydr%v_aroot_layer(2)) + allocate(ccohort_hydr%v_aroot_layer_init(2)) + + ! Current volume inputs: v_ag(1) = 5.0 m3, v_troot = 10.0 m3, v_aroot = [1.0 m3, 2.0 m3] + ccohort_hydr%v_ag(1) = 5.0_r8 + ccohort_hydr%v_troot = 10.0_r8 + ccohort_hydr%v_aroot_layer(1) = 1.0_r8 + ccohort_hydr%v_aroot_layer(2) = 2.0_r8 + + ! Initial target arrays set to zero + ccohort_hydr%v_ag_init(:) = 0.0_r8 + ccohort_hydr%v_troot_init = 0.0_r8 + ccohort_hydr%v_aroot_layer_init(:) = 0.0_r8 + + ! Execute save routine + call SavePreviousCompartmentVolumes(ccohort_hydr) + + ! Verify exact copy of current volumes to _init variables + @assertEqual(5.0_r8, ccohort_hydr%v_ag_init(1)) + @assertEqual(10.0_r8, ccohort_hydr%v_troot_init) + @assertEqual(1.0_r8, ccohort_hydr%v_aroot_layer_init(1)) + @assertEqual(2.0_r8, ccohort_hydr%v_aroot_layer_init(2)) + + deallocate(ccohort_hydr%v_aroot_layer) + deallocate(ccohort_hydr%v_aroot_layer_init) + end subroutine test_SavePreviousCompartmentVolumes + + ! Summary: Tests SumBetweenDepths depth-weighted integration across non-uniform soil layer interfaces, partial overlaps, and boundary limits. + @Test + subroutine test_SumBetweenDepths(this) + class(TestPlantHydraulicsGeom), intent(inout) :: this + type(ed_site_hydr_type) :: csite_hydr + ! Layer profile input quantities + real(r8) :: array_in(3) + real(r8) :: res_full, res_partial, res_out_of_bounds + real(r8), parameter :: tol = 1.e-6_r8 + + ! 3 rhizosphere soil layers setup with non-uniform thickness + csite_hydr%nlevrhiz = 3 + allocate(csite_hydr%zi_rhiz(3)) + allocate(csite_hydr%dz_rhiz(3)) + ! Soil depth interfaces [m] (0.25m top, 0.75m middle, 1.75m bottom interface) + csite_hydr%zi_rhiz(1) = 0.25_r8 + csite_hydr%zi_rhiz(2) = 0.75_r8 + csite_hydr%zi_rhiz(3) = 1.75_r8 + ! Soil layer thicknesses [m] (0.25m, 0.50m, 1.00m) + csite_hydr%dz_rhiz(1) = 0.25_r8 + csite_hydr%dz_rhiz(2) = 0.50_r8 + csite_hydr%dz_rhiz(3) = 1.00_r8 + + ! Quantity array per soil layer [units/layer] + array_in(1) = 10.0_r8 + array_in(2) = 20.0_r8 + array_in(3) = 30.0_r8 + + ! Case 1: Integration over full layer 2 (0.25m to 0.75m) + ! Rationale: Verifies encasement sum (20.0 / 0.5 = 40.0 depth-weighted rate) + res_full = SumBetweenDepths(csite_hydr, 0.25_r8, 0.75_r8, array_in) + @assertEqual(40.0_r8, res_full, tolerance=tol) + + ! Case 2: Integration over layers 1 and 2 (0.0m to 0.75m) + ! Rationale: Encase layers 1 and 2 (sum 10 + 20 = 30 over 0.75m interval => 40.0) + res_partial = SumBetweenDepths(csite_hydr, 0.0_r8, 0.75_r8, array_in) + @assertEqual(40.0_r8, res_partial, tolerance=tol) + + ! Case 3: Integration depth beyond soil profile (2.0m to 3.0m) + ! Rationale: Deeper than soil column returns 0.0 + res_out_of_bounds = SumBetweenDepths(csite_hydr, 2.0_r8, 3.0_r8, array_in) + @assertEqual(0.0_r8, res_out_of_bounds, tolerance=tol) + + ! Invariants: Non-negative integration density + @assertTrue(res_full >= 0.0_r8) + @assertTrue(res_partial >= 0.0_r8) + + deallocate(csite_hydr%zi_rhiz) + deallocate(csite_hydr%dz_rhiz) + end subroutine test_SumBetweenDepths + +end module FatesPlantHydraulicsMod_geom_test diff --git a/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_init_test.pf b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_init_test.pf new file mode 100644 index 0000000000..881a8951d6 --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_init_test.pf @@ -0,0 +1 @@ +! This file was deleted. diff --git a/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_sizedep_test.pf b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_sizedep_test.pf new file mode 100644 index 0000000000..07c78be9d5 --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_sizedep_test.pf @@ -0,0 +1,185 @@ +! Summary: Unit tests for FatesPlantHydraulicsMod size-dependent state updates when plant geometry changes. +!========================================================================================= +! Tests for FatesPlantHydraulicsMod size-dependent state updates when plant geometry changes +!========================================================================================= +module FatesPlantHydraulicsMod_sizedep_test + use FatesPlantHydraulicsMod, only : UpdateSizeDepPlantHydProps, UpdateSizeDepPlantHydStates, UpdateSizeDepRhizHydProps, InitHydroGlobals + use funit + + use FatesHydraulicsMemMod, only : ed_cohort_hydr_type, ed_site_hydr_type, n_hypool_ag, n_hypool_leaf, & + aroot_p_media, leaf_p_media, stem_p_media, troot_p_media, n_plant_media + use EDPftvarcon, only : EDPftvarcon_inst + use FatesConstantsMod, only : r8 => fates_r8, denh2o => dens_fresh_liquid_water + use EDTypesMod, only : ed_site_type + use FatesCohortMod, only : fates_cohort_type + use FatesPatchMod, only : fates_patch_type + use FatesInterfaceTypesMod, only : bc_in_type + + implicit none + + @TestCase + type, extends(TestCase) :: TestPlantHydraulicsSizedep + contains + end type TestPlantHydraulicsSizedep + +contains + + + ! Summary: Tests size-dependent hydraulic property recalculation (node lengths, volumes, conductances) when plant geometry changes. + @Test + subroutine test_UpdateSizeDepPlantHydProps(this) + use FatesUnitTestPRTParametersMod, only : init_mock_prt_params + use EDPftvarconMockMod, only : init_mock_edpftvarcon + use EDTypesMockMod, only : init_mock_site_patch_cohort + + class(TestPlantHydraulicsSizedep), intent(inout) :: this + type(ed_site_type), pointer :: csite => null() + type(fates_patch_type), pointer :: cpatch => null() + type(fates_cohort_type), pointer :: ccohort => null() + integer :: nlevrhiz = 1 + real(r8) :: expected_kmax + real(r8), parameter :: tol = 1.e-6_r8 + + ! Initialize mock global structures and hydraulics parameters + call init_mock_prt_params() + call init_mock_edpftvarcon() + call InitHydroGlobals() + + call init_mock_site_patch_cohort(csite, cpatch, ccohort) + ccohort%pft = 1 + + ! Configure 1 soil layer and initialize hydraulics site and cohort structures + csite%si_hydr%nlevrhiz = nlevrhiz + call csite%si_hydr%InitHydrSite(1, 1, 1, 1) + call ccohort%co_hydr%AllocateHydrCohortArrays(nlevrhiz) + + ! Soil geometry: thickness dz [m], interface zi [m] + csite%si_hydr%zi_rhiz(1) = 0.5_r8 + csite%si_hydr%dz_rhiz(1) = 0.5_r8 + csite%si_hydr%l_aroot_layer(1) = 1.0_r8 + csite%si_hydr%map_r2s(1,1) = 1 + csite%si_hydr%map_r2s(1,2) = 1 + + ! Execute size-dependent hydraulic property update routine + call UpdateSizeDepPlantHydProps(csite, ccohort) + + ! Expected analytical kmax [kg/m2/MPa/s] derived from sapwood cross-sectional area and organ geometry + expected_kmax = 0.006733541266373205_r8 + + @assertEqual(expected_kmax, ccohort%co_hydr%kmax_stem_lower(1), tolerance=tol) + + ! Clean up allocated cohort structures + call ccohort%co_hydr%DeallocateHydrCohortArrays() + end subroutine test_UpdateSizeDepPlantHydProps + + ! Summary: Tests relative water content (theta) scaling to conserve total water mass when compartment volume changes. + @Test + subroutine test_UpdateSizeDepPlantHydStates(this) + class(TestPlantHydraulicsSizedep), intent(inout) :: this + type(ed_site_type), target :: currentSite + type(fates_cohort_type), target :: ccohort + integer :: nlevrhiz = 2 + integer :: ft = 1 + + ! Setup node saturation and residual moisture thresholds for mock testing + if (.not. allocated(EDPftvarcon_inst%hydr_thetas_node)) then + allocate(EDPftvarcon_inst%hydr_thetas_node(1, n_plant_media)) + allocate(EDPftvarcon_inst%hydr_resid_node(1, n_plant_media)) + end if + ! Saturation moisture threshold [m3/m3] + EDPftvarcon_inst%hydr_thetas_node(1, :) = 0.5_r8 + ! Residual moisture threshold [m3/m3] + EDPftvarcon_inst%hydr_resid_node(1, :) = 0.05_r8 + + allocate(currentSite%si_hydr) + currentSite%si_hydr%nlevrhiz = nlevrhiz + currentSite%si_hydr%h2oveg_growturn_err = 0.0_r8 + + ccohort%pft = ft + ! Cohort density n [indiv/m2] + ccohort%n = 1.0_r8 + allocate(ccohort%co_hydr) + call ccohort%co_hydr%AllocateHydrCohortArrays(nlevrhiz) + + ! Test inputs: initial volume v_init [m3], current volume v_curr [m3], initial theta [m3/m3] + ccohort%co_hydr%v_ag = 2.0_r8 + ccohort%co_hydr%v_ag_init = 1.0_r8 + ccohort%co_hydr%th_ag = 0.3_r8 + + ccohort%co_hydr%v_troot = 2.0_r8 + ccohort%co_hydr%v_troot_init = 1.0_r8 + ccohort%co_hydr%th_troot = 0.3_r8 + + ccohort%co_hydr%v_aroot_layer(:) = 2.0_r8 + ccohort%co_hydr%v_aroot_layer_init(:) = 1.0_r8 + ccohort%co_hydr%th_aroot(:) = 0.3_r8 + + ! Execute state update routine + call UpdateSizeDepPlantHydStates(currentSite, ccohort) + + ! Analytical solution derivation: + ! Water mass conservation: theta_new * v_curr = theta_old * v_init + ! theta_new = 0.3 * (1.0 / 2.0) = 0.15 m3/m3 (uncapped since 0.05 <= 0.15 <= 0.5) + @assertEqual(0.15_r8, ccohort%co_hydr%th_ag(1), tolerance=1.e-8_r8) + @assertEqual(0.15_r8, ccohort%co_hydr%th_troot, tolerance=1.e-8_r8) + @assertEqual(0.15_r8, ccohort%co_hydr%th_aroot(1), tolerance=1.e-8_r8) + + ! Invariant D: Volume scaling water mass conservation (theta_new * v_curr == theta_old * v_init) + @assertEqual(0.3_r8 * 1.0_r8, ccohort%co_hydr%th_ag(1) * ccohort%co_hydr%v_ag(1), tolerance=1.e-8_r8) + @assertEqual(0.3_r8 * 1.0_r8, ccohort%co_hydr%th_troot * ccohort%co_hydr%v_troot, tolerance=1.e-8_r8) + @assertEqual(0.3_r8 * 1.0_r8, ccohort%co_hydr%th_aroot(1) * ccohort%co_hydr%v_aroot_layer(1), tolerance=1.e-8_r8) + + ! Teardown + call ccohort%co_hydr%DeallocateHydrCohortArrays() + deallocate(ccohort%co_hydr) + deallocate(currentSite%si_hydr) + end subroutine test_UpdateSizeDepPlantHydStates + + ! Summary: Tests rhizosphere shell volume and root surface area update routines when root lengths change. + @Test + subroutine test_UpdateSizeDepRhizHydProps(this) + class(TestPlantHydraulicsSizedep), intent(inout) :: this + type(ed_site_type), target :: currentSite + type(bc_in_type) :: bc_in + integer :: nlevrhiz = 1 + + allocate(currentSite%si_hydr) + currentSite%si_hydr%nlevrhiz = nlevrhiz + call currentSite%si_hydr%InitHydrSite(1, 1, 1, 1) + + allocate(bc_in%hksat_sisl(1)) + allocate(bc_in%dz_sisl(1)) + bc_in%hksat_sisl(1) = 0.1_r8 + bc_in%dz_sisl(1) = 0.5_r8 + + ! Initial values: l_aroot = 10.0 m, r_node_shell = 2.0 m, v_shell = 3.0 m3 + currentSite%si_hydr%l_aroot_layer(:) = 10.0_r8 + currentSite%si_hydr%r_node_shell(:,:) = 2.0_r8 + currentSite%si_hydr%v_shell(:,:) = 3.0_r8 + currentSite%si_hydr%dz_rhiz(:) = 0.5_r8 + currentSite%si_hydr%rs1(:) = 0.01_r8 + currentSite%si_hydr%map_r2s(1,1) = 1 + currentSite%si_hydr%map_r2s(1,2) = 1 + + ! Empty patch topology: youngest_patch = null + nullify(currentSite%youngest_patch) + + ! Execute rhizosphere property update routine + call UpdateSizeDepRhizHydProps(currentSite, bc_in) + + ! SavePreviousRhizVolumes copies current properties to _init variables before recalculation + @assertEqual(10.0_r8, currentSite%si_hydr%l_aroot_layer_init(1), tolerance=1.e-8_r8) + @assertEqual(2.0_r8, currentSite%si_hydr%r_node_shell_init(1,1), tolerance=1.e-8_r8) + @assertEqual(3.0_r8, currentSite%si_hydr%v_shell_init(1,1), tolerance=1.e-8_r8) + + ! With no patches present, l_aroot_layer is re-summed to 0.0 m + @assertEqual(0.0_r8, currentSite%si_hydr%l_aroot_layer(1), tolerance=1.e-8_r8) + + ! Teardown + deallocate(bc_in%hksat_sisl) + deallocate(bc_in%dz_sisl) + deallocate(currentSite%si_hydr) + end subroutine test_UpdateSizeDepRhizHydProps + +end module FatesPlantHydraulicsMod_sizedep_test +