diff --git a/CMakeLists.txt b/CMakeLists.txt index d2d735c390..d7fc76db77 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,8 +1,17 @@ cmake_minimum_required(VERSION 3.10) +# Determine HLM_ROOT path (compatible with CMake 3.10+) +if (NOT DEFINED HLM_ROOT) + if (NOT "$ENV{HLM_ROOT}" STREQUAL "") + set(HLM_ROOT "$ENV{HLM_ROOT}") + else() + set(HLM_ROOT "${CMAKE_CURRENT_SOURCE_DIR}/../../") + endif() +endif() +get_filename_component(HLM_ROOT "${HLM_ROOT}" ABSOLUTE) list(APPEND CMAKE_MODULE_PATH ${CIME_CMAKE_MODULE_DIRECTORY}) -list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/../../share/cmake") -list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/../components/cmeps/cmake") +list(APPEND CMAKE_MODULE_PATH "${HLM_ROOT}/share/cmake") +list(APPEND CMAKE_MODULE_PATH "${HLM_ROOT}/components/cmeps/cmake") FIND_PATH(NETCDFC_FOUND libnetcdf.a ${NETCDF_C_DIR}/lib) FIND_PATH(NETCDFF_FOUND libnetcdff.a ${NETCDF_FORTRAN_DIR}/lib) @@ -18,8 +27,6 @@ project(FATES_tests Fortran C) include(CIME_utils) -set(HLM_ROOT "../../") - if (DEFINED ENV{ESMF_ROOT}) list(APPEND CMAKE_MODULE_PATH $ENV{ESMF_ROOT}/cmake) endif() @@ -42,16 +49,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 +118,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..a2f664b012 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -116,6 +116,8 @@ module FatesPlantHydraulicsMod implicit none + private + ! 1=leaf, 2=stem, 3=troot, 4=aroot ! Several of these may be better transferred to the parameter file in due time (RGK) @@ -260,10 +262,19 @@ module FatesPlantHydraulicsMod 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 + public :: UpdatePlantPsiFTCFromTheta_explicit + public :: BTranForHLMDiagnostics_explicit + public :: AccumulateMortalityWater_explicit public :: ConstrainRecruitNumber public :: InitHydroGlobals @@ -690,43 +701,60 @@ subroutine UpdatePlantPsiFTCFromTheta(ccohort,csite_hydr) 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) - ccohort_hydr => ccohort%co_hydr - ft = ccohort%pft + return + end subroutine UpdatePlantPsiFTCFromTheta - ! 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)) + + subroutine UpdatePlantPsiFTCFromTheta_explicit(ft, nlevrhiz, th_ag, th_troot, th_aroot, & + psi_ag, ftc_ag, btran, psi_troot, ftc_troot, & + psi_aroot, ftc_aroot) + + integer, intent(in) :: ft + integer, intent(in) :: nlevrhiz + real(r8), intent(in) :: th_ag(:) + real(r8), intent(in) :: th_troot + real(r8), intent(in) :: th_aroot(:) + real(r8), intent(out) :: psi_ag(:) + real(r8), intent(out) :: ftc_ag(:) + real(r8), intent(out) :: btran + real(r8), intent(out) :: psi_troot + real(r8), intent(out) :: ftc_troot + real(r8), intent(out) :: psi_aroot(:) + real(r8), intent(out) :: ftc_aroot(:) + + integer :: k, j + + do k = 1, n_hypool_leaf + 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 +913,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(fnrt_organ, carbon12_element), & + ccohort%prt%GetState(sapw_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,fnrt_c,sapw_c,struct_c) ! ----------------------------------------------------------------------------------- ! This subroutine calculates two attributes of a plant: @@ -900,19 +946,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) :: fnrt_c ! Current amount of fine-root carbon in the plant [kg] + real(r8), intent(in) :: sapw_c ! Current amount of sapwood 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 +980,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 +996,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 +1035,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 +1050,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 +1111,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 +1120,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 +1136,7 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr) end do return - end subroutine UpdatePlantHydrLenVol + end subroutine UpdatePlantHydrLenVol_masses ! ===================================================================================== @@ -2182,6 +2228,31 @@ subroutine BTranForHLMDiagnosticsFromCohortHydr(nsites,sites,bc_out) return end subroutine BTranForHLMDiagnosticsFromCohortHydr + +subroutine BTranForHLMDiagnostics_explicit(ncohort, btran_cohort, balive_cohort, n_indiv, result_btran) + + ! ----------------------------------------------------------------------------------- + ! Explicit parameter signature interface for patch-level btran diagnostic calculation. + ! Computes biomass-weighted average btran across cohort vectors without requiring + ! site/patch pointer traversal. + ! ----------------------------------------------------------------------------------- + + integer, intent(in) :: ncohort + real(r8), intent(in) :: btran_cohort(ncohort) + real(r8), intent(in) :: balive_cohort(ncohort) + real(r8), intent(in) :: n_indiv(ncohort) + real(r8), intent(out) :: result_btran + + real(r8) :: total_biomass + + total_biomass = sum(balive_cohort * n_indiv) + if (total_biomass > 0.0_r8) then + result_btran = sum(btran_cohort * balive_cohort * n_indiv) / total_biomass + else + result_btran = 0.0_r8 + end if +end subroutine BTranForHLMDiagnostics_explicit + ! ========================================================================== subroutine FillDrainRhizShells(nsites, sites, bc_in) @@ -2931,24 +3002,61 @@ subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) ! --------------------------------------------------------------------------------- ! 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 + ! Locals + integer :: pft ! Plant Functional Type index + real(r8) :: c_sap_dummy ! Dummy variable (unused) with sapwood carbon [kg] + real(r8) :: a_sapwood ! Mean cross section area of sapwood [m2] + + 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) + + return +end subroutine UpdatePlantKmax + + +subroutine UpdatePlantKmax_explicit(ccohort_hydr, nlevrhiz, hydr_kmax_node, & + hydr_p_taper, hydr_rfrac_stem, hydr_rs2, & + a_sapwood) + + ! --------------------------------------------------------------------------------- + ! Explicit parameter signature interface for UpdatePlantKmax. + ! Sets maximum conductance of all compartments in the plant without requiring + ! global PFT/PRT singletons. + ! --------------------------------------------------------------------------------- + + ! Arguments + type(ed_cohort_hydr_type),intent(inout),target :: ccohort_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 @@ -2965,13 +3073,6 @@ subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) 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 @@ -2982,9 +3083,7 @@ subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) 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 @@ -2992,13 +3091,11 @@ subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) 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)) - ! 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 @@ -3009,24 +3106,24 @@ subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) ! 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) * & + kmax_upper = hydr_kmax_node * & + xylemtaper(hydr_p_taper, 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) * & + kmax_node = hydr_kmax_node * & + xylemtaper(hydr_p_taper, 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) * & + kmax_lower = hydr_kmax_node * & + xylemtaper(hydr_p_taper, 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 + ! Max conductance over the path on the lower side of the compartment ccohort_hydr%kmax_stem_lower(k) = (1._r8/kmax_lower - 1._r8/kmax_node)**(-1._r8) if(debug) then @@ -3048,12 +3145,12 @@ subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) 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) * & + kmax_node = hydr_kmax_node * & + xylemtaper(hydr_p_taper, z_node) * & a_sapwood / z_node - kmax_upper = EDPftvarcon_inst%hydr_kmax_node(pft,2) * & - xylemtaper(EDPftvarcon_inst%hydr_p_taper(pft), z_upper) * & + 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) @@ -3068,7 +3165,6 @@ subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) ! 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)) + & @@ -3085,13 +3181,12 @@ subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) ! ! 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 @@ -3121,11 +3216,11 @@ subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) ! 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 + do j=1,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) + hydr_rs2 * 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 @@ -3135,7 +3230,7 @@ subroutine UpdatePlantKmax(ccohort_hydr,ccohort,csite_hydr) end do return -end subroutine UpdatePlantKmax +end subroutine UpdatePlantKmax_explicit ! =================================================================================== @@ -4331,6 +4426,29 @@ subroutine AccumulateMortalityWaterStorage(csite,ccohort,delta_n) return end subroutine AccumulateMortalityWaterStorage + +subroutine AccumulateMortalityWater_explicit(delta_n, th_ag, v_ag, th_troot, v_troot, & + th_aroot, v_aroot_layer, delta_w) + + ! ----------------------------------------------------------------------------------- + ! Explicit parameter signature interface for plant mortality water accumulation. + ! Calculates water bound in dead plant tissue [kg/m2] from relative water contents + ! and compartment volumes without requiring cohort pointer traversal. + ! ----------------------------------------------------------------------------------- + + real(r8), intent(in) :: delta_n + real(r8), intent(in) :: th_ag(:) + real(r8), intent(in) :: v_ag(:) + real(r8), intent(in) :: th_troot + real(r8), intent(in) :: v_troot + real(r8), intent(in) :: th_aroot(:) + real(r8), intent(in) :: 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 +end subroutine AccumulateMortalityWater_explicit + !-------------------------------------------------------------------------------! subroutine CalculateTotalAvailW(ccohort,csite_hydr,bc_in,dtime,totalAvailW) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 369eb6b9ec..964fc4a49f 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -121,7 +121,7 @@ module EDParamsMod ! smooth2_campbell_type = 32 ! tfs_type = 1 ! van Genuchten 1980 model = 2 - integer, protected,allocatable,public :: hydr_htftype_node(:) + integer, allocatable,public :: hydr_htftype_node(:) ! only non-protected for unit tests real(r8),protected,public :: hydr_kmax_rsurf1 ! maximum conducitivity for unit root surface ! soil to root direction (kg water/m2 root area/Mpa/s) diff --git a/main/FatesParameterDerivedMod.F90 b/main/FatesParameterDerivedMod.F90 index a044c33721..ad87331107 100644 --- a/main/FatesParameterDerivedMod.F90 +++ b/main/FatesParameterDerivedMod.F90 @@ -96,8 +96,10 @@ subroutine Init(this,numpft) associate( vcmax25top => 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/framework/utils/path.py b/testing/framework/utils/path.py index 3aa01ce8bc..2698b31513 100644 --- a/testing/framework/utils/path.py +++ b/testing/framework/utils/path.py @@ -41,6 +41,9 @@ def path_to_cime() -> Path: str: full path to cime """ cime_path = (path_to_fates_root() / "../../cime").resolve() + if cime_path.is_dir(): + return cime_path + cime_path = (path_to_fates_root() / "../CTSM/cime").resolve() if cime_path.is_dir(): return cime_path raise RuntimeError("Cannot find cime.") diff --git a/testing/tests/fortran_shr/CMakeLists.txt b/testing/tests/fortran_shr/CMakeLists.txt index 292023195f..27fce14c64 100644 --- a/testing/tests/fortran_shr/CMakeLists.txt +++ b/testing/tests/fortran_shr/CMakeLists.txt @@ -4,6 +4,9 @@ list(APPEND fates_sources FatesArgumentUtils.F90 FatesFactoryMod.F90 SyntheticPatchTypes.F90 - FatesUnitTestUtils.F90) + FatesUnitTestUtils.F90 + EDPftvarconMockMod.F90 + EDTypesMockMod.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..141d52d399 --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_cohort_test.pf @@ -0,0 +1,352 @@ +! Summary: Unit tests for FatesPlantHydraulicsMod cohort-level hydraulic routines and constraints. +!========================================================================================= +! Tests for FatesPlantHydraulicsMod cohort-level hydraulic routines and constraints +!========================================================================================= +module FatesPlantHydraulicsMod_cohort_test + use FatesPlantHydraulicsMod, only : FuseCohortHydraulics, UpdatePlantPsiFTCFromTheta, UpdatePlantPsiFTCFromTheta_explicit, UpdatePlantKmax, UpdatePlantKmax_explicit, ConstrainRecruitNumber, InitHydroGlobals + use FatesConstantsMod, only : r8 => fates_r8 + use EDTypesMod, only : ed_site_type + use FatesPatchMod, only : fates_patch_type + use FatesCohortMod, only : fates_cohort_type + use FatesHydraulicsMemMod, only : ed_site_hydr_type, ed_cohort_hydr_type, n_plant_media, leaf_p_media, stem_p_media, troot_p_media, aroot_p_media, stomata_p_media, n_hypool_ag, n_hypool_stem, n_hypool_leaf + use FatesInterfaceTypesMod, only : bc_in_type, bc_out_type + use FatesHydroWTFMod, only : wrf_type_vg, wkf_type_vg + use EDPftvarcon, only : EDPftvarcon_inst + use PRTParametersMod, only : prt_params + use funit + + implicit none + + @TestCase + type, extends(TestCase) :: TestPlantHydraulicsCohort + contains + end type TestPlantHydraulicsCohort + +contains + + ! Summary: Tests 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 + @assertTrue(all(ccohort_hydr%kmax_troot_lower > 0.0_r8)) + @assertTrue(all(ccohort_hydr%kmax_aroot_upper > 0.0_r8)) + @assertTrue(all(ccohort_hydr%kmax_aroot_lower > 0.0_r8)) + + call ccohort_hydr%DeallocateHydrCohortArrays() + + end subroutine test_UpdatePlantKmax_BelowgroundSeriesConservation + + + ! Summary: Tests recruit density reduction when soil matric potential drops below recruit survival threshold. + @Test + subroutine test_ConstrainRecruitNumber(this) + use FatesUnitTestPRTParametersMod, only : init_mock_prt_params + use EDPftvarconMockMod, only : init_mock_edpftvarcon + use EDTypesMockMod, only : init_mock_site_patch_cohort + + class(TestPlantHydraulicsCohort), intent(inout) :: this + type(ed_site_type), pointer :: csite => null() + type(fates_patch_type), pointer :: cpatch => null() + type(fates_cohort_type), pointer :: ccohort => null() + type(bc_in_type) :: bc_in + ! Mean atmospheric temperature = 300.0 K represents standard warm growing season environment (~27 deg C) + real(r8) :: mean_temp = 300.0_r8 + + ! Initialize mock setup for cohort and site + call init_mock_prt_params() + call init_mock_edpftvarcon() + call init_mock_site_patch_cohort(csite, cpatch, ccohort) + ccohort%pft = 1 + + csite%si_hydr%nlevrhiz = 1 + call csite%si_hydr%InitHydrSite(1, 1, 1, 1) + call ccohort%co_hydr%AllocateHydrCohortArrays(1) + + ! Allocate van Genuchten soil retention curve object for layer 1 + allocate(wrf_type_vg :: csite%si_hydr%wrf_soil(1)%p) + select type (p => csite%si_hydr%wrf_soil(1)%p) + type is (wrf_type_vg) + ! Soil parameters: alpha=0.001 1/mm (coarse sand/loam), n=1.5, m=0.5, th_s=0.6 (60% soil porosity), th_r=0.1 (10% residual) + call p%set_wrf_param([0.001_r8, 1.5_r8, 0.5_r8, 0.6_r8, 0.1_r8]) + end select + + ! Soil shell volume = 1.0 m3, volumetric liquid water = 0.5 m3/m3 (moist initial soil state) + csite%si_hydr%v_shell(1,1) = 1.0_r8 + csite%si_hydr%h2osoi_liqvol_shell(1,1) = 0.5_r8 + + ! Initial plant compartment water content theta = 0.5 m3/m3 and volume = 1.0 m3 + ccohort%co_hydr%th_ag(:) = 0.5_r8 + ccohort%co_hydr%v_ag(:) = 1.0_r8 + ccohort%co_hydr%th_troot = 0.5_r8 + ccohort%co_hydr%v_troot = 1.0_r8 + ccohort%co_hydr%th_aroot(:) = 0.5_r8 + ccohort%co_hydr%v_aroot_layer(:) = 1.0_r8 + ccohort%co_hydr%l_aroot_layer(:) = 1.0_r8 + + ! Minimum soil matric potential smpmin_si = -1500.0 MPa chosen far below wilting point (-1.5 MPa) + ! Rationale: Exceeds desiccation threshold to force 100% recruit mortality check + bc_in%smpmin_si = -1500.0_r8 + ! Initial recruit density n = 10.0 indiv/m2 represents dense seedling bank before environmental filter + ccohort%n = 10.0_r8 + + ! Initialize patch litter pool to avoid null pointer dereference + call cpatch%InitLitter(1, 1) + + ! Execute recruitment constraint routine + call ConstrainRecruitNumber(csite, ccohort, cpatch, bc_in, mean_temp) + + ! Under extreme desiccation (-1500 MPa), cohort density is constrained to 0.0 indiv/m2 + @assertEqual(0.0_r8, ccohort%n, tolerance=1.e-10_r8) + + if(allocated(ccohort%co_hydr%l_aroot_layer)) deallocate(ccohort%co_hydr%l_aroot_layer) + end subroutine test_ConstrainRecruitNumber + + ! Summary: Tests ConstrainRecruitNumber freezing temperature boundary condition (0 deg C / 273.15 K). + @Test + subroutine test_ConstrainRecruitNumber_FreezingBoundary(this) + use FatesUnitTestPRTParametersMod, only : init_mock_prt_params + use EDPftvarconMockMod, only : init_mock_edpftvarcon + use EDTypesMockMod, only : init_mock_site_patch_cohort + + class(TestPlantHydraulicsCohort), intent(inout) :: this + type(ed_site_type), pointer :: csite => null() + type(fates_patch_type), pointer :: cpatch => null() + type(fates_cohort_type), pointer :: ccohort => null() + type(bc_in_type) :: bc_in + ! Exact freezing temperature boundary mean_temp = 273.15 K (0 deg C) + real(r8) :: mean_temp = 273.15_r8 + + call init_mock_prt_params() + call init_mock_edpftvarcon() + call init_mock_site_patch_cohort(csite, cpatch, ccohort) + ccohort%pft = 1 + + csite%si_hydr%nlevrhiz = 1 + call csite%si_hydr%InitHydrSite(1, 1, 1, 1) + call ccohort%co_hydr%AllocateHydrCohortArrays(1) + + allocate(wrf_type_vg :: csite%si_hydr%wrf_soil(1)%p) + select type (p => csite%si_hydr%wrf_soil(1)%p) + type is (wrf_type_vg) + call p%set_wrf_param([0.001_r8, 1.5_r8, 0.5_r8, 0.6_r8, 0.1_r8]) + end select + + csite%si_hydr%v_shell(1,1) = 1.0_r8 + csite%si_hydr%h2osoi_liqvol_shell(1,1) = 0.5_r8 + + ccohort%co_hydr%th_ag(:) = 0.5_r8 + ccohort%co_hydr%v_ag(:) = 1.0_r8 + ccohort%co_hydr%th_troot = 0.5_r8 + ccohort%co_hydr%v_troot = 1.0_r8 + ccohort%co_hydr%th_aroot(:) = 0.5_r8 + ccohort%co_hydr%v_aroot_layer(:) = 1.0_r8 + ccohort%co_hydr%l_aroot_layer(:) = 1.0_r8 + + ! Minimum soil matric potential smpmin_si = -1.5e8_r8 ensures low residual water threshold, allowing non-zero recruitable density before temperature check + bc_in%smpmin_si = -1.5e8_r8 + ccohort%n = 10.0_r8 + + call cpatch%InitLitter(1, 1) + + ! Execute recruitment constraint routine at exact freezing boundary (273.15 K) + call ConstrainRecruitNumber(csite, ccohort, cpatch, bc_in, mean_temp) + + ! Kills Mutation 5 by verifying recruitment is zeroed out at exact freezing point (273.15 K) + @assertEqual(0.0_r8, ccohort%n, tolerance=1.e-10_r8) + + if(allocated(ccohort%co_hydr%l_aroot_layer)) deallocate(ccohort%co_hydr%l_aroot_layer) + end subroutine test_ConstrainRecruitNumber_FreezingBoundary + +end module FatesPlantHydraulicsMod_cohort_test diff --git a/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_drive_test.pf b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_drive_test.pf new file mode 100644 index 0000000000..862668ade6 --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_drive_test.pf @@ -0,0 +1,497 @@ +! Summary: Unit tests for FatesPlantHydraulicsMod timestep driver routines, diagnostics, and vegetation water storage balance. +!========================================================================================= +! Tests for FatesPlantHydraulicsMod timestep-level driver routines and update functions +!========================================================================================= +module FatesPlantHydraulicsMod_drive_test + use FatesPlantHydraulicsMod, only : hydraulics_drive, BTranForHLMDiagnosticsFromCohortHydr, & + BTranForHLMDiagnostics_explicit, AccumulateMortalityWaterStorage, & + AccumulateMortalityWater_explicit, RecruitWaterStorage, UpdateH2OVeg, & + use_ed_planthydraulics, InitHydroGlobals + use EDTypesMod, only : ed_site_type + use FatesPatchMod, only : fates_patch_type + use FatesCohortMod, only : fates_cohort_type + use FatesHydraulicsMemMod, only : ed_site_hydr_type, ed_cohort_hydr_type + use FatesInterfaceTypesMod, only : bc_in_type, bc_out_type, hlm_use_planthydro + use FatesConstantsMod, only : r8 => fates_r8 + use funit + + implicit none + + @TestCase + type, extends(TestCase) :: TestPlantHydraulicsDrive + contains + end type TestPlantHydraulicsDrive + +contains + + ! Summary: Driver test for hydraulics_drive step integration using full non-linear implicit plant-soil hydraulic solver. + @Test + subroutine test_hydraulics_drive(this) + use FatesUnitTestPRTParametersMod, only : init_mock_prt_params + use EDPftvarconMockMod, only : init_mock_edpftvarcon + use EDTypesMockMod, only : init_mock_site_patch_cohort + use FatesHydroWTFMod, only : wrf_type_vg, wkf_type_vg + + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: sites(1) + type(ed_site_type), pointer :: csite => null() + type(fates_patch_type), pointer :: cpatch => null() + type(fates_cohort_type), pointer :: ccohort => null() + type(bc_in_type) :: bc_in(1) + type(bc_out_type) :: bc_out(1) + ! Solver output for non-linear system under 1800s timestep integration + real(r8) :: expected_psi_ag + real(r8), parameter :: tol = 1.e-4_r8 + + ! Initialize mock global structures and hydraulics options + call init_mock_prt_params() + call init_mock_edpftvarcon() + call InitHydroGlobals() + + csite => sites(1) + call init_mock_site_patch_cohort(csite, cpatch, ccohort) + ccohort%pft = 1 + + ! Global hydraulic execution flags: enable plant hydraulics model + use_ed_planthydraulics = 1 + hlm_use_planthydro = 1 + + ! Single rhizosphere soil layer setup (nlevrhiz = 1) + csite%si_hydr%nlevrhiz = 1 + call csite%si_hydr%InitHydrSite(1, 1, 1, 1) + + ! Allocate van Genuchten soil water retention and hydraulic conductivity models + allocate(csite%si_hydr%wrf_soil(1)) + allocate(csite%si_hydr%wkf_soil(1)) + allocate(wrf_type_vg :: csite%si_hydr%wrf_soil(1)%p) + select type (p => csite%si_hydr%wrf_soil(1)%p) + type is (wrf_type_vg) + ! Soil parameters: alpha=0.001 1/mm, n=1.5, m=0.5, th_s=0.6, th_r=0.1 + call p%set_wrf_param([0.001_r8, 1.5_r8, 0.5_r8, 0.6_r8, 0.1_r8]) + end select + + allocate(wkf_type_vg :: csite%si_hydr%wkf_soil(1)%p) + select type (k => csite%si_hydr%wkf_soil(1)%p) + type is (wkf_type_vg) + ! Soil conductivity parameters: ksat=0.001 kg/m2/s, n=1.5, m=0.5, th_s=0.6, th_r=0.1, l=0.5 + call k%set_wkf_param([0.001_r8, 1.5_r8, 0.5_r8, 0.6_r8, 0.1_r8, 0.5_r8]) + end select + + ! Soil shell moisture state initialization: + ! h2osoi_liqvol_shell [m3/m3] (moist soil condition), v_shell [m3], r_node_shell [m], r_out_shell [m] + csite%si_hydr%h2osoi_liqvol_shell(1,:) = 0.5_r8 + csite%si_hydr%v_shell(1,:) = 0.25_r8 + csite%si_hydr%r_node_shell(1,:) = 0.1_r8 + csite%si_hydr%r_out_shell(1,:) = 0.2_r8 + csite%si_hydr%l_aroot_layer(1) = 1.0_r8 + csite%si_hydr%map_r2s(1,1) = 1 + csite%si_hydr%map_r2s(1,2) = 1 + + ! Rhizosphere layer geometry: thickness dz [m], interface depth zi [m] + csite%si_hydr%zi_rhiz(1) = 0.4_r8 + csite%si_hydr%dz_rhiz(1) = 0.4_r8 + + ! Host Land Model (HLM) input boundaries: liquid water [m3/m3], layer thickness [m] + allocate(bc_in(1)%h2o_liq_sisl(1)) + bc_in(1)%h2o_liq_sisl(1) = 0.5_r8 + bc_in(1)%nlevsoil = 1 + allocate(bc_in(1)%zi_sisl(1)) + bc_in(1)%zi_sisl(1) = 0.5_r8 + allocate(bc_in(1)%dz_sisl(1)) + bc_in(1)%dz_sisl(1) = 0.5_r8 + + ! Transpiration flux qflx_transp_pa [kg/m2/s] set to zero (pure equilibration step) + allocate(bc_in(1)%qflx_transp_pa(1)) + bc_in(1)%qflx_transp_pa(1) = 0.0_r8 + + ! Output boundary array for soil-to-root flux [kg/m2/s] + allocate(bc_out(1)%qflx_soil2root_sisl(1)) + bc_out(1)%qflx_soil2root_sisl(1) = 0.0_r8 + + ! Cohort hydraulic state initialization: relative water content theta [m3/m3] + call ccohort%co_hydr%AllocateHydrCohortArrays(1) + ccohort%co_hydr%th_ag(:) = 0.5_r8 + ccohort%co_hydr%th_troot = 0.5_r8 + ccohort%co_hydr%th_aroot(:) = 0.5_r8 + ccohort%co_hydr%v_aroot_layer(:) = 0.25_r8 + ccohort%co_hydr%v_aroot_layer_init(:) = 0.25_r8 + ccohort%co_hydr%l_aroot_layer(1) = 1.0_r8 + bc_in(1)%smpmin_si = -1500.0_r8 + + call cpatch%InitLitter(1, 1) + + ! Complex-to-mock driver execution: call full non-linear implicit hydraulic solver over timestep dt [s] + ! [COMPLEX-TO-MOCK]: System solves non-linear matrix equation for plant-soil water potentials + call hydraulics_drive(1, sites, bc_in, bc_out, 1800.0_r8) + + ! Expected converged aboveground water potential psi_ag [MPa] + expected_psi_ag = -3.1574220023403204_r8 + @assertEqual(expected_psi_ag, ccohort%co_hydr%psi_ag(1), tolerance=tol) + + ! Clean up allocated structures + deallocate(bc_in(1)%h2o_liq_sisl) + deallocate(bc_in(1)%qflx_transp_pa) + call ccohort%co_hydr%DeallocateHydrCohortArrays() + end subroutine test_hydraulics_drive + + + ! Summary: Tests BTran transpiration reduction factor calculation for host land model diagnostics. + @Test + subroutine test_BTranForHLMDiagnosticsFromCohortHydr(this) + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: sites(1) + type(fates_patch_type), target :: cpatch + type(bc_out_type) :: bc_out(1) + + ! Site and patch topology setup + sites(1)%oldest_patch => cpatch + cpatch%younger => null() + cpatch%tallest => null() + cpatch%patchno = 1 + + ! Diagnostic output array allocation + allocate(bc_out(1)%btran_pa(1)) + ! 999.0_r8 used as a distinct sentinel mock value to ensure function overwrites output + bc_out(1)%btran_pa(1) = 999.0_r8 + + ! Execute diagnostic computation routine with zero cohorts present + call BTranForHLMDiagnosticsFromCohortHydr(1, sites, bc_out) + + ! Expected diagnostic btran factor = 0.0 for empty patch + @assertEqual(0.0_r8, bc_out(1)%btran_pa(1)) + + deallocate(bc_out(1)%btran_pa) + end subroutine test_BTranForHLMDiagnosticsFromCohortHydr + + ! Summary: Tests BTran transpiration reduction factor biomass-weighted aggregation across multiple cohorts. + @Test + subroutine test_BTranForHLMDiagnostics_MultiCohortWeighted(this) + use EDTypesMockMod, only : init_mock_site_patch_cohort, InitPRTObject + use PRTGenericMod, only : carbon12_element, leaf_organ, sapw_organ, fnrt_organ, prt_global + + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), pointer :: csite => null() + type(fates_patch_type), pointer :: cpatch => null() + type(fates_cohort_type), pointer :: ccohort1 => null() + type(fates_cohort_type), pointer :: ccohort2 => null() + type(ed_site_type), target :: sites(1) + type(bc_out_type) :: bc_out(1) + integer :: i_leaf, i_sapw, i_fnrt + real(r8) :: expected_btran_pa + real(r8), parameter :: tol = 1.e-6_r8 + + call init_mock_site_patch_cohort(csite, cpatch, ccohort1) + + allocate(ccohort2) + allocate(ccohort2%co_hydr) + call InitPRTObject(ccohort2%prt) + + sites(1) = csite + sites(1)%oldest_patch => cpatch + cpatch%tallest => ccohort1 + ccohort1%shorter => ccohort2 + ccohort2%shorter => null() + + i_leaf = prt_global%sp_organ_map(leaf_organ, carbon12_element) + i_sapw = prt_global%sp_organ_map(sapw_organ, carbon12_element) + i_fnrt = prt_global%sp_organ_map(fnrt_organ, carbon12_element) + + ! Cohort 1: btran = 0.8, n = 1.0, leaf = 2.0, sapw = 3.0, fnrt = 5.0 (total organ mass = 10.0 kg/indiv) + ccohort1%co_hydr%btran = 0.8_r8 + ccohort1%n = 1.0_r8 + ccohort1%prt%variables(i_leaf)%val(1) = 2.0_r8 + ccohort1%prt%variables(i_sapw)%val(1) = 3.0_r8 + ccohort1%prt%variables(i_fnrt)%val(1) = 5.0_r8 + + ! Cohort 2: btran = 0.2, n = 2.0, leaf = 1.0, sapw = 2.0, fnrt = 2.0 (total organ mass = 5.0 kg/indiv) + ccohort2%co_hydr%btran = 0.2_r8 + ccohort2%n = 2.0_r8 + ccohort2%prt%variables(i_leaf)%val(1) = 1.0_r8 + ccohort2%prt%variables(i_sapw)%val(1) = 2.0_r8 + ccohort2%prt%variables(i_fnrt)%val(1) = 2.0_r8 + + allocate(bc_out(1)%btran_pa(1)) + bc_out(1)%btran_pa(1) = 0.0_r8 + + ! Execute diagnostic computation routine with two active cohorts + call BTranForHLMDiagnosticsFromCohortHydr(1, sites, bc_out) + + ! Analytical solution including fine-root carbon (fnrt): + ! Cohort 1 biomass = 10.0 kg, weighted btran = 0.8 * 10.0 * 1.0 = 8.0 + ! Cohort 2 biomass = 5.0 kg, weighted btran = 0.2 * 5.0 * 2.0 = 2.0 + ! Total live patch biomass = 10.0*1.0 + 5.0*2.0 = 20.0 kg/m2 + ! Expected btran_pa = (8.0 + 2.0) / 20.0 = 0.50 + expected_btran_pa = 0.50_r8 + @assertEqual(expected_btran_pa, bc_out(1)%btran_pa(1), tolerance=tol) + + deallocate(bc_out(1)%btran_pa) + deallocate(ccohort2%co_hydr) + deallocate(ccohort2) + end subroutine test_BTranForHLMDiagnostics_MultiCohortWeighted + + ! Summary: Tests BTranForHLMDiagnostics_explicit pure calculation function with explicit array inputs. + @Test + subroutine test_BTranForHLMDiagnostics_explicit(this) + class(TestPlantHydraulicsDrive), intent(inout) :: this + real(r8) :: btran_cohort(2) + real(r8) :: balive_cohort(2) + real(r8) :: n_indiv(2) + real(r8) :: result_btran + real(r8), parameter :: tol = 1.e-6_r8 + + ! Cohort 1: btran = 0.8, balive = 10.0, n = 1.0 => weighted btran = 8.0 + btran_cohort(1) = 0.8_r8 + balive_cohort(1) = 10.0_r8 + n_indiv(1) = 1.0_r8 + + ! Cohort 2: btran = 0.2, balive = 5.0, n = 2.0 => weighted btran = 2.0 + btran_cohort(2) = 0.2_r8 + balive_cohort(2) = 5.0_r8 + n_indiv(2) = 2.0_r8 + + ! Execute explicit subroutine (total biomass = 10*1 + 5*2 = 20, expected btran = (8 + 2)/20 = 0.50) + call BTranForHLMDiagnostics_explicit(2, btran_cohort, balive_cohort, n_indiv, result_btran) + @assertEqual(0.50_r8, result_btran, tolerance=tol) + + ! Zero biomass edge case returns 0.0 + balive_cohort(:) = 0.0_r8 + call BTranForHLMDiagnostics_explicit(2, btran_cohort, balive_cohort, n_indiv, result_btran) + @assertEqual(0.0_r8, result_btran, tolerance=tol) + end subroutine test_BTranForHLMDiagnostics_explicit + + ! Summary: Tests water mass balance conservation when accumulating water storage from cohort mortality. + @Test + subroutine test_AccumulateMortalityWaterStorage(this) + use FatesConstantsMod, only: denh2o => dens_fresh_liquid_water + use EDTypesMod, only: AREA_INV + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: csite + type(fates_cohort_type), target :: ccohort + type(ed_site_hydr_type), target :: csite_hydr + type(ed_cohort_hydr_type), target :: co_hydr + ! Cohort mortality density delta_n = 1.0 indiv/m2 + real(r8) :: delta_n = 1.0_r8 + real(r8) :: expected_delta_w + + csite%si_hydr => csite_hydr + ccohort%co_hydr => co_hydr + + ! Compartment volumes: v_ag = 1.0 m3, v_troot = 2.0 m3, v_aroot = 3.0 m3 + ! Compartment water contents: theta_ag = 0.5, theta_troot = 0.5, theta_aroot = 0.5 m3/m3 + allocate(co_hydr%th_aroot(1), co_hydr%v_aroot_layer(1)) + co_hydr%th_ag(:) = 0.5_r8 + co_hydr%v_ag(:) = 1.0_r8 + co_hydr%th_troot = 0.5_r8 + co_hydr%v_troot = 2.0_r8 + co_hydr%th_aroot(1) = 0.5_r8 + co_hydr%v_aroot_layer(1) = 3.0_r8 + + ! Initial dead vegetation water pool = 0.0 kg/m2, initial live vegetation water pool = 100.0 kg/m2 + csite_hydr%h2oveg_dead = 0.0_r8 + csite_hydr%h2oveg = 100.0_r8 + + ! Execute mortality water accumulation routine + call AccumulateMortalityWaterStorage(csite, ccohort, delta_n) + + ! Mass conservation analytical calculation: + ! delta_w = sum(theta_i * V_i) * denh2o * delta_n * AREA_INV + ! delta_w = (0.5*1.0 + 0.5*2.0 + 0.5*3.0) * 1000.0 * 1.0 * 1.0 = 3.0 * 1000.0 = 3000.0 kg/m2 + expected_delta_w = (sum(co_hydr%th_ag(:)*co_hydr%v_ag(:)) + & + co_hydr%th_troot*co_hydr%v_troot + & + sum(co_hydr%th_aroot(:)*co_hydr%v_aroot_layer(:))) * & + denh2o * delta_n * AREA_INV + + ! Verify exact mass balance transfer from live to dead pool + @assertEqual(expected_delta_w, csite_hydr%h2oveg_dead) + @assertEqual(100.0_r8 - expected_delta_w, csite_hydr%h2oveg) + + ! Invariant A: Total water mass conservation during mortality transfer (initial_total == final_total) + @assertEqual(100.0_r8, csite_hydr%h2oveg + csite_hydr%h2oveg_dead, tolerance=1.e-8_r8) + + deallocate(co_hydr%th_aroot, co_hydr%v_aroot_layer) + end subroutine test_AccumulateMortalityWaterStorage + + ! Summary: Tests AccumulateMortalityWaterStorage absorbing root vector dot product across multiple soil layers. + @Test + subroutine test_AccumulateMortalityWaterStorage_MultiLayer(this) + use FatesConstantsMod, only: denh2o => dens_fresh_liquid_water + use EDTypesMod, only: AREA_INV + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: csite + type(fates_cohort_type), target :: ccohort + type(ed_site_hydr_type), target :: csite_hydr + type(ed_cohort_hydr_type), target :: co_hydr + real(r8) :: delta_n = 1.0_r8 + real(r8) :: expected_delta_w + real(r8), parameter :: tol = 1.e-6_r8 + + csite%si_hydr => csite_hydr + ccohort%co_hydr => co_hydr + + ! Multi-layer soil setup (nlevrhiz = 2): + ! th_aroot = [0.2, 0.8] m3/m3, v_aroot_layer = [3.0, 1.0] m3 + ! Dot product = 0.2*3.0 + 0.8*1.0 = 1.4 m3 + ! Product of sums = (1.0)*(4.0) = 4.0 m3 (Mutation 7 fails) + allocate(co_hydr%th_aroot(2), co_hydr%v_aroot_layer(2)) + co_hydr%th_ag(:) = 0.5_r8 + co_hydr%v_ag(:) = 1.0_r8 + co_hydr%th_troot = 0.5_r8 + co_hydr%v_troot = 2.0_r8 + + co_hydr%th_aroot(1) = 0.2_r8 + co_hydr%th_aroot(2) = 0.8_r8 + co_hydr%v_aroot_layer(1) = 3.0_r8 + co_hydr%v_aroot_layer(2) = 1.0_r8 + + csite_hydr%h2oveg_dead = 0.0_r8 + csite_hydr%h2oveg = 10000.0_r8 + + call AccumulateMortalityWaterStorage(csite, ccohort, delta_n) + + ! Analytical expected delta_w: + ! (0.5*1.0 + 0.5*2.0 + (0.2*3.0 + 0.8*1.0)) * 1000.0 * 1.0 * 0.0001 = 3.4 * 0.1 = 0.34 kg/m2 + expected_delta_w = 0.34_r8 + @assertEqual(expected_delta_w, csite_hydr%h2oveg_dead, tolerance=tol) + + deallocate(co_hydr%th_aroot, co_hydr%v_aroot_layer) + end subroutine test_AccumulateMortalityWaterStorage_MultiLayer + + ! Summary: Tests AccumulateMortalityWater_explicit pure calculation function with explicit array inputs. + @Test + subroutine test_AccumulateMortalityWater_explicit(this) + use FatesConstantsMod, only: denh2o => dens_fresh_liquid_water + use EDTypesMod, only: AREA_INV + class(TestPlantHydraulicsDrive), intent(inout) :: this + real(r8) :: delta_n + real(r8) :: th_ag(1), v_ag(1) + real(r8) :: th_troot, v_troot + real(r8) :: th_aroot(2), v_aroot_layer(2) + real(r8) :: delta_w + real(r8), parameter :: tol = 1.e-6_r8 + + delta_n = 1.0_r8 + th_ag(1) = 0.5_r8 + v_ag(1) = 1.0_r8 + th_troot = 0.5_r8 + v_troot = 2.0_r8 + th_aroot(1) = 0.2_r8; v_aroot_layer(1) = 3.0_r8 + th_aroot(2) = 0.8_r8; v_aroot_layer(2) = 1.0_r8 + + call AccumulateMortalityWater_explicit(delta_n, th_ag, v_ag, th_troot, v_troot, & + th_aroot, v_aroot_layer, delta_w) + ! Analytical expected delta_w: + ! (0.5*1.0 + 0.5*2.0 + (0.2*3.0 + 0.8*1.0)) * 1000.0 * 1.0 * 0.0001 = 2.9 * 0.1 = 0.29 kg/m2 + @assertEqual(0.29_r8, delta_w, tolerance=tol) + end subroutine test_AccumulateMortalityWater_explicit + + + ! Summary: Tests recruit water storage pool initialization during recruitment events. + @Test + subroutine test_RecruitWaterStorage(this) + use FatesConstantsMod, only: denh2o => dens_fresh_liquid_water + use EDTypesMod, only: AREA_INV + use FatesInterfaceTypesMod, only: hlm_use_planthydro + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: sites(1) + type(fates_patch_type), target :: cpatch + type(fates_cohort_type), target :: ccohort + type(ed_site_hydr_type), target :: csite_hydr + type(ed_cohort_hydr_type), target :: co_hydr + real(r8) :: expected_recruit_w + + hlm_use_planthydro = 1 + sites(1)%si_hydr => csite_hydr + sites(1)%oldest_patch => cpatch + cpatch%tallest => ccohort + cpatch%younger => null() + ccohort%shorter => null() + ccohort%co_hydr => co_hydr + + ! Cohort recruited status = true, recruit density n = 2.0 indiv/m2 + co_hydr%is_newly_recruited = .true. + ccohort%n = 2.0_r8 + + ! Compartment volumes: v_ag = 1.0 m3, v_troot = 2.0 m3, v_aroot = 3.0 m3 + allocate(co_hydr%th_aroot(1), co_hydr%v_aroot_layer(1)) + co_hydr%th_ag(:) = 0.5_r8 + co_hydr%v_ag(:) = 1.0_r8 + co_hydr%th_troot = 0.5_r8 + co_hydr%v_troot = 2.0_r8 + co_hydr%th_aroot(1) = 0.5_r8 + co_hydr%v_aroot_layer(1) = 3.0_r8 + + csite_hydr%h2oveg_recruit = 0.0_r8 + + ! Execute recruit water storage update + call RecruitWaterStorage(1, sites) + + ! Analytical solution: recruit_w = (0.5*1 + 0.5*2 + 0.5*3) * 1000 * 2.0 * 1.0 = 6000.0 kg/m2 + expected_recruit_w = (sum(co_hydr%th_ag(:)*co_hydr%v_ag(:)) + & + co_hydr%th_troot*co_hydr%v_troot + & + sum(co_hydr%th_aroot(:)*co_hydr%v_aroot_layer(:))) * & + denh2o * ccohort%n * AREA_INV + + ! Verify exact recruit water storage value + @assertEqual(expected_recruit_w, csite_hydr%h2oveg_recruit) + + deallocate(co_hydr%th_aroot, co_hydr%v_aroot_layer) + end subroutine test_RecruitWaterStorage + + ! Summary: Tests total vegetation water storage summation across all live plant cohorts. + @Test + subroutine test_UpdateH2OVeg(this) + use FatesConstantsMod, only: denh2o => dens_fresh_liquid_water + use EDTypesMod, only: AREA_INV + use FatesInterfaceTypesMod, only: hlm_use_planthydro + class(TestPlantHydraulicsDrive), intent(inout) :: this + type(ed_site_type), target :: csite + type(fates_patch_type), target :: cpatch + type(fates_cohort_type), target :: ccohort + type(ed_site_hydr_type), target :: csite_hydr + type(ed_cohort_hydr_type), target :: co_hydr + type(bc_out_type) :: bc_out + real(r8) :: expected_h2oveg + + hlm_use_planthydro = 1 + csite%si_hydr => csite_hydr + csite%oldest_patch => cpatch + cpatch%tallest => ccohort + cpatch%younger => null() + ccohort%shorter => null() + ccohort%co_hydr => co_hydr + + ! Cohort properties: density n = 10.0 indiv/m2 + co_hydr%is_newly_recruited = .false. + allocate(co_hydr%th_aroot(1), co_hydr%v_aroot_layer(1)) + co_hydr%th_ag(:) = 0.5_r8 + co_hydr%v_ag(:) = 1.0_r8 + co_hydr%th_troot = 0.5_r8 + co_hydr%v_troot = 2.0_r8 + co_hydr%th_aroot(1) = 0.5_r8 + co_hydr%v_aroot_layer(1) = 3.0_r8 + ccohort%n = 10.0_r8 + + ! Zero error terms before diagnostic summation + csite_hydr%h2oveg_dead = 0.0_r8 + csite_hydr%h2oveg_growturn_err = 0.0_r8 + csite_hydr%h2oveg_hydro_err = 0.0_r8 + csite_hydr%trans_err = 0.0_r8 + csite_hydr%h2oveg = 0.0_r8 + + ! Execute vegetation water diagnostic update + call UpdateH2OVeg(csite, bc_out) + + ! Analytical expected total water storage = 30,000 kg/m2 + expected_h2oveg = (sum(co_hydr%th_ag(:)*co_hydr%v_ag(:)) + & + co_hydr%th_troot*co_hydr%v_troot + & + sum(co_hydr%th_aroot(:)*co_hydr%v_aroot_layer(:))) * & + denh2o * ccohort%n * AREA_INV + + ! Verify calculated plant stored water matches analytical total + @assertEqual(expected_h2oveg, bc_out%plant_stored_h2o_si) + + deallocate(co_hydr%th_aroot, co_hydr%v_aroot_layer) + end subroutine test_UpdateH2OVeg + +end module FatesPlantHydraulicsMod_drive_test + diff --git a/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_geom_test.pf b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_geom_test.pf new file mode 100644 index 0000000000..3440bfdf97 --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_geom_test.pf @@ -0,0 +1,204 @@ +! 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 + 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, fnrt_c = 0.5 kgC, sapw_c = 5.0 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, & + 0.5_r8, & + 5.0_r8, & + 10.0_r8) + + ! Analytical / allometric target values: + ! Leaf volume (v_ag(1)): 0.00329127901571 m3 (calculated from denleaf and c2b) + ! Stem sapwood volume (v_ag(2)): 0.0000675 m3 (a_sapwood * z_stem) + ! Transport root volume (v_troot): 0.02770796326795 m3 + @assertEqual(0.0032912790157100975_r8, ccohort%co_hydr%v_ag(1), tolerance=tol) + @assertEqual(0.0000675_r8, ccohort%co_hydr%v_ag(2), tolerance=tol) + @assertEqual(0.02770796326795_r8, ccohort%co_hydr%v_troot, tolerance=tol) + + ! Total fine root length = 10,000 m distributed across soil layers + @assertEqual(10000.0_r8, sum(ccohort%co_hydr%l_aroot_layer), tolerance=1.0_r8) + ! Total fine root volume = pi * r^2 * L_rhiz = pi * (0.001)^2 * 8820.0 = 0.02770796326795 m3 + @assertEqual(0.02770796326795_r8, sum(ccohort%co_hydr%v_aroot_layer), tolerance=1.e-4_r8) + + call ccohort%co_hydr%DeallocateHydrCohortArrays() + deallocate(ccohort%co_hydr) + deallocate(csite_hydr%zi_rhiz) + deallocate(csite_hydr%dz_rhiz) + end subroutine test_UpdatePlantHydrLenVol + + ! Summary: Tests copying current compartment volumes into previous timestep volume arrays for size-dependent updates. + @Test + subroutine test_SavePreviousCompartmentVolumes(this) + class(TestPlantHydraulicsGeom), intent(inout) :: this + type(ed_cohort_hydr_type) :: ccohort_hydr + + allocate(ccohort_hydr%v_aroot_layer(2)) + allocate(ccohort_hydr%v_aroot_layer_init(2)) + + ! Current volume inputs: v_ag(1) = 5.0 m3, v_troot = 10.0 m3, v_aroot = [1.0 m3, 2.0 m3] + ccohort_hydr%v_ag(1) = 5.0_r8 + ccohort_hydr%v_troot = 10.0_r8 + ccohort_hydr%v_aroot_layer(1) = 1.0_r8 + ccohort_hydr%v_aroot_layer(2) = 2.0_r8 + + ! Initial target arrays set to zero + ccohort_hydr%v_ag_init(:) = 0.0_r8 + ccohort_hydr%v_troot_init = 0.0_r8 + ccohort_hydr%v_aroot_layer_init(:) = 0.0_r8 + + ! Execute save routine + call SavePreviousCompartmentVolumes(ccohort_hydr) + + ! Verify exact copy of current volumes to _init variables + @assertEqual(5.0_r8, ccohort_hydr%v_ag_init(1)) + @assertEqual(10.0_r8, ccohort_hydr%v_troot_init) + @assertEqual(1.0_r8, ccohort_hydr%v_aroot_layer_init(1)) + @assertEqual(2.0_r8, ccohort_hydr%v_aroot_layer_init(2)) + + deallocate(ccohort_hydr%v_aroot_layer) + deallocate(ccohort_hydr%v_aroot_layer_init) + end subroutine test_SavePreviousCompartmentVolumes + +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..1cae6dec39 --- /dev/null +++ b/testing/tests/unit/fates_plant_hydraulics_test/FatesPlantHydraulicsMod_sizedep_test.pf @@ -0,0 +1,243 @@ +! 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 residual moisture floor clamping in UpdateSizeDepPlantHydStates when organ volume expands rapidly. + @Test + subroutine test_UpdateSizeDepPlantHydStates_ResidualClamping(this) + class(TestPlantHydraulicsSizedep), intent(inout) :: this + type(ed_site_type), target :: currentSite + type(fates_cohort_type), target :: ccohort + integer :: nlevrhiz = 1 + integer :: ft = 1 + real(r8) :: expected_clamped_theta + ! Small theta offset parameter used in constrain_water_contents to prevent exact boundary singularities [m3/m3] + real(r8), parameter :: small_theta_num = 1.e-7_r8 + + 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 + ccohort%n = 1.0_r8 + allocate(ccohort%co_hydr) + call ccohort%co_hydr%AllocateHydrCohortArrays(nlevrhiz) + + ! Large volume expansion: v_init = 1.0 m3, v_curr = 100.0 m3, th_init = 0.10 m3/m3 + ! Uncorrected theta = 0.10 * (1.0 / 100.0) = 0.001 m3/m3, which drops below 0.05 residual floor + ccohort%co_hydr%v_ag = 100.0_r8 + ccohort%co_hydr%v_ag_init = 1.0_r8 + ccohort%co_hydr%th_ag = 0.10_r8 + + ccohort%co_hydr%v_troot = 100.0_r8 + ccohort%co_hydr%v_troot_init = 1.0_r8 + ccohort%co_hydr%th_troot = 0.10_r8 + + ccohort%co_hydr%v_aroot_layer(:) = 100.0_r8 + ccohort%co_hydr%v_aroot_layer_init(:) = 1.0_r8 + ccohort%co_hydr%th_aroot(:) = 0.10_r8 + + call UpdateSizeDepPlantHydStates(currentSite, ccohort) + + ! Expected moisture content clamped to residual threshold + small_theta_num offset: 0.05 + 1.e-7 = 0.0500001 m3/m3 + expected_clamped_theta = 0.05_r8 + small_theta_num + + @assertEqual(expected_clamped_theta, ccohort%co_hydr%th_ag(1), tolerance=1.e-8_r8) + @assertEqual(expected_clamped_theta, ccohort%co_hydr%th_troot, tolerance=1.e-8_r8) + @assertEqual(expected_clamped_theta, ccohort%co_hydr%th_aroot(1), tolerance=1.e-8_r8) + + call ccohort%co_hydr%DeallocateHydrCohortArrays() + deallocate(ccohort%co_hydr) + deallocate(currentSite%si_hydr) + end subroutine test_UpdateSizeDepPlantHydStates_ResidualClamping + + ! 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 +