From 33dc1db3edf7978832e5715868928b6269a2a9b6 Mon Sep 17 00:00:00 2001 From: Chonggang Xu Date: Wed, 6 Apr 2022 20:39:07 -0600 Subject: [PATCH 1/7] put a transpiration flux check to increase numerical stability --- biogeophys/FatesPlantHydraulicsMod.F90 | 93 ++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 5571adc089..cd8b577761 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -2475,6 +2475,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 +2592,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 numberical 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 CaculateTotalAvailW(ccohort,site_hydr,bc_in(s),dtime, availWater); + if(qflx_tran_veg_indiv*dtime > 0.8_r8*availWater)then + site_hydr%trans_err = site_hydr%trans_err + dtime*ccohort%n*(qflx_tran_veg_indiv & + - 0.8_r8*availWater) + qflx_tran_veg_indiv=0.8_r8*availWater + endif + + endif ! Save the transpiration flux for diagnostics (currently its a constant boundary condition) ccohort_hydr%qtop = qflx_tran_veg_indiv*dtime @@ -4316,6 +4335,80 @@ subroutine AccumulateMortalityWaterStorage(csite,ccohort,delta_n) return end subroutine AccumulateMortalityWaterStorage +!-------------------------------------------------------------------------------! + +subroutine CaculateTotalAvailW(cohort,site_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 :: cohort + type(ed_site_hydr_type), intent(inout),target :: site_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 :: cohort_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 + + + cohort_hydr => cohort%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 = cohort_hydr%v_ag(i) + th_node = cohort_hydr%th_ag(i) + thr = EDPftvarcon_inst%hydr_resid_node(cohort%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 = cohort_hydr%v_troot + th_node = cohort_hydr%th_troot + thr = EDPftvarcon_inst%hydr_resid_node(cohort%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,site_hydr%nlevrhiz + v_node = cohort_hydr%v_aroot_layer(ilayer) + th_node = cohort_hydr%th_aroot(ilayer) + thr = EDPftvarcon_inst%hydr_resid_node(cohort%pft,i) + totalAvailW = totalAvailW +max(th_node-thr,0._r8)*v_node + enddo + else + do ilayer=1,site_hydr%nlevrhiz + thr=site_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(cohort_hydr%l_aroot_layer(ilayer)>nearzero)then + aroot_frac_plant = cohort_hydr%l_aroot_layer(ilayer)/site_hydr%l_aroot_layer(ilayer) + else + aroot_frac_plant = 0._r8 + end if + ! The volume of the Rhizosphere for a single plant + v_node = site_hydr%v_shell(ilayer,ishell)*aroot_frac_plant + th_node = site_hydr%h2osoi_liqvol_shell(ilayer,ishell) + shell_sum_v = sum(site_hydr%h2osoi_liqvol_shell(ilayer,:)) + recruit_water_demand = 0._r8 + if(shell_sum_v>0._r8.and. & + (.not.isnan(site_hydr%recruit_w_uptake(ilayer))))then + recruit_water_demand = site_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 CaculateTotalAvailW + + !-------------------------------------------------------------------------------! subroutine RecruitWaterStorage(nsites,sites,bc_out) From aa290b3b8ad95dc86214cb3ac0ff7772128514cf Mon Sep 17 00:00:00 2001 From: Chonggang Xu Date: Wed, 6 Apr 2022 20:40:00 -0600 Subject: [PATCH 2/7] add a transpiration error term --- main/FatesHydraulicsMemMod.F90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) 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)) From a80641d65fcf16030f73c57b27c1868845dd1243 Mon Sep 17 00:00:00 2001 From: Chonggang Xu Date: Thu, 7 Apr 2022 14:06:26 -0600 Subject: [PATCH 3/7] add the transpiration error term for mass balance check --- biogeophys/FatesPlantHydraulicsMod.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index cd8b577761..cc587f0aa2 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 From eec7cf713fe5cb83d86ab1e338eacfd4d2c6ec23 Mon Sep 17 00:00:00 2001 From: Marcos Longo <5891904+mpaiao@users.noreply.github.com> Date: Mon, 7 Nov 2022 11:30:19 -0800 Subject: [PATCH 4/7] Update FatesPlantHydraulicsMod.F90 Fix a few typos in FatesPlantHydraulicsMod (no change in the code functionality). --- biogeophys/FatesPlantHydraulicsMod.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index cc587f0aa2..130d612cbc 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -2595,7 +2595,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) end if !Hydro uses the transpiration flux from last time step - !For numberical stability, make sure that transpiration flux is less than 80% of total water + !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 @@ -2603,7 +2603,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) !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 CaculateTotalAvailW(ccohort,site_hydr,bc_in(s),dtime, availWater); + call CalculateTotalAvailW(ccohort,site_hydr,bc_in(s),dtime, availWater); if(qflx_tran_veg_indiv*dtime > 0.8_r8*availWater)then site_hydr%trans_err = site_hydr%trans_err + dtime*ccohort%n*(qflx_tran_veg_indiv & - 0.8_r8*availWater) @@ -4338,7 +4338,7 @@ end subroutine AccumulateMortalityWaterStorage !-------------------------------------------------------------------------------! -subroutine CaculateTotalAvailW(cohort,site_hydr,bc_in,dtime,totalAvailW) +subroutine CalculateTotalAvailW(cohort,site_hydr,bc_in,dtime,totalAvailW) ! --------------------------------------------------------------------------- ! This subroutine estimates the total available water for transpiration ! for an individual plant @@ -4407,7 +4407,7 @@ subroutine CaculateTotalAvailW(cohort,site_hydr,bc_in,dtime,totalAvailW) end if end do -end subroutine CaculateTotalAvailW +end subroutine CalculateTotalAvailW !-------------------------------------------------------------------------------! From 2ae1f8005cb9a5918e189eb400859171ebcf9f3d Mon Sep 17 00:00:00 2001 From: Marcos Longo <5891904+mpaiao@users.noreply.github.com> Date: Fri, 13 Jan 2023 11:49:31 -0300 Subject: [PATCH 5/7] Unit conversion fix: FatesPlantHydraulicsMod.F90 Revised the transpiration error quantification. The previous version seemed to have the wrong units, and BalanceCheckMod.F90 was often failing because errh2o was several orders of magnitude larger than the tolerance. With these edits, all my test simulations seem to be going through, although I have not made extensive tests. @xuchongang, could you please check and see if these changes make sense? Thanks! --- biogeophys/FatesPlantHydraulicsMod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 130d612cbc..c19e402573 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -2605,9 +2605,9 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) if (qflx_tran_veg_indiv>0.0_r8)then call CalculateTotalAvailW(ccohort,site_hydr,bc_in(s),dtime, availWater); if(qflx_tran_veg_indiv*dtime > 0.8_r8*availWater)then - site_hydr%trans_err = site_hydr%trans_err + dtime*ccohort%n*(qflx_tran_veg_indiv & - - 0.8_r8*availWater) - qflx_tran_veg_indiv=0.8_r8*availWater + site_hydr%trans_err = site_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 From e45ae4461c8330c76a45c0076df4f01d97740c98 Mon Sep 17 00:00:00 2001 From: Marcos Longo <5891904+mpaiao@users.noreply.github.com> Date: Fri, 27 Jan 2023 05:53:49 -0800 Subject: [PATCH 6/7] Update FatesPlantHydraulicsMod.F90 Fix minor typos to my previous commit. site_hydr was renamed to csite_hydr at some point, and my previous commit was using the old notation. I fixed that and also edited CalculateTotalAvailW to use csite and ccohort to be consistent with the other subroutines. --- biogeophys/FatesPlantHydraulicsMod.F90 | 52 +++++++++++++------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index c19e402573..8fffd4fb70 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -2605,8 +2605,8 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) if (qflx_tran_veg_indiv>0.0_r8)then call CalculateTotalAvailW(ccohort,site_hydr,bc_in(s),dtime, availWater); if(qflx_tran_veg_indiv*dtime > 0.8_r8*availWater)then - site_hydr%trans_err = site_hydr%trans_err & - + ccohort%n*AREA_INV*(dtime*qflx_tran_veg_indiv - 0.8_r8*availWater) + 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 @@ -4338,17 +4338,17 @@ end subroutine AccumulateMortalityWaterStorage !-------------------------------------------------------------------------------! -subroutine CalculateTotalAvailW(cohort,site_hydr,bc_in,dtime,totalAvailW) +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 :: cohort - type(ed_site_hydr_type), intent(inout),target :: site_hydr ! ED site_hydr structure + 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 :: cohort_hydr + 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) @@ -4358,46 +4358,46 @@ subroutine CalculateTotalAvailW(cohort,site_hydr,bc_in,dtime,totalAvailW) real(r8)::availW,recruit_water_demand - cohort_hydr => cohort%co_hydr + 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 = cohort_hydr%v_ag(i) - th_node = cohort_hydr%th_ag(i) - thr = EDPftvarcon_inst%hydr_resid_node(cohort%pft,i) + 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 = cohort_hydr%v_troot - th_node = cohort_hydr%th_troot - thr = EDPftvarcon_inst%hydr_resid_node(cohort%pft,i) + 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,site_hydr%nlevrhiz - v_node = cohort_hydr%v_aroot_layer(ilayer) - th_node = cohort_hydr%th_aroot(ilayer) - thr = EDPftvarcon_inst%hydr_resid_node(cohort%pft,i) + 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,site_hydr%nlevrhiz - thr=site_hydr%wrf_soil(ilayer)%p%th_from_psi(& + 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(cohort_hydr%l_aroot_layer(ilayer)>nearzero)then - aroot_frac_plant = cohort_hydr%l_aroot_layer(ilayer)/site_hydr%l_aroot_layer(ilayer) + 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 = site_hydr%v_shell(ilayer,ishell)*aroot_frac_plant - th_node = site_hydr%h2osoi_liqvol_shell(ilayer,ishell) - shell_sum_v = sum(site_hydr%h2osoi_liqvol_shell(ilayer,:)) + 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(site_hydr%recruit_w_uptake(ilayer))))then - recruit_water_demand = site_hydr%recruit_w_uptake(ilayer)* & + (.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 From 17f2b859124ffefb46b08dcc7a05c1e0880be352 Mon Sep 17 00:00:00 2001 From: Marcos Longo <5891904+mpaiao@users.noreply.github.com> Date: Fri, 27 Jan 2023 05:58:14 -0800 Subject: [PATCH 7/7] Update FatesPlantHydraulicsMod.F90 Additional update of site_hydr name to csite_hydr. --- biogeophys/FatesPlantHydraulicsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 8fffd4fb70..7efea5a403 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -2603,7 +2603,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) !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,site_hydr,bc_in(s),dtime, availWater); + 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)