fix(hydraulics): fix SumBetweenDepths layer integration logic and add pFUnit tests

biogeophys/FatesPlantHydraulicsMod.F90

@@ -116,6 +116,8 @@ module FatesPlantHydraulicsMod
 
   implicit none
 
+  private
+
 
   ! 1=leaf, 2=stem, 3=troot, 4=aroot
   ! Several of these may be better transferred to the parameter file in due time (RGK)
@@ -260,9 +262,15 @@ module FatesPlantHydraulicsMod
   public :: UpdateSizeDepRhizHydProps
   public :: RestartHydrStates
   public :: SavePreviousCompartmentVolumes
-  public :: SavePreviousRhizVolumes
+  public :: SumBetweenDepths ! only public for unit tests
+
   public :: UpdatePlantHydrNodes
   public :: UpdatePlantHydrLenVol
+
+  interface UpdatePlantHydrLenVol
+     module procedure UpdatePlantHydrLenVol_cohort
+     module procedure UpdatePlantHydrLenVol_masses
+  end interface UpdatePlantHydrLenVol
   public :: UpdatePlantKmax
   public :: ConstrainRecruitNumber
   public :: InitHydroGlobals
@@ -885,7 +893,25 @@ end subroutine UpdateSizeDepPlantHydProps
 
   ! =====================================================================================
 
-  subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
+  subroutine UpdatePlantHydrLenVol_cohort(ccohort, csite_hydr)
+    type(fates_cohort_type),intent(inout)  :: ccohort
+    type(ed_site_hydr_type),intent(in)  :: csite_hydr
+
+    call UpdatePlantHydrLenVol_masses(ccohort%co_hydr, csite_hydr, &
+         ccohort%pft, ccohort%dbh, ccohort%height, &
+         ccohort%crowndamage, ccohort%canopy_trim, ccohort%efstem_coh, &
+         ccohort%size_class, &
+         ccohort%prt%GetState(leaf_organ, carbon12_element), &
+         ccohort%prt%GetState(sapw_organ, carbon12_element), &
+         ccohort%prt%GetState(fnrt_organ, carbon12_element), &
+         ccohort%prt%GetState(struct_organ, carbon12_element))
+
+  end subroutine UpdatePlantHydrLenVol_cohort
+
+
+  subroutine UpdatePlantHydrLenVol_masses(ccohort_hydr,csite_hydr,pft,dbh,height, &
+                                          crowndamage,canopy_trim,efstem_coh,size_class, &
+                                          leaf_c,sapw_c,fnrt_c,struct_c)
 
     ! -----------------------------------------------------------------------------------
     ! This subroutine calculates two attributes of a plant:
@@ -900,19 +926,25 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
     ! -----------------------------------------------------------------------------------
 
     ! Arguments
-    type(fates_cohort_type),intent(inout)  :: ccohort
+    type(ed_cohort_hydr_type),intent(inout) :: ccohort_hydr
     type(ed_site_hydr_type),intent(in)  :: csite_hydr
+    integer, intent(in) :: pft
+    real(r8), intent(in) :: dbh
+    real(r8), intent(in) :: height
+    integer, intent(in) :: crowndamage
+    real(r8), intent(in) :: canopy_trim
+    real(r8), intent(in) :: efstem_coh
+    integer, intent(in) :: size_class
+    real(r8), intent(in) :: leaf_c                       ! Current amount of leaf carbon in the plant                            [kg]
+    real(r8), intent(in) :: sapw_c                       ! Current amount of sapwood carbon in the plant                         [kg]
+    real(r8), intent(in) :: fnrt_c                       ! Current amount of fine-root carbon in the plant                       [kg]
+    real(r8), intent(in) :: struct_c                     ! Current amount of structural carbon in the plant                      [kg]
 
-    type(ed_cohort_hydr_type),pointer :: ccohort_hydr     ! Plant hydraulics structure
     integer  :: j,k
     integer  :: ft                           ! Plant functional type index
     real(r8) :: roota                        ! root profile parameter a zeng2001_crootfr
     real(r8) :: rootb                        ! root profile parameter b zeng2001_crootfr
-    real(r8) :: leaf_c                       ! Current amount of leaf carbon in the plant                            [kg]
     real(r8) :: leaf_c_target                ! Target leaf carbon (with some conditions) [kgC]
-    real(r8) :: fnrt_c                       ! Current amount of fine-root carbon in the plant                       [kg]
-    real(r8) :: sapw_c                       ! Current amount of sapwood carbon in the plant                         [kg]
-    real(r8) :: struct_c                     ! Current amount of structural carbon in the plant                      [kg]
     real(r8) :: woody_bg_c                   ! belowground woody biomass in carbon units                             [kgC/indiv]
     real(r8) :: z_stem                       ! the height of the plants stem below crown [m]
     real(r8) :: sla                          ! specific leaf area                                                    [cm2/g]
@@ -928,7 +960,6 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
     real(r8) :: crown_depth                  ! Depth of the plant's crown [m]
     real(r8) :: norm                         ! total root fraction used <1
     integer  :: nlevrhiz                     ! number of rhizosphere levels
-    real(r8) :: dbh                          ! the dbh of current cohort                                             [cm]   
     real(r8) :: z_fr                         ! rooting depth of a cohort                                             [cm]
     real(r8) :: v_leaf_donate(1:n_hypool_leaf)   ! the volume that leaf will donate to xylem     
 
@@ -945,13 +976,8 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
     ! to estimate maximum leaf carbon
 
 
-    ccohort_hydr => ccohort%co_hydr
-    ft           = ccohort%pft
+    ft           = pft
     nlevrhiz     = csite_hydr%nlevrhiz
-    leaf_c       = ccohort%prt%GetState(leaf_organ, carbon12_element)
-    sapw_c       = ccohort%prt%GetState(sapw_organ, carbon12_element)
-    fnrt_c       = ccohort%prt%GetState(fnrt_organ, carbon12_element)
-    struct_c     = ccohort%prt%GetState(struct_organ, carbon12_element)
     roota        = prt_params%fnrt_prof_a(ft)
     rootb        = prt_params%fnrt_prof_b(ft)
 
@@ -989,8 +1015,8 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
     ! Lets also avoid super-low targets that have very low trimming functions
 
     ! efleaf_coh hard-coded to 1 in the call below to avoid zero leaf volume
-    call bleaf(ccohort%dbh,ccohort%pft,ccohort%crowndamage, &
-         max(ccohort%canopy_trim,min_trim),1.0_r8, leaf_c_target)
+    call bleaf(dbh,ft,crowndamage, &
+         max(canopy_trim,min_trim),1.0_r8, leaf_c_target)
 
     ccohort_hydr%v_ag(1:n_hypool_leaf) = max(leaf_c,min_leaf_frac*leaf_c_target) * &
          prt_params%c2b(ft) / denleaf/ real(n_hypool_leaf,r8)
@@ -1004,19 +1030,19 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
     ! v_stem       = c_stem_biom / (prt_params%wood_density(ft) * kg_per_g * cm3_per_m3 )
 
     ! calculate the sapwood cross-sectional area
-    call bsap_allom(ccohort%dbh,ccohort%pft,ccohort%crowndamage, &
-         ccohort%canopy_trim, ccohort%efstem_coh, a_sapwood_target,sapw_c_target)
+    call bsap_allom(dbh,ft,crowndamage, &
+         canopy_trim, efstem_coh, a_sapwood_target,sapw_c_target)
 
     ! uncomment this if you want to use
     ! the actual sapwood, which may be lower than target due to branchfall.
     a_sapwood = a_sapwood_target  ! * sapw_c / sapw_c_target
 
     ! alternative cross section calculation
-    ! a_sapwood    = a_leaf_tot / ( 0.001_r8 + 0.025_r8 * ccohort%height ) * 1.e-4_r8
+    ! a_sapwood    = a_leaf_tot / ( 0.001_r8 + 0.025_r8 * height ) * 1.e-4_r8
 
-    !call CrownDepth(ccohort%height,ft,crown_depth)
-    crown_depth  = min(ccohort%height,0.1_r8)
-    z_stem       = ccohort%height - crown_depth
+    !call CrownDepth(height,ft,crown_depth)
+    crown_depth  = min(height,0.1_r8)
+    z_stem       = height - crown_depth
     v_sapwood    = a_sapwood * z_stem    ! + 0.333_r8*a_sapwood*crown_depth
 
     ! Junyan changed the following code to calculate the above ground node volume
@@ -1065,7 +1091,7 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
     ! calculations.
 
 
-    call MaximumRootingDepth(ccohort%dbh,ft,csite_hydr%zi_rhiz(nlevrhiz),z_fr)
+    call MaximumRootingDepth(dbh,ft,csite_hydr%zi_rhiz(nlevrhiz),z_fr)
 
     do j=1,nlevrhiz
 
@@ -1074,7 +1100,7 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
 
        if(debug)then
           write(fates_log(),*) 'check rooting depth of cohort '
-          write(fates_log(),*) 'dbh: ',ccohort%dbh,' sice class: ',ccohort%size_class
+          write(fates_log(),*) 'dbh: ',dbh,' sice class: ',size_class
           write(fates_log(),*) 'csite_hydr%dz_rhiz(j) is: ', csite_hydr%dz_rhiz(j)
           write(fates_log(),*) 'z_max cohort: ',z_fr
           write(fates_log(),*) 'layer:  ',j,' bottom depth (m): ',csite_hydr%zi_rhiz(j),' rooting fraction:',rootfr
@@ -1090,7 +1116,7 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
     end do
 
     return
-  end subroutine UpdatePlantHydrLenVol
+  end subroutine UpdatePlantHydrLenVol_masses
 
   ! =====================================================================================
 
@@ -5603,16 +5629,26 @@ function SumBetweenDepths(csite_hydr,depth_t,depth_b,array_in) result(depth_sum)
    depth_sum = depth_sum + sum(array_in(i_rhiz_t:i_rhiz_b))
    end if
 
-   ! Find fraction contribution from top partial layer (if any)
-   if(i_rhiz_t>1) then
-   frac = (csite_hydr%zi_rhiz(i_rhiz_t-1)-depth_t)/csite_hydr%dz_rhiz(i_rhiz_t-1)
-   depth_sum = depth_sum + frac*array_in(i_rhiz_t-1)
-   end if
+   if (i_rhiz_t - 1 == i_rhiz_b + 1) then
+       ! Top and bottom depths fall within the same layer
+       frac = (depth_b - depth_t) / csite_hydr%dz_rhiz(i_rhiz_t-1)
+       depth_sum = depth_sum + frac*array_in(i_rhiz_t-1)
+   else
+       ! Find fraction contribution from top partial layer (if any)
+       if(i_rhiz_t>1) then
+           frac = (csite_hydr%zi_rhiz(i_rhiz_t-1)-depth_t)/csite_hydr%dz_rhiz(i_rhiz_t-1)
+           depth_sum = depth_sum + frac*array_in(i_rhiz_t-1)
+       end if
 
-   ! Find fraction contribution from bottom partial layer (if any)
-   if(i_rhiz_b<nlevrhiz) then
-   frac = (depth_b-csite_hydr%zi_rhiz(i_rhiz_b))/csite_hydr%dz_rhiz(i_rhiz_b+1)
-   depth_sum = depth_sum + frac*array_in(i_rhiz_b+1)
+       ! Find fraction contribution from bottom partial layer (if any)
+       if(i_rhiz_b<nlevrhiz) then
+           if(i_rhiz_b == 0) then
+               frac = depth_b / csite_hydr%dz_rhiz(1)
+           else
+               frac = (depth_b-csite_hydr%zi_rhiz(i_rhiz_b))/csite_hydr%dz_rhiz(i_rhiz_b+1)
+           end if
+           depth_sum = depth_sum + frac*array_in(i_rhiz_b+1)
+       end if
    end if
 
    depth_sum = depth_sum/(min(depth_b,csite_hydr%zi_rhiz(nlevrhiz))-depth_t)

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)

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
 

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

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)

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

testing/tests/fortran_shr/CMakeLists.txt

@@ -4,6 +4,8 @@ list(APPEND fates_sources
   FatesArgumentUtils.F90
   FatesFactoryMod.F90
   SyntheticPatchTypes.F90
-  FatesUnitTestUtils.F90)
+  FatesUnitTestUtils.F90
+  EDPftvarconMockMod.F90
+  FatesUnitTestPRTParametersMod.F90)
 
-sourcelist_to_parent(fates_sources)
+sourcelist_to_parent(fates_sources)

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