diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 5571adc089..7efea5a403 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/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) diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90 index 65163e1c8c..b871132538 100644 --- a/main/FatesHydraulicsMemMod.F90 +++ b/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))