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

Adding Medlyn stomatal conductance model into CTSM-FATES #608

Merged
merged 17 commits into from
Jun 21, 2020
Merged

Adding Medlyn stomatal conductance model into CTSM-FATES #608

Show file tree
Hide file tree
Changes from 11 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
3fa3292
Added Medlyn stomatal model, first pass through.
Feb 5, 2020
b1ddca4
Added Medlyn stomatal model, first pass through.
Feb 5, 2020
786c764
Added Medlyn stomatal model, first pass through.
Feb 9, 2020
0957319
Add switch for choosing between Medlyn stomatal model and Ball-Berry …
May 8, 2020
565fd71
update parameter file
May 10, 2020
9ce3bc5
Update parameter file
May 11, 2020
927170e
Revert the parameter file
May 13, 2020
1639e58
Update the parameter file
May 13, 2020
00b06c4
Merge branch 'master' of https://github.com/NGEET/fates into fates_ma…
May 14, 2020
b51ffe0
Added the logging of Medlyn model error
May 15, 2020
bf8d5f6
Small fix about the logging of Medlyn model error
May 15, 2020
b636984
small edits on parameter description
Jun 2, 2020
dc7c257
Conflict Resolutions between medlyn PR and parameter file updates
rgknox Jun 5, 2020
a57c17b
Resolved conflicts with photosynthesis code
rgknox Jun 5, 2020
15d8efd
removed redundant declaration of stomatal_intercept
rgknox Jun 5, 2020
5e9e1ee
bug fix on setting stomatal model as uninitialized
rgknox Jun 5, 2020
7f84a4b
Merge pull request #1 from rgknox/fates_master_medlyn
Qianyuxuan Jun 8, 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
106 changes: 68 additions & 38 deletions biogeophys/FatesPlantRespPhotosynthMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ module FATESPlantRespPhotosynthMod
use PRTGenericMod, only : store_organ
use PRTGenericMod, only : repro_organ
use PRTGenericMod, only : struct_organ
use EDParamsMod, only : ED_val_bbopt_c3, ED_val_bbopt_c4, ED_val_base_mr_20
use EDParamsMod, only : ED_val_base_mr_20, stomatal_model

! CIME Globals
use shr_log_mod , only : errMsg => shr_log_errMsg
Expand All @@ -64,6 +64,13 @@ module FATESPlantRespPhotosynthMod
real(r8),parameter :: rsmax0 = 2.e8_r8

logical :: debug = .false.
!-------------------------------------------------------------------------------------

! Ratio of H2O/CO2 gas diffusion in stomatal airspace (approximate)
real(r8),parameter :: h2o_co2_stoma_diffuse_ratio = 1.6_r8

! Ratio of H2O/CO2 gass diffusion in the leaf boundary layer (approximate)
real(r8),parameter :: h2o_co2_bl_diffuse_ratio = 1.4_r8

contains

Expand Down Expand Up @@ -162,7 +169,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
real(r8) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa)
real(r8) :: co2_cpoint ! CO2 compensation point (Pa)
real(r8) :: btran_eff ! effective transpiration wetness factor (0 to 1)
real(r8) :: bbb ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s)
real(r8) :: stomatal_intercept_btran ! water-stressed minimum stomatal conductance (umol H2O/m**2/s)
real(r8) :: kn ! leaf nitrogen decay coefficient
real(r8) :: cf ! s m**2/umol -> s/m (ideal gas conversion) [umol/m3]
real(r8) :: gb_mol ! leaf boundary layer conductance (molar form: [umol /m**2/s])
Expand Down Expand Up @@ -243,17 +250,17 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
! Ball-Berry minimum leaf conductance, unstressed (umol H2O/m**2/s)
! For C3 and C4 plants
! -----------------------------------------------------------------------------------
real(r8), dimension(0:1) :: bbbopt


associate( &
c3psn => EDPftvarcon_inst%c3psn , &
slatop => EDPftvarcon_inst%slatop , & ! specific leaf area at top of canopy,
! projected area basis [m^2/gC]
woody => EDPftvarcon_inst%woody) ! Is vegetation woody or not?
woody => EDPftvarcon_inst%woody , & ! Is vegetation woody or not?
stomatal_intercept => EDPftvarcon_inst%stomatal_intercept ) !Unstressed minimum stomatal conductance



bbbopt(0) = ED_val_bbopt_c4
bbbopt(1) = ED_val_bbopt_c3

do s = 1,nsites

Expand Down Expand Up @@ -396,8 +403,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)

if (hlm_use_planthydro.eq.itrue ) then

bbb = max( cf/rsmax0, bbbopt(nint(c3psn(ft)))*currentCohort%co_hydr%btran )
btran_eff = currentCohort%co_hydr%btran
stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentCohort%co_hydr%btran )
btran_eff = currentCohort%co_hydr%btran

! dinc_ed is the total vegetation area index of each "leaf" layer
! we convert to the leaf only portion of the increment
Expand All @@ -413,7 +420,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)

else

bbb = max( cf/rsmax0, bbbopt(nint(c3psn(ft)))*currentPatch%btran_ft(ft) )
stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentPatch%btran_ft(ft) )
btran_eff = currentPatch%btran_ft(ft)
! For consistency sake, we use total LAI here, and not exposed
! if the plant is under-snow, it will be effectively dormant for
Expand Down Expand Up @@ -518,7 +525,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
bc_in(s)%cair_pa(ifp), & ! in
bc_in(s)%oair_pa(ifp), & ! in
btran_eff, & ! in
bbb, & ! in
stomatal_intercept_btran, & ! in
cf, & ! in
gb_mol, & ! in
ceair, & ! in
Expand Down Expand Up @@ -822,7 +829,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
can_co2_ppress, & ! in
can_o2_ppress, & ! in
btran, & ! in
bbb, & ! in
stomatal_intercept_btran, & ! in
cf, & ! in
gb_mol, & ! in
ceair, & ! in
Expand Down Expand Up @@ -877,7 +884,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
real(r8), intent(in) :: can_co2_ppress ! Partial pressure of CO2 NEAR the leaf surface (Pa)
real(r8), intent(in) :: can_o2_ppress ! Partial pressure of O2 NEAR the leaf surface (Pa)
real(r8), intent(in) :: btran ! transpiration wetness factor (0 to 1)
real(r8), intent(in) :: bbb ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s)
real(r8), intent(in) :: stomatal_intercept_btran !water-stressed minimum stomatal conductance (umol H2O/m**2/s)
real(r8), intent(in) :: cf ! s m**2/umol -> s/m (ideal gas conversion) [umol/m3]
real(r8), intent(in) :: gb_mol ! leaf boundary layer conductance (umol /m**2/s)
real(r8), intent(in) :: ceair ! vapor pressure of air, constrained (Pa)
Expand Down Expand Up @@ -920,9 +927,8 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
real(r8) :: leaf_co2_ppress ! CO2 partial pressure at leaf surface (Pa)
real(r8) :: init_co2_inter_c ! First guess intercellular co2 specific to C path

real(r8), dimension(0:1) :: bbbopt ! Cuticular conductance at full water potential (umol H2O /m2/s)


real(r8) :: term ! intermediate variable in Medlyn stomatal conductance model
real(r8) :: vpd ! water vapor deficit in Medlyn stomatal model (KPa)

! Parameters
! ------------------------------------------------------------------------
Expand Down Expand Up @@ -951,15 +957,17 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in

! empirical curvature parameter for ap photosynthesis co-limitation
real(r8),parameter :: theta_ip = 0.999_r8


associate( bb_slope => EDPftvarcon_inst%BB_slope) ! slope of BB relationship
associate( bb_slope => EDPftvarcon_inst%BB_slope ,& ! slope of BB relationship, unitless
medlyn_slope=> EDPftvarcon_inst%medlyn_slope , & ! Slope for Medlyn stomatal conductance model method, the unit is KPa^0.5
stomatal_intercept=> EDPftvarcon_inst%stomatal_intercept ) !Unstressed minimum stomatal conductance, the unit is umol/m**2/s



! photosynthetic pathway: 0. = c4, 1. = c3
c3c4_path_index = nint(EDPftvarcon_inst%c3psn(ft))

bbbopt(0) = ED_val_bbopt_c4
bbbopt(1) = ED_val_bbopt_c3

if (c3c4_path_index == 1) then
init_co2_inter_c = init_a2l_co2_c3 * can_co2_ppress
Expand All @@ -978,7 +986,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
! The cuticular conductance already factored in maximum resistance as a bound
! no need to re-bound it

rstoma_out = cf/bbb
rstoma_out = cf/stomatal_intercept_btran

c13disc_z = 0.0_r8 !carbon 13 discrimination in night time carbon flux, note value of 1.0 is used in CLM

Expand Down Expand Up @@ -1094,21 +1102,37 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
end if

! Quadratic gs_mol calculation with an known. Valid for an >= 0.
! With an <= 0, then gs_mol = bbb

leaf_co2_ppress = can_co2_ppress- 1.4_r8/gb_mol * anet * can_press
! With an <= 0, then gs_mol = stomatal_intercept_btran
leaf_co2_ppress = can_co2_ppress- h2o_co2_bl_diffuse_ratio/gb_mol * anet * can_press
leaf_co2_ppress = max(leaf_co2_ppress,1.e-06_r8)
aquad = leaf_co2_ppress
bquad = leaf_co2_ppress*(gb_mol - bbb) - bb_slope(ft) * anet * can_press
cquad = -gb_mol*(leaf_co2_ppress*bbb + &
if ( stomatal_model == 2 ) then
!stomatal conductance calculated from Medlyn et al. (2011), the numerical &
!implementation was adapted from the equations in CLM5.0
vpd = max((veg_esat - ceair), 50._r8) * 0.001_r8 !addapted from CLM5. Put some constraint on VPD
!when Medlyn stomatal conductance is being used, the unit is KPa. Ignoring the constraint will cause errors when model runs.
term = h2o_co2_stoma_diffuse_ratio * anet / (leaf_co2_ppress / can_press)
aquad = 1.0_r8
bquad = -(2.0 * (stomatal_intercept_btran+ term) + (medlyn_slope(ft) * term)**2 / &
(gb_mol * vpd ))
cquad = stomatal_intercept_btran*stomatal_intercept_btran + &
(2.0*stomatal_intercept_btran + term * &
(1.0 - medlyn_slope(ft)* medlyn_slope(ft) / vpd)) * term

call quadratic_f (aquad, bquad, cquad, r1, r2)
gs_mol = max(r1,r2)

else if ( stomatal_model == 1 ) then !stomatal conductance calculated from Ball et al. (1987)
aquad = leaf_co2_ppress
bquad = leaf_co2_ppress*(gb_mol - stomatal_intercept_btran) - bb_slope(ft) * anet * can_press
cquad = -gb_mol*(leaf_co2_ppress*stomatal_intercept_btran + &
bb_slope(ft)*anet*can_press * ceair/ veg_esat )

call quadratic_f (aquad, bquad, cquad, r1, r2)
gs_mol = max(r1,r2)

end if
! Derive new estimate for co2_inter_c
co2_inter_c = can_co2_ppress - anet * can_press * &
(1.4_r8*gs_mol+1.6_r8*gb_mol) / (gb_mol*gs_mol)
(h2o_co2_bl_diffuse_ratio*gs_mol+h2o_co2_stoma_diffuse_ratio*gb_mol) / (gb_mol*gs_mol)

! Check for co2_inter_c convergence. Delta co2_inter_c/pair = mol/mol.
! Multiply by 10**6 to convert to umol/mol (ppm). Exit iteration if
Expand All @@ -1121,17 +1145,18 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
end if
end do !iteration loop

! End of co2_inter_c iteration. Check for an < 0, in which case gs_mol = bbb
! End of co2_inter_c iteration. Check for an < 0, in which case
! gs_mol =stomatal_intercept_btran
if (anet < 0._r8) then
gs_mol = bbb
gs_mol = stomatal_intercept_btran
end if

! Final estimates for leaf_co2_ppress and co2_inter_c
! (needed for early exit of co2_inter_c iteration when an < 0)
leaf_co2_ppress = can_co2_ppress - 1.4_r8/gb_mol * anet * can_press
leaf_co2_ppress = can_co2_ppress - h2o_co2_bl_diffuse_ratio/gb_mol * anet * can_press
leaf_co2_ppress = max(leaf_co2_ppress,1.e-06_r8)
co2_inter_c = can_co2_ppress - anet * can_press * &
(1.4_r8*gs_mol+1.6_r8*gb_mol) / (gb_mol*gs_mol)
(h2o_co2_bl_diffuse_ratio*gs_mol+h2o_co2_stoma_diffuse_ratio*gb_mol) / (gb_mol*gs_mol)

! Convert gs_mol (umol /m**2/s) to gs (m/s) and then to rs (s/m)
gs = gs_mol / cf
Expand Down Expand Up @@ -1166,13 +1191,18 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
write (fates_log(),*)'gs_mol= ',gs_mol
call endrun(msg=errMsg(sourcefile, __LINE__))
end if

! Compare with Ball-Berry model: gs_mol = m * an * hs/leaf_co2_ppress p + b
hs = (gb_mol*ceair + gs_mol* veg_esat ) / ((gb_mol+gs_mol)*veg_esat )
gs_mol_err = bb_slope(ft)*max(anet, 0._r8)*hs/leaf_co2_ppress*can_press + bbb


! Compare with Medlyn model: gs_mol = 1.6*(1+m/sqrt(vpd)) * an/leaf_co2_ppress*p + b
if ( stomatal_model == 2 ) then
gs_mol_err = h2o_co2_stoma_diffuse_ratio*(1 + medlyn_slope(ft)/sqrt(vpd))*max(anet,0._r8)/leaf_co2_ppress*can_press + stomatal_intercept_btran
! Compare with Ball-Berry model: gs_mol = m * an * hs/leaf_co2_ppress*p + b
else if ( stomatal_model == 1 ) then
hs = (gb_mol*ceair + gs_mol* veg_esat ) / ((gb_mol+gs_mol)*veg_esat )
gs_mol_err = bb_slope(ft)*max(anet, 0._r8)*hs/leaf_co2_ppress*can_press + stomatal_intercept_btran
end if

if (abs(gs_mol-gs_mol_err) > 1.e-01_r8) then
write (fates_log(),*) 'CF: Ball-Berry error check - stomatal conductance error:'
write (fates_log(),*) 'Stomatal model error check - stomatal conductance error:'
write (fates_log(),*) gs_mol, gs_mol_err
end if

Expand All @@ -1190,7 +1220,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in

psn_out = 0._r8
anet_av_out = 0._r8
rstoma_out = min(rsmax0, cf/(stem_cuticle_loss_frac*bbbopt(c3c4_path_index)))
rstoma_out = min(rsmax0,cf/(stem_cuticle_loss_frac*stomatal_intercept(ft)))
c13disc_z = 0.0_r8

end if !is there leaf area?
Expand Down
Loading