Merge pull request #525 from rgknox/xuchongang/HYDRO_UPDATE

fixes to hydraulics

Author: rgknox

biogeophys/FatesPlantHydraulicsMod.F90

@@ -40,8 +40,9 @@ module FatesPlantHydraulicsMod
    use FatesConstantsMod, only : cm2_per_m2
    use FatesConstantsMod, only : g_per_kg
 
-   use EDParamsMod       , only : hydr_kmax_rsurf
-
+   use EDParamsMod       , only : hydr_kmax_rsurf1
+   !   use EDParamsMod       , only : hydr_kmax_rsurf2
+   
    use EDTypesMod        , only : ed_site_type
    use EDTypesMod        , only : ed_patch_type
    use EDTypesMod        , only : ed_cohort_type
@@ -141,6 +142,7 @@ module FatesPlantHydraulicsMod
    public :: CopyCohortHydraulics
    public :: FuseCohortHydraulics
    public :: updateSizeDepTreeHydProps
+   public :: updateWaterDepTreeHydProps
    public :: updateSizeDepTreeHydStates
    public :: initTreeHydStates
    public :: updateSizeDepRhizHydProps
@@ -150,6 +152,7 @@ module FatesPlantHydraulicsMod
    public :: SavePreviousRhizVolumes
    public :: UpdateTreeHydrNodes
    public :: UpdateTreeHydrLenVolCond
+   public :: UpdateWaterDepTreeHydrCond
    public :: ConstrainRecruitNumber
 
    !------------------------------------------------------------------------------
@@ -325,7 +328,8 @@ subroutine initTreeHydStates(site_p, cc_p, bc_in)
        !   watsat(c,j), watres(c,j), alpha_VG(c,j), n_VG(c,j), m_VG(c,j), l_VG(c,j), &
        !   smp)
        !ccohort_hydr%psi_aroot(j) = smp
-       ccohort_hydr%psi_aroot(j) = csite%si_hydr%psisoi_liq_innershell(j)
+       !ccohort_hydr%psi_aroot(j) = csite%si_hydr%psisoi_liq_innershell(j)
+       ccohort_hydr%psi_aroot(j) = -0.2_r8 !do not assume the equalibrium between soil and root
 
        call th_from_psi(ft, 4, ccohort_hydr%psi_aroot(j), ccohort_hydr%th_aroot(j), csite%si_hydr, bc_in )
     end do
@@ -346,6 +350,7 @@ subroutine initTreeHydStates(site_p, cc_p, bc_in)
     !hydrostatic equilibrium with the water potential immediately below
     dz = ccohort_hydr%z_node_ag(n_hypool_ag) - ccohort_hydr%z_node_troot(1)
     ccohort_hydr%psi_ag(n_hypool_ag) = ccohort_hydr%psi_troot(1) - 1.e-6_r8*denh2o*grav*dz
+    if (ccohort_hydr%psi_ag(n_hypool_ag)>0.0_r8) ccohort_hydr%psi_ag(n_hypool_ag) = -0.01_r8
     call th_from_psi(ft, 2, ccohort_hydr%psi_ag(n_hypool_ag), ccohort_hydr%th_ag(n_hypool_ag), csite%si_hydr, bc_in)
     do k=n_hypool_ag-1, 1, -1
        dz = ccohort_hydr%z_node_ag(k) - ccohort_hydr%z_node_ag(k+1)
@@ -543,9 +548,44 @@ subroutine updateSizeDepTreeHydProps(currentSite,ccohort,bc_in)
     ! volumes, and UpdateTreeHydrNodes is called prior to this.
 
     call UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in)
+      
+  end subroutine updateSizeDepTreeHydProps
+  
+  ! =====================================================================================
+  
+  subroutine updateWaterDepTreeHydProps(currentSite,ccohort,bc_in)
+
+
+    ! DESCRIPTION: Updates absorbing root length (total and its vertical distribution)
+    ! as well as the consequential change in the size of the 'representative' rhizosphere
+    ! shell radii, volumes, and compartment volumes of plant tissues
+
+    ! !USES:
+    use shr_sys_mod        , only : shr_sys_abort
+
+    ! ARGUMENTS:
+    type(ed_site_type)     , intent(in)             :: currentSite ! Site stuff
+    type(ed_cohort_type)   , intent(inout)          :: ccohort     ! current cohort pointer
+    type(bc_in_type)       , intent(in)             :: bc_in       ! Boundary Conditions
+
+    ! Locals
+    integer                            :: nlevsoi_hyd             ! Number of total soil layers
+    type(ed_cohort_hydr_type), pointer :: ccohort_hydr
+    integer                            :: ft
+
+    nlevsoi_hyd                =  currentSite%si_hydr%nlevsoi_hyd
+    ccohort_hydr               => ccohort%co_hydr
+    ft                         =  ccohort%pft
 
+    ! This updates plant compartment volumes, lengths and 
+    ! maximum conductances. Make sure for already
+    ! initialized vegetation, that SavePreviousCompartment
+    ! volumes, and UpdateTreeHydrNodes is called prior to this.
 
-  end subroutine updateSizeDepTreeHydProps
+    call UpdateWaterDepTreeHydrCond(currentSite,ccohort,nlevsoi_hyd,bc_in)
+    
+    
+  end subroutine updateWaterDepTreeHydProps
 
   ! =====================================================================================
 
@@ -616,7 +656,6 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in)
                                                      ! hydraulic conductance        [kg s-1 MPa-1]
       real(r8) :: kmax_tot                           ! total tree (leaf to root tip) 
                                                      ! hydraulic conductance        [kg s-1 MPa-1]
-      
       real(r8),parameter :: taper_exponent = 1._r8/3._r8 ! Savage et al. (2010) xylem taper exponent [-]
 
       ccohort_hydr => ccohort%co_hydr
@@ -810,10 +849,84 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in)
                ccohort_hydr%kmax_treebg_layer(j) = rootfr*ccohort_hydr%kmax_treebg_tot
             end do
          end if
+
       end if !check for bleaf
 
     end subroutine UpdateTreeHydrLenVolCond
+    
+    
+    
+   !=====================================================================================
 
+  subroutine UpdateWaterDepTreeHydrCond(currentSite,ccohort,nlevsoi_hyd,bc_in)
+    
+      ! -----------------------------------------------------------------------------------
+      ! This subroutine calculates update the conductivity for the soil-root interface,
+      ! depending on the plant water uptake/loss. 
+      ! we assume that the conductivitity for water uptake is larger than 
+      ! water loss due to composite regulation of resistance the roots
+      ! hydraulic vs osmostic with and without transpiration
+      ! Steudle, E. Water uptake by roots: effects of water deficit. 
+      ! J Exp Bot 51, 1531-1542, doi:DOI 10.1093/jexbot/51.350.1531 (2000).
+      ! -----------------------------------------------------------------------------------
+
+      ! Arguments
+      type(ed_site_type)    , intent(in)  :: currentSite   ! Site target
+      type(ed_cohort_type),intent(inout)  :: ccohort       ! cohort target
+      integer,intent(in)                  :: nlevsoi_hyd   ! number of soil hydro layers
+      type(bc_in_type)       , intent(in) :: bc_in         ! Boundary Conditions
+      
+      type(ed_cohort_hydr_type),pointer :: ccohort_hydr     ! Plant hydraulics structure
+      type(ed_site_hydr_type),pointer :: csite_hydr
+    
+      integer  :: j,k
+      real(r8) :: hksat_s                            ! hksat converted to units of 10^6sec 
+      real(r8) :: kmax_root_surf_total               ! maximum conducitivity for total root surface(kg water/Mpa/s)
+      real(r8) :: kmax_soil_total                    ! maximum conducitivity for from root surface to soil shell(kg water/Mpa/s) 
+                                                     ! which is equiv to   [kg m-1 s-1 MPa-1]    
+      real(r8) :: kmax_root_surf                     ! maximum conducitivity for unit root surface (kg water/m2 root area/Mpa/s) 
+
+      ! (RGK 4-2019) THE FOLLOWING SHOULD BE ADDED TO THE PARAMETER FILE IN NEXT GROUPING
+      real(r8), parameter :: hydr_kmax_rsurf2 = 0.0001_r8  ! kg water/m2 root area/Mpa/s
+
+
+      ccohort_hydr => ccohort%co_hydr
+      csite_hydr => currentSite%si_hydr
+      k = 1   !only for the first soil shell
+      do j=1, nlevsoi_hyd
+	 
+         hksat_s = bc_in%hksat_sisl(j) * 1.e-3_r8 * 1/grav * 1.e6_r8
+	 if(ccohort_hydr%psi_aroot(j)<csite_hydr%psisoi_liq_innershell(j))then
+	     kmax_root_surf = hydr_kmax_rsurf1
+	 else
+	     kmax_root_surf = hydr_kmax_rsurf2
+	 endif
+         kmax_root_surf_total = kmax_root_surf*2._r8*pi_const *csite_hydr%rs1(j)* &
+            csite_hydr%l_aroot_layer(j)	 
+         if(csite_hydr%r_node_shell(j,1) <= csite_hydr%rs1(j)) then
+	     !csite_hydr%kmax_upper_shell(j,k)  = large_kmax_bound
+             !csite_hydr%kmax_bound_shell(j,k)  = large_kmax_bound
+             !csite_hydr%kmax_lower_shell(j,k)  = large_kmax_bound
+             ccohort_hydr%kmax_innershell(j)  = kmax_root_surf_total
+
+          else
+
+	     kmax_soil_total = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
+                   log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
+
+             !csite_hydr%kmax_upper_shell(j,k)  = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
+	     !      log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
+             !csite_hydr%kmax_bound_shell(j,k)  = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
+	     !      log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
+             !csite_hydr%kmax_lower_shell(j,k)  = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
+	     !      log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
+
+             ccohort_hydr%kmax_innershell(j)  = (1._r8/kmax_root_surf_total + &
+		                       1._r8/kmax_soil_total)**(-1._r8) 
+          end if		 
+       end do
+  
+  end subroutine UpdateWaterDepTreeHydrCond
 
   ! =====================================================================================
   subroutine updateSizeDepTreeHydStates(currentSite,ccohort)
@@ -951,6 +1064,7 @@ subroutine CopyCohortHydraulics(newCohort, oldCohort)
      ! quantities indexed by soil layer
      ncohort_hydr%z_node_aroot       = ocohort_hydr%z_node_aroot
      ncohort_hydr%kmax_treebg_layer  = ocohort_hydr%kmax_treebg_layer
+     ncohort_hydr%kmax_innershell    = ocohort_hydr%kmax_innershell     
      ncohort_hydr%v_aroot_layer_init = ocohort_hydr%v_aroot_layer_init
      ncohort_hydr%v_aroot_layer      = ocohort_hydr%v_aroot_layer
      ncohort_hydr%l_aroot_layer      = ocohort_hydr%l_aroot_layer
@@ -1609,11 +1723,7 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in)
     real(r8)                       :: large_kmax_bound = 1.e4_r8   ! for replacing kmax_bound_shell wherever the 
                                                                    ! innermost shell radius is less than the assumed 
                                                                    ! absorbing root radius rs1
-    real(r8)                       :: kmax_root_surf               ! maximum conducitivity for unit root surface
-                                                                   ! (kg water/m2 root area/Mpa/s)
                                                                    ! 1.e-5_r8 from Rudinger et al 1994             	
-    real(r8)                       :: kmax_root_surf_total         !maximum conducitivity for total root surface(kg water/Mpa/s)
-    real(r8)                       :: kmax_soil_total              !maximum conducitivity for total root surface(kg water/Mpa/s)
     integer                        :: nlevsoi_hyd	  
 
     !-----------------------------------------------------------------------
@@ -1633,7 +1743,7 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in)
        enddo !cohort
        cPatch => cPatch%older
     enddo !patch
-    kmax_root_surf = hydr_kmax_rsurf
+   
     csite_hydr%l_aroot_1D = sum( csite_hydr%l_aroot_layer(:))
 
     ! update outer radii of column-level rhizosphere shells (same across patches and cohorts)
@@ -1646,6 +1756,9 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in)
     enddo
     call shellGeom( csite_hydr%l_aroot_1D, csite_hydr%rs1(1), AREA, sum(bc_in%dz_sisl(1:nlevsoi_hyd)), &
           csite_hydr%r_out_shell_1D(:), csite_hydr%r_node_shell_1D(:), csite_hydr%v_shell_1D(:))
+	  
+    !update the conductitivity for first soil shell is done at subroutine UpdateWaterDepTreeHydrCond
+    !which is dependant on whether it is water uptake or loss for every 30 minutes
 
     do j = 1,csite_hydr%nlevsoi_hyd
 
@@ -1654,50 +1767,8 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in)
        ! proceed only if the total absorbing root length (site-level) has changed in this layer
        if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then
 
-          do k = 1,nshell
-	     if(k == 1) then
-	        kmax_root_surf_total = kmax_root_surf*2._r8*pi_const *csite_hydr%rs1(j)* &
-		                       csite_hydr%l_aroot_layer(j)
-                if(csite_hydr%r_node_shell(j,k) <= csite_hydr%rs1(j)) then
-		   !csite_hydr%kmax_upper_shell(j,k)  = large_kmax_bound
-                   !csite_hydr%kmax_bound_shell(j,k)  = large_kmax_bound
-                   !csite_hydr%kmax_lower_shell(j,k)  = large_kmax_bound
-                   csite_hydr%kmax_upper_shell(j,k)  = kmax_root_surf_total
-                   csite_hydr%kmax_bound_shell(j,k)  = kmax_root_surf_total
-                   csite_hydr%kmax_lower_shell(j,k)  = kmax_root_surf_total
-
-                else
-
-		   kmax_soil_total = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
-                         log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
-
-                   !csite_hydr%kmax_upper_shell(j,k)  = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
-		   !      log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
-                   !csite_hydr%kmax_bound_shell(j,k)  = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
-		   !      log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
-                   !csite_hydr%kmax_lower_shell(j,k)  = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
-		   !      log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s
-
-                   csite_hydr%kmax_upper_shell(j,k)  = (1._r8/kmax_root_surf_total + &
-		                       1._r8/kmax_soil_total)**(-1._r8)    
-                   csite_hydr%kmax_bound_shell(j,k)  = (1._r8/kmax_root_surf_total + &
-		                       1._r8/kmax_soil_total)**(-1._r8) 
-                   csite_hydr%kmax_lower_shell(j,k)  = (1._r8/kmax_root_surf_total + &
-		                       1._r8/kmax_soil_total)**(-1._r8)
-                end if
-		if(j == 1) then
-                   if(csite_hydr%r_node_shell(j,k) <= csite_hydr%rs1(j)) then
-                     csite_hydr%kmax_upper_shell_1D(k)  = csite_hydr%kmax_upper_shell(1,k)
-                     csite_hydr%kmax_bound_shell_1D(k)  = csite_hydr%kmax_bound_shell(1,k)
-                     csite_hydr%kmax_lower_shell_1D(k)  = csite_hydr%kmax_lower_shell(1,k)
-                   else
-                      csite_hydr%kmax_upper_shell_1D(k) = csite_hydr%kmax_upper_shell(1,k)
-                      csite_hydr%kmax_bound_shell_1D(k) = csite_hydr%kmax_bound_shell(1,k)
-                      csite_hydr%kmax_lower_shell_1D(k) = csite_hydr%kmax_lower_shell(1,k)
-                   end if
-                end if
-             else
-                csite_hydr%kmax_upper_shell(j,k)        = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
+          do k = 2,nshell
+	        csite_hydr%kmax_upper_shell(j,k)        = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
                       log(csite_hydr%r_node_shell(j,k)/csite_hydr%r_out_shell(j,k-1))*hksat_s
                 csite_hydr%kmax_bound_shell(j,k)        = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / &
                       log(csite_hydr%r_node_shell(j,k)/csite_hydr%r_node_shell(j,k-1))*hksat_s
@@ -1711,16 +1782,14 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in)
                    csite_hydr%kmax_lower_shell_1D(k)    = 2._r8*pi_const*csite_hydr%l_aroot_1D / &
                          log(csite_hydr%r_out_shell_1D( k)/csite_hydr%r_node_shell_1D(k  ))*hksat_s
                 end if
-             end if
           enddo ! loop over rhizosphere shells
        end if !has l_aroot_layer changed?
     enddo ! loop over soil layers
 
 
-
     return
   end subroutine UpdateSizeDepRhizVolLenCon
-
+    
 
   ! =====================================================================================
 
@@ -2445,7 +2514,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime )
               ! [mm H2O/cohort/s] = [mm H2O / patch / s] / [cohort/patch]
 !!              qflx_tran_veg_patch_coh      = qflx_trans_patch_vol * qflx_rel_tran_coh
 
-
+	      call updateWaterDepTreeHydProps(sites(s),ccohort,bc_in(s))
 
               if(site_hydr%nlevsoi_hyd > 1) then
                  ! BUCKET APPROXIMATION OF THE SOIL-ROOT HYDRAULIC GRADIENT (weighted average across layers)
@@ -2516,6 +2585,9 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime )
 		   kmax_lower(                  1 : n_hypool_ag    ) = ccohort_hydr%kmax_lower(:)
 		   kmax_upper((        n_hypool_ag+1)              ) = ccohort_hydr%kmax_upper_troot
                    if(site_hydr%nlevsoi_hyd == 1) then
+		      site_hydr%kmax_upper_shell_1D(1) =  ccohort_hydr%kmax_innershell(1)
+		      site_hydr%kmax_lower_shell_1D(1) =  ccohort_hydr%kmax_innershell(1)
+		      site_hydr%kmax_bound_shell_1D(1) =  ccohort_hydr%kmax_innershell(1)
                       !! estimate troot-aroot and aroot-radial components as a residual:
                       !! 25% each of total (surface of aroots to leaves) resistance
                       kmax_bound((        n_hypool_ag+1):(n_hypool_ag+2 )) = 2._r8 * ccohort_hydr%kmax_treebg_tot
@@ -2645,11 +2717,14 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime )
                          v_node_1l((n_hypool_ag+n_hypool_troot+1)           ) = ccohort_hydr%v_aroot_layer(j)
                          v_node_1l((n_hypool_tot-nshell+1):(n_hypool_tot)) = site_hydr%v_shell(j,:) * &
                                ccohort_hydr%l_aroot_layer(j)/bc_in(s)%dz_sisl(j)
+			 site_hydr%kmax_bound_shell(j,1)=ccohort_hydr%kmax_innershell(j)
+			 site_hydr%kmax_upper_shell(j,1)=ccohort_hydr%kmax_innershell(j)
+			 site_hydr%kmax_lower_shell(j,1)=ccohort_hydr%kmax_innershell(j)
                          kmax_bound_1l(:) = 0._r8 
                          kmax_bound_shell_1l(:) = site_hydr%kmax_bound_shell(j,:) * &
                                                   ccohort_hydr%l_aroot_layer(j) / site_hydr%l_aroot_layer(j)
 
-
+                           
                          ! transporting-to-absorbing root conductance: factor of 2 means one-half of the total 
                          ! belowground resistance in layer j      
                          kmax_bound_1l((n_hypool_ag+1)) = 2._r8 * ccohort_hydr%kmax_treebg_layer(j)