Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fixes to hydraulics #525

Merged
merged 17 commits into from
May 30, 2019
Merged

fixes to hydraulics #525

Show file tree
Hide file tree
Changes from 10 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
59c7c31
change the tree hydro initial condition so that we do not assume equi…
xuchongang Apr 12, 2019
e5c422b
update the hydro codes so that it the soil-root conductance is higher…
xuchongang Apr 12, 2019
e5cb467
update soil-root interface hydraulic conductivity for water uptake vs…
xuchongang Apr 12, 2019
00dfc52
add the hydraulic conductivity from absorbing roots to inner soil she…
xuchongang Apr 12, 2019
f9e0dbe
update soil-root conductivity parameters so that the conductivity is …
xuchongang Apr 12, 2019
b87c41f
resolve a bug for the patch loop to calculate the soil-root interface…
xuchongang Apr 19, 2019
def6d20
update the code of inner soil shell hydraulic conductivity depending …
xuchongang Apr 22, 2019
e360ec2
update the comments on 1st soil shell hydraulic conductivity
xuchongang Apr 22, 2019
3b15841
update the hardwired threshold for stomata condutance to avoid too mu…
xuchongang Apr 25, 2019
4bf50a4
Trivial, changes some indents.
rgknox Apr 25, 2019
113a3ca
Recalculated bbb (cuticular conductance) so that the same value is us…
rgknox Apr 25, 2019
3e71640
Removed condition on bbb per discussion with C. Xu. Added a parameter…
rgknox Apr 25, 2019
f311ef5
Temporarily setting hydr_kmax_rsurf2 to be a hard-wired parameter con…
rgknox Apr 26, 2019
b4d643e
Small fix on declaring bbbopt as a local
rgknox Apr 26, 2019
ce61db1
Merge pull request #2 from rgknox/xuchongang/HYDRO_UPDATE_v2
rgknox Apr 29, 2019
7042533
Reverted parameter file to not-include kmax_rsurf2, will add in next …
rgknox May 29, 2019
5dbff93
Merge resolution with master, resolved kmax_rsurf1 name change in EDP…
rgknox May 29, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
189 changes: 130 additions & 59 deletions biogeophys/FatesPlantHydraulicsMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,8 +40,9 @@ module FatesPlantHydraulicsMod
use FatesConstantsMod, only : cm2_per_m2
use FatesConstantsMod, only : g_per_kg

use EDParamsMod , only : hydr_kmax_rsurf

use EDParamsMod , only : hydr_kmax_rsurf1
use EDParamsMod , only : hydr_kmax_rsurf2

use EDTypesMod , only : ed_site_type
use EDTypesMod , only : ed_patch_type
use EDTypesMod , only : ed_cohort_type
Expand Down Expand Up @@ -141,6 +142,7 @@ module FatesPlantHydraulicsMod
public :: CopyCohortHydraulics
public :: FuseCohortHydraulics
public :: updateSizeDepTreeHydProps
public :: updateWaterDepTreeHydProps
public :: updateSizeDepTreeHydStates
public :: initTreeHydStates
public :: updateSizeDepRhizHydProps
Expand All @@ -150,6 +152,7 @@ module FatesPlantHydraulicsMod
public :: SavePreviousRhizVolumes
public :: UpdateTreeHydrNodes
public :: UpdateTreeHydrLenVolCond
public :: UpdateWaterDepTreeHydrCond
public :: ConstrainRecruitNumber

!------------------------------------------------------------------------------
Expand Down Expand Up @@ -325,7 +328,8 @@ subroutine initTreeHydStates(site_p, cc_p, bc_in)
! watsat(c,j), watres(c,j), alpha_VG(c,j), n_VG(c,j), m_VG(c,j), l_VG(c,j), &
! smp)
!ccohort_hydr%psi_aroot(j) = smp
ccohort_hydr%psi_aroot(j) = csite%si_hydr%psisoi_liq_innershell(j)
!ccohort_hydr%psi_aroot(j) = csite%si_hydr%psisoi_liq_innershell(j)
ccohort_hydr%psi_aroot(j) = -0.2_r8 !do not assume the equalibrium between soil and root

call th_from_psi(ft, 4, ccohort_hydr%psi_aroot(j), ccohort_hydr%th_aroot(j), csite%si_hydr, bc_in )
end do
Expand All @@ -346,6 +350,7 @@ subroutine initTreeHydStates(site_p, cc_p, bc_in)
!hydrostatic equilibrium with the water potential immediately below
dz = ccohort_hydr%z_node_ag(n_hypool_ag) - ccohort_hydr%z_node_troot(1)
ccohort_hydr%psi_ag(n_hypool_ag) = ccohort_hydr%psi_troot(1) - 1.e-6_r8*denh2o*grav*dz
if (ccohort_hydr%psi_ag(n_hypool_ag)>0.0_r8) ccohort_hydr%psi_ag(n_hypool_ag) = -0.01_r8
call th_from_psi(ft, 2, ccohort_hydr%psi_ag(n_hypool_ag), ccohort_hydr%th_ag(n_hypool_ag), csite%si_hydr, bc_in)
do k=n_hypool_ag-1, 1, -1
dz = ccohort_hydr%z_node_ag(k) - ccohort_hydr%z_node_ag(k+1)
Expand Down Expand Up @@ -543,9 +548,44 @@ subroutine updateSizeDepTreeHydProps(currentSite,ccohort,bc_in)
! volumes, and UpdateTreeHydrNodes is called prior to this.

call UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in)

end subroutine updateSizeDepTreeHydProps

! =====================================================================================

subroutine updateWaterDepTreeHydProps(currentSite,ccohort,bc_in)


! DESCRIPTION: Updates absorbing root length (total and its vertical distribution)
! as well as the consequential change in the size of the 'representative' rhizosphere
! shell radii, volumes, and compartment volumes of plant tissues

! !USES:
use shr_sys_mod , only : shr_sys_abort

! ARGUMENTS:
type(ed_site_type) , intent(in) :: currentSite ! Site stuff
type(ed_cohort_type) , intent(inout) :: ccohort ! current cohort pointer
type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions

! Locals
integer :: nlevsoi_hyd ! Number of total soil layers
type(ed_cohort_hydr_type), pointer :: ccohort_hydr
integer :: ft

nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd
ccohort_hydr => ccohort%co_hydr
ft = ccohort%pft

! This updates plant compartment volumes, lengths and
! maximum conductances. Make sure for already
! initialized vegetation, that SavePreviousCompartment
! volumes, and UpdateTreeHydrNodes is called prior to this.

call UpdateWaterDepTreeHydrCond(currentSite,ccohort,nlevsoi_hyd,bc_in)

end subroutine updateSizeDepTreeHydProps

end subroutine updateWaterDepTreeHydProps

! =====================================================================================

Expand Down Expand Up @@ -616,7 +656,6 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in)
! hydraulic conductance [kg s-1 MPa-1]
real(r8) :: kmax_tot ! total tree (leaf to root tip)
! hydraulic conductance [kg s-1 MPa-1]

real(r8),parameter :: taper_exponent = 1._r8/3._r8 ! Savage et al. (2010) xylem taper exponent [-]

ccohort_hydr => ccohort%co_hydr
Expand Down Expand Up @@ -810,10 +849,80 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in)
ccohort_hydr%kmax_treebg_layer(j) = rootfr*ccohort_hydr%kmax_treebg_tot
end do
end if

end if !check for bleaf

end subroutine UpdateTreeHydrLenVolCond



!=====================================================================================

subroutine UpdateWaterDepTreeHydrCond(currentSite,ccohort,nlevsoi_hyd,bc_in)

! -----------------------------------------------------------------------------------
! This subroutine calculates update the conductivity for the soil-root interface,
! depending on the plant water uptake/loss.
! we assume that the conductivitity for water uptake is larger than
! water loss due to composite regulation of resistance the roots
! hydraulic vs osmostic with and without transpiration
! Steudle, E. Water uptake by roots: effects of water deficit.
! J Exp Bot 51, 1531-1542, doi:DOI 10.1093/jexbot/51.350.1531 (2000).
! -----------------------------------------------------------------------------------

! Arguments
type(ed_site_type) , intent(in) :: currentSite ! Site target
type(ed_cohort_type),intent(inout) :: ccohort ! cohort target
integer,intent(in) :: nlevsoi_hyd ! number of soil hydro layers
type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions

type(ed_cohort_hydr_type),pointer :: ccohort_hydr ! Plant hydraulics structure
type(ed_site_hydr_type),pointer :: csite_hydr

integer :: j,k
real(r8) :: hksat_s ! hksat converted to units of 10^6sec
real(r8) :: kmax_root_surf_total ! maximum conducitivity for total root surface(kg water/Mpa/s)
real(r8) :: kmax_soil_total ! maximum conducitivity for from root surface to soil shell(kg water/Mpa/s)
! which is equiv to [kg m-1 s-1 MPa-1]
real(r8) :: kmax_root_surf ! maximum conducitivity for unit root surface (kg water/m2 root area/Mpa/s)

ccohort_hydr => ccohort%co_hydr
csite_hydr => currentSite%si_hydr
k = 1 !only for the first soil shell
do j=1, nlevsoi_hyd

hksat_s = bc_in%hksat_sisl(j) * 1.e-3_r8 * 1/grav * 1.e6_r8
if(ccohort_hydr%psi_aroot(j)<csite_hydr%psisoi_liq_innershell(j))then
kmax_root_surf = hydr_kmax_rsurf1
else
kmax_root_surf = hydr_kmax_rsurf2
endif
kmax_root_surf_total = kmax_root_surf*2._r8*pi_const *csite_hydr%rs1(j)* &
csite_hydr%l_aroot_layer(j)
if(csite_hydr%r_node_shell(j,1) <= csite_hydr%rs1(j)) then
!csite_hydr%kmax_upper_shell(j,k) = large_kmax_bound
!csite_hydr%kmax_bound_shell(j,k) = large_kmax_bound
!csite_hydr%kmax_lower_shell(j,k) = large_kmax_bound
ccohort_hydr%kmax_innershell(j) = kmax_root_surf_total

else

kmax_soil_total = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s

!csite_hydr%kmax_upper_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
!csite_hydr%kmax_bound_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
!csite_hydr%kmax_lower_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s

ccohort_hydr%kmax_innershell(j) = (1._r8/kmax_root_surf_total + &
1._r8/kmax_soil_total)**(-1._r8)
end if
end do

end subroutine UpdateWaterDepTreeHydrCond

! =====================================================================================
subroutine updateSizeDepTreeHydStates(currentSite,ccohort)
Expand Down Expand Up @@ -951,6 +1060,7 @@ subroutine CopyCohortHydraulics(newCohort, oldCohort)
! quantities indexed by soil layer
ncohort_hydr%z_node_aroot = ocohort_hydr%z_node_aroot
ncohort_hydr%kmax_treebg_layer = ocohort_hydr%kmax_treebg_layer
ncohort_hydr%kmax_innershell = ocohort_hydr%kmax_innershell
ncohort_hydr%v_aroot_layer_init = ocohort_hydr%v_aroot_layer_init
ncohort_hydr%v_aroot_layer = ocohort_hydr%v_aroot_layer
ncohort_hydr%l_aroot_layer = ocohort_hydr%l_aroot_layer
Expand Down Expand Up @@ -1609,11 +1719,7 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in)
real(r8) :: large_kmax_bound = 1.e4_r8 ! for replacing kmax_bound_shell wherever the
! innermost shell radius is less than the assumed
! absorbing root radius rs1
real(r8) :: kmax_root_surf ! maximum conducitivity for unit root surface
! (kg water/m2 root area/Mpa/s)
! 1.e-5_r8 from Rudinger et al 1994
real(r8) :: kmax_root_surf_total !maximum conducitivity for total root surface(kg water/Mpa/s)
real(r8) :: kmax_soil_total !maximum conducitivity for total root surface(kg water/Mpa/s)
integer :: nlevsoi_hyd

!-----------------------------------------------------------------------
Expand All @@ -1633,7 +1739,7 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in)
enddo !cohort
cPatch => cPatch%older
enddo !patch
kmax_root_surf = hydr_kmax_rsurf

csite_hydr%l_aroot_1D = sum( csite_hydr%l_aroot_layer(:))

! update outer radii of column-level rhizosphere shells (same across patches and cohorts)
Expand All @@ -1646,6 +1752,9 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in)
enddo
call shellGeom( csite_hydr%l_aroot_1D, csite_hydr%rs1(1), AREA, sum(bc_in%dz_sisl(1:nlevsoi_hyd)), &
csite_hydr%r_out_shell_1D(:), csite_hydr%r_node_shell_1D(:), csite_hydr%v_shell_1D(:))

!update the conductitivity for first soil shell is done at subroutine UpdateWaterDepTreeHydrCond
!which is dependant on whether it is water uptake or loss for every 30 minutes

do j = 1,csite_hydr%nlevsoi_hyd

Expand All @@ -1654,50 +1763,8 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in)
! proceed only if the total absorbing root length (site-level) has changed in this layer
if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then

do k = 1,nshell
if(k == 1) then
kmax_root_surf_total = kmax_root_surf*2._r8*pi_const *csite_hydr%rs1(j)* &
csite_hydr%l_aroot_layer(j)
if(csite_hydr%r_node_shell(j,k) <= csite_hydr%rs1(j)) then
!csite_hydr%kmax_upper_shell(j,k) = large_kmax_bound
!csite_hydr%kmax_bound_shell(j,k) = large_kmax_bound
!csite_hydr%kmax_lower_shell(j,k) = large_kmax_bound
csite_hydr%kmax_upper_shell(j,k) = kmax_root_surf_total
csite_hydr%kmax_bound_shell(j,k) = kmax_root_surf_total
csite_hydr%kmax_lower_shell(j,k) = kmax_root_surf_total

else

kmax_soil_total = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s

!csite_hydr%kmax_upper_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
!csite_hydr%kmax_bound_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
!csite_hydr%kmax_lower_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s

csite_hydr%kmax_upper_shell(j,k) = (1._r8/kmax_root_surf_total + &
1._r8/kmax_soil_total)**(-1._r8)
csite_hydr%kmax_bound_shell(j,k) = (1._r8/kmax_root_surf_total + &
1._r8/kmax_soil_total)**(-1._r8)
csite_hydr%kmax_lower_shell(j,k) = (1._r8/kmax_root_surf_total + &
1._r8/kmax_soil_total)**(-1._r8)
end if
if(j == 1) then
if(csite_hydr%r_node_shell(j,k) <= csite_hydr%rs1(j)) then
csite_hydr%kmax_upper_shell_1D(k) = csite_hydr%kmax_upper_shell(1,k)
csite_hydr%kmax_bound_shell_1D(k) = csite_hydr%kmax_bound_shell(1,k)
csite_hydr%kmax_lower_shell_1D(k) = csite_hydr%kmax_lower_shell(1,k)
else
csite_hydr%kmax_upper_shell_1D(k) = csite_hydr%kmax_upper_shell(1,k)
csite_hydr%kmax_bound_shell_1D(k) = csite_hydr%kmax_bound_shell(1,k)
csite_hydr%kmax_lower_shell_1D(k) = csite_hydr%kmax_lower_shell(1,k)
end if
end if
else
csite_hydr%kmax_upper_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
do k = 2,nshell
csite_hydr%kmax_upper_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
log(csite_hydr%r_node_shell(j,k)/csite_hydr%r_out_shell(j,k-1))*hksat_s
csite_hydr%kmax_bound_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
log(csite_hydr%r_node_shell(j,k)/csite_hydr%r_node_shell(j,k-1))*hksat_s
Expand All @@ -1711,16 +1778,14 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in)
csite_hydr%kmax_lower_shell_1D(k) = 2._r8*pi_const*csite_hydr%l_aroot_1D / &
log(csite_hydr%r_out_shell_1D( k)/csite_hydr%r_node_shell_1D(k ))*hksat_s
end if
end if
enddo ! loop over rhizosphere shells
end if !has l_aroot_layer changed?
enddo ! loop over soil layers



return
end subroutine UpdateSizeDepRhizVolLenCon


! =====================================================================================

Expand Down Expand Up @@ -2445,7 +2510,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime )
! [mm H2O/cohort/s] = [mm H2O / patch / s] / [cohort/patch]
!! qflx_tran_veg_patch_coh = qflx_trans_patch_vol * qflx_rel_tran_coh


call updateWaterDepTreeHydProps(sites(s),ccohort,bc_in(s))

if(site_hydr%nlevsoi_hyd > 1) then
! BUCKET APPROXIMATION OF THE SOIL-ROOT HYDRAULIC GRADIENT (weighted average across layers)
Expand Down Expand Up @@ -2516,6 +2581,9 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime )
kmax_lower( 1 : n_hypool_ag ) = ccohort_hydr%kmax_lower(:)
kmax_upper(( n_hypool_ag+1) ) = ccohort_hydr%kmax_upper_troot
if(site_hydr%nlevsoi_hyd == 1) then
site_hydr%kmax_upper_shell_1D(1) = ccohort_hydr%kmax_innershell(1)
site_hydr%kmax_lower_shell_1D(1) = ccohort_hydr%kmax_innershell(1)
site_hydr%kmax_bound_shell_1D(1) = ccohort_hydr%kmax_innershell(1)
!! estimate troot-aroot and aroot-radial components as a residual:
!! 25% each of total (surface of aroots to leaves) resistance
kmax_bound(( n_hypool_ag+1):(n_hypool_ag+2 )) = 2._r8 * ccohort_hydr%kmax_treebg_tot
Expand Down Expand Up @@ -2645,11 +2713,14 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime )
v_node_1l((n_hypool_ag+n_hypool_troot+1) ) = ccohort_hydr%v_aroot_layer(j)
v_node_1l((n_hypool_tot-nshell+1):(n_hypool_tot)) = site_hydr%v_shell(j,:) * &
ccohort_hydr%l_aroot_layer(j)/bc_in(s)%dz_sisl(j)
site_hydr%kmax_bound_shell(j,1)=ccohort_hydr%kmax_innershell(j)
site_hydr%kmax_upper_shell(j,1)=ccohort_hydr%kmax_innershell(j)
site_hydr%kmax_lower_shell(j,1)=ccohort_hydr%kmax_innershell(j)
kmax_bound_1l(:) = 0._r8
kmax_bound_shell_1l(:) = site_hydr%kmax_bound_shell(j,:) * &
ccohort_hydr%l_aroot_layer(j) / site_hydr%l_aroot_layer(j)


! transporting-to-absorbing root conductance: factor of 2 means one-half of the total
! belowground resistance in layer j
kmax_bound_1l((n_hypool_ag+1)) = 2._r8 * ccohort_hydr%kmax_treebg_layer(j)
Expand Down
Loading