improve hydro stability

biogeophys/FatesPlantHydraulicsMod.F90

@@ -1809,7 +1809,8 @@ subroutine UpdateH2OVeg(csite,bc_out,prev_site_h2o,icall)
     ! 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
+          csite_hydr%h2oveg_hydro_err-&
+          csite_hydr%trans_err
 
 
     ! Perform a conservation check if desired
@@ -2475,6 +2476,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime)
 
   integer, parameter :: soilz_disagg = 0   ! disaggregate rhizosphere layers based on depth
   integer, parameter :: soilk_disagg = 1   ! disaggregate rhizosphere layers based on conductance
+  real(r8) :: availWater          ! available water for transpriation [kg water/plant]
 
   integer, parameter :: rootflux_disagg = soilk_disagg
 
@@ -2591,6 +2593,24 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime)
               else
                  qflx_tran_veg_indiv     = 0._r8
               end if
+
+	      !Hydro uses the transpiration flux from last time step
+              !For numerical stability, make sure that transpiration flux is less than 80% of total water
+              !availability. If it is larger than that, we will adjust the
+              !transpiration flux and record the excessive demand to the plant
+              !water storage pool, which will be dumped later in HLM
+              !This should not affect the science of hydrodynamics as this
+              !basically delays the water limition by a few steps with a benefit
+              !of numerical stability. 
+              if (qflx_tran_veg_indiv>0.0_r8)then
+                  call CalculateTotalAvailW(ccohort,csite_hydr,bc_in(s),dtime, availWater);
+                  if(qflx_tran_veg_indiv*dtime > 0.8_r8*availWater)then
+                     csite_hydr%trans_err = csite_hydr%trans_err &
+                                          + ccohort%n*AREA_INV*(dtime*qflx_tran_veg_indiv - 0.8_r8*availWater)
+                     qflx_tran_veg_indiv = 0.8_r8 * availWater / dtime
+                  endif
+
+              endif
 
               ! Save the transpiration flux for diagnostics (currently its a constant boundary condition)
               ccohort_hydr%qtop = qflx_tran_veg_indiv*dtime
@@ -4316,6 +4336,80 @@ subroutine AccumulateMortalityWaterStorage(csite,ccohort,delta_n)
    return
 end subroutine AccumulateMortalityWaterStorage
 
+!-------------------------------------------------------------------------------!
+
+subroutine CalculateTotalAvailW(ccohort,csite_hydr,bc_in,dtime,totalAvailW)
+  ! ---------------------------------------------------------------------------
+  ! This subroutine estimates the total available water for transpiration
+  ! for an individual plant
+  ! ---------------------------------------------------------------------------
+  type(ed_cohort_type) , intent(inout), target  :: ccohort  
+  type(ed_site_hydr_type), intent(inout),target :: csite_hydr        ! ED site_hydr structure
+  type(bc_in_type),intent(in) :: bc_in
+  real(r8), intent(in)::dtime !time step (seconds)
+  real(r8), intent(out)::totalAvailW
+  type(ed_cohort_hydr_type), pointer :: ccohort_hydr
+  integer :: i, ilayer,ishell
+  real (r8) :: aroot_frac_plant !fraction of absorbing rooting in the layer
+  real(r8)::thr   !residual water content (m3/m3)
+  real(r8)::th_node   !residual water content for the node (m3/m3)
+  real(r8)::v_node   !volume of the node (m3)
+  real(r8)::shell_sum_v   !total volume of the shells(m3)
+  real(r8)::availW,recruit_water_demand
+
+
+  ccohort_hydr => ccohort%co_hydr
+
+  totalAvailW = 0._r8
+  do i = 1,n_hypool_tot
+
+        if (i<=n_hypool_ag) then                    ! leaf and stem, n_hypool_ag= 2
+           v_node  = ccohort_hydr%v_ag(i)
+           th_node = ccohort_hydr%th_ag(i) 
+           thr     = EDPftvarcon_inst%hydr_resid_node(ccohort%pft,i) 
+           totalAvailW = totalAvailW +max(th_node-thr,0._r8)*v_node
+        elseif (i==n_hypool_ag+1) then              ! i=3, transport root
+           v_node  = ccohort_hydr%v_troot
+           th_node = ccohort_hydr%th_troot
+           thr     = EDPftvarcon_inst%hydr_resid_node(ccohort%pft,i)
+           totalAvailW = totalAvailW +max(th_node-thr,0._r8)*v_node
+        elseif (i==n_hypool_ag+2) then              ! i=4, fine roots
+           do ilayer=1,csite_hydr%nlevrhiz
+              v_node  = ccohort_hydr%v_aroot_layer(ilayer)
+              th_node = ccohort_hydr%th_aroot(ilayer)
+              thr     = EDPftvarcon_inst%hydr_resid_node(ccohort%pft,i)
+              totalAvailW = totalAvailW +max(th_node-thr,0._r8)*v_node
+           enddo
+        else
+           do ilayer=1,csite_hydr%nlevrhiz
+              thr=csite_hydr%wrf_soil(ilayer)%p%th_from_psi(&
+                        bc_in%smpmin_si*denh2o*grav_earth*m_per_mm*mpa_per_pa)
+              ishell  = i-(n_hypool_ag+2)              ! i>=5, rhizosphere
+	      if(ccohort_hydr%l_aroot_layer(ilayer)>nearzero)then
+                aroot_frac_plant = ccohort_hydr%l_aroot_layer(ilayer)/csite_hydr%l_aroot_layer(ilayer)
+              else
+                aroot_frac_plant = 0._r8
+              end if
+              ! The volume of the Rhizosphere for a single plant
+              v_node  = csite_hydr%v_shell(ilayer,ishell)*aroot_frac_plant
+              th_node = csite_hydr%h2osoi_liqvol_shell(ilayer,ishell)
+	      shell_sum_v = sum(csite_hydr%h2osoi_liqvol_shell(ilayer,:))
+              recruit_water_demand = 0._r8
+              if(shell_sum_v>0._r8.and. &
+                    (.not.isnan(csite_hydr%recruit_w_uptake(ilayer))))then
+	           recruit_water_demand = csite_hydr%recruit_w_uptake(ilayer)* &
+                      v_node/shell_sum_v*dtime*AREA
+              endif
+	      availW = max(th_node-thr,0._r8)*v_node
+              totalAvailW = totalAvailW +max(availw-recruit_water_demand, 0._r8)
+
+          enddo
+        end if
+     end do
+
+end subroutine CalculateTotalAvailW
+
+
 !-------------------------------------------------------------------------------!
 
 subroutine RecruitWaterStorage(nsites,sites,bc_out)

main/FatesHydraulicsMemMod.F90

@@ -120,7 +120,10 @@ module FatesHydraulicsMemMod
                                                     !  insufficient plant water available to 
                                                     !  support transpiration
 
-
+     real(r8) :: trans_err                          !error water pool for cases
+                                                    !where there is not enough
+                                                    !water in the plant or soil
+                                                    !for transpiration
      ! Useful diagnostics
      ! ----------------------------------------------------------------------------------
 
@@ -410,7 +413,8 @@ subroutine InitHydrSite(this,numpft,numlevsclass,hydr_solver_type,nlevsoil)
 
          this%h2oveg_growturn_err = 0.0_r8
          this%h2oveg_hydro_err    = 0.0_r8
-
+         this%trans_err    = 0.0_r8
+
          ! We have separate water transfer functions and parameters
          ! for each soil layer, and each plant compartment type
          allocate(this%wrf_soil(1:nlevrhiz))