Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Hydraulic redistribution #1126

Merged
merged 5 commits into from
Oct 21, 2020
Merged

Hydraulic redistribution #1126

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
d621803
working patch HR
djk2120 Aug 26, 2020
664eaba
add VPD_CAN and predawn and local noon vegwp
djk2120 Aug 27, 2020
ba8bd74
add vpd_can with PHS
djk2120 Aug 27, 2020
be4af0e
rework implementation of vegwp_pd
djk2120 Aug 27, 2020
241b788
enforce that vegwp_pd is at night and also before noon
djk2120 Aug 27, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 12 additions & 1 deletion src/biogeophys/CanopyStateType.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,8 @@ module CanopyStateType
real(r8) , pointer :: rscanopy_patch (:) ! patch canopy stomatal resistance (s/m) (ED specific)

real(r8) , pointer :: vegwp_patch (:,:) ! patch vegetation water matric potential (mm)
real(r8) , pointer :: vegwp_ln_patch (:,:) ! patch vegetation water matric potential at local noon (mm)
real(r8) , pointer :: vegwp_pd_patch (:,:) ! patch predawn vegetation water matric potential (mm)

real(r8) :: leaf_mr_vcm = spval ! Scalar constant of leaf respiration with Vcmax

Expand Down Expand Up @@ -131,7 +133,8 @@ subroutine InitAllocate(this, bounds)
allocate(this%rscanopy_patch (begp:endp)) ; this%rscanopy_patch (:) = nan
! allocate(this%gccanopy_patch (begp:endp)) ; this%gccanopy_patch (:) = 0.0_r8
allocate(this%vegwp_patch (begp:endp,1:nvegwcs)) ; this%vegwp_patch (:,:) = nan

allocate(this%vegwp_ln_patch (begp:endp,1:nvegwcs)) ; this%vegwp_ln_patch (:,:) = nan
allocate(this%vegwp_pd_patch (begp:endp,1:nvegwcs)) ; this%vegwp_pd_patch (:,:) = nan
end subroutine InitAllocate

!-----------------------------------------------------------------------
Expand Down Expand Up @@ -243,6 +246,14 @@ subroutine InitHistory(this, bounds)
call hist_addfld2d (fname='VEGWP', units='mm', type2d='nvegwcs', &
avgflag='A', long_name='vegetation water matric potential for sun/sha canopy,xyl,root segments', &
ptr_patch=this%vegwp_patch)
this%vegwp_ln_patch(begp:endp,:) = spval
call hist_addfld2d (fname='VEGWPLN', units='mm', type2d='nvegwcs', &
avgflag='A', long_name='vegetation water matric potential for sun/sha canopy,xyl,root at local noon', &
ptr_patch=this%vegwp_ln_patch)
this%vegwp_pd_patch(begp:endp,:) = spval
call hist_addfld2d (fname='VEGWPPD', units='mm', type2d='nvegwcs', avgflag='A', &
long_name='predawn vegetation water matric potential for sun/sha canopy,xyl,root', &
ptr_patch=this%vegwp_pd_patch)
end if

end subroutine InitHistory
Expand Down
46 changes: 40 additions & 6 deletions src/biogeophys/PhotosynthesisMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ module PhotosynthesisMod
use clm_varctl , only : use_c13, use_c14, use_cn, use_cndv, use_fates, use_luna, use_hydrstress
use clm_varctl , only : iulog
use clm_varpar , only : nlevcan, nvegwcs, mxpft
use clm_varcon , only : namep, c14ratio, spval
use clm_varcon , only : namep, c14ratio, spval, isecspday
use decompMod , only : bounds_type
use QuadraticMod , only : quadratic
use pftconMod , only : pftcon
Expand Down Expand Up @@ -148,7 +148,7 @@ module PhotosynthesisMod
real(r8), pointer, private :: gs_mol_patch (:,:) ! patch leaf stomatal conductance (umol H2O/m**2/s)
real(r8), pointer, private :: gb_mol_patch (:) ! patch leaf boundary layer conductance (umol H2O/m**2/s)
real(r8), pointer, private :: rh_leaf_patch (:) ! patch fractional humidity at leaf surface (dimensionless)

real(r8), pointer, private :: vpd_can_patch (:) ! patch canopy vapor pressure deficit (kPa)
real(r8), pointer, private :: alphapsnsun_patch (:) ! patch sunlit 13c fractionation ([])
real(r8), pointer, private :: alphapsnsha_patch (:) ! patch shaded 13c fractionation ([])

Expand Down Expand Up @@ -293,7 +293,7 @@ subroutine InitAllocate(this, bounds)
allocate(this%mbb_patch (begp:endp)) ; this%mbb_patch (:) = nan
allocate(this%gb_mol_patch (begp:endp)) ; this%gb_mol_patch (:) = nan
allocate(this%rh_leaf_patch (begp:endp)) ; this%rh_leaf_patch (:) = nan

allocate(this%vpd_can_patch (begp:endp)) ; this%vpd_can_patch (:) = nan
allocate(this%psnsun_patch (begp:endp)) ; this%psnsun_patch (:) = nan
allocate(this%psnsha_patch (begp:endp)) ; this%psnsha_patch (:) = nan
allocate(this%c13_psnsun_patch (begp:endp)) ; this%c13_psnsun_patch (:) = nan
Expand Down Expand Up @@ -376,6 +376,14 @@ subroutine InitHistory(this, bounds)
call hist_addfld1d (fname='RH_LEAF', units='fraction', &
avgflag='A', long_name='fractional humidity at leaf surface', &
ptr_patch=this%rh_leaf_patch, set_spec=spval, default='inactive')

this%vpd_can_patch(begp:endp) = spval
call hist_addfld1d (fname='VPD_CAN', units='kPa', &
avgflag='A', long_name='canopy vapor pressure deficit', &
ptr_patch=this%vpd_can_patch, set_spec=spval, default='inactive')



this%lnca_patch(begp:endp) = spval
call hist_addfld1d (fname='LNC', units='gN leaf/m^2', &
avgflag='A', long_name='leaf N concentration', &
Expand Down Expand Up @@ -1224,6 +1232,7 @@ subroutine Photosynthesis ( bounds, fn, filterp, &
bbb => photosyns_inst%bbb_patch , & ! Output: [real(r8) (:) ] Ball-Berry minimum leaf conductance (umol H2O/m**2/s)
mbb => photosyns_inst%mbb_patch , & ! Output: [real(r8) (:) ] Ball-Berry slope of conductance-photosynthesis relationship
rh_leaf => photosyns_inst%rh_leaf_patch , & ! Output: [real(r8) (:) ] fractional humidity at leaf surface (dimensionless)
vpd_can => photosyns_inst%vpd_can_patch , & ! Output: [real(r8) (:) ] canopy vapor pressure deficit (kPa)
lnc => photosyns_inst%lnca_patch , & ! Output: [real(r8) (:) ] top leaf layer leaf N concentration (gN leaf/m^2)
light_inhibit=> photosyns_inst%light_inhibit , & ! Input: [logical ] flag if light should inhibit respiration
leafresp_method=> photosyns_inst%leafresp_method , & ! Input: [integer ] method type to use for leaf-maint.-respiration at 25C canopy top
Expand Down Expand Up @@ -1646,6 +1655,7 @@ subroutine Photosynthesis ( bounds, fn, filterp, &
else if ( stomatalcond_mtd == stomatalcond_mtd_medlyn2011 )then
! Put some constraints on RH in the canopy when Medlyn stomatal conductance is being used
rh_can = max((esat_tv(p) - ceair), 50._r8) * 0.001_r8
vpd_can(p) = rh_can
end if

! Electron transport rate for C3 plants. Convert par from W/m2 to
Expand Down Expand Up @@ -2730,6 +2740,7 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, &
bbb => photosyns_inst%bbb_patch , & ! Output: [real(r8) (:) ] Ball-Berry minimum leaf conductance (umol H2O/m**2/s)
mbb => photosyns_inst%mbb_patch , & ! Output: [real(r8) (:) ] Ball-Berry slope of conductance-photosynthesis relationship
rh_leaf => photosyns_inst%rh_leaf_patch , & ! Output: [real(r8) (:) ] fractional humidity at leaf surface (dimensionless)
vpd_can => photosyns_inst%vpd_can_patch , & ! Output: [real(r8) (:) ] canopy vapor pressure deficit (kPa)
lnc => photosyns_inst%lnca_patch , & ! Output: [real(r8) (:) ] top leaf layer leaf N concentration (gN leaf/m^2)
light_inhibit=> photosyns_inst%light_inhibit , & ! Input: [logical ] flag if light should inhibit respiration
leafresp_method=> photosyns_inst%leafresp_method , & ! Input: [integer ] method type to use for leaf-maint.-respiration at 25C canopy top
Expand Down Expand Up @@ -3278,6 +3289,7 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, &
else if ( stomatalcond_mtd == stomatalcond_mtd_medlyn2011 )then
! Put some constraints on RH in the canopy when Medlyn stomatal conductance is being used
rh_can = max((esat_tv(p) - ceair), 50._r8) * 0.001_r8
vpd_can(p) = rh_can
end if

! Electron transport rate for C3 plants. Convert par from W/m2 to
Expand Down Expand Up @@ -3310,7 +3322,7 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, &
end if

!find ci and stomatal conductance
call hybrid_PHS(ci_z_sun(p,iv), ci_z_sha(p,iv), p, iv, c, gb_mol(p), bsun(p),bsha(p), je_sun, &
call hybrid_PHS(ci_z_sun(p,iv), ci_z_sha(p,iv), p, iv, c, g, gb_mol(p), bsun(p),bsha(p), je_sun, &
je_sha, cair(p), oair(p), lmr_z_sun(p,iv), lmr_z_sha(p,iv), &
par_z_sun(p,iv), par_z_sha(p,iv), rh_can, gs_mol_sun(p,iv), gs_mol_sha(p,iv), &
qsatl(p), qaf(p), iter1, iter2, atm2lnd_inst, photosyns_inst, &
Expand Down Expand Up @@ -3534,7 +3546,7 @@ end subroutine PhotosynthesisHydraulicStress
!------------------------------------------------------------------------------

!--------------------------------------------------------------------------------
subroutine hybrid_PHS(x0sun, x0sha, p, iv, c, gb_mol, bsun, bsha, jesun, jesha, &
subroutine hybrid_PHS(x0sun, x0sha, p, iv, c, g, gb_mol, bsun, bsha, jesun, jesha, &
cair, oair, lmr_z_sun, lmr_z_sha, par_z_sun, par_z_sha, rh_can, &
gs_mol_sun, gs_mol_sha, qsatl, qaf, iter1, iter2, atm2lnd_inst, photosyns_inst, &
canopystate_inst, waterdiagnosticbulk_inst, soilstate_inst, temperature_inst, waterfluxbulk_inst)
Expand All @@ -3551,13 +3563,15 @@ subroutine hybrid_PHS(x0sun, x0sha, p, iv, c, gb_mol, bsun, bsha, jesun, jesha,
!
!
!!USES:
use clm_time_manager , only : is_near_local_noon
!
!! ARGUMENTS:
implicit none
real(r8), intent(inout) :: x0sun,x0sha ! initial guess and final value of the solution for cisun/cisha
integer , intent(in) :: p ! pft index
integer , intent(in) :: iv ! radiation canopy layer index
integer , intent(in) :: c ! column index
integer , intent(in) :: g ! gridcell index
real(r8), intent(in) :: gb_mol ! leaf boundary layer conductance (umol H2O/m**2/s)
real(r8), intent(out) :: bsun ! sunlit canopy transpiration wetness factor (0 to 1)
real(r8), intent(out) :: bsha ! shaded canopy transpiration wetness factor (0 to 1)
Expand Down Expand Up @@ -3618,7 +3632,8 @@ subroutine hybrid_PHS(x0sun, x0sha, p, iv, c, gb_mol, bsun, bsha, jesun, jesha,

associate( &
qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Input: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm)
vegwp => canopystate_inst%vegwp_patch & ! Input/Output: [real(r8) (:,:) ] vegetation water matric potential (mm)
vegwp => canopystate_inst%vegwp_patch ,& ! Input/Output: [real(r8) (:,:) ] vegetation water matric potential (mm)
vegwp_ln => canopystate_inst%vegwp_ln_patch & ! Output: [real(r8) (:,:) ] vegetation water matric potential (mm) at local noon
)


Expand Down Expand Up @@ -3767,6 +3782,14 @@ subroutine hybrid_PHS(x0sun, x0sha, p, iv, c, gb_mol, bsun, bsha, jesun, jesha,
call getvegwp(p, c, x, gb_mol, gs_mol_sun, gs_mol_sha, qsatl, qaf, soilflux, &
atm2lnd_inst, canopystate_inst, waterdiagnosticbulk_inst, soilstate_inst, temperature_inst)
vegwp(p,:)=x

!write out local noon vwp (within +/- 1hr)
if ( is_near_local_noon( grc%londeg(g), deltasec=3600 ) )then
vegwp_ln(p,:) = vegwp(p,:)
else
vegwp_ln(p,:) = spval
end if

if (soilflux<0._r8) soilflux = 0._r8
qflx_tran_veg(p) = soilflux

Expand Down Expand Up @@ -4186,6 +4209,7 @@ subroutine calcstress(p,c,x,bsun,bsha,gb_mol,gs_mol_sun,gs_mol_sha,qsatl,qaf, &
! USES
use clm_varpar , only : nlevsoi
use clm_varcon , only : rgas
use clm_time_manager , only : get_local_time
!!
! !ARGUMENTS:
integer , intent(in) :: p ! pft index
Expand Down Expand Up @@ -4218,6 +4242,7 @@ subroutine calcstress(p,c,x,bsun,bsha,gb_mol,gs_mol_sun,gs_mol_sha,qsatl,qaf, &
real(r8) :: gs0sun,gs0sha ! local gs_mol copies
real(r8) :: qsun,qsha ! attenuated transpiration fluxes
integer :: j ! index
integer :: g ! gridcell index
real(r8) :: cf ! s m**2/umol -> s/m
integer :: iter ! newton's method iteration number
logical :: flag ! signal that matrix was not invertible
Expand All @@ -4241,6 +4266,7 @@ subroutine calcstress(p,c,x,bsun,bsha,gb_mol,gs_mol_sun,gs_mol_sha,qsatl,qaf, &
tgcm => temperature_inst%thm_patch , & ! Input: [real(r8) (:) ] air temperature at agcm reference height (kelvin)
bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b"
qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Input: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm)
vegwp_pd => canopystate_inst%vegwp_pd_patch , & ! Output: [real(r8) (:,:) ] vegetation water matric potential (mm) predawn
sucsat => soilstate_inst%sucsat_col & ! Input: [real(r8) (:,:) ] minimum soil suction (mm)
)

Expand Down Expand Up @@ -4380,6 +4406,14 @@ subroutine calcstress(p,c,x,bsun,bsha,gb_mol,gs_mol_sun,gs_mol_sha,qsatl,qaf, &
if (soilflux<0._r8) soilflux = 0._r8
qflx_tran_veg(p) = soilflux
endif

!save predawn vegwp
g = patch%gridcell(p)
if (night .and. get_local_time(grc%londeg(g))<(isecspday/2)) then
vegwp_pd(p,:) = x
else
vegwp_pd(p,:) = spval
end if


end associate
Expand Down
33 changes: 21 additions & 12 deletions src/biogeophys/SoilWaterPlantSinkMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -274,27 +274,30 @@ subroutine Compute_EffecRootFrac_And_VertTranSink_HydStress( bounds, &
integer :: pi ! patch index
real(r8) :: temp(bounds%begc:bounds%endc) ! accumulator for rootr weighting
real(r8) :: grav2 ! soil layer gravitational potential relative to surface (mm H2O)
real(r8) :: patchflux ! patch level soil-to-plant water flux (mm/s)
integer , parameter :: soil=1,root=4 ! index values
!-----------------------------------------------------------------------

associate(&
k_soil_root => soilstate_inst%k_soil_root_patch , & ! Input: [real(r8) (:,:) ]
k_soil_root => soilstate_inst%k_soil_root_patch , & ! Input: [real(r8) (:,:) ]
! soil-root interface conductance (mm/s)
qflx_phs_neg_col => waterfluxbulk_inst%qflx_phs_neg_col , & ! Input: [real(r8) (:) ] n
! net neg hydraulic redistribution flux(mm H2O/s)
qflx_tran_veg_col => waterfluxbulk_inst%qflx_tran_veg_col , & ! Input: [real(r8) (:) ]
qflx_phs_neg_col => waterfluxbulk_inst%qflx_phs_neg_col , & ! Output: [real(r8) (:) ]
! net neg hydraulic redistribution flux col (mm H2O/s)
qflx_hydr_redist_patch => waterfluxbulk_inst%qflx_hydr_redist_patch, & ! Output: [real(r8) (:) ]
! net neg hydraulic redistribution flux patch (mm H2O/s)
qflx_tran_veg_col => waterfluxbulk_inst%qflx_tran_veg_col , & ! Input: [real(r8) (:) ]
! vegetation transpiration (mm H2O/s) (+ = to atm)
qflx_tran_veg_patch => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Input: [real(r8) (:) ]
qflx_tran_veg_patch => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Input: [real(r8) (:) ]
! vegetation transpiration (mm H2O/s) (+ = to atm)
qflx_rootsoi_col => waterfluxbulk_inst%qflx_rootsoi_col , & ! Output: [real(r8) (:) ]
qflx_rootsoi_col => waterfluxbulk_inst%qflx_rootsoi_col , & ! Output: [real(r8) (:) ]
! col root and soil water
! exchange [mm H2O/s] [+ into root]
smp => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix pot. [mm]
frac_veg_nosno => canopystate_inst%frac_veg_nosno_patch , & ! Input: [integer (:) ]
smp => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix pot. [mm]
frac_veg_nosno => canopystate_inst%frac_veg_nosno_patch , & ! Input: [integer (:) ]
! fraction of vegetation not
! covered by snow (0 OR 1) [-]
z => col%z , & ! Input: [real(r8) (:,:) ] layer node depth (m)
vegwp => canopystate_inst%vegwp_patch & ! Input: [real(r8) (:,:) ] vegetation water
z => col%z , & ! Input: [real(r8) (:,:) ] layer node depth (m)
vegwp => canopystate_inst%vegwp_patch & ! Input: [real(r8) (:,:) ] vegetation water
! matric potential (mm)
)

Expand All @@ -308,10 +311,16 @@ subroutine Compute_EffecRootFrac_And_VertTranSink_HydStress( bounds, &
do pi = 1,max_patch_per_col
if (pi <= col%npatches(c)) then
p = col%patchi(c) + pi - 1
if (j == 1) then
qflx_hydr_redist_patch(p) = 0._r8
end if
if (patch%active(p).and.frac_veg_nosno(p)>0) then
if (patch%wtcol(p) > 0._r8) then
temp(c) = temp(c) + k_soil_root(p,j) &
* (smp(c,j) - vegwp(p,4) - grav2)* patch%wtcol(p)
patchflux = k_soil_root(p,j) * (smp(c,j) - vegwp(p,4) - grav2)
if (patchflux <0) then
qflx_hydr_redist_patch(p) = qflx_hydr_redist_patch(p) + patchflux
end if
temp(c) = temp(c) + patchflux * patch%wtcol(p)
endif
end if
end if
Expand Down
Loading