diff --git a/biogeophys/FatesHydroWTFMod.F90 b/biogeophys/FatesHydroWTFMod.F90
index 8356b2b165..69c7fce9aa 100644
--- a/biogeophys/FatesHydroWTFMod.F90
+++ b/biogeophys/FatesHydroWTFMod.F90
@@ -28,8 +28,14 @@ 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
 
+  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,6 +50,9 @@ module FatesHydroWTFMod
                                             ! elastic-caviation region
 
 
+  !real(r8), parameter :: min_theta_cch = 0.05_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
 
@@ -79,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
@@ -89,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
@@ -240,7 +250,8 @@ module FatesHydroWTFMod
      procedure :: set_wkf_param     => set_wkf_param_tfs
   end type wkf_type_tfs
 
-
+  public :: get_min_ftc_weight
+  
 contains
 
   ! =====================================================================================
@@ -256,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
 
   ! ============================================================================
 
@@ -341,7 +354,36 @@ 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
+
+    ! 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))
+
+    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(:)
@@ -433,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
@@ -611,7 +653,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
 
@@ -624,8 +670,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
 
@@ -633,6 +689,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)
@@ -654,7 +713,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
 
@@ -664,8 +727,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)<nearzero ) then
           dftcdpsi = 0._r8   ! We cap ftc, so derivative is zero
        else
 
@@ -684,6 +747,17 @@ function dftcdpsi_from_psi_vg(this,psi) result(dftcdpsi)
           
           dftcdpsi = 2._r8*(1._r8-t1*t2)*(t1*dt2 + t2*dt1)/t3 - &
                t3**(-2._r8)*dt3*(1._r8-t1*t2)**2._r8
+          
+          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 if
@@ -708,12 +782,17 @@ 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      = fates_unset_r8
-    this%psi_min     = fates_unset_r8
-    this%dpsidth_min = fates_unset_r8
+    this%psi_min     = min_psi_cch
+    ! 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
   end subroutine set_wrf_param_cch
 
@@ -751,8 +830,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(psi<this%psi_min) then
+           th = this%th_sat*(this%psi_min/this%psi_sat)**(-1.0_r8/this%beta)
+        else
+           th = this%th_sat*(psi/this%psi_sat)**(-1.0_r8/this%beta)
+        end if
     end if
     return
   end function th_from_psi_cch
@@ -767,8 +850,12 @@ function psi_from_th_cch(this,th) result(psi)
 
     if(th>this%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(th<this%th_min) then
+           psi = this%psi_sat*(this%th_min/this%th_sat)**(-this%beta)
+        else
+           psi = this%psi_sat*(th/this%th_sat)**(-this%beta)
+        end if
     end if
 
   end function psi_from_th_cch
@@ -783,9 +870,15 @@ function dpsidth_from_th_cch(this,th) result(dpsidth)
 
     ! Differentiate:
     if(th>this%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<this%th_min) then
+          ! We shouldn't really actually encounter this
+          ! scenario because we cap th at th_min, but have this here for edge cases
+          dpsidth = -this%beta*this%psi_sat/this%th_sat * (this%th_min/this%th_sat)**(-this%beta-1._r8)
+       else
+          dpsidth = -this%beta*this%psi_sat/this%th_sat * (th/this%th_sat)**(-this%beta-1._r8)
+       end if
     end if
 
   end function dpsidth_from_th_cch
@@ -798,7 +891,11 @@ function ftc_from_psi_cch(this,psi) result(ftc)
     real(r8),intent(in) :: psi
     real(r8)            :: psi_eff
     real(r8)            :: ftc
+    real(r8)            :: min_ftc_weight,dmin_ftc_weight_dpsi
+    real(r8), pointer   :: psi_min
 
+    psi_min => this%wrf%psi_min
+    
     ! th = th_sat * (psi/psi_sat)^(-1/b)
 
     ! ftc = (th/th_sat)^(2*b+3)
@@ -806,11 +903,19 @@ function ftc_from_psi_cch(this,psi) result(ftc)
     !     = ((psi/psi_sat)^(-1/b))^(2*b+3)
     !     = (psi/psi_sat)^(-2-3/b)
 
-
     psi_eff = min(psi,this%psi_sat)
 
     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
 
   ! ====================================================================================
@@ -820,15 +925,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%psi_sat)then
-       dftcdpsi = (-2._r8-3._r8/this%beta) / this%psi_sat * &
-            (psi/this%psi_sat)**(-3._r8-3._r8/this%beta)
+
+       ftc = this%ftc_from_psi(psi)
+       if ( abs(ftc-min_ftc)<nearzero ) then
+
+          dftcdpsi = 0._r8
+          
+       else
+          
+          dftcdpsi = (-2._r8-3._r8/this%beta) / this%psi_sat * &
+               (psi/this%psi_sat)**(-3._r8-3._r8/this%beta)
+          
+          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
     else
        dftcdpsi = 0._r8
     end if
@@ -888,7 +1015,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
@@ -915,7 +1043,7 @@ 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%th_min      = 1.e-8_r8
+    this%th_min      = fates_unset_r8
     this%psi_min     = fates_unset_r8
     this%dpsidth_min = fates_unset_r8
 
@@ -965,9 +1093,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
@@ -1037,8 +1167,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
 
@@ -1195,7 +1324,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.
@@ -1231,21 +1360,22 @@ end function dpsidth_from_th_smooth_cch
   function ftc_from_psi_smooth_cch(this,psi) result(ftc)
 
     class(wkf_type_smooth_cch) :: this
-    real(r8),intent(in) :: psi
-    real(r8)            :: ftc
-
-    real(r8)                     :: pc
-    real(r8)                     :: kr
-    real(r8)                     :: dkr_dP
-    !
-    real(r8)                     :: sat_res
-    real(r8)                     :: alpha
-    real(r8)                     :: lambda
-    real(r8)                     :: Se
-    real(r8)                     :: deltaPc
-    real(r8)                     :: dSe_dpc
-    real(r8)                     :: dkr_dSe
-
+    real(r8),intent(in)        :: psi
+    real(r8) :: ftc
+    real(r8) :: pc
+    real(r8) :: kr
+    real(r8) :: dkr_dP
+    real(r8) :: sat_res
+    real(r8) :: alpha
+    real(r8) :: lambda
+    real(r8) :: Se
+    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
+
+    psi_min => this%wrf%psi_min
     pc = psi
     sat_res = 0._r8
     alpha   = -1._r8/this%psi_sat
@@ -1269,8 +1399,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
 
@@ -1281,11 +1421,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
@@ -1293,46 +1432,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)<nearzero)then
 
-    if( pc <= this%scch_pu ) then
-       ! Unsaturated full Brooks-Corey regime.
-       ! Here, `pc <= pu < 0`.
-       Se      = (-alpha*pc)**(-lambda)
+       dftcdpsi = 0._r8
 
-       dSe_dpc = -lambda*Se/pc
+    else
 
-       kr      = Se ** (3.d0 + 2.d0/lambda)
+       pc = psi
+       sat_res = 0._r8
+       alpha   = -1._r8/this%psi_sat
+       lambda  = 1._r8/this%beta
 
-       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)
+       if( pc <= this%scch_pu ) then
+          ! Unsaturated full Brooks-Corey regime.
+          ! Here, `pc <= pu < 0`.
+          Se      = (-alpha*pc)**(-lambda)
 
-       dSe_dpc = deltaPc*(2*this%scch_b2 + 3*deltaPc*this%scch_b3)
+          dSe_dpc = -lambda*Se/pc
 
-       kr      = Se ** (2.5d0 + 2.d0/lambda)
+          kr      = Se ** (3.d0 + 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
-    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
 
@@ -1489,7 +1647,7 @@ subroutine set_wrf_param_tfs(this,params_in)
     this%cap_slp  = params_in(8)
     this%pmedia   = int(params_in(9))
 
-    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_tfs
@@ -1730,11 +1888,27 @@ 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 = max(psi_min,min(-nearzero,psi))
+    
+    ftc = 1._r8/(1._r8 + (psi_eff/this%p50)**this%avuln)
 
-    psi_eff = min(-nearzero,psi)
+    if(ftc<=min_ftc) then
+       ftc = min_ftc
+    else
+       ! Add protections and ensure no conductance at incredibly
+       ! low suction
 
-    ftc = max(min_ftc,1._r8/(1._r8 + (psi_eff/this%p50)**this%avuln))
+       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
 
   ! ====================================================================================
@@ -1747,20 +1921,40 @@ 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
+    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 = 1._r8/(1._r8 + (psi/this%p50)**this%avuln)
-       if(ftc<min_ftc) then
-          dftcdpsi = 0._r8
+
+       ftc = this%ftc_from_psi(psi_eff)
+       
+       if ( abs(ftc-min_ftc)<nearzero ) then
+          dftcdpsi = 0._r8   ! We cap ftc, so derivative is zero
        else
-          fx  = 1._r8 + (psi/this%p50)**this%avuln
-          dfx = this%avuln*(psi/this%p50)**(this%avuln-1._r8) * (1._r8/this%p50)
+
+          fx  = 1._r8 + (psi_eff/this%p50)**this%avuln
+          dfx = this%avuln*(psi_eff/this%p50)**(this%avuln-1._r8) * (1._r8/this%p50)
           dftcdpsi = -fx**(-2._r8)*dfx
+          
+          call get_min_ftc_weight(psi_min,psi_eff,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 if
 
diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90
index 03977fa0ec..454bf5d28d 100644
--- a/biogeophys/FatesPlantHydraulicsMod.F90
+++ b/biogeophys/FatesPlantHydraulicsMod.F90
@@ -498,6 +498,15 @@ subroutine RestartHydrStates(sites,nsites,bc_in,bc_out)
           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
+
+       
        ! Update static quantities related to the rhizosphere
        call UpdateSizeDepRhizVolLenCon(sites(s), bc_in(s))
 
@@ -654,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. &
@@ -702,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))
@@ -1656,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
 
@@ -1720,7 +1737,7 @@ subroutine UpdateH2OVeg(csite,bc_out,prev_site_h2o,icall)
   if( hlm_use_planthydro.eq.ifalse ) return
 
     csite_hydr => 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
@@ -1740,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
 
@@ -2439,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
 
@@ -2460,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
@@ -2482,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
 
@@ -2513,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
@@ -2651,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
@@ -2744,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 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))
@@ -2767,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. Lets create a simple weighting
+           !       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)
+                 sumweight     = sumweight + weight_sl(j_bc)
+              end do
+           end if
+           
+           
            do j_bc = j_t,j_b
               
               ! Fill the output array to the HLM
@@ -2799,13 +2830,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
 
@@ -3111,7 +3150,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
@@ -3147,7 +3186,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
@@ -3163,8 +3204,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<nshell) r_bg = r_bg + 1._r8/(kmax_lo*ftc_shell )
+           r_bg = r_bg + 1._r8/(kmax_up*max(ftc_shell,neglibible_cond) )
+           if(k<nshell) r_bg = r_bg + 1._r8/(kmax_lo*max(ftc_shell,neglibible_cond) )
         end do
 
         !! upper bound limited to size()-1 b/c of zero-flux outer boundary condition
@@ -3206,7 +3247,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
@@ -3239,8 +3280,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])
 
@@ -3329,8 +3370,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
@@ -3891,7 +3932,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
@@ -3940,7 +3981,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
@@ -3994,28 +4035,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<nshell) then
+         write(fates_log(),*) '  mean soil-to-soil kmax with more outer shell: ',&
+              1._r8/(1._r8/(csite_hydr%kmax_lower_shell(ilayer,i)*aroot_frac_plant) + 1._r8/(csite_hydr%kmax_upper_shell(ilayer,i+1)*aroot_frac_plant))
+      end if
+   end do
+   
    write(fates_log(),*) 'kmax_aroot_radial_out: ',cohort_hydr%kmax_aroot_radial_out(ilayer)
    write(fates_log(),*) 'surf area of root: ',2._r8 * pi_const * EDPftvarcon_inst%hydr_rs2(ft) * cohort_hydr%l_aroot_layer(ilayer)
    write(fates_log(),*) 'aroot_frac_plant: ',aroot_frac_plant,cohort_hydr%l_aroot_layer(ilayer),csite_hydr%l_aroot_layer(ilayer)
@@ -4025,17 +4077,8 @@ subroutine Report1DError(cohort, csite_hydr, ilayer, z_node, v_node, &
    write(fates_log(),*) 'tree lai: ',cohort%treelai,' m2/m2 crown'
    write(fates_log(),*) 'area and area to volume ratios'
    write(fates_log(),*) ''
-   write(fates_log(),*) 'a:',v_node(4)
    write(fates_log(),*) '                      ',2._r8 * pi_const * EDPftvarcon_inst%hydr_rs2(ft) * cohort_hydr%l_aroot_layer(ilayer)
-   write(fates_log(),*) 'r1:',v_node(5)
    write(fates_log(),*) '                      ',2._r8 * pi_const * csite_hydr%r_out_shell(ilayer,1) * cohort_hydr%l_aroot_layer(ilayer)
-   write(fates_log(),*) 'r2:',v_node(6)
-   write(fates_log(),*) '                      '
-   write(fates_log(),*) 'r3:',v_node(7)
-   write(fates_log(),*) '                      '
-   write(fates_log(),*) 'r4:',v_node(8)
-   write(fates_log(),*) '                      '
-   write(fates_log(),*) 'r5:',v_node(9)
    write(fates_log(),*) 'inner shell kmaxs: ',csite_hydr%kmax_lower_shell(:,1)*aroot_frac_plant
 
    deallocate(psi_node)
@@ -4110,18 +4153,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
@@ -5977,7 +6032,6 @@ subroutine PicardSolve2D(csite_hydr,cohort,cohort_hydr, &
           if(abs(wb_error) < max_allowed_residual .or. maxval(abs(residual(:))) < 1.e-3_r8 .or. maxval(abs(th_node(:) - th_prev(:))) < 1.e-3) exit picardloop
 
           if(icnv == 1 ) then
-             print *,'dtime-',dtime,tm
              exit picardloop !explicit integration with small time step
           end if
 
@@ -6218,6 +6272,7 @@ subroutine InitHydroGlobals()
 
    integer :: ft            ! PFT index
    integer :: pm            ! plant media index
+   integer :: node_type     
    integer :: inode         ! compartment node index
    real(r8) :: cap_corr     ! correction for nonzero psi0x (TFS)
    real(r8) :: cap_slp      ! slope of capillary region of curve
@@ -6235,6 +6290,7 @@ subroutine InitHydroGlobals()
    ! -----------------------------------------------------------------------------------
 
    do pm = 1, n_plant_media
+
       select case(hydr_htftype_node(pm))
       case(van_genuchten_type)
          do ft = 1,numpft
@@ -6303,6 +6359,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
@@ -6316,6 +6380,15 @@ 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
+   ! (for stomatal conductance, point to the internal leaf retention structure)
+   do ft = 1,numpft
+      wkf_plant(stomata_p_media,ft)%p%wrf =>  wrf_plant(1,ft)%p
+   end do
+
+   
 
    return
 end subroutine InitHydroGlobals
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/functional_unit_testing/hydro/HydroUTestDriver.py b/functional_unit_testing/hydro/HydroUTestDriver.py
index 5c2b026361..0d210d119c 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
@@ -150,18 +150,65 @@ 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, 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
+    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 +225,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,16 +261,23 @@ 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')
 
 
+    #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)
 
-    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,10 +416,120 @@ 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 = 5
+    for icl in range(ntex):
+        clay_frac = float(icl)/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 
+
+    #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 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
diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90
index 61e97173c7..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
@@ -168,19 +170,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
@@ -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
 
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 )