diff --git a/CMakeLists.txt b/CMakeLists.txt index d2d735c390..7aaa2fde0b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -42,16 +42,16 @@ add_subdirectory(${HLM_ROOT}/share/src csm_share) add_subdirectory(${HLM_ROOT}/share/unit_test_stubs/util csm_share_stubs) # Add FATES source directories -add_subdirectory(${HLM_ROOT}/src/fates/main fates_main) -add_subdirectory(${HLM_ROOT}/src/fates/biogeochem fates_biogeochem) -add_subdirectory(${HLM_ROOT}/src/fates/biogeophys fates_biogeophys) -add_subdirectory(${HLM_ROOT}/src/fates/parteh fates_parteh) -add_subdirectory(${HLM_ROOT}/src/fates/fire fates_fire) -add_subdirectory(${HLM_ROOT}/src/fates/radiation fates_radiation) +add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/main fates_main) +add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/biogeochem fates_biogeochem) +add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/biogeophys fates_biogeophys) +add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/parteh fates_parteh) +add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/fire fates_fire) +add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/radiation fates_radiation) # Testing directories -add_subdirectory(${HLM_ROOT}/src/fates/testing/tests/fortran_shr test_share) -add_subdirectory(${HLM_ROOT}/src/fates/testing/tests/functional/fire/shr fire_share) +add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/testing/tests/fortran_shr test_share) +add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/testing/tests/functional/fire/shr fire_share) # Remove shr_mpi_mod from share_sources. # This is needed because we want to use the mock shr_mpi_mod in place of the real one @@ -111,4 +111,4 @@ link_directories(${CMAKE_CURRENT_BINARY_DIR}) link_libraries(esmf) # Add the main test directory -add_subdirectory(${HLM_ROOT}/src/fates/testing) +add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/testing) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 5d15bf1dd5..71127f3c94 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) @@ -244,11 +246,13 @@ module FatesPlantHydraulicsMod ! ! !PUBLIC MEMBER FUNCTIONS: public :: AccumulateMortalityWaterStorage + public :: AccumulateMortalityWater_explicit ! only non-protected for unit tests public :: RecruitWaterStorage public :: hydraulics_drive public :: InitHydrSites public :: HydrSiteColdStart public :: BTranForHLMDiagnosticsFromCohortHydr + public :: BTranForHLMDiagnostics_explicit ! only non-protected for unit tests public :: InitHydrCohort public :: DeallocateHydrCohort public :: UpdateH2OVeg @@ -256,16 +260,38 @@ module FatesPlantHydraulicsMod public :: UpdateSizeDepPlantHydProps public :: UpdateSizeDepPlantHydStates public :: UpdatePlantPsiFTCFromTheta + public :: UpdatePlantPsiFTCFromTheta_explicit ! only non-protected for unit tests + + interface UpdatePlantPsiFTCFromTheta + module procedure UpdatePlantPsiFTCFromTheta_cohort + module procedure UpdatePlantPsiFTCFromTheta_explicit + end interface UpdatePlantPsiFTCFromTheta + public :: InitPlantHydStates public :: UpdateSizeDepRhizHydProps public :: RestartHydrStates public :: SavePreviousCompartmentVolumes - public :: SavePreviousRhizVolumes + public :: UpdatePlantHydrNodes public :: UpdatePlantHydrLenVol + + interface UpdatePlantHydrLenVol + module procedure UpdatePlantHydrLenVol_cohort + module procedure UpdatePlantHydrLenVol_masses + end interface UpdatePlantHydrLenVol + public :: UpdatePlantKmax + public :: UpdatePlantKmax_explicit ! only non-protected for unit tests + + interface UpdatePlantKmax + module procedure UpdatePlantKmax_cohort + module procedure UpdatePlantKmax_explicit + end interface UpdatePlantKmax public :: ConstrainRecruitNumber public :: InitHydroGlobals + public :: OrderLayersForSolve1D, wrf_plant, wkf_plant + public :: shellGeom + public :: SumBetweenDepths ! only non-protected for unit tests ! RGK 12-2021: UpdateSizeDepRhizHydStates was removed ! this code can be found in tags prior to @@ -681,52 +707,63 @@ end subroutine InitPlantHydStates ! ===================================================================================== - subroutine UpdatePlantPsiFTCFromTheta(ccohort,csite_hydr) + subroutine UpdatePlantPsiFTCFromTheta_cohort(ccohort,csite_hydr) - ! This subroutine updates the potential and the fractional - ! of total conductivity based on the relative water - ! content - ! Arguments + ! Unpacking wrapper around UpdatePlantPsiFTCFromTheta_explicit type(fates_cohort_type),intent(inout), target :: ccohort type(ed_site_hydr_type),intent(in), target :: csite_hydr - - ! Locals - integer :: ft ! Plant functional type - integer :: k ! loop index for compartments - integer :: j ! Loop index for soil layers type(ed_cohort_hydr_type), pointer :: ccohort_hydr + ccohort_hydr => ccohort%co_hydr + call UpdatePlantPsiFTCFromTheta_explicit( ccohort%pft, csite_hydr%nlevrhiz, & + ccohort_hydr%th_ag, ccohort_hydr%th_troot, ccohort_hydr%th_aroot, & + ccohort_hydr%psi_ag, ccohort_hydr%ftc_ag, ccohort_hydr%btran, & + ccohort_hydr%psi_troot, ccohort_hydr%ftc_troot, & + ccohort_hydr%psi_aroot, ccohort_hydr%ftc_aroot ) + end subroutine UpdatePlantPsiFTCFromTheta_cohort + subroutine UpdatePlantPsiFTCFromTheta_explicit( ft, nlevrhiz, th_ag, th_troot, th_aroot, & + psi_ag, ftc_ag, btran, psi_troot, ftc_troot, & + psi_aroot, ftc_aroot ) - ccohort_hydr => ccohort%co_hydr - ft = ccohort%pft + ! This subroutine updates the potential and the fractional + ! of total conductivity based on explicit relative water content arrays + integer, intent(in) :: ft, nlevrhiz + real(r8), intent(in) :: th_ag(:), th_troot, th_aroot(:) + real(r8), intent(out) :: psi_ag(:), ftc_ag(:), btran + real(r8), intent(out) :: psi_troot, ftc_troot + real(r8), intent(out) :: psi_aroot(:), ftc_aroot(:) + + ! Locals + integer :: k ! loop index for compartments + integer :: j ! Loop index for soil layers ! Update Psi and FTC in above-ground compartments ! ----------------------------------------------------------------------------------- do k = 1,n_hypool_leaf - ccohort_hydr%psi_ag(k) = wrf_plant(leaf_p_media,ft)%p%psi_from_th(ccohort_hydr%th_ag(k)) - ccohort_hydr%ftc_ag(k) = wkf_plant(leaf_p_media,ft)%p%ftc_from_psi(ccohort_hydr%psi_ag(k)) + psi_ag(k) = wrf_plant(leaf_p_media,ft)%p%psi_from_th(th_ag(k)) + ftc_ag(k) = wkf_plant(leaf_p_media,ft)%p%ftc_from_psi(psi_ag(k)) end do - ccohort_hydr%btran = wkf_plant(stomata_p_media,ft)%p%ftc_from_psi(ccohort_hydr%psi_ag(1)) + btran = wkf_plant(stomata_p_media,ft)%p%ftc_from_psi(psi_ag(1)) do k = n_hypool_leaf+1, n_hypool_ag - ccohort_hydr%psi_ag(k) = wrf_plant(stem_p_media,ft)%p%psi_from_th(ccohort_hydr%th_ag(k)) - ccohort_hydr%ftc_ag(k) = wkf_plant(stem_p_media,ft)%p%ftc_from_psi(ccohort_hydr%psi_ag(k)) + psi_ag(k) = wrf_plant(stem_p_media,ft)%p%psi_from_th(th_ag(k)) + ftc_ag(k) = wkf_plant(stem_p_media,ft)%p%ftc_from_psi(psi_ag(k)) end do ! Update the Psi and FTC for the transporting root compartment - ccohort_hydr%psi_troot = wrf_plant(troot_p_media,ft)%p%psi_from_th(ccohort_hydr%th_troot) - ccohort_hydr%ftc_troot = wkf_plant(troot_p_media,ft)%p%ftc_from_psi(ccohort_hydr%psi_troot) + psi_troot = wrf_plant(troot_p_media,ft)%p%psi_from_th(th_troot) + ftc_troot = wkf_plant(troot_p_media,ft)%p%ftc_from_psi(psi_troot) ! Update the Psi and FTC for the absorbing roots - do j = 1, csite_hydr%nlevrhiz - ccohort_hydr%psi_aroot(j) = wrf_plant(aroot_p_media,ft)%p%psi_from_th(ccohort_hydr%th_aroot(j)) - ccohort_hydr%ftc_aroot(j) = wkf_plant(aroot_p_media,ft)%p%ftc_from_psi(ccohort_hydr%psi_aroot(j)) + do j = 1, nlevrhiz + psi_aroot(j) = wrf_plant(aroot_p_media,ft)%p%psi_from_th(th_aroot(j)) + ftc_aroot(j) = wkf_plant(aroot_p_media,ft)%p%ftc_from_psi(psi_aroot(j)) end do return - end subroutine UpdatePlantPsiFTCFromTheta + end subroutine UpdatePlantPsiFTCFromTheta_explicit ! ===================================================================================== @@ -885,7 +922,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 +955,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 +989,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 +1005,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 +1044,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 +1059,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 +1120,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 +1129,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 +1145,7 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) end do return - end subroutine UpdatePlantHydrLenVol + end subroutine UpdatePlantHydrLenVol_masses ! ===================================================================================== @@ -2184,6 +2239,24 @@ end subroutine BTranForHLMDiagnosticsFromCohortHydr ! ========================================================================== +subroutine BTranForHLMDiagnostics_explicit(n_cohorts, btran_cohort, balive_cohort, n_indiv, result_btran) + integer, intent(in) :: n_cohorts + real(r8), intent(in) :: btran_cohort(n_cohorts) + real(r8), intent(in) :: balive_cohort(n_cohorts) + real(r8), intent(in) :: n_indiv(n_cohorts) + real(r8), intent(out) :: result_btran + + real(r8) :: total_balive + + total_balive = sum(balive_cohort(:) * n_indiv(:)) + if (total_balive > 0.0_r8) then + result_btran = sum(btran_cohort(:) * balive_cohort(:) * n_indiv(:)) / total_balive + else + result_btran = 0.0_r8 + end if + return +end subroutine BTranForHLMDiagnostics_explicit + subroutine FillDrainRhizShells(nsites, sites, bc_in) ! ! Created by Brad Christoffersen, Jan 2016 @@ -2905,237 +2978,129 @@ end subroutine Hydraulics_BC -subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) +subroutine UpdatePlantKmax_cohort(ccohort_hydr,ccohort,csite_hydr) + + ! Unpacking wrapper around UpdatePlantKmax_explicit + type(ed_cohort_hydr_type),intent(inout),target :: ccohort_hydr + type(fates_cohort_type),intent(in),target :: ccohort + type(ed_site_hydr_type),intent(in),target :: csite_hydr + + integer :: pft + real(r8) :: a_sapwood, c_sap_dummy + + pft = ccohort%pft + + ! Get the cross-section of the plant's sapwood area [m2] + call bsap_allom(ccohort%dbh,pft,ccohort%crowndamage, & + ccohort%canopy_trim, ccohort%efstem_coh, a_sapwood,c_sap_dummy) + + call UpdatePlantKmax_explicit( ccohort_hydr, csite_hydr%nlevrhiz, & + EDPftvarcon_inst%hydr_kmax_node(pft,2), & + EDPftvarcon_inst%hydr_p_taper(pft), & + EDPftvarcon_inst%hydr_rfrac_stem(pft), & + EDPftvarcon_inst%hydr_rs2(pft), & + a_sapwood ) + +end subroutine UpdatePlantKmax_cohort + + +subroutine UpdatePlantKmax_explicit( ccohort_hydr, nlevrhiz, hydr_kmax_node, & + hydr_p_taper, hydr_rfrac_stem, hydr_rs2, & + a_sapwood ) ! --------------------------------------------------------------------------------- - ! ! This routine sets the maximum conductance of all compartments in the plant, from - ! leaves, to stem, to transporting root, to the absorbing roots. - ! These properties are dependent only on the materials (conductivity) and the - ! geometry of the compartments. - ! The units of all K_max values are [kg H2O s-1 MPa-1] - ! - ! There are some different ways to represent overall conductance from node-to-node - ! throughout the hydraulic system. Universally, all can make use of a system - ! where we separate the hydraulic compartments of the nodes into the upper (closer - ! to the sky) and lower (away from the sky) portions of the compartment. It is - ! possible that due to things like xylem taper, the two portions may have different - ! conductivity, and therefore differnet conductances. - ! - ! Assumption 0. This routine calculates maximum conductivity for 1 plant. - ! Assumption 1. The compartment volumes, heights and lengths have all been - ! determined, probably called just before this routine. - ! - ! Steudle, E. Water uptake by roots: effects of water deficit. - ! J Exp Bot 51, 1531-1542, doi:DOI 10.1093/jexbot/51.350.1531 (2000). + ! leaves, to stem, to transporting root, to the absorbing roots, based on explicit + ! conductivity parameters and geometry. ! --------------------------------------------------------------------------------- ! Arguments - type(ed_cohort_hydr_type),intent(inout),target :: ccohort_hydr - type(fates_cohort_type),intent(in),target :: ccohort - type(ed_site_hydr_type),intent(in),target :: csite_hydr + integer, intent(in) :: nlevrhiz + real(r8), intent(in) :: hydr_kmax_node + real(r8), intent(in) :: hydr_p_taper + real(r8), intent(in) :: hydr_rfrac_stem + real(r8), intent(in) :: hydr_rs2 + real(r8), intent(in) :: a_sapwood ! Locals integer :: k ! Compartment (node) index integer :: j ! Soil layer index integer :: k_ag ! Compartment index for above-ground indexed array - integer :: pft ! Plant Functional Type index - real(r8) :: c_sap_dummy ! Dummy variable (unused) with sapwood carbon [kg] real(r8) :: z_lower ! distance between lower edge and mean petiole height [m] real(r8) :: z_upper ! distance between upper edge and mean petiole height [m] real(r8) :: z_node ! distance between compartment center and mph [m] real(r8) :: kmax_lower ! Max conductance from compartment edge to mph [kg s-1 Mpa-1] real(r8) :: kmax_node ! Max conductance from compartment edge to mph [kg s-1 Mpa-1] real(r8) :: kmax_upper ! Max conductance from compartment edge to mph [kg s-1 Mpa-1] - real(r8) :: a_sapwood ! Mean cross section area of sapwood [m2] real(r8) :: rmin_ag ! Minimum total resistance of all above ground pathways - ! [kg-1 s MPa] real(r8) :: kmax_bg ! Total maximum conductance of all below-ground pathways - ! from the absorbing roots center nodes to the - ! transporting root center node - real(r8) :: rootfr ! fraction of absorbing root in each soil layer - ! assumes propotion of absorbing root is equal - ! to proportion of total root real(r8) :: kmax_layer ! max conductance between transporting root node - ! and absorbing root node in each layer [kg s-1 MPa-1] real(r8) :: surfarea_aroot_layer ! Surface area of absorbing roots in each - ! soil layer [m2] real(r8) :: sum_l_aroot ! sum of plant's total root length real(r8),parameter :: min_pet_stem_dz = 0.00001_r8 ! Force at least a small difference - ! in the top of stem and petiole - - - pft = ccohort%pft - - ! Get the cross-section of the plant's sapwood area [m2] - call bsap_allom(ccohort%dbh,pft,ccohort%crowndamage, & - ccohort%canopy_trim, ccohort%efstem_coh, a_sapwood,c_sap_dummy) ! Leaf Maximum Hydraulic Conductance - ! The starting hypothesis is that there is no resistance inside the - ! leaf, between the petiole and the center of storage. To override - ! this, make provisions by changing the kmax to a not-absurdly high - ! value. It is assumed that the conductance in this default case, - ! is regulated completely by the stem conductance from the stem's - ! center of storage, to the petiole. - ccohort_hydr%kmax_petiole_to_leaf = 1.e8_r8 - ! Stem Maximum Hydraulic Conductance - do k=1, n_hypool_stem - - ! index for "above-ground" arrays, that contain stem and leaf - ! in one vector k_ag = k+n_hypool_leaf - ! Depth from the petiole to the lower, node and upper compartment edges - z_lower = ccohort_hydr%z_node_ag(n_hypool_leaf) - ccohort_hydr%z_lower_ag(k_ag) z_node = ccohort_hydr%z_node_ag(n_hypool_leaf) - ccohort_hydr%z_node_ag(k_ag) z_upper = max( min_pet_stem_dz,ccohort_hydr%z_node_ag(n_hypool_leaf) - & ccohort_hydr%z_upper_ag(k_ag)) + kmax_upper = hydr_kmax_node * xylemtaper(hydr_p_taper, z_upper) * a_sapwood / z_upper + kmax_node = hydr_kmax_node * xylemtaper(hydr_p_taper, z_node) * a_sapwood / z_node + kmax_lower = hydr_kmax_node * xylemtaper(hydr_p_taper, z_lower) * a_sapwood / z_lower - ! Then we calculate the maximum conductance from each the lower, node and upper - ! edges of the compartment to the petiole. The xylem taper factor requires - ! that the kmax it is scaling is from the point of interest to the mean height - ! of the petioles. Then we can back out the conductance over just the path - ! of the upper and lower compartments, but subtracting them as resistors in - ! series. - - ! max conductance from upper edge to mean petiole height - ! If there is no height difference between the upper compartment edge and - ! the petiole, at least give it some nominal amount to void FPE's - kmax_upper = EDPftvarcon_inst%hydr_kmax_node(pft,2) * & - xylemtaper(EDPftvarcon_inst%hydr_p_taper(pft), z_upper) * & - a_sapwood / z_upper - - ! max conductance from node to mean petiole height - kmax_node = EDPftvarcon_inst%hydr_kmax_node(pft,2) * & - xylemtaper(EDPftvarcon_inst%hydr_p_taper(pft), z_node) * & - a_sapwood / z_node - - ! max conductance from lower edge to mean petiole height - kmax_lower = EDPftvarcon_inst%hydr_kmax_node(pft,2) * & - xylemtaper(EDPftvarcon_inst%hydr_p_taper(pft), z_lower) * & - a_sapwood / z_lower - - ! Max conductance over the path of the upper side of the compartment ccohort_hydr%kmax_stem_upper(k) = (1._r8/kmax_node - 1._r8/kmax_upper)**(-1._r8) - - ! Max conductance over the path on the loewr side of the compartment ccohort_hydr%kmax_stem_lower(k) = (1._r8/kmax_lower - 1._r8/kmax_node)**(-1._r8) if(debug) then - ! The following clauses should never be true: - if( (z_lower < z_node) .or. & - (z_node < z_upper) ) then + if( (z_lower < z_node) .or. (z_node < z_upper) ) then write(fates_log(),*) 'Problem calculating stem Kmax' write(fates_log(),*) z_lower, z_node, z_upper write(fates_log(),*) kmax_lower*z_lower, kmax_node*z_node, kmax_upper*z_upper call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if - enddo ! Maximum conductance of the upper compartment in the transporting root - ! that connects to the lowest stem (btw: z_lower_ag(n_hypool_ag) == 0) - z_upper = ccohort_hydr%z_lower_ag(n_hypool_leaf) z_node = ccohort_hydr%z_lower_ag(n_hypool_leaf)-ccohort_hydr%z_node_troot - kmax_node = EDPftvarcon_inst%hydr_kmax_node(pft,2) * & - xylemtaper(EDPftvarcon_inst%hydr_p_taper(pft), z_node) * & - a_sapwood / z_node - - kmax_upper = EDPftvarcon_inst%hydr_kmax_node(pft,2) * & - xylemtaper(EDPftvarcon_inst%hydr_p_taper(pft), z_upper) * & - a_sapwood / z_upper + kmax_node = hydr_kmax_node * xylemtaper(hydr_p_taper, z_node) * a_sapwood / z_node + kmax_upper = hydr_kmax_node * xylemtaper(hydr_p_taper, z_upper) * a_sapwood / z_upper ccohort_hydr%kmax_troot_upper = (1._r8/kmax_node - 1._r8/kmax_upper)**(-1._r8) - ! The maximum conductance between the center node of the transporting root - ! compartment, and the center node of the absorbing root compartment, is calculated - ! as a residual. Specifically, we look at the total resistance the plant has in - ! the stem so far, by adding those resistances in series. - ! Then we use a parameter to specify what fraction of the resistance - ! should be below-ground between the transporting root node and the absorbing roots. - ! After that total is calculated, we then convert to a conductance, and split the - ! conductance in parallel between root layers, based on the root fraction. - ! Note* The inverse of max conductance (KMax) is minimum resistance: - - rmin_ag = 1._r8/ccohort_hydr%kmax_petiole_to_leaf + & sum(1._r8/ccohort_hydr%kmax_stem_upper(1:n_hypool_stem)) + & sum(1._r8/ccohort_hydr%kmax_stem_lower(1:n_hypool_stem)) + & 1._r8/ccohort_hydr%kmax_troot_upper - ! Calculate the residual resistance below ground, as a resistor - ! in series with the existing above ground - ! Invert to find below-ground kmax - ! (rmin_ag+rmin_bg)*fr = rmin_ag - ! rmin_ag + rmin_bg = rmin_ag/fr - ! rmin_bg = (1/fr-1) * rmin_ag - ! - ! if kmax_bg = 1/rmin_bg : - ! - ! kmax_bg = 1/((1/fr-1) * rmin_ag) - - kmax_bg = 1._r8/(rmin_ag*(1._r8/EDPftvarcon_inst%hydr_rfrac_stem(pft) - 1._r8)) - + kmax_bg = 1._r8/(rmin_ag*(1._r8/hydr_rfrac_stem - 1._r8)) - ! The max conductance of each layer is in parallel, therefore - ! the kmax terms of each layer, should sum to kmax_bg sum_l_aroot = sum(ccohort_hydr%l_aroot_layer(:)) - do j=1,csite_hydr%nlevrhiz - + do j=1,nlevrhiz kmax_layer = kmax_bg*ccohort_hydr%l_aroot_layer(j)/sum_l_aroot - - ! Two transport pathways, in two compartments exist in each layer. - ! These pathways are connected in serial. - ! For simplicity, we simply split the resistance between the two. - ! Mathematically, this results in simply doubling the conductance - ! and applying to both paths. Here are the two paths: - ! 1) is the path between the transporting root's center node, to - ! the boundary of the transporting root with the boundary of - ! the absorbing root (kmax_troot_lower) - ! 2) is the path between the boundary of the absorbing root and - ! transporting root, with the absorbing root's center node - ! (kmax_aroot_upper) - ccohort_hydr%kmax_troot_lower(j) = 3.0_r8 * kmax_layer ccohort_hydr%kmax_aroot_upper(j) = 3.0_r8 * kmax_layer ccohort_hydr%kmax_aroot_lower(j) = 3.0_r8 * kmax_layer - end do - ! Finally, we calculate maximum radial conductance from the root - ! surface to its center node. This transport is not a xylem transport - ! like the calculations prior to this. This transport is through the - ! exodermis, cortex, casparian strip and endodermis. The actual conductance - ! will possibly depend on the potential gradient (whether out-of the root, - ! or in-to the root). So we calculate the kmax's for both cases, - ! and save them for the final conductance calculation. - - do j=1,csite_hydr%nlevrhiz - - ! Surface area of the absorbing roots for a single plant in this layer [m2] - surfarea_aroot_layer = 2._r8 * pi_const * & - EDPftvarcon_inst%hydr_rs2(ccohort%pft) * ccohort_hydr%l_aroot_layer(j) - - ! Convert from surface conductivity [kg H2O m-2 s-1 MPa-1] to [kg H2O s-1 MPa-1] - ccohort_hydr%kmax_aroot_radial_in(j) = hydr_kmax_rsurf1 * surfarea_aroot_layer - + do j=1,nlevrhiz + surfarea_aroot_layer = 2._r8 * pi_const * hydr_rs2 * ccohort_hydr%l_aroot_layer(j) + ccohort_hydr%kmax_aroot_radial_in(j) = hydr_kmax_rsurf1 * surfarea_aroot_layer ccohort_hydr%kmax_aroot_radial_out(j) = hydr_kmax_rsurf2 * surfarea_aroot_layer - end do return -end subroutine UpdatePlantKmax +end subroutine UpdatePlantKmax_explicit ! =================================================================================== @@ -4318,14 +4283,11 @@ subroutine AccumulateMortalityWaterStorage(csite,ccohort,delta_n) ccohort_hydr => ccohort%co_hydr csite_hydr => csite%si_hydr - delta_w = (sum(ccohort_hydr%th_ag(:)*ccohort_hydr%v_ag(:)) + & - ccohort_hydr%th_troot*ccohort_hydr%v_troot + & - sum(ccohort_hydr%th_aroot(:)*ccohort_hydr%v_aroot_layer(:)))* & - denh2o*delta_n*AREA_INV + call AccumulateMortalityWater_explicit(delta_n, ccohort_hydr%th_ag, ccohort_hydr%v_ag, & + ccohort_hydr%th_troot, ccohort_hydr%v_troot, & + ccohort_hydr%th_aroot, ccohort_hydr%v_aroot_layer, delta_w) csite_hydr%h2oveg_dead = csite_hydr%h2oveg_dead + delta_w - - csite_hydr%h2oveg = csite_hydr%h2oveg - delta_w return @@ -4333,6 +4295,22 @@ end subroutine AccumulateMortalityWaterStorage !-------------------------------------------------------------------------------! +subroutine AccumulateMortalityWater_explicit(delta_n, th_ag, v_ag, th_troot, v_troot, & + th_aroot, v_aroot_layer, delta_w) + real(r8), intent(in) :: delta_n + real(r8), intent(in) :: th_ag(:), v_ag(:) + real(r8), intent(in) :: th_troot, v_troot + real(r8), intent(in) :: th_aroot(:), v_aroot_layer(:) + real(r8), intent(out) :: delta_w + + delta_w = (sum(th_ag(:) * v_ag(:)) + th_troot * v_troot + & + sum(th_aroot(:) * v_aroot_layer(:))) * & + denh2o * delta_n * AREA_INV + return +end subroutine AccumulateMortalityWater_explicit + +!-------------------------------------------------------------------------------! + subroutine CalculateTotalAvailW(ccohort,csite_hydr,bc_in,dtime,totalAvailW) ! --------------------------------------------------------------------------- @@ -4856,6 +4834,30 @@ subroutine MatSolve2D(csite_hydr,cohort,cohort_hydr, & sapflow,rootuptake,wb_err_plant , dwat_plant, & dth_layershell_site) + type(ed_site_hydr_type), intent(inout),target :: csite_hydr + type(ed_cohort_hydr_type), target :: cohort_hydr + type(fates_cohort_type) , intent(inout), target :: cohort + real(r8),intent(in) :: tmx + real(r8),intent(in) :: qtop + real(r8),intent(out) :: sapflow + real(r8),intent(out) :: rootuptake(:) + real(r8),intent(out) :: wb_err_plant + real(r8),intent(out) :: dwat_plant + real(r8),intent(inout) :: dth_layershell_site(:,:) + + call MatSolve2D_explicit(csite_hydr, cohort_hydr, cohort%pft, cohort%n, & + tmx, qtop, sapflow, rootuptake, wb_err_plant, & + dwat_plant, dth_layershell_site) + +end subroutine MatSolve2D + +! ===================================================================================== + +subroutine MatSolve2D_explicit(csite_hydr,cohort_hydr,ft,cohort_n, & + tmx,qtop, & + sapflow,rootuptake,wb_err_plant , dwat_plant, & + dth_layershell_site) + ! --------------------------------------------------------------------------------- ! This solution to the plant water flux equations casts all the fluxes through a @@ -4891,7 +4893,8 @@ subroutine MatSolve2D(csite_hydr,cohort,cohort_hydr, & ! ----------------------------------------------------------------------------------- type(ed_site_hydr_type), intent(inout),target :: csite_hydr ! ED csite_hydr structure type(ed_cohort_hydr_type), target :: cohort_hydr - type(fates_cohort_type) , intent(inout), target :: cohort + integer, intent(in) :: ft + real(r8), intent(in) :: cohort_n real(r8),intent(in) :: tmx ! time interval to integrate over [s] real(r8),intent(in) :: qtop real(r8),intent(out) :: sapflow ! time integrated mass flux between transp-root and stem [kg] @@ -5042,8 +5045,7 @@ subroutine MatSolve2D(csite_hydr,cohort,cohort_hydr, & dth_node => csite_hydr%dth_node, & node_layer => csite_hydr%node_layer, & h_node => csite_hydr%h_node, & - dftc_dpsi_node => csite_hydr%dftc_dpsi_node, & - ft => cohort%pft) + dftc_dpsi_node => csite_hydr%dftc_dpsi_node) ! This NaN's the scratch arrays @@ -5157,7 +5159,7 @@ subroutine MatSolve2D(csite_hydr,cohort,cohort_hydr, & write(fates_log(),*) 'could not converge on a solution.' write(fates_log(),*) 'Perhaps try increasing iteration cap,' write(fates_log(),*) 'and decreasing relaxation factors.' - write(fates_log(),*) 'pft: ',ft,' dbh: ',cohort%dbh + write(fates_log(),*) 'pft: ',ft call endrun(msg=errMsg(sourcefile, __LINE__)) endif @@ -5549,7 +5551,7 @@ subroutine MatSolve2D(csite_hydr,cohort,cohort_hydr, & ishell = k-1 dth_layershell_site(j,ishell) = dth_layershell_site(j,ishell) + & dth_node(inode) * cohort_hydr%l_aroot_layer(j) * & - cohort%n / csite_hydr%l_aroot_layer(j) + cohort_n / csite_hydr%l_aroot_layer(j) endif enddo @@ -5611,8 +5613,12 @@ function SumBetweenDepths(csite_hydr,depth_t,depth_b,array_in) result(depth_sum) ! 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..313cf964ff --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_cohort_test.pf @@ -0,0 +1,539 @@ +! 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, MatSolve2D_explicit + 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 explicitly passed arrays for relative water content (theta) to matric potential (psi) and FTC mapping over a full gradient. + @Test + subroutine test_UpdatePlantPsiFTCFromTheta_sweep(this) + use FatesUnitTestPRTParametersMod, only : init_mock_prt_params + use EDPftvarconMockMod, only : init_mock_edpftvarcon + use FatesHydraulicsMemMod, only : n_plant_media, leaf_p_media + use EDPftvarcon, only : EDPftvarcon_inst + + class(TestPlantHydraulicsCohort), intent(inout) :: this + integer :: pft, nlevrhiz, i + real(r8) :: th_ag(n_hypool_ag), th_troot, th_aroot(1) + 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 + real(r8) :: alpha, m_vg, n_vg, th_sat, th_res, p50, avuln, se, th_test + real(r8), parameter :: tol = 1.e-6_r8 + + call init_mock_prt_params() + call init_mock_edpftvarcon() + call InitHydroGlobals() + + pft = 1 + nlevrhiz = 1 + + ! Fetch the mathematical parameters for the analytical oracles directly from the initialized mock state + alpha = EDPftvarcon_inst%hydr_vg_alpha_node(pft, leaf_p_media) + m_vg = EDPftvarcon_inst%hydr_vg_m_node(pft, leaf_p_media) + n_vg = EDPftvarcon_inst%hydr_vg_n_node(pft, leaf_p_media) + th_sat = EDPftvarcon_inst%hydr_thetas_node(pft, leaf_p_media) + th_res = EDPftvarcon_inst%hydr_resid_node(pft, leaf_p_media) + p50 = EDPftvarcon_inst%hydr_p50_node(pft, leaf_p_media) + avuln = EDPftvarcon_inst%hydr_avuln_node(pft, leaf_p_media) + + ! Parameter Sweep: Evaluate from 98% saturation down to 10% (near residual) + do i = 1, 10 + ! Generate a test theta + th_test = th_res + (th_sat - th_res) * (1.0_r8 - (real(i, r8) / 11.0_r8)) + + th_ag(:) = th_test + th_troot = th_test + th_aroot(1) = th_test + + ! 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 Oracle 1: Van Genuchten Psi + se = (th_test - th_res) / (th_sat - th_res) + expected_psi = -(1.0_r8 / alpha) * ((se**(-1.0_r8 / m_vg) - 1.0_r8)**(1.0_r8 / n_vg)) + + ! Analytical Oracle 2: Weibull FTC + expected_ftc = 1.0_r8 / (1.0_r8 + (expected_psi / p50)**avuln) + + @assertEqual(expected_psi, psi_ag(1), tol, message="Failed VG Psi curve") + @assertEqual(expected_ftc, ftc_ag(1), tol, message="Failed Weibull FTC curve") + @assertEqual(expected_psi, psi_troot, tol) + @assertEqual(expected_ftc, ftc_troot, tol) + @assertEqual(expected_psi, psi_aroot(1), tol) + @assertEqual(expected_ftc, ftc_aroot(1), tol) + + ! Invariants + @assertTrue(all(psi_ag <= 0.0_r8)) + @assertTrue(all(ftc_ag >= 0.0_r8 .and. ftc_ag <= 1.0_r8)) + end do + + end subroutine test_UpdatePlantPsiFTCFromTheta_sweep + + + ! 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 + @assertEqual(4.1827368176100528_r8, ccohort_hydr%kmax_troot_lower(1), tolerance=1.e-6_r8) + @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 + + ! Summary: Level 1 Smoke Test - Verifies non-NaN outputs and basic solver convergence under standard conditions. + @Test + subroutine test_MatSolve2D_explicit_smoke(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() + integer :: nlevrhiz, nshell + real(r8) :: tmx, qtop, sapflow, wb_err_plant, dwat_plant + real(r8), allocatable :: rootuptake(:) + real(r8), allocatable :: dth_layershell_site(:,:) + + call init_mock_prt_params() + call init_mock_edpftvarcon() + call InitHydroGlobals() + + call init_mock_site_patch_cohort(csite, cpatch, ccohort) + ccohort%pft = 1 + ccohort%n = 1.0_r8 + + nlevrhiz = 1 + nshell = 1 + csite%si_hydr%nlevrhiz = nlevrhiz + call csite%si_hydr%InitHydrSite(1, 1, 1, 1) + call ccohort%co_hydr%AllocateHydrCohortArrays(nlevrhiz) + + allocate(rootuptake(nlevrhiz)) + allocate(dth_layershell_site(nlevrhiz, nshell)) + rootuptake(:) = 0.0_r8 + dth_layershell_site(:,:) = 0.0_r8 + + 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.33_r8, 0.5_r8, 0.05_r8]) + end select + + tmx = 1800.0_r8 + qtop = 0.0_r8 + + call MatSolve2D_explicit(csite%si_hydr, ccohort%co_hydr, ccohort%pft, ccohort%n, & + tmx, qtop, sapflow, rootuptake, wb_err_plant, dwat_plant, & + dth_layershell_site) + + ! Level 1 Invariants: Non-NaN and basic mass conservation + @assertFalse(ieee_is_nan(sapflow)) + @assertFalse(ieee_is_nan(dwat_plant)) + @assertFalse(ieee_is_nan(wb_err_plant)) + @assertEqual(0.0_r8, wb_err_plant, tolerance=1.e-6_r8) + + deallocate(rootuptake) + deallocate(dth_layershell_site) + call ccohort%co_hydr%DeallocateHydrCohortArrays() + end subroutine test_MatSolve2D_explicit_smoke + + ! Summary: Level 2 Boundary & Edge Case Test - Verifies stability under extreme transpiration demand and multi-layer soil topology. + @Test + subroutine test_MatSolve2D_explicit_boundary_cases(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() + integer :: nlevrhiz, nshell, j + real(r8) :: tmx, qtop, sapflow, wb_err_plant, dwat_plant + real(r8), allocatable :: rootuptake(:) + real(r8), allocatable :: dth_layershell_site(:,:) + + call init_mock_prt_params() + call init_mock_edpftvarcon() + call InitHydroGlobals() + + call init_mock_site_patch_cohort(csite, cpatch, ccohort) + ccohort%pft = 1 + ccohort%n = 1.0_r8 + + ! Multi-layer (2 soil layers) and 2 rhizosphere shells setup + nlevrhiz = 2 + nshell = 2 + csite%si_hydr%nlevrhiz = nlevrhiz + call csite%si_hydr%InitHydrSite(1, 1, 1, 1) + call ccohort%co_hydr%AllocateHydrCohortArrays(nlevrhiz) + + allocate(rootuptake(nlevrhiz)) + allocate(dth_layershell_site(nlevrhiz, nshell)) + rootuptake(:) = 0.0_r8 + dth_layershell_site(:,:) = 0.0_r8 + + do j = 1, nlevrhiz + allocate(wrf_type_vg :: csite%si_hydr%wrf_soil(j)%p) + select type (p => csite%si_hydr%wrf_soil(j)%p) + type is (wrf_type_vg) + call p%set_wrf_param([0.001_r8, 1.5_r8, 0.33_r8, 0.5_r8, 0.05_r8]) + end select + end do + + tmx = 1800.0_r8 + ! Pathological extreme transpiration demand [100.0 kg/s] + qtop = 100.0_r8 + + call MatSolve2D_explicit(csite%si_hydr, ccohort%co_hydr, ccohort%pft, ccohort%n, & + tmx, qtop, sapflow, rootuptake, wb_err_plant, dwat_plant, & + dth_layershell_site) + + ! Level 2 Invariants: Physical bounds and solver resilience under extreme desiccation + @assertFalse(ieee_is_nan(sapflow)) + @assertFalse(ieee_is_nan(dwat_plant)) + @assertFalse(ieee_is_nan(wb_err_plant)) + + deallocate(rootuptake) + deallocate(dth_layershell_site) + call ccohort%co_hydr%DeallocateHydrCohortArrays() + end subroutine test_MatSolve2D_explicit_boundary_cases + + ! Summary: Level 3 Parameter Sweep & Oracle Test - Verifies mass conservation across a continuous gradient of transpiration demands. + @Test + subroutine test_MatSolve2D_explicit_mass_conservation_sweep(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() + integer :: nlevrhiz, nshell, i + real(r8) :: tmx, qtop, sapflow, wb_err_plant, dwat_plant + real(r8), allocatable :: rootuptake(:) + real(r8), allocatable :: dth_layershell_site(:,:) + + ! Transpiration demand gradient [kg/s] + real(r8), parameter :: qtop_sweep(5) = [0.0_r8, 1.e-6_r8, 5.e-6_r8, 1.e-5_r8, 5.e-5_r8] + + call init_mock_prt_params() + call init_mock_edpftvarcon() + call InitHydroGlobals() + + call init_mock_site_patch_cohort(csite, cpatch, ccohort) + ccohort%pft = 1 + ccohort%n = 1.0_r8 + + nlevrhiz = 1 + nshell = 1 + csite%si_hydr%nlevrhiz = nlevrhiz + call csite%si_hydr%InitHydrSite(1, 1, 1, 1) + call ccohort%co_hydr%AllocateHydrCohortArrays(nlevrhiz) + + allocate(rootuptake(nlevrhiz)) + allocate(dth_layershell_site(nlevrhiz, nshell)) + + 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.33_r8, 0.5_r8, 0.05_r8]) + end select + + tmx = 1800.0_r8 + + ! Sweep transpiration demand across 5 continuous points + do i = 1, size(qtop_sweep) + qtop = qtop_sweep(i) + rootuptake(:) = 0.0_r8 + dth_layershell_site(:,:) = 0.0_r8 + + call MatSolve2D_explicit(csite%si_hydr, ccohort%co_hydr, ccohort%pft, ccohort%n, & + tmx, qtop, sapflow, rootuptake, wb_err_plant, dwat_plant, & + dth_layershell_site) + + ! Level 3 Mathematical Oracle Assertion: + ! Mass balance error MUST remain strictly zero (< 1e-5 kg) across the full gradient + @assertEqual(0.0_r8, wb_err_plant, tolerance=1.e-5_r8, & + message="Mass conservation violated during transpiration sweep") + end do + + deallocate(rootuptake) + deallocate(dth_layershell_site) + call ccohort%co_hydr%DeallocateHydrCohortArrays() + end subroutine test_MatSolve2D_explicit_mass_conservation_sweep + +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..e261386bf1 --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_geom_test.pf @@ -0,0 +1,327 @@ +! 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) + ! Summary: Tests SumBetweenDepths integration when requested bottom depth exceeds total soil column depth. + @Test + subroutine test_SumBetweenDepths_ExceedingSoilColumn(this) + use FatesPlantHydraulicsMod, only : SumBetweenDepths + class(TestPlantHydraulicsGeom), intent(inout) :: this + type(ed_site_hydr_type) :: csite_hydr + real(r8) :: array_in(2) + real(r8) :: result_sum + real(r8), parameter :: tol = 1.e-6_r8 + + ! Soil column with 2 layers: layer 1 [0.0 - 0.5m, dz=0.5], layer 2 [0.5 - 1.0m, dz=0.5] + csite_hydr%nlevrhiz = 2 + allocate(csite_hydr%zi_rhiz(2)) + allocate(csite_hydr%dz_rhiz(2)) + csite_hydr%zi_rhiz(1) = 0.5_r8 + csite_hydr%zi_rhiz(2) = 1.0_r8 + csite_hydr%dz_rhiz(1) = 0.5_r8 + csite_hydr%dz_rhiz(2) = 0.5_r8 + + ! Quantity array_in: 10.0 in layer 1, 20.0 in layer 2 + ! Total integral across full 1.0m soil column = 10*0.5 + 20*0.5 = 15.0 + ! Average flux per meter over 1.0m soil column = 15.0 / 1.0 = 15.0 + array_in(1) = 10.0_r8 + array_in(2) = 20.0_r8 + + ! Query integration depth depth_t = 0.0m, depth_b = 3.0m (exceeds soil column 1.0m) + ! Proper soil column capping: integral is 30.0 / (min(3.0, 1.0) - 0.0) = 30.0 / 1.0 = 30.0 + ! Mutated code (uncapped): integral is 30.0 / (3.0 - 0.0) = 10.0 (FAILS!) + result_sum = SumBetweenDepths(csite_hydr, 0.0_r8, 3.0_r8, array_in) + + @assertEqual(30.0_r8, result_sum, tolerance=tol) + + deallocate(csite_hydr%zi_rhiz) + deallocate(csite_hydr%dz_rhiz) + end subroutine test_SumBetweenDepths_ExceedingSoilColumn + + ! Summary: Tests SumBetweenDepths depth-weighted integration across non-uniform soil layer interfaces, partial overlaps, and boundary limits. + @Test + subroutine test_SumBetweenDepths(this) + use FatesPlantHydraulicsMod, only : SumBetweenDepths + 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) + end subroutine test_SumBetweenDepths + + + ! Summary: Tests SumBetweenDepths for shallow root depth integration (i_rhiz_b == 0) on non-uniform vertical grid spacing. + @Test + subroutine test_SumBetweenDepths_ShallowPartialLayer(this) + class(TestPlantHydraulicsGeom), intent(inout) :: this + type(ed_site_hydr_type) :: csite_hydr + real(r8) :: array_in(3) + real(r8) :: res_shallow + ! Numerical tolerance for depth-weighted sum + real(r8), parameter :: tol = 1.e-6_r8 + + ! 3 rhizosphere soil layers setup with non-uniform thickness (5cm, 15cm, 80cm) + csite_hydr%nlevrhiz = 3 + allocate(csite_hydr%zi_rhiz(3)) + allocate(csite_hydr%dz_rhiz(3)) + csite_hydr%zi_rhiz(1) = 0.05_r8 + csite_hydr%zi_rhiz(2) = 0.20_r8 + csite_hydr%zi_rhiz(3) = 1.00_r8 + csite_hydr%dz_rhiz(1) = 0.05_r8 + csite_hydr%dz_rhiz(2) = 0.15_r8 + csite_hydr%dz_rhiz(3) = 0.80_r8 + + array_in(1) = 10.0_r8 + array_in(2) = 20.0_r8 + array_in(3) = 30.0_r8 + + ! Integration from surface 0.0m down to shallow depth 0.02m (within layer 1 where i_rhiz_b == 0) + ! Analytical derivation: + ! Fraction of layer 1 = 0.02 / 0.05 = 0.40 + ! Integrated quantity = 0.40 * 10.0 = 4.0 + ! Normalized depth-weighted sum = 4.0 / (0.02 - 0.00) = 200.0 + res_shallow = SumBetweenDepths(csite_hydr, 0.0_r8, 0.02_r8, array_in) + @assertEqual(200.0_r8, res_shallow, tolerance=tol) + + deallocate(csite_hydr%zi_rhiz) + deallocate(csite_hydr%dz_rhiz) + end subroutine test_SumBetweenDepths_ShallowPartialLayer + +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 +