Skip to content

Commit

Permalink
Merge remote branch 'djk/hydraulic_redistribution' for PR#1126
Browse files Browse the repository at this point in the history
  • Loading branch information
@negin513
negin513 committed Oct 14, 2020
2 parents 09dc996 + 241b788 commit 8a76816
Show file tree
Hide file tree
Showing 4 changed files with 84 additions and 21 deletions.
13 changes: 12 additions & 1 deletion src/biogeophys/CanopyStateType.F90
View file
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
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
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

0 comments on commit 8a76816

Please sign in to comment.