From 2e8608cfa24fb006eef8b9112eb91d3265e687e4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 22 Feb 2024 16:04:58 -0500 Subject: [PATCH 01/13] Added protections to Campbell Clapp-Hornberger water transfer functions so that it will tolerate incredibly small water contents --- biogeophys/FatesHydroWTFMod.F90 | 43 ++++++++++++++++++++++++--------- 1 file changed, 31 insertions(+), 12 deletions(-) diff --git a/biogeophys/FatesHydroWTFMod.F90 b/biogeophys/FatesHydroWTFMod.F90 index 8356b2b165..ea5ac7aba0 100644 --- a/biogeophys/FatesHydroWTFMod.F90 +++ b/biogeophys/FatesHydroWTFMod.F90 @@ -44,6 +44,8 @@ module FatesHydroWTFMod ! elastic-caviation region + real(r8), parameter :: min_psi_cch = -9._r8 ! Minimum psi we are willing to track in cch + ! Generic class that can be extended to describe ! specific water retention functions @@ -146,6 +148,7 @@ module FatesHydroWTFMod ! Water Retention Function type, public, extends(wrf_type) :: wrf_type_cch real(r8) :: th_sat ! Saturation volumetric water content [m3/m3] + real(r8) :: th_res ! Residual volumetric water content [m3/m3] real(r8) :: psi_sat ! Bubbling pressure (potential at saturation) [Mpa] real(r8) :: beta ! Clapp-Hornberger "beta" parameter [-] contains @@ -174,6 +177,7 @@ module FatesHydroWTFMod ! Water Retention Function type, public, extends(wrf_type) :: wrf_type_smooth_cch real(r8) :: th_sat ! Saturation volumetric water content [m3/m3] + real(r8) :: th_res ! Residual volumetric water content [m3/m3] real(r8) :: psi_sat ! Bubbling pressure (potential at saturation) [Mpa] real(r8) :: beta ! Clapp-Hornberger "beta" parameter [-] real(r8) :: scch_pu ! An estimated breakpoint capillary pressure, below which the specified water retention curve is applied. It is also the lower limit when the smoothing function is applied. [Mpa] @@ -710,10 +714,11 @@ subroutine set_wrf_param_cch(this,params_in) this%th_max = max_sf_interp*this%th_sat this%psi_max = this%psi_from_th(this%th_max-tiny(this%th_max)) this%dpsidth_max = this%dpsidth_from_th(this%th_max-tiny(this%th_max)) - this%th_min = fates_unset_r8 - this%psi_min = fates_unset_r8 - this%dpsidth_min = fates_unset_r8 + this%psi_min = min_psi_cch + this%th_min = this%th_from_psi(min_psi_cch+tiny(this%th_max)) + this%dpsidth_min = this%dpsidth_from_th(this%th_min+tiny(this%th_max)) + return end subroutine set_wrf_param_cch @@ -751,8 +756,12 @@ function th_from_psi_cch(this,psi) result(th) if(psi>this%psi_max) then ! Linear range for extreme values th = this%th_max + (psi-this%psi_max)/this%dpsidth_max - else - th = this%th_sat*(psi/this%psi_sat)**(-1.0_r8/this%beta) + else + if(psithis%th_max) then psi = this%psi_max + this%dpsidth_max*(th-max_sf_interp*this%th_sat) - else - psi = this%psi_sat*(th/this%th_sat)**(-this%beta) + else + if(ththis%th_max) then - dpsidth = this%dpsidth_max + dpsidth = this%dpsidth_max else - dpsidth = -this%beta*this%psi_sat/this%th_sat * (th/this%th_sat)**(-this%beta-1._r8) + if(th Date: Thu, 22 Feb 2024 16:23:08 -0500 Subject: [PATCH 02/13] reverting the cap on the smoothed cch function, it already has a protection --- biogeophys/FatesHydroWTFMod.F90 | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/biogeophys/FatesHydroWTFMod.F90 b/biogeophys/FatesHydroWTFMod.F90 index ea5ac7aba0..e1f1a40df2 100644 --- a/biogeophys/FatesHydroWTFMod.F90 +++ b/biogeophys/FatesHydroWTFMod.F90 @@ -148,7 +148,6 @@ module FatesHydroWTFMod ! Water Retention Function type, public, extends(wrf_type) :: wrf_type_cch real(r8) :: th_sat ! Saturation volumetric water content [m3/m3] - real(r8) :: th_res ! Residual volumetric water content [m3/m3] real(r8) :: psi_sat ! Bubbling pressure (potential at saturation) [Mpa] real(r8) :: beta ! Clapp-Hornberger "beta" parameter [-] contains @@ -177,7 +176,6 @@ module FatesHydroWTFMod ! Water Retention Function type, public, extends(wrf_type) :: wrf_type_smooth_cch real(r8) :: th_sat ! Saturation volumetric water content [m3/m3] - real(r8) :: th_res ! Residual volumetric water content [m3/m3] real(r8) :: psi_sat ! Bubbling pressure (potential at saturation) [Mpa] real(r8) :: beta ! Clapp-Hornberger "beta" parameter [-] real(r8) :: scch_pu ! An estimated breakpoint capillary pressure, below which the specified water retention curve is applied. It is also the lower limit when the smoothing function is applied. [Mpa] @@ -932,11 +930,9 @@ subroutine set_wrf_param_smooth_cch(this,params_in) this%th_max = max_sf_interp*this%th_sat this%psi_max = this%psi_from_th(this%th_max-tiny(this%th_max)) this%dpsidth_max = this%dpsidth_from_th(this%th_max-tiny(this%th_max)) - - this%psi_min = min_psi_cch - this%th_min = this%th_from_psi(min_psi_cch+tiny(this%th_max)) - this%dpsidth_min = this%dpsidth_from_th(this%th_min+tiny(this%th_max)) - + this%th_min = fates_unset_r8 + this%psi_min = fates_unset_r8 + this%dpsidth_min = fates_unset_r8 return end subroutine set_wrf_param_smooth_cch @@ -1056,8 +1052,7 @@ function th_from_psi_smooth_cch(this,psi) result(th) sat = 1.d0 endif th = sat * this%th_sat - - + return end function th_from_psi_smooth_cch @@ -1214,7 +1209,7 @@ function dpsidth_from_th_smooth_cch(this,th) result(dpsidth) sat_res = 0._r8 alpha = -1._r8/this%psi_sat lambda = 1._r8/this%beta - + pc = 1._r8 * this%psi_from_th(th) if( pc <= this%scch_pu ) then ! Unsaturated full Brooks-Corey regime. From fce4287b45b738a6438780d65a77c3d718bbc4c4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 13 Mar 2024 14:20:02 -0400 Subject: [PATCH 03/13] Updated error reporting in 1d taylor solve for hydro --- biogeophys/FatesPlantHydraulicsMod.F90 | 57 +++++++++++++------------- 1 file changed, 29 insertions(+), 28 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index e38e042252..0fe4829f97 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -3993,28 +3993,39 @@ subroutine Report1DError(cohort, csite_hydr, ilayer, z_node, v_node, & write(fates_log(),*) 'psi_z: ',h_node(:)-psi_node(:) write(fates_log(),*) 'vol, theta, H, Psi, kmax-' write(fates_log(),*) 'flux: ', q_top_eff*dt_step - write(fates_log(),*) 'l:',v_node(1),th_node(1),h_node(1),psi_node(1) - write(fates_log(),*) ' ',cohort_hydr%kmax_stem_upper(1)*rootfr_scaler - write(fates_log(),*) 's:',v_node(2),th_node(2),h_node(2),psi_node(2) - write(fates_log(),*) ' ',1._r8/(1._r8/(cohort_hydr%kmax_stem_lower(1)*rootfr_scaler) + 1._r8/(cohort_hydr%kmax_troot_upper*rootfr_scaler)) - write(fates_log(),*) 't:',v_node(3),th_node(3),h_node(3) - write(fates_log(),*) ' ',1._r8/(1._r8/cohort_hydr%kmax_troot_lower(ilayer)+ 1._r8/cohort_hydr%kmax_aroot_upper(ilayer)) - write(fates_log(),*) 'a:',v_node(4),th_node(4),h_node(4) - write(fates_log(),*) ' in:',1._r8/(1._r8/cohort_hydr%kmax_aroot_radial_in(ilayer) + & + + do i = 1,n_hypool_leaf + k = i + write(fates_log(),*) 'leaf node ',k,v_node(k),th_node(k),h_node(k),psi_node(k) + end do + do i = 1,n_hypool_stem + k = i+n_hypool_leaf + write(fates_log(),*) 'stem node',k,v_node(k),th_node(k),h_node(k),psi_node(k) + write(fates_log(),*) ' ',cohort_hydr%kmax_stem_upper(k)*rootfr_scaler + end do + write(fates_log(),*) 'troot to stem kmax: ', & + 1._r8/(1._r8/(cohort_hydr%kmax_stem_lower(n_hypool_stem)*rootfr_scaler) + 1._r8/(cohort_hydr%kmax_troot_upper*rootfr_scaler)) + k = n_hypool_leaf + n_hypool_stem + 1 + write(fates_log(),*) 'troot node:',k,v_node(k),th_node(k),h_node(k) + write(fates_log(),*) 'aroot to troot kmax: ', & + 1._r8/(1._r8/cohort_hydr%kmax_troot_lower(ilayer)+ 1._r8/cohort_hydr%kmax_aroot_upper(ilayer)) + k = n_hypool_leaf + n_hypool_stem + 2 + write(fates_log(),*) 'aroot node:',k,v_node(k),th_node(k),h_node(k) + write(fates_log(),*) ' kmax soil-root in:',1._r8/(1._r8/cohort_hydr%kmax_aroot_radial_in(ilayer) + & 1._r8/(csite_hydr%kmax_upper_shell(ilayer,1)*aroot_frac_plant) + & 1._r8/cohort_hydr%kmax_aroot_upper(ilayer)) - write(fates_log(),*) ' out:',1._r8/(1._r8/cohort_hydr%kmax_aroot_radial_out(ilayer) + & + write(fates_log(),*) ' kmax soil-root out:',1._r8/(1._r8/cohort_hydr%kmax_aroot_radial_out(ilayer) + & 1._r8/(csite_hydr%kmax_upper_shell(ilayer,1)*aroot_frac_plant) + & 1._r8/cohort_hydr%kmax_aroot_upper(ilayer)) - write(fates_log(),*) 'r1:',v_node(5),th_node(5),h_node(5) - write(fates_log(),*) ' ',1._r8/(1._r8/(csite_hydr%kmax_lower_shell(ilayer,1)*aroot_frac_plant) + 1._r8/(csite_hydr%kmax_upper_shell(ilayer,2)*aroot_frac_plant)) - write(fates_log(),*) 'r2:',v_node(6),th_node(6),h_node(6) - write(fates_log(),*) ' ',1._r8/(1._r8/(csite_hydr%kmax_lower_shell(ilayer,2)*aroot_frac_plant) + 1._r8/(csite_hydr%kmax_upper_shell(ilayer,3)*aroot_frac_plant)) - write(fates_log(),*) 'r3:',v_node(7),th_node(7),h_node(7) - write(fates_log(),*) ' ',1._r8/(1._r8/(csite_hydr%kmax_lower_shell(ilayer,3)*aroot_frac_plant) + 1._r8/(csite_hydr%kmax_upper_shell(ilayer,4)*aroot_frac_plant)) - write(fates_log(),*) 'r4:',v_node(8),th_node(8),h_node(8) - write(fates_log(),*) ' ',1._r8/(1._r8/(csite_hydr%kmax_lower_shell(ilayer,4)*aroot_frac_plant) + 1._r8/(csite_hydr%kmax_upper_shell(ilayer,5)*aroot_frac_plant)) - write(fates_log(),*) 'r5:',v_node(9),th_node(9),h_node(9) + do i = 1,nshell + k = = n_hypool_leaf + n_hypool_stem + 2 + i + write(fates_log(),*) 'rhizo shell k:',k,v_node(k),th_node(k),h_node(k) + if(i Date: Sun, 17 Mar 2024 17:49:26 -0600 Subject: [PATCH 04/13] correct minor build errors --- biogeophys/FatesPlantHydraulicsMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 0fe4829f97..73c4f442e7 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -3939,7 +3939,7 @@ subroutine Report1DError(cohort, csite_hydr, ilayer, z_node, v_node, & logical, intent(in) :: recruitflag type(ed_cohort_hydr_type),pointer :: cohort_hydr - integer :: i + integer :: i, k integer :: ft real(r8) :: leaf_water real(r8) :: stem_water @@ -4018,7 +4018,7 @@ subroutine Report1DError(cohort, csite_hydr, ilayer, z_node, v_node, & 1._r8/(csite_hydr%kmax_upper_shell(ilayer,1)*aroot_frac_plant) + & 1._r8/cohort_hydr%kmax_aroot_upper(ilayer)) do i = 1,nshell - k = = n_hypool_leaf + n_hypool_stem + 2 + i + k = n_hypool_leaf + n_hypool_stem + 2 + i write(fates_log(),*) 'rhizo shell k:',k,v_node(k),th_node(k),h_node(k) if(i Date: Mon, 25 Mar 2024 13:05:26 -0400 Subject: [PATCH 05/13] Adding ftc protections (ie 0 flux) for perilously low wcs --- biogeophys/FatesHydroWTFMod.F90 | 12 +++++------ biogeophys/FatesPlantHydraulicsMod.F90 | 30 ++++++++++++++++++-------- 2 files changed, 27 insertions(+), 15 deletions(-) diff --git a/biogeophys/FatesHydroWTFMod.F90 b/biogeophys/FatesHydroWTFMod.F90 index e1f1a40df2..c67ab65a6e 100644 --- a/biogeophys/FatesHydroWTFMod.F90 +++ b/biogeophys/FatesHydroWTFMod.F90 @@ -28,7 +28,7 @@ module FatesHydroWTFMod __FILE__ - real(r8), parameter :: min_ftc = 0.00001_r8 ! Minimum allowed fraction of total conductance + real(r8), parameter :: min_ftc = 0.0_r8 ! Minimum allowed fraction of total conductance ! Bounds on saturated fraction, outside of which we use linear PV or stop flow ! In this context, the saturated fraction is defined by the volumetric WC "th" @@ -44,7 +44,7 @@ module FatesHydroWTFMod ! elastic-caviation region - real(r8), parameter :: min_psi_cch = -9._r8 ! Minimum psi we are willing to track in cch + real(r8), parameter :: min_theta_cch = 0.01_r8 ! Minimum theta (matches ctsm) ! Generic class that can be extended to describe ! specific water retention functions @@ -627,7 +627,7 @@ function ftc_from_psi_vg(this,psi) result(ftc) ! Make sure this is well behaved ftc = min(1._r8,max(min_ftc,num/den)) - + else ftc = 1._r8 @@ -713,9 +713,9 @@ subroutine set_wrf_param_cch(this,params_in) this%psi_max = this%psi_from_th(this%th_max-tiny(this%th_max)) this%dpsidth_max = this%dpsidth_from_th(this%th_max-tiny(this%th_max)) - this%psi_min = min_psi_cch - this%th_min = this%th_from_psi(min_psi_cch+tiny(this%th_max)) - this%dpsidth_min = this%dpsidth_from_th(this%th_min+tiny(this%th_max)) + this%th_min = min_theta_cch + this%psi_min = this%psi_from_th(min_theta_cch+tiny(this%th_max)) + this%dpsidth_min = this%dpsidth_from_th(min_theta_cch+tiny(this%th_max)) return end subroutine set_wrf_param_cch diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 73c4f442e7..f666b3409b 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -4111,18 +4111,30 @@ subroutine GetImTaylorKAB(kmax_up,kmax_dn, & end if - ! Calculate total effective conductance over path [kg s-1 MPa-1] - k_eff = 1._r8/(1._r8/(ftc_up*kmax_up)+1._r8/(ftc_dn*kmax_dn)) - - ! "A" term, which operates on the downstream node (closer to atm) - a_term = k_eff**2.0_r8 * h_diff * kmax_dn**(-1.0_r8) * ftc_dn**(-2.0_r8) & - * dftc_dtheta_dn - k_eff * dpsi_dtheta_dn + ! Prevent issue with zero flux + if( (ftc_up*kmax_up)>nearzero .and. (ftc_dn*kmax_dn)>nearzero ) then + + ! Calculate total effective conductance over path [kg s-1 MPa-1] + k_eff = 1._r8/(1._r8/(ftc_up*kmax_up)+1._r8/(ftc_dn*kmax_dn)) + + ! "A" term, which operates on the downstream node (closer to atm) + a_term = k_eff**2.0_r8 * h_diff * kmax_dn**(-1.0_r8) * ftc_dn**(-2.0_r8) & + * dftc_dtheta_dn - k_eff * dpsi_dtheta_dn + + + ! "B" term, which operates on the upstream node (further from atm) + b_term = k_eff**2.0_r8 * h_diff * kmax_up**(-1.0_r8) * ftc_up**(-2.0_r8) & + * dftc_dtheta_up + k_eff * dpsi_dtheta_up + else - ! "B" term, which operates on the upstream node (further from atm) - b_term = k_eff**2.0_r8 * h_diff * kmax_up**(-1.0_r8) * ftc_up**(-2.0_r8) & - * dftc_dtheta_up + k_eff * dpsi_dtheta_up + k_eff = 0._r8 + a_term = 0._r8 + b_term = 0._r8 + + end if + ! Restore ftc ftc_dn = ftc_dn_tmp ftc_up = ftc_up_tmp From 1abba2fdaf72f9c70b54a674f2d632a1358ee671 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 26 Mar 2024 09:56:24 -0400 Subject: [PATCH 06/13] hydro safemath, looking into pinning zero ftc to min psi --- biogeophys/FatesHydroWTFMod.F90 | 9 +- .../hydro/HydroUTestDriver.py | 187 +++++++++++++++--- main/FatesHydraulicsMemMod.F90 | 24 +-- 3 files changed, 179 insertions(+), 41 deletions(-) diff --git a/biogeophys/FatesHydroWTFMod.F90 b/biogeophys/FatesHydroWTFMod.F90 index c67ab65a6e..9c0ef60a89 100644 --- a/biogeophys/FatesHydroWTFMod.F90 +++ b/biogeophys/FatesHydroWTFMod.F90 @@ -903,7 +903,8 @@ subroutine set_wrf_param_smooth_cch(this,params_in) this%scch_b2 = 0.d0 this%scch_b3 = (2.d0 - bcAtPu*(2.d0+lambdaDeltaPuOnPu)) * oneOnDeltaPu * oneOnDeltaPu * oneOnDeltaPu if( this%scch_b3 <= 0.d0 ) then - write(fates_log(),*) 'set_wrf_param_smooth_cch b3 <=0',pu,ps,alpha,lambda,oneOnDeltaPu,lambdaDeltaPuOnPu,bcAtPu,this%psi_sat + write(fates_log(),*) 'set_wrf_param_smooth_cch b3 <=0',pu,ps,alpha,lambda,oneOnDeltaPu, & + lambdaDeltaPuOnPu,bcAtPu,this%psi_sat call endrun(msg=errMsg(sourcefile, __LINE__)) endif else @@ -980,9 +981,11 @@ subroutine set_wkf_param_smooth_cch(this,params_in) ! Store coefficients for cubic function. this%scch_b2 = 0.d0 - this%scch_b3 = (2.d0 - bcAtPu*(2.d0+lambdaDeltaPuOnPu)) * oneOnDeltaPu * oneOnDeltaPu * oneOnDeltaPu + this%scch_b3 = (2.d0 - bcAtPu*(2.d0+lambdaDeltaPuOnPu)) * & + oneOnDeltaPu * oneOnDeltaPu * oneOnDeltaPu if( this%scch_b3 <= 0.d0 ) then - write(fates_log(),*) 'set_wrf_param_smooth_cch b3 <=0',pu,ps,alpha,lambda,oneOnDeltaPu,lambdaDeltaPuOnPu,bcAtPu,this%psi_sat + write(fates_log(),*) 'set_wrf_param_smooth_cch b3 <=0', & + pu,ps,alpha,lambda,oneOnDeltaPu,lambdaDeltaPuOnPu,bcAtPu,this%psi_sat call endrun(msg=errMsg(sourcefile, __LINE__)) endif else diff --git a/functional_unit_testing/hydro/HydroUTestDriver.py b/functional_unit_testing/hydro/HydroUTestDriver.py index 5c2b026361..933c2aa837 100644 --- a/functional_unit_testing/hydro/HydroUTestDriver.py +++ b/functional_unit_testing/hydro/HydroUTestDriver.py @@ -150,18 +150,60 @@ def __init__(self,index,p50,avuln): iret = setwkf(ci(index),ci(tfs_type),ci(len(init_wkf_args)),c8_arr(init_wkf_args)) +def OMParams(zsoi): + + zsapric= 0.5 + om_watsat = max(0.93 - 0.1 *(zsoi/zsapric), 0.83) + om_bsw = min(2.7 + 9.3 *(zsoi/zsapric), 12.0) + om_sucsat = min(10.3 - 0.2 *(zsoi/zsapric), 10.1) + + #om_sucsat = SuctionMMtoMPa(om_sucsat) + + return(om_watsat,om_sucsat,om_bsw) + +def CCHParmsCosby84T5(zsoi,om_frac,sand_frac,clay_frac): + + # cosby_1984_table5 + + # Get pedotransfer for soil matrix + watsat = 0.489 - 0.00126*sand_frac + bsw = 2.91 + 0.159*clay_frac + sucsat = 10. * ( 10.**(1.88-0.0131*sand_frac) ) + + [om_watsat,om_sucsat,om_bsw] = OMParams(zsoi) + + # Update pedotransfer to include organic material + watsat = (1.0 - om_frac) * watsat + om_watsat * om_frac + bsw = (1.0 - om_frac) * bsw + om_bsw * om_frac + sucsat = (1.0 - om_frac) * sucsat + om_sucsat * om_frac + + # Convert from mm to MPa + sucsat = SuctionMMtoMPa(sucsat) + + + # psi = psi_sat*(th/th_sat)**(-beta) + # th = th_sat*(psi / psi_sat)^(-1/beta) + + return(watsat,sucsat,bsw) + + +def SuctionMMtoMPa(suction_mm): + + denh2o = 1.0e3 # kg/m3 + grav_earth = 9.8 # m/s2 + mpa_per_pa = 1.0e-6 # MPa per Pa + m_per_mm = 1.0e-3 # meters per millimeter + + suction_mpa = (-1.0)*suction_mm*denh2o*grav_earth*mpa_per_pa*m_per_mm + + return(suction_mpa) + def main(argv): - # First check to make sure python 2.7 is being used + version = platform.python_version() verlist = version.split('.') - if( not ((verlist[0] == '2') & (verlist[1] == '7') & (int(verlist[2])>=15) ) ): - print("The PARTEH driver mus be run with python 2.7") - print(" with tertiary version >=15.") - print(" your version is {}".format(version)) - print(" exiting...") - sys.exit(2) # Read in the arguments # ======================================================================================= @@ -178,21 +220,22 @@ def main(argv): # min_theta = np.full(shape=(2),dtype=np.float64,fill_value=np.nan) - -# wrf_type = [vg_type, vg_type, cch_type, cch_type] -# wkf_type = [vg_type, tfs_type, cch_type, tfs_type] - -# th_ress = [0.01, 0.10, -9, -9] -# th_sats = [0.55, 0.55, 0.65, 0.65] -# alphas = [1.0, 1.0, 1.0, 1.0] -# psds = [2.7, 2.7, 2.7, 2.7] -# tort = [0.5, 0.5, 0.5, 0.5] -# beta = [-9, -9, 6, 9] -# avuln = [2.0, 2.0, 2.5, 2.5] -# p50 = [-1.5, -1.5, -2.25, -2.25] - - ncomp= 4 - + # wrf_type = [vg_type, vg_type, cch_type, cch_type] + # wkf_type = [vg_type, tfs_type, cch_type, tfs_type] + + # th_ress = [0.01, 0.10, -9, -9] + # th_sats = [0.55, 0.55, 0.65, 0.65] + # alphas = [1.0, 1.0, 1.0, 1.0] + # psds = [2.7, 2.7, 2.7, 2.7] + # tort = [0.5, 0.5, 0.5, 0.5] + # beta = [-9, -9, 6, 9] + # avuln = [2.0, 2.0, 2.5, 2.5] + # p50 = [-1.5, -1.5, -2.25, -2.25] + + + ncomp = 4 + ncomp_tot = 20 + rwc_fd = [1.0,0.958,0.958,0.958] rwccap = [1.0,0.947,0.947,0.947] cap_slp = [] @@ -213,7 +256,7 @@ def main(argv): # Allocate memory to our objective classes - iret = initalloc_wtfs(ci(ncomp),ci(ncomp)) + iret = initalloc_wtfs(ci(ncomp_tot),ci(ncomp_tot)) print('Allocated') @@ -221,8 +264,8 @@ def main(argv): # vg_wrf(1,alpha=1.0,psd=2.7,th_sat=0.55,th_res=0.1) # vg_wkf(1,alpha=1.0,psd=2.7,th_sat=0.55,th_res=0.1,tort=0.5) - cch_wrf(1,th_sat=0.55, psi_sat=-1.56e-3, beta=6) - cch_wkf(1,th_sat=0.55, psi_sat=-1.56e-3, beta=6) + cch_wrf(1,th_sat=0.55, psi_sat=-1.56e-3, beta=3) + cch_wkf(1,th_sat=0.55, psi_sat=-1.56e-3, beta=3) # cch_wrf(3,th_sat=0.55, psi_sat=-1.56e-3, beta=6) # tfs_wkf(3,p50=-2.25, avuln=2.0) @@ -361,11 +404,103 @@ def main(argv): ax1.set_ylim([0,2]) # ax1.set_ylim([0,10]) ax1.legend(loc='upper right') - plt.show() + ndepth = 1000 + zdepth = np.linspace(0.01,10.0,num = ndepth) + + om_watsat_z = np.zeros(shape=(ndepth,1)) + om_sucsat_z = np.zeros(shape=(ndepth,1)) + om_bsw_z = np.zeros(shape=(ndepth,1)) + + for icl in range(ndepth): + [om_watsat_z[icl],om_sucsat_z[icl],om_bsw_z[icl]] = OMParams(zdepth[icl]) + + fig4,((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2,figsize=(9,6)) + + ax1.plot(om_watsat_z,-zdepth) + ax2.plot(om_sucsat_z,-zdepth) + ax3.plot(om_bsw_z,-zdepth) + ax4.axis('off') + + + ax1.set_ylabel('soil depth') + ax1.set_xlabel('Saturated WC [m3/m3]') + ax2.set_xlabel('Saturated Suction [mm]') + ax3.set_ylabel('soil depth') + ax3.set_xlabel('beta ') + # ax1.set_xlim([-10,3]) + # ax1.set_ylim([0,2]) + # ax1.set_ylim([0,10]) + # ax1.legend(loc='upper right') + plt.tight_layout() + + + # Soil texture distributions + + #clay_frac = np.zeros(shape=(100,1)) + #sand_frac = np.zeros(shape=(100,1)) + #om_frac = np.zeros(shape=(100,1)) + #clay_frac = np.linspace(0.0, 0.7, num=npts) + + watsat_v = [] + sucsat_v = [] + bsw_v = [] + + ntex = 100 + for icl in range(ntex): + clay_frac = float(ic)/float(ntex) + for isa in range(ntex): + sand_frac = float(isa)/float(ntex) + + if( (clay_frac+sand_frac)<1.0 ): + for om_frac in [0.0,1.0]: + for zsoi in [0.01,10.0]: + [watsat,sucsat,bsw] = CCHParmsCosby84T5(zsoi,om_frac,sand_frac,clay_frac) + watsat_v.append(watsat) + sucsat_v.append(sucsat) + bsw_v.append(bsw) + + + fig5,((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2,figsize=(9,6)) + + ax1.hist(watsat_v,bins=50) + ax2.hist(sucsat_v,bins=50) + ax3.hist(bsw_v,bins=50) + ax4.axis('off') + + + ax1.set_xlabel('Sat. WC [m3/m3]') + ax2.set_xlabel('Sat. Sucation [MPa]') + ax3.set_xlabel('Beta') + + plt.tight_layout() + + ind = [] + ind.append(np.argmin(bsw_v)) + ind.append(np.argmax(bsw_v)) + ind.append(np.argmin(sucsat_v)) + ind.append(np.argmax(sucsat_v)) + ind.append(np.argmin(watsat_v)) + ind.append(np.argmax(watsat_v)) + + # Now lets test the lowest possible water contents using these + + ii=4 + + for i in ind: + ii = ii+1 + cch_wrf(ii,th_sat=watsat_v[i], psi_sat=sucsat_v[i], beta=bsw_v[i]) + cch_wkf(ii,th_sat=watsat_v[i], psi_sat=sucsat_v[i], beta=bsw_v[i]) + print('---th sat: {}, psi sat: {}, beat: {}'.format(watsat_v[i],sucsat_v[i],bsw_v[i])) + for psi in -np.linspace(1,24,24): + th = th_from_psi(ci(ii),c8(psi)) + print('{}, psi: {}, th: {}'.format(ii,psi,th)) + + plt.show() + # code.interact(local=dict(globals(), **locals())) # Helper code to plot negative logs diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90 index 61e97173c7..f31ffe7e25 100644 --- a/main/FatesHydraulicsMemMod.F90 +++ b/main/FatesHydraulicsMemMod.F90 @@ -168,19 +168,19 @@ module FatesHydraulicsMemMod real(r8), allocatable :: residual(:) real(r8), allocatable :: ajac(:,:) ! Jacobian (N terms, N equations) - real(r8), allocatable :: th_node_init(:) + real(r8), allocatable :: th_node_init(:) real(r8), allocatable :: th_node_prev(:) - real(r8), allocatable :: th_node(:) - real(r8), allocatable :: dth_node(:) - real(r8), allocatable :: h_node(:) - real(r8), allocatable :: v_node(:) - real(r8), allocatable :: z_node(:) - real(r8), allocatable :: psi_node(:) - real(r8), allocatable :: q_flux(:) - real(r8), allocatable :: dftc_dpsi_node(:) - real(r8), allocatable :: ftc_node(:) - real(r8), allocatable :: kmax_up(:) - real(r8), allocatable :: kmax_dn(:) + real(r8), allocatable :: th_node(:) ! Relative water content (theta) of node [m3/m3] + real(r8), allocatable :: dth_node(:) ! Change (time derivative) in water content of node + real(r8), allocatable :: h_node(:) ! + real(r8), allocatable :: v_node(:) ! Volume of the node [m3] + real(r8), allocatable :: z_node(:) ! Eleveation potential of the node (datum 0 is surface) + real(r8), allocatable :: psi_node(:) ! Suction of the node [MPa] + real(r8), allocatable :: q_flux(:) ! Mass flux of pathways between nodes [] + real(r8), allocatable :: dftc_dpsi_node(:) ! Differential of fraction total conductivity with suction + real(r8), allocatable :: ftc_node(:) ! fraction of total conductivity [-] + real(r8), allocatable :: kmax_up(:) ! Maximum conductivity for upstream side of compartment + real(r8), allocatable :: kmax_dn(:) ! Maximum conductivity for downstream side of compartment ! Scratch arrays real(r8) :: cohort_recruit_water_layer(nlevsoi_hyd_max) ! the recruit water requirement for a cohort From 35b396dacd30b87928f5dec9c232e2749f0860cd Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 29 Mar 2024 16:23:56 -0400 Subject: [PATCH 07/13] added hydro code to force zero conductance a psi_min --- biogeophys/FatesHydroWTFMod.F90 | 301 ++++++++++++++---- biogeophys/FatesPlantHydraulicsMod.F90 | 34 ++ .../hydro/HydroUTestDriver.py | 45 ++- .../hydro/f90_src/HydroUnitWrapMod.F90 | 15 +- 4 files changed, 317 insertions(+), 78 deletions(-) diff --git a/biogeophys/FatesHydroWTFMod.F90 b/biogeophys/FatesHydroWTFMod.F90 index 9c0ef60a89..a7c47efb7a 100644 --- a/biogeophys/FatesHydroWTFMod.F90 +++ b/biogeophys/FatesHydroWTFMod.F90 @@ -30,6 +30,12 @@ module FatesHydroWTFMod real(r8), parameter :: min_ftc = 0.0_r8 ! Minimum allowed fraction of total conductance + real(r8), parameter :: min_ftc_scalar=2.0_r8 ! This scalar is used to define the attenuation + ! of the weighting that we use for imposing + ! ftc min at the minimum allowable psi + ! A value of two yielded a wieghting factor of + ! about 0.175 after 1 MPa, and 0.025 after 2 MPA + ! Bounds on saturated fraction, outside of which we use linear PV or stop flow ! In this context, the saturated fraction is defined by the volumetric WC "th" ! and the volumetric residual and saturation "th_res" and "th_sat": (th-th_r)/(th_sat-th_res) @@ -44,7 +50,8 @@ module FatesHydroWTFMod ! elastic-caviation region - real(r8), parameter :: min_theta_cch = 0.01_r8 ! Minimum theta (matches ctsm) + !real(r8), parameter :: min_theta_cch = 0.05_r8 ! Minimum theta (matches ctsm) + real(r8), parameter :: min_psi_cch = -20.0_r8 ! Minimum theta (matches ctsm) ! Generic class that can be extended to describe ! specific water retention functions @@ -91,6 +98,7 @@ module FatesHydroWTFMod ! water conductance functions type, public :: wkf_type + type(wrf_type), pointer :: wrf ! Pointer to the matching water retention function contains procedure :: ftc_from_psi => ftc_from_psi_base procedure :: dftcdpsi_from_psi => dftcdpsi_from_psi_base @@ -242,7 +250,8 @@ module FatesHydroWTFMod procedure :: set_wkf_param => set_wkf_param_tfs end type wkf_type_tfs - + public :: get_min_ftc_weight + contains ! ===================================================================================== @@ -343,7 +352,34 @@ function th_linear_res(this,psi) result(th) end function th_linear_res ! =========================================================================== + + subroutine get_min_ftc_weight(psi_min,psi,min_ftc_weight,dmin_ftc_weight_dpsi) + + ! This routine determines the weighting factor used + ! to generate a smooth curve to impose that min_ftc happens at min_psi + ! This method is to be used with any and all water transfer functions + + real(r8),intent(in) :: psi ! MPa suction + real(r8),intent(in) :: psi_min ! Minimum value of psi we track, value that pins to ftc_min + real(r8),intent(out):: min_ftc_weight ! weighting factor for min_ftc + real(r8),intent(out):: dmin_ftc_weight_dpsi ! derivative of weighting factor wrt psi + min_ftc_weight = exp(min_ftc_scalar*(psi_min-psi)) + dmin_ftc_weight_dpsi = -min_ftc_scalar*exp(min_ftc_scalar*(psi_min-psi)) + + if(min_ftc_weight>=1.)then + min_ftc_weight = 1._r8 + dmin_ftc_weight_dpsi = 0._r8 + elseif(min_ftc_weight<=0._r8)then + min_ftc_weight = 0._r8 + dmin_ftc_weight_dpsi = 0._r8 + end if + + return + end subroutine get_min_ftc_weight + + ! ============================================================================= + subroutine set_wrf_param_base(this,params_in) class(wrf_type) :: this real(r8),intent(in) :: params_in(:) @@ -613,7 +649,11 @@ function ftc_from_psi_vg(this,psi) result(ftc) real(r8) :: psi_eff real(r8) :: m ! pore size distribution param () real(r8) :: n ! pore size distribution param (psd) + real(r8) :: min_ftc_weight + real(r8) :: dmin_ftc_weight_dpsi + real(r8), pointer :: psi_min + psi_min => this%wrf%psi_min n = this%n_vg m = this%m_vg @@ -626,8 +666,18 @@ function ftc_from_psi_vg(this,psi) result(ftc) den = (1._r8 + (this%alpha*psi_eff)**n)**(this%tort*(m)) ! Make sure this is well behaved - ftc = min(1._r8,max(min_ftc,num/den)) - + ftc = min(1._r8,num/den) + + if(ftc<=min_ftc) then + ftc = min_ftc + else + ! Add protections and ensure no conductance at incredibly + ! low suction + call get_min_ftc_weight(psi_min,psi,min_ftc_weight,dmin_ftc_weight_dpsi) + + ftc = ftc*(1._r8 - min_ftc_weight) + min_ftc*min_ftc_weight + end if + else ftc = 1._r8 @@ -635,6 +685,9 @@ function ftc_from_psi_vg(this,psi) result(ftc) end function ftc_from_psi_vg + + + ! ==================================================================================== function dftcdpsi_from_psi_vg(this,psi) result(dftcdpsi) @@ -656,7 +709,11 @@ function dftcdpsi_from_psi_vg(this,psi) result(dftcdpsi) real(r8) :: dftcdpsi ! change in frac total cond wrt psi real(r8) :: m ! pore size distribution param (1/psd) real(r8) :: n + real(r8) :: min_ftc_weight + real(r8) :: dmin_ftc_weight_dpsi + real(r8), pointer :: psi_min + psi_min => this%wrf%psi_min n =this%n_vg m =this%m_vg @@ -666,8 +723,8 @@ function dftcdpsi_from_psi_vg(this,psi) result(dftcdpsi) psi_eff = -psi ! switch VG 1980 convention ftc = this%ftc_from_psi(psi) - - if(ftc<=min_ftc) then + + if ( abs(ftc-min_ftc) this%wrf%psi_min + ! th = th_sat * (psi/psi_sat)^(-1/b) ! ftc = (th/th_sat)^(2*b+3) @@ -826,6 +900,15 @@ function ftc_from_psi_cch(this,psi) result(ftc) ftc = (psi_eff/this%psi_sat)**(-2._r8-3._r8/this%beta) + if(ftc<=min_ftc) then + ftc = min_ftc + else + call get_min_ftc_weight(psi_min,psi,min_ftc_weight,dmin_ftc_weight_dpsi) + + ftc = ftc*(1._r8 - min_ftc_weight) + min_ftc*min_ftc_weight + + end if + end function ftc_from_psi_cch ! ==================================================================================== @@ -835,15 +918,37 @@ function dftcdpsi_from_psi_cch(this,psi) result(dftcdpsi) class(wkf_type_cch) :: this real(r8),intent(in) :: psi real(r8) :: dftcdpsi ! change in frac total cond wrt psi + real(r8) :: min_ftc_weight,dmin_ftc_weight_dpsi + real(r8) :: ftc + real(r8), pointer :: psi_min - ! Differentiate: - ! ftc = (psi/this%psi_sat)**(-2._r8-3._r8/this%beta) - + psi_min => this%wrf%psi_min + ! Note that if we assume a constant, capped FTC=1.0 ! at saturation, then the derivative is zero there if(psi this%wrf%psi_min pc = psi sat_res = 0._r8 alpha = -1._r8/this%psi_sat @@ -1286,8 +1392,18 @@ function ftc_from_psi_smooth_cch(this,psi) result(ftc) ! Here, `pc >= ps`. kr = 1.d0 endif - ftc = max(kr, min_ftc) + !ftc = max(kr, min_ftc) + + if(kr<=min_ftc) then + ftc = min_ftc + else + call get_min_ftc_weight(psi_min,psi,min_ftc_weight,dmin_ftc_weight_dpsi) + + ftc = kr*(1._r8 - min_ftc_weight) + min_ftc*min_ftc_weight + + end if + end function ftc_from_psi_smooth_cch @@ -1298,11 +1414,10 @@ function dftcdpsi_from_psi_smooth_cch(this,psi) result(dftcdpsi) class(wkf_type_smooth_cch) :: this real(r8),intent(in) :: psi real(r8) :: dftcdpsi ! change in frac total cond wrt psi - + real(r8) :: ftc real(r8) :: pc real(r8) :: kr real(r8) :: dkr_dP - ! real(r8) :: sat_res real(r8) :: alpha real(r8) :: lambda @@ -1310,46 +1425,65 @@ function dftcdpsi_from_psi_smooth_cch(this,psi) result(dftcdpsi) real(r8) :: deltaPc real(r8) :: dSe_dpc real(r8) :: dkr_dSe + real(r8) :: min_ftc_weight,dmin_ftc_weight_dpsi + real(r8), pointer :: psi_min - pc = psi - sat_res = 0._r8 - alpha = -1._r8/this%psi_sat - lambda = 1._r8/this%beta + psi_min => this%wrf%psi_min + ftc = this%ftc_from_psi(psi) + + if( abs(ftc-min_ftc)= ps`. - kr = 1.d0 - dkr_dP = 0.d0 - endif - dftcdpsi = dkr_dP - if(kr<=min_ftc) then - dftcdpsi = 0._r8 - endif + dkr_dSe = (3.d0 + 2.d0/lambda)*kr/Se + dkr_dp = dkr_dSe*dSe_dpc + elseif( pc < this%scch_ps ) then + ! Cubic smoothing regime. + ! Here, `pu < pc < ps <= 0`. + deltaPc = pc - this%scch_ps + Se = 1.d0 + deltaPc*deltaPc*(this%scch_b2 + deltaPc*this%scch_b3) + + dSe_dpc = deltaPc*(2*this%scch_b2 + 3*deltaPc*this%scch_b3) + + kr = Se ** (2.5d0 + 2.d0/lambda) + + dkr_dSe = (2.5d0 + 2.d0/lambda)*kr/Se + dkr_dp = dkr_dSe*dSe_dpc + else + ! Saturated regime. + ! Here, `pc >= ps`. + kr = 1.d0 + dkr_dP = 0.d0 + endif + dftcdpsi = dkr_dP + call get_min_ftc_weight(psi_min,psi,min_ftc_weight,dmin_ftc_weight_dpsi) + + ! differentiate: + ! ftc = ftc*(1._r8 - min_ftc_weight) + min_ftc*min_ftc_weight + ! ftc = ftc - ftc*min_ftc_weight + min_ftc*min_ftc_weight + + dftcdpsi = dftcdpsi - & + (dftcdpsi*min_ftc_weight + ftc*dmin_ftc_weight_dpsi) + & + min_ftc*dmin_ftc_weight_dpsi + + end if end function dftcdpsi_from_psi_smooth_cch @@ -1747,11 +1881,26 @@ function ftc_from_psi_tfs(this,psi) result(ftc) real(r8),intent(in) :: psi ! real(r8) :: ftc real(r8) :: psi_eff + real(r8) :: min_ftc_weight,dmin_ftc_weight_dpsi + real(r8), pointer :: psi_min + psi_min => this%wrf%psi_min psi_eff = min(-nearzero,psi) - ftc = max(min_ftc,1._r8/(1._r8 + (psi_eff/this%p50)**this%avuln)) + ftc = 1._r8/(1._r8 + (psi_eff/this%p50)**this%avuln) + + if(ftc<=min_ftc) then + ftc = min_ftc + else + ! Add protections and ensure no conductance at incredibly + ! low suction + call get_min_ftc_weight(psi_min,psi_eff,min_ftc_weight,dmin_ftc_weight_dpsi) + + ftc = ftc*(1._r8 - min_ftc_weight) + min_ftc*min_ftc_weight + end if + + end function ftc_from_psi_tfs ! ==================================================================================== @@ -1764,20 +1913,38 @@ function dftcdpsi_from_psi_tfs(this,psi) result(dftcdpsi) real(r8) :: fx ! portion of ftc function real(r8) :: dfx ! differentiation of portion of func real(r8) :: dftcdpsi ! change in frac total cond wrt psi + real(r8) :: min_ftc_weight,dmin_ftc_weight_dpsi + real(r8), pointer :: psi_min + psi_min => this%wrf%psi_min + ! Differentiate ! ftc = 1._r8/(1._r8 + (psi/this%p50(ft))**this%avuln(ft)) if(psi>0._r8)then dftcdpsi = 0._r8 else - ftc = 1._r8/(1._r8 + (psi/this%p50)**this%avuln) - if(ftc sites(s)%si_hydr%wrf_soil(j)%p + end do + + ! Update static quantities related to the rhizosphere call UpdateSizeDepRhizVolLenCon(sites(s), bc_in(s)) @@ -1655,6 +1664,15 @@ subroutine HydrSiteColdStart(sites, bc_in ) write(fates_log(),*) 'TFS conductance not used in soil' call endrun(msg=errMsg(sourcefile, __LINE__)) end select + + ! The fraction of total conductance functions need to know psi_min + ! to handle very very low conductances, therefore we need to construct + ! a pointer in conductance structure to the water retention structure + + do j=1,sites(s)%si_hydr%nlevrhiz + sites(s)%si_hydr%wkf_soil(j)%p%wrf => sites(s)%si_hydr%wrf_soil(j)%p + end do + end do @@ -6315,6 +6333,14 @@ subroutine InitHydroGlobals() write(fates_log(),*) 'undefined water conductance type for plants, pm:',pm,'type: ',hydr_htftype_node(pm) call endrun(msg=errMsg(sourcefile, __LINE__)) end select + + ! The fraction of total conductance functions need to know psi_min + ! to handle very very low conductances, therefore we need to construct + ! a pointer in conductance structure to the water retention structure + do ft = 1,numpft + wkf_plant(pm,ft)%p%wrf => wrf_plant(pm,ft)%p + end do + end do ! There is only 1 stomata conductance hypothesis which uses the p50 and @@ -6328,6 +6354,14 @@ subroutine InitHydroGlobals() EDPftvarcon_inst%hydr_avuln_gs(ft)]) end do + ! The fraction of total conductance functions need to know psi_min + ! to handle very very low conductances, therefore we need to construct + ! a pointer in conductance structure to the water retention structure + do ft = 1,numpft + wkf_plant(stomata_p_media,ft)%p%wrf => wrf_plant(stomata_p_media,ft)%p + end do + + return end subroutine InitHydroGlobals diff --git a/functional_unit_testing/hydro/HydroUTestDriver.py b/functional_unit_testing/hydro/HydroUTestDriver.py index 933c2aa837..98c38a63c2 100644 --- a/functional_unit_testing/hydro/HydroUTestDriver.py +++ b/functional_unit_testing/hydro/HydroUTestDriver.py @@ -58,10 +58,10 @@ ftc_from_psi.restype = c_double dftcdpsi_from_psi = f90_hydrounitwrap_obj.__hydrounitwrapmod_MOD_wrapdftcdpsi dftcdpsi_from_psi.restype = c_double - +ftcminwt_from_psi = f90_hydrounitwrap_obj.__hydrounitwrapmod_MOD_wrapftcminweightfrompsi +ftcminwt_from_psi.restype = c_double # Some constants -rwcft = [1.0,0.958,0.958,0.958] rwccap = [1.0,0.947,0.947,0.947] pm_leaf = 1 pm_stem = 2 @@ -260,6 +260,13 @@ def main(argv): print('Allocated') + #min_psi = -10. + #min_psi_falloff = 1 + #test_psi = np.linspace(min_psi, 0, num=1000) + #weight = y = e^(2*(-10-x)) + + + # Define the funcions and their parameters # vg_wrf(1,alpha=1.0,psd=2.7,th_sat=0.55,th_res=0.1) # vg_wkf(1,alpha=1.0,psd=2.7,th_sat=0.55,th_res=0.1,tort=0.5) @@ -447,9 +454,9 @@ def main(argv): sucsat_v = [] bsw_v = [] - ntex = 100 + ntex = 5 for icl in range(ntex): - clay_frac = float(ic)/float(ntex) + clay_frac = float(icl)/float(ntex) for isa in range(ntex): sand_frac = float(isa)/float(ntex) @@ -486,19 +493,37 @@ def main(argv): # Now lets test the lowest possible water contents using these + #psi_v = -np.linspace(0.01,24,24) + psi_v = -np.logspace(-3,1.5,60) + fig10,(ax1,ax2,ax3) = plt.subplots(3,1,figsize=(6,9)) + + ii=4 - + th_v = [] + ftc_v = [] + ftc_min_wt_v = [] for i in ind: ii = ii+1 cch_wrf(ii,th_sat=watsat_v[i], psi_sat=sucsat_v[i], beta=bsw_v[i]) cch_wkf(ii,th_sat=watsat_v[i], psi_sat=sucsat_v[i], beta=bsw_v[i]) print('---th sat: {}, psi sat: {}, beat: {}'.format(watsat_v[i],sucsat_v[i],bsw_v[i])) - for psi in -np.linspace(1,24,24): - th = th_from_psi(ci(ii),c8(psi)) - print('{}, psi: {}, th: {}'.format(ii,psi,th)) + for psi in psi_v: + th_v.append(th_from_psi(ci(ii),c8(psi))) + ftc_v.append(ftc_from_psi(ci(ii),c8(psi))) + ftc_min_wt_v.append(ftcminwt_from_psi(ci(ii),c8(psi))) + + ax1.plot(psi_v,th_v) + ax2.plot(psi_v,ftc_v) + ax3.plot(psi_v,ftc_min_wt_v) + th_v = [] + ftc_v = [] + ftc_min_wt_v = [] + + ax2.set_xlabel('Suction MPa') + ax1.set_ylabel('Theta') + ax2.set_ylabel('FTC') + ax3.set_ylabel('FTC Min Weight') - - plt.show() # code.interact(local=dict(globals(), **locals())) diff --git a/functional_unit_testing/hydro/f90_src/HydroUnitWrapMod.F90 b/functional_unit_testing/hydro/f90_src/HydroUnitWrapMod.F90 index 03e95a6a32..02167209ef 100644 --- a/functional_unit_testing/hydro/f90_src/HydroUnitWrapMod.F90 +++ b/functional_unit_testing/hydro/f90_src/HydroUnitWrapMod.F90 @@ -14,7 +14,8 @@ module HydroUnitWrapMod use FatesHydroWTFMod, only : wrf_type,wrf_type_vg,wrf_type_cch use FatesHydroWTFMod, only : wkf_type,wkf_type_vg,wkf_type_cch,wkf_type_tfs use FatesHydroWTFMod, only : wrf_arr_type,wkf_arr_type,wrf_type_tfs - + use FatesHydroWTFMod, only : get_min_ftc_weight + implicit none public save @@ -107,6 +108,8 @@ subroutine SetWKF(index,itype,npvals,pvals) stop end if + wkfs(index)%p%wrf => wrfs(index)%p + return end subroutine SetWKF @@ -167,5 +170,15 @@ function WrapFTCFromPSI(index,psi) result(ftc) return end function WrapFTCFromPSI + function WrapFTCMinWeightFromPsi(index,psi) result(min_ftc_weight) + integer, intent(in) :: index + real(r8),intent(in) :: psi + real(r8) :: min_ftc_weight,dmin_ftc_weight_dpsi + + print*,"PSI_MIN: ",wkfs(index)%p%wrf%psi_min + + call get_min_ftc_weight(wkfs(index)%p%wrf%psi_min,psi,min_ftc_weight,dmin_ftc_weight_dpsi) + + end function WrapFTCMinWeightFromPsi end module HydroUnitWrapMod From 57b10d32409aa3ff9f3789eabc717c3c89de26ea Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 2 Apr 2024 10:51:04 -0400 Subject: [PATCH 08/13] Getting zero flux at psi_min in fates hydro --- biogeophys/FatesHydroWTFMod.F90 | 25 ++++++++++++++----------- biogeophys/FatesPlantHydraulicsMod.F90 | 17 +++++++++++------ main/FatesInterfaceMod.F90 | 3 ++- 3 files changed, 27 insertions(+), 18 deletions(-) diff --git a/biogeophys/FatesHydroWTFMod.F90 b/biogeophys/FatesHydroWTFMod.F90 index a7c47efb7a..54833f979b 100644 --- a/biogeophys/FatesHydroWTFMod.F90 +++ b/biogeophys/FatesHydroWTFMod.F90 @@ -364,9 +364,9 @@ subroutine get_min_ftc_weight(psi_min,psi,min_ftc_weight,dmin_ftc_weight_dpsi) real(r8),intent(out):: min_ftc_weight ! weighting factor for min_ftc real(r8),intent(out):: dmin_ftc_weight_dpsi ! derivative of weighting factor wrt psi - min_ftc_weight = exp(min_ftc_scalar*(psi_min-psi)) - dmin_ftc_weight_dpsi = -min_ftc_scalar*exp(min_ftc_scalar*(psi_min-psi)) - + min_ftc_weight = exp(min_ftc_scalar*max(psi_min-psi,-10._r8)) + dmin_ftc_weight_dpsi = -min_ftc_scalar*exp(min_ftc_scalar*max(psi_min-psi,-10._r8)) + if(min_ftc_weight>=1.)then min_ftc_weight = 1._r8 dmin_ftc_weight_dpsi = 0._r8 @@ -785,7 +785,7 @@ subroutine set_wrf_param_cch(this,params_in) !this%psi_min = this%psi_from_th(min_theta_cch+tiny(this%th_max)) this%psi_min = min_psi_cch this%th_min = this%th_from_psi(min_psi_cch+tiny(this%th_max)) - this%dpsidth_min = this%dpsidth_from_th(this%th_min+tiny(this%th_max)) + this%dpsidth_min = this%dpsidth_from_th(this%th_min) return end subroutine set_wrf_param_cch @@ -1885,8 +1885,8 @@ function ftc_from_psi_tfs(this,psi) result(ftc) real(r8), pointer :: psi_min psi_min => this%wrf%psi_min - psi_eff = min(-nearzero,psi) - + psi_eff = max(psi_min,min(-nearzero,psi)) + ftc = 1._r8/(1._r8 + (psi_eff/this%p50)**this%avuln) if(ftc<=min_ftc) then @@ -1894,6 +1894,7 @@ function ftc_from_psi_tfs(this,psi) result(ftc) else ! Add protections and ensure no conductance at incredibly ! low suction + call get_min_ftc_weight(psi_min,psi_eff,min_ftc_weight,dmin_ftc_weight_dpsi) ftc = ftc*(1._r8 - min_ftc_weight) + min_ftc*min_ftc_weight @@ -1915,27 +1916,29 @@ function dftcdpsi_from_psi_tfs(this,psi) result(dftcdpsi) real(r8) :: dftcdpsi ! change in frac total cond wrt psi real(r8) :: min_ftc_weight,dmin_ftc_weight_dpsi real(r8), pointer :: psi_min + real(r8) :: psi_eff psi_min => this%wrf%psi_min + psi_eff = max(psi_min,min(-nearzero,psi)) ! Differentiate ! ftc = 1._r8/(1._r8 + (psi/this%p50(ft))**this%avuln(ft)) - if(psi>0._r8)then + if(psi_eff>0._r8)then dftcdpsi = 0._r8 else - ftc = this%ftc_from_psi(psi) + ftc = this%ftc_from_psi(psi_eff) if ( abs(ftc-min_ftc) sites(s)%si_hydr%wrf_soil(j)%p end do @@ -3128,7 +3128,7 @@ subroutine OrderLayersForSolve1D(csite_hydr,cohort,cohort_hydr,ordered, kbg_laye integer :: tmp ! temporarily holds a soil layer index integer :: ft ! functional type index of plant integer :: j,jj,k ! layer and shell indices - + real(r8), parameter :: neglibible_cond = 1.e-10_r8 kbg_tot = 0._r8 kbg_layer(:) = 0._r8 @@ -3164,7 +3164,9 @@ subroutine OrderLayersForSolve1D(csite_hydr,cohort,cohort_hydr,ordered, kbg_laye ! Calculate total effective conductance over path [kg s-1 MPa-1] ! from absorbing root node to 1st rhizosphere shell - r_bg = 1._r8/(kmax_aroot*ftc_aroot) + ! (since this is just about ordering, its ok to create a pseudo resistance + ! for nodes with zero conductance..) + r_bg = 1._r8/(kmax_aroot*max(ftc_aroot,neglibible_cond)) ! Path is across the upper an lower rhizosphere comparment ! on each side of the nodes. Since there is no flow across the outer @@ -3180,8 +3182,8 @@ subroutine OrderLayersForSolve1D(csite_hydr,cohort,cohort_hydr,ordered, kbg_laye ftc_shell = csite_hydr%wkf_soil(j)%p%ftc_from_psi(psi_shell) - r_bg = r_bg + 1._r8/(kmax_up*ftc_shell) - if(k wrf_plant(stomata_p_media,ft)%p + wkf_plant(stomata_p_media,ft)%p%wrf => wrf_plant(1,ft)%p end do diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 3c2093fa37..99e13fa5f4 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -41,6 +41,7 @@ module FatesInterfaceMod use FatesConstantsMod , only : n_landuse_cats use FatesConstantsMod , only : primaryland use FatesConstantsMod , only : secondaryland + use FatesConstantsMod , only : n_term_mort_types use FatesGlobals , only : fates_global_verbose use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun @@ -877,7 +878,7 @@ subroutine SetFatesGlobalElements2(use_fates) end if fates_maxElementsPerSite = max(fates_maxPatchesPerSite * fates_maxElementsPerPatch, & - numWatermem*numpft, num_vegtemp_mem, num_elements, nlevsclass*numpft) + numWatermem*numpft, num_vegtemp_mem, num_elements, nlevsclass*numpft*n_term_mort_types) if(hlm_use_planthydro==itrue)then fates_maxElementsPerSite = max(fates_maxElementsPerSite, nshell*nlevsoi_hyd_max ) From a269473b505ac731e5c4fc0be8205b2a89f3f767 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 3 Apr 2024 16:18:02 -0700 Subject: [PATCH 09/13] adjusted ordering of initializing bounds on cch functions to avoid undefined variables, reduce minimum suction of cch and fixed an error check --- biogeophys/FatesHydroWTFMod.F90 | 51 +++++++++++++++----------- biogeophys/FatesPlantHydraulicsMod.F90 | 31 ++++++++++------ 2 files changed, 49 insertions(+), 33 deletions(-) diff --git a/biogeophys/FatesHydroWTFMod.F90 b/biogeophys/FatesHydroWTFMod.F90 index 54833f979b..f83aa4d5d9 100644 --- a/biogeophys/FatesHydroWTFMod.F90 +++ b/biogeophys/FatesHydroWTFMod.F90 @@ -51,7 +51,7 @@ module FatesHydroWTFMod !real(r8), parameter :: min_theta_cch = 0.05_r8 ! Minimum theta (matches ctsm) - real(r8), parameter :: min_psi_cch = -20.0_r8 ! Minimum theta (matches ctsm) + real(r8), parameter :: min_psi_cch = -15.0_r8 ! Minimum theta (matches ctsm) ! Generic class that can be extended to describe ! specific water retention functions @@ -88,7 +88,7 @@ module FatesHydroWTFMod procedure, non_overridable :: psi_linear_res procedure, non_overridable :: th_linear_sat procedure, non_overridable :: th_linear_res - procedure, non_overridable :: set_min_max + procedure, non_overridable :: set_min_max_from_satres procedure, non_overridable :: get_thmin end type wrf_type @@ -267,23 +267,25 @@ module FatesHydroWTFMod ! of numerical integration. ! ============================================================================ - subroutine set_min_max(this,th_res,th_sat) + subroutine set_min_max_from_satres(this,th_res,th_sat) ! This routine uses max_sf_interp and min_sft_interp ! to define the bounds of where the linear ranges start and stop - + ! This only works for functions defined by a saturation and + ! a residual value + class(wrf_type) :: this real(r8),intent(in) :: th_res real(r8),intent(in) :: th_sat this%th_max = max_sf_interp*(th_sat-th_res)+th_res this%th_min = min_sf_interp*(th_sat-th_res)+th_res - this%psi_max = this%psi_from_th(this%th_max-tiny(this%th_max)) - this%dpsidth_max = this%dpsidth_from_th(this%th_max-tiny(this%th_max)) - this%psi_min = this%psi_from_th(this%th_min+tiny(this%th_min)) - this%dpsidth_min = this%dpsidth_from_th(this%th_min+tiny(this%th_min)) + this%psi_max = this%psi_from_th(this%th_max) + this%dpsidth_max = this%dpsidth_from_th(this%th_max) + this%psi_min = this%psi_from_th(this%th_min) + this%dpsidth_min = this%dpsidth_from_th(this%th_min) - end subroutine set_min_max + end subroutine set_min_max_from_satres ! ============================================================================ @@ -364,6 +366,8 @@ subroutine get_min_ftc_weight(psi_min,psi,min_ftc_weight,dmin_ftc_weight_dpsi) real(r8),intent(out):: min_ftc_weight ! weighting factor for min_ftc real(r8),intent(out):: dmin_ftc_weight_dpsi ! derivative of weighting factor wrt psi + ! If the difference between psi and psi-min is greater than 10MPa + ! just assume there is no effect of the minimum function (ie weight 0) min_ftc_weight = exp(min_ftc_scalar*max(psi_min-psi,-10._r8)) dmin_ftc_weight_dpsi = -min_ftc_scalar*exp(min_ftc_scalar*max(psi_min-psi,-10._r8)) @@ -471,7 +475,7 @@ subroutine set_wrf_param_vg(this,params_in) this%th_sat = params_in(4) this%th_res = params_in(5) - call this%set_min_max(this%th_res,this%th_sat) + call this%set_min_max_from_satres(this%th_res,this%th_sat) return end subroutine set_wrf_param_vg @@ -778,13 +782,15 @@ subroutine set_wrf_param_cch(this,params_in) ! Set DERIVED constants ! used for interpolating in extreme ranges this%th_max = max_sf_interp*this%th_sat - this%psi_max = this%psi_from_th(this%th_max-tiny(this%th_max)) - this%dpsidth_max = this%dpsidth_from_th(this%th_max-tiny(this%th_max)) - - !this%th_min = min_theta_cch - !this%psi_min = this%psi_from_th(min_theta_cch+tiny(this%th_max)) this%psi_min = min_psi_cch - this%th_min = this%th_from_psi(min_psi_cch+tiny(this%th_max)) + ! Need a temporary th_min(this can't be uninitialized while calculating psi_max) + this%th_min = 0.001_r8 + + this%psi_max = this%psi_from_th(this%th_max) + this%dpsidth_max = this%dpsidth_from_th(this%th_max) + + ! Get the actual th_min which is equivalent to psi_min + this%th_min = this%th_from_psi(min_psi_cch) this%dpsidth_min = this%dpsidth_from_th(this%th_min) return @@ -826,7 +832,7 @@ function th_from_psi_cch(this,psi) result(th) th = this%th_max + (psi-this%psi_max)/this%dpsidth_max else if(psi csite%si_hydr - csite_hydr%h2oveg = 0.0_r8 + csite_hydr%h2oveg = 0.0_r8 currentPatch => csite%oldest_patch do while(associated(currentPatch)) currentCohort=>currentPatch%tallest @@ -1757,13 +1757,14 @@ subroutine UpdateH2OVeg(csite,bc_out,prev_site_h2o,icall) currentPatch => currentPatch%younger enddo !end patch loop + ! convert from kg/site to kg/m2 csite_hydr%h2oveg = csite_hydr%h2oveg*AREA_INV ! Note that h2oveg_dead is incremented wherever we have litter fluxes ! and it will be reduced via an evaporation term - ! growturn_err is a term to accomodate error in growth or - ! turnover. need to be improved for future(CX) - bc_out%plant_stored_h2o_si = csite_hydr%h2oveg + csite_hydr%h2oveg_dead - & + ! growturn_err is a term to accomodate error in growth or + ! turnover. need to be improved for future(CX) + bc_out%plant_stored_h2o_si = csite_hydr%h2oveg + csite_hydr%h2oveg_dead - & csite_hydr%h2oveg_growturn_err - & csite_hydr%h2oveg_hydro_err @@ -2816,13 +2817,21 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) delta_soil_storage = sum(csite_hydr%h2osoi_liqvol_shell(:,:) * & csite_hydr%v_shell(:,:)) * denh2o * AREA_INV - prev_h2osoil - if(abs(delta_plant_storage - (root_flux - transp_flux)) > error_thresh ) then + ! This is to check closure and include the known error + ! The error is essentially the overestimate transpiration + ! versus change in state (q_top_eff*dt_substep) - (w_tot_beg-w_tot_end) + ! That is why we remove the error from the transpiration in this check + if(abs(delta_plant_storage - (root_flux + csite_hydr%errh2o_hyd - transp_flux)) > error_thresh ) then write(fates_log(),*) 'Site plant water balance does not close' + write(fates_log(),*) 'Allowable, actual error (kg/m2): ',error_thresh, & + abs(delta_plant_storage - (root_flux + csite_hydr%dwat_veg - transp_flux)) write(fates_log(),*) 'delta plant storage: ',delta_plant_storage,' [kg/m2]' write(fates_log(),*) 'integrated root flux: ',root_flux,' [kg/m2]' write(fates_log(),*) 'transpiration flux: ',transp_flux,' [kg/m2]' write(fates_log(),*) 'end storage: ',csite_hydr%h2oveg write(fates_log(),*) 'pre_h2oveg', prev_h2oveg + write(fates_log(),*) 'csite_hydr%errh2o_hyd:',csite_hydr%errh2o_hyd + write(fates_log(),*) 'csite_hydr%dwat_veg:',csite_hydr%dwat_veg call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -3225,7 +3234,7 @@ end subroutine OrderLayersForSolve1D subroutine ImTaylorSolve1D(slat, slon,recruitflag,csite_hydr,cohort,cohort_hydr,dtime,q_top, & ordered,kbg_layer, sapflow,rootuptake,& - wb_err_plant,dwat_plant,dth_layershell_col) + wb_err_ps,dwat_plant,dth_layershell_col) ! ------------------------------------------------------------------------------- ! Calculate the hydraulic conductances across a list of paths. The list is a 1D vector, and @@ -3258,8 +3267,8 @@ subroutine ImTaylorSolve1D(slat, slon,recruitflag,csite_hydr,cohort,cohort_hydr, real(r8),intent(out) :: sapflow ! time integrated mass flux between transp-root and stem [kg] real(r8),intent(out) :: rootuptake(:) ! time integrated mass flux between rhizosphere and aroot [kg] - real(r8),intent(out) :: wb_err_plant ! total error from the plant, transpiration - ! should match change in storage [kg] + real(r8),intent(out) :: wb_err_ps ! total error from the plant-soil system, transpiration + ! should match change in storage [kg] real(r8),intent(out) :: dwat_plant ! Change in plant stored water [kg] real(r8),intent(inout) :: dth_layershell_col(:,:) ! accumulated water content change over all cohorts in a column [m3 m-3]) @@ -3348,8 +3357,8 @@ subroutine ImTaylorSolve1D(slat, slon,recruitflag,csite_hydr,cohort,cohort_hydr, cohort_hydr%iterh1 = 0 cohort_hydr%iterh2 = 0 - ! Initialize plant water error (integrated flux-storage) - wb_err_plant = 0._r8 + ! Initialize plant-soil water error (integrated flux-storage) + wb_err_ps = 0._r8 ! Initialize integrated change in total plant water dwat_plant = 0._r8 @@ -3910,7 +3919,7 @@ subroutine ImTaylorSolve1D(slat, slon,recruitflag,csite_hydr,cohort,cohort_hydr, dth_node(n_hypool_ag+2)*cohort_hydr%v_aroot_layer(ilayer))*denh2o ! Remember the error for the cohort - wb_err_plant = wb_err_plant + wb_err_layer + wb_err_ps = wb_err_ps + wb_err_layer ! Save the change in water mass in the rhizosphere. Note that we did ! not immediately update the state variables upon completing each From cf054753fc1000276e24b5343efb8d3ceb00945e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 11 Apr 2024 13:33:47 -0700 Subject: [PATCH 10/13] Added btran to hydro restart, and protecting downscaling rhiz layer fluxes to soil layer fluxes when no conductance (trivial bug fix). --- biogeophys/FatesPlantHydraulicsMod.F90 | 48 ++++++++++++++-------- biogeophys/FatesPlantRespPhotosynthMod.F90 | 16 ++++++-- main/FatesHydraulicsMemMod.F90 | 11 +++-- 3 files changed, 51 insertions(+), 24 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index c4061a0788..874ab80a94 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -663,9 +663,6 @@ subroutine InitPlantHydStates(site, cohort) cohort_hydr%btran = wkf_plant(stomata_p_media,ft)%p%ftc_from_psi(cohort_hydr%psi_ag(1)) - - !flc_gs_from_psi(cohort_hydr%psi_ag(1),cohort%pft) - ! We do allow for positive pressures. ! But starting off with positive pressures is something we try to avoid if ( (cohort_hydr%psi_troot>0.0_r8) .or. & @@ -711,6 +708,8 @@ subroutine UpdatePlantPsiFTCFromTheta(ccohort,csite_hydr) ccohort_hydr%ftc_ag(k) = wkf_plant(leaf_p_media,ft)%p%ftc_from_psi(ccohort_hydr%psi_ag(k)) end do + ccohort_hydr%btran = wkf_plant(stomata_p_media,ft)%p%ftc_from_psi(ccohort_hydr%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)) @@ -2458,8 +2457,15 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) csite_hydr => sites(s)%si_hydr + bc_out(s)%qflx_soil2root_sisl(:) = 0._r8 + csite_hydr%sapflow_scpf(:,:) = 0._r8 + csite_hydr%rootuptake_sl(:) = 0._r8 + csite_hydr%rootuptake0_scpf(:,:) = 0._r8 + csite_hydr%rootuptake10_scpf(:,:) = 0._r8 + csite_hydr%rootuptake50_scpf(:,:) = 0._r8 + csite_hydr%rootuptake100_scpf(:,:) = 0._r8 + if( sum(csite_hydr%l_aroot_layer) == 0._r8 ) then - bc_out(s)%qflx_soil2root_sisl(:) = 0._r8 cycle end if @@ -2479,14 +2485,6 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) bc_out(s)%qflx_ro_sisl(:) = 0._r8 - ! Zero out diagnotsics that rely on accumulation - csite_hydr%sapflow_scpf(:,:) = 0._r8 - csite_hydr%rootuptake_sl(:) = 0._r8 - csite_hydr%rootuptake0_scpf(:,:) = 0._r8 - csite_hydr%rootuptake10_scpf(:,:) = 0._r8 - csite_hydr%rootuptake50_scpf(:,:) = 0._r8 - csite_hydr%rootuptake100_scpf(:,:) = 0._r8 - ! Initialize water mass balancing terms [kg h2o / m2] ! -------------------------------------------------------------------------------- transp_flux = 0._r8 @@ -2501,6 +2499,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ifp = 0 cpatch => sites(s)%oldest_patch do while (associated(cpatch)) + if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then ifp = ifp + 1 @@ -2532,7 +2531,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) end if ccohort=>cpatch%tallest - do while(associated(ccohort)) + co_loop1: do while(associated(ccohort)) ccohort_hydr => ccohort%co_hydr ft = ccohort%pft @@ -2670,9 +2669,8 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ccohort_hydr%btran = wkf_plant(stomata_p_media,ft)%p%ftc_from_psi(ccohort_hydr%psi_ag(1)) - ccohort => ccohort%shorter - enddo !cohort + enddo co_loop1 !cohort endif ! not bareground patch cpatch => cpatch%younger enddo !patch @@ -2763,9 +2761,9 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) sumweight = 0._r8 do j_bc = j_t,j_b if(rootflux_disagg == soilk_disagg)then - ! Weight disaggregation by K*dz, but only for flux - ! into the root, othersize weight by depth - if(qflx_soil2root_rhiz>0._r8)then + if(qflx_soil2root_rhiz>0._r8 )then + ! Weight disaggregation by K*dz, but only for flux + ! into the root, othersize weight by depth ! h2osoi_liqvol: [kg/m2] / [m] / [kg/m3] = [m3/m3] eff_por = bc_in(s)%eff_porosity_sl(j_bc) h2osoi_liqvol = min(eff_por, bc_in(s)%h2o_liq_sisl(j_bc)/(bc_in(s)%dz_sisl(j_bc)*denh2o)) @@ -2786,6 +2784,20 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) end do ! Second pass, apply normalized weighting factors for fluxes + + ! Note: It is possible that the soil is so dry there is no conductance + ! In these cases, solver error may create some non-zero, yet + ! trivially small fluxes. The conductances are exactly zero + ! so we need to renormalize the weighting function + if(sumweight < nearzero)then + sumweight = 0._r8 + do j_bc = j_t,j_b + weight_sl(j_bc) = csite_hydr%rootl_sl(j_bc)*bc_in(s)%h2o_liq_sisl(j_bc) + sumweight = sumweight + weight_sl(j_bc) + end do + end if + + do j_bc = j_t,j_b ! Fill the output array to the HLM diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index df0450016b..379baaff54 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -1313,6 +1313,9 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! empirical curvature parameter for ap photosynthesis co-limitation real(r8),parameter :: theta_ip = 0.999_r8 + ! minimum Leaf area to solve, too little has shown instability + real(r8), parameter :: min_la_to_solve = 0.0000000001_r8 + associate( bb_slope => EDPftvarcon_inst%bb_slope ,& ! slope of BB relationship, unitless medlyn_slope=> EDPftvarcon_inst%medlyn_slope , & ! Slope for Medlyn stomatal conductance model method, the unit is KPa^0.5 stomatal_intercept=> EDPftvarcon_inst%stomatal_intercept ) !Unstressed minimum stomatal conductance, the unit is umol/m**2/s @@ -1360,7 +1363,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! absorbed per unit leaf area. if(sunsha == 1)then !sunlit - if(( laisun_lsl * canopy_area_lsl) > 0.0000000001_r8)then + if(( laisun_lsl * canopy_area_lsl) > min_la_to_solve)then qabs = parsun_lsl / (laisun_lsl * canopy_area_lsl ) qabs = qabs * 0.5_r8 * (1._r8 - fnps) * 4.6_r8 @@ -1370,9 +1373,16 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in end if else - qabs = parsha_lsl / (laisha_lsl * canopy_area_lsl) - qabs = qabs * 0.5_r8 * (1._r8 - fnps) * 4.6_r8 + if( (parsha_lsl>nearzero) .and. (laisha_lsl * canopy_area_lsl) > min_la_to_solve ) then + qabs = parsha_lsl / (laisha_lsl * canopy_area_lsl) + qabs = qabs * 0.5_r8 * (1._r8 - fnps) * 4.6_r8 + else + ! The radiative transfer schemes are imperfect + ! they can sometimes generate negative values here + qabs = 0._r8 + end if + end if !convert the absorbed par into absorbed par per m2 of leaf, diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90 index f31ffe7e25..2baaa9d240 100644 --- a/main/FatesHydraulicsMemMod.F90 +++ b/main/FatesHydraulicsMemMod.F90 @@ -3,10 +3,12 @@ module FatesHydraulicsMemMod use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : fates_unset_r8 use FatesGlobals, only : fates_log + use FatesGlobals, only : endrun => fates_endrun use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use FatesConstantsMod, only : itrue,ifalse use FatesHydroWTFMod, only : wrf_arr_type use FatesHydroWTFMod, only : wkf_arr_type + use shr_log_mod , only : errMsg => shr_log_errMsg implicit none private @@ -319,7 +321,11 @@ module FatesHydraulicsMemMod procedure :: Dump end type ed_cohort_hydr_type - + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + + contains subroutine CopyCohortHydraulics(ncohort_hydr, ocohort_hydr) @@ -372,8 +378,7 @@ subroutine CopyCohortHydraulics(ncohort_hydr, ocohort_hydr) ncohort_hydr%iterh2 = ocohort_hydr%iterh2 ncohort_hydr%iterlayer = ocohort_hydr%iterlayer ncohort_hydr%errh2o = ocohort_hydr%errh2o - - + ! BC PLANT HYDRAULICS - flux terms ncohort_hydr%qtop = ocohort_hydr%qtop From 25fb73aa1f2b46f7fc806260198463f2b17f5066 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 29 Apr 2024 11:23:52 -0400 Subject: [PATCH 11/13] Updated some comments in low water potential protections --- biogeophys/FatesHydroWTFMod.F90 | 2 +- functional_unit_testing/hydro/HydroUTestDriver.py | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/biogeophys/FatesHydroWTFMod.F90 b/biogeophys/FatesHydroWTFMod.F90 index f83aa4d5d9..69c7fce9aa 100644 --- a/biogeophys/FatesHydroWTFMod.F90 +++ b/biogeophys/FatesHydroWTFMod.F90 @@ -51,7 +51,7 @@ module FatesHydroWTFMod !real(r8), parameter :: min_theta_cch = 0.05_r8 ! Minimum theta (matches ctsm) - real(r8), parameter :: min_psi_cch = -15.0_r8 ! Minimum theta (matches ctsm) + real(r8), parameter :: min_psi_cch = -15.0_r8 ! Minimum suction (MPa) ! Generic class that can be extended to describe ! specific water retention functions diff --git a/functional_unit_testing/hydro/HydroUTestDriver.py b/functional_unit_testing/hydro/HydroUTestDriver.py index 98c38a63c2..0d210d119c 100644 --- a/functional_unit_testing/hydro/HydroUTestDriver.py +++ b/functional_unit_testing/hydro/HydroUTestDriver.py @@ -162,7 +162,12 @@ def OMParams(zsoi): return(om_watsat,om_sucsat,om_bsw) def CCHParmsCosby84T5(zsoi,om_frac,sand_frac,clay_frac): - + + # Cosby, B.J., Hornberger, G.M., Clapp, R.B., and Ginn, T.R. 1984. + # A statistical exploration of the relationships of soil moisture + # characteristics to the physical properties of soils. Water Resour. + # Res. 20:682-690. + # cosby_1984_table5 # Get pedotransfer for soil matrix From a5d9fd953e0b76dd98c56d5874d662b3089b3256 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 29 Apr 2024 11:34:30 -0400 Subject: [PATCH 12/13] update to comments in hydro flux disaggregation code --- biogeophys/FatesPlantHydraulicsMod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 874ab80a94..8263b98937 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -2763,7 +2763,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) if(rootflux_disagg == soilk_disagg)then if(qflx_soil2root_rhiz>0._r8 )then ! Weight disaggregation by K*dz, but only for flux - ! into the root, othersize weight by depth + ! into the root, othersize weight by root length ! h2osoi_liqvol: [kg/m2] / [m] / [kg/m3] = [m3/m3] eff_por = bc_in(s)%eff_porosity_sl(j_bc) h2osoi_liqvol = min(eff_por, bc_in(s)%h2o_liq_sisl(j_bc)/(bc_in(s)%dz_sisl(j_bc)*denh2o)) @@ -2787,8 +2787,8 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ! Note: It is possible that the soil is so dry there is no conductance ! In these cases, solver error may create some non-zero, yet - ! trivially small fluxes. The conductances are exactly zero - ! so we need to renormalize the weighting function + ! trivially small fluxes. Lets create a simple weighting + ! function based on root length and available water. if(sumweight < nearzero)then sumweight = 0._r8 do j_bc = j_t,j_b From ec2720740693bf2507a8319aff460d726994a7b3 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 1 May 2024 15:03:20 -0600 Subject: [PATCH 13/13] disaggregating fluxes based on root length only when sumweights is zero, to enable first time-step issues in ctsm --- biogeophys/FatesPlantHydraulicsMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 8263b98937..454bf5d28d 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -2788,11 +2788,11 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ! Note: It is possible that the soil is so dry there is no conductance ! In these cases, solver error may create some non-zero, yet ! trivially small fluxes. Lets create a simple weighting - ! function based on root length and available water. + ! function based on root length. if(sumweight < nearzero)then sumweight = 0._r8 do j_bc = j_t,j_b - weight_sl(j_bc) = csite_hydr%rootl_sl(j_bc)*bc_in(s)%h2o_liq_sisl(j_bc) + weight_sl(j_bc) = csite_hydr%rootl_sl(j_bc) sumweight = sumweight + weight_sl(j_bc) end do end if