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.

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

parameter file updates #659

Merged
merged 11 commits into from
Jun 3, 2020
56 changes: 23 additions & 33 deletions biogeophys/FatesPlantRespPhotosynthMod.F90
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

! CIME Globals
use shr_log_mod , only : errMsg => shr_log_errMsg
Expand Down Expand Up @@ -161,8 +161,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
real(r8) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa)
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) :: btran_eff ! effective transpiration wetness factor (0 to 1)
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,18 +243,14 @@ 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( &
stomatal_intercept => EDPftvarcon_inst%stomatal_intercept, &
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?


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

do s = 1,nsites

! Multi-layer parameters scaled by leaf nitrogen profile.
Expand Down Expand Up @@ -393,10 +389,10 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
(hlm_use_planthydro.eq.itrue) .or. &
(nleafage > 1) .or. &
(hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then
if (hlm_use_planthydro.eq.itrue ) then
bbb = max( cf/rsmax0, bbbopt(nint(c3psn(ft)))*currentCohort%co_hydr%btran )

if (hlm_use_planthydro.eq.itrue ) then

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
Expand All @@ -412,8 +408,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
cumulative_lai = lai_canopy_above + lai_layers_above + 0.5*lai_current

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 +514,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 +818,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 +873,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 +916,6 @@ 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)



! Parameters
! ------------------------------------------------------------------------
Expand Down Expand Up @@ -952,15 +945,12 @@ 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
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
else
Expand All @@ -978,7 +968,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,13 +1084,13 @@ 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
! With an <= 0, then gs_mol = stomatal_intercept_btran

leaf_co2_ppress = can_co2_ppress- 1.4_r8/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 + &
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)
Expand All @@ -1121,9 +1111,9 @@ 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
Expand Down Expand Up @@ -1169,7 +1159,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in

! 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
gs_mol_err = bb_slope(ft)*max(anet, 0._r8)*hs/leaf_co2_ppress*can_press + stomatal_intercept_btran

if (abs(gs_mol-gs_mol_err) > 1.e-01_r8) then
write (fates_log(),*) 'CF: Ball-Berry error check - stomatal conductance error:'
Expand All @@ -1190,7 +1180,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
44 changes: 18 additions & 26 deletions main/EDParamsMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,6 @@ module EDParamsMod
real(r8),protected, public :: ED_val_understorey_death
real(r8),protected, public :: ED_val_cwd_fcel
real(r8),protected, public :: ED_val_cwd_flig
real(r8),protected, public :: ED_val_bbopt_c3
real(r8),protected, public :: ED_val_bbopt_c4
real(r8),protected, public :: ED_val_base_mr_20
real(r8),protected, public :: ED_val_phen_drought_threshold
real(r8),protected, public :: ED_val_phen_doff_time
Expand All @@ -44,6 +42,7 @@ module EDParamsMod
real(r8),protected, public :: ED_val_cohort_age_fusion_tol
real(r8),protected, public :: ED_val_patch_fusion_tol
real(r8),protected, public :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry
integer,protected, public :: stomatal_model !switch for choosing between stomatal conductance models, 1 for Ball-Berry, 2 for Medlyn


logical,protected, public :: active_crown_fire ! flag, 1=active crown fire 0=no active crown fire
Expand All @@ -68,8 +67,6 @@ module EDParamsMod
character(len=param_string_length),parameter,public :: ED_name_understorey_death = "fates_mort_understorey_death"
character(len=param_string_length),parameter,public :: ED_name_cwd_fcel= "fates_cwd_fcel"
character(len=param_string_length),parameter,public :: ED_name_cwd_flig= "fates_cwd_flig"
character(len=param_string_length),parameter,public :: ED_name_bbopt_c3= "fates_bbopt_c3"
character(len=param_string_length),parameter,public :: ED_name_bbopt_c4= "fates_bbopt_c4"
character(len=param_string_length),parameter,public :: ED_name_base_mr_20= "fates_base_mr_20"
character(len=param_string_length),parameter,public :: ED_name_phen_drought_threshold= "fates_phen_drought_threshold"
character(len=param_string_length),parameter,public :: ED_name_phen_doff_time= "fates_phen_doff_time"
Expand All @@ -84,6 +81,7 @@ module EDParamsMod
character(len=param_string_length),parameter,public :: ED_name_cohort_age_fusion_tol = "fates_cohort_age_fusion_tol"
character(len=param_string_length),parameter,public :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol"
character(len=param_string_length),parameter,public :: ED_name_canopy_closure_thresh= "fates_canopy_closure_thresh"
character(len=param_string_length),parameter,public :: ED_name_stomatal_model= "fates_leaf_stomatal_model"

! Resistance to active crown fire

Expand Down Expand Up @@ -171,8 +169,6 @@ subroutine FatesParamsInit()
ED_val_understorey_death = nan
ED_val_cwd_fcel = nan
ED_val_cwd_flig = nan
ED_val_bbopt_c3 = nan
ED_val_bbopt_c4 = nan
ED_val_base_mr_20 = nan
ED_val_phen_drought_threshold = nan
ED_val_phen_doff_time = nan
Expand All @@ -186,7 +182,8 @@ subroutine FatesParamsInit()
ED_val_cohort_size_fusion_tol = nan
ED_val_cohort_age_fusion_tol = nan
ED_val_patch_fusion_tol = nan
ED_val_canopy_closure_thresh = nan
ED_val_canopy_closure_thresh = nan
stomatal_model = -9
hydr_kmax_rsurf1 = nan
hydr_kmax_rsurf2 = nan

Expand Down Expand Up @@ -230,6 +227,7 @@ subroutine FatesRegisterParams(fates_params)
character(len=param_string_length), parameter :: dim_names_height(1) = (/dimension_name_history_height_bins/)
character(len=param_string_length), parameter :: dim_names_coageclass(1) = (/dimension_name_history_coage_bins/)


call FatesParamsInit()

call fates_params%RegisterParameter(name=ED_name_mort_disturb_frac, dimension_shape=dimension_shape_scalar, &
Expand All @@ -253,12 +251,6 @@ subroutine FatesRegisterParams(fates_params)
call fates_params%RegisterParameter(name=ED_name_cwd_flig, dimension_shape=dimension_shape_scalar, &
dimension_names=dim_names_scalar)

call fates_params%RegisterParameter(name=ED_name_bbopt_c3, dimension_shape=dimension_shape_scalar, &
dimension_names=dim_names_scalar)

call fates_params%RegisterParameter(name=ED_name_bbopt_c4, dimension_shape=dimension_shape_scalar, &
dimension_names=dim_names_scalar)

call fates_params%RegisterParameter(name=ED_name_base_mr_20, dimension_shape=dimension_shape_scalar, &
dimension_names=dim_names_scalar)

Expand Down Expand Up @@ -300,6 +292,9 @@ subroutine FatesRegisterParams(fates_params)

call fates_params%RegisterParameter(name=ED_name_canopy_closure_thresh, dimension_shape=dimension_shape_scalar, &
dimension_names=dim_names_scalar)

call fates_params%RegisterParameter(name=ED_name_stomatal_model, dimension_shape=dimension_shape_scalar, &
dimension_names=dim_names_scalar)

call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf1, dimension_shape=dimension_shape_scalar, &
dimension_names=dim_names_scalar)
Expand Down Expand Up @@ -378,8 +373,8 @@ subroutine FatesReceiveParams(fates_params)

class(fates_parameters_type), intent(inout) :: fates_params

real(r8) :: active_crown_fire_real !Local temp to transfer real data in file

real(r8) :: tmpreal ! local real variable for changing type on read
call fates_params%RetreiveParameter(name=ED_name_mort_disturb_frac, &
data=fates_mortality_disturbance_fraction)

Expand All @@ -401,12 +396,6 @@ subroutine FatesReceiveParams(fates_params)
call fates_params%RetreiveParameter(name=ED_name_cwd_flig, &
data=ED_val_cwd_flig)

call fates_params%RetreiveParameter(name=ED_name_bbopt_c3, &
data=ED_val_bbopt_c3)

call fates_params%RetreiveParameter(name=ED_name_bbopt_c4, &
data=ED_val_bbopt_c4)

call fates_params%RetreiveParameter(name=ED_name_base_mr_20, &
data=ED_val_base_mr_20)

Expand Down Expand Up @@ -449,6 +438,10 @@ subroutine FatesReceiveParams(fates_params)
call fates_params%RetreiveParameter(name=ED_name_canopy_closure_thresh, &
data=ED_val_canopy_closure_thresh)

call fates_params%RetreiveParameter(name=ED_name_stomatal_model, &
data=tmpreal)
stomatal_model = nint(tmpreal)

call fates_params%RetreiveParameter(name=hydr_name_kmax_rsurf1, &
data=hydr_kmax_rsurf1)

Expand Down Expand Up @@ -495,8 +488,8 @@ subroutine FatesReceiveParams(fates_params)
data=q10_froz)

call fates_params%RetreiveParameter(name=fates_name_active_crown_fire, &
data=active_crown_fire_real)
active_crown_fire = (abs(active_crown_fire_real-1.0_r8)<nearzero)
data=tmpreal)
active_crown_fire = (abs(tmpreal-1.0_r8)<nearzero)

call fates_params%RetreiveParameter(name=fates_name_cg_strikes, &
data=cg_strikes)
Expand Down Expand Up @@ -536,8 +529,6 @@ subroutine FatesReportParams(is_master)
write(fates_log(),fmt0) 'ED_val_understorey_death = ',ED_val_understorey_death
write(fates_log(),fmt0) 'ED_val_cwd_fcel = ',ED_val_cwd_fcel
write(fates_log(),fmt0) 'ED_val_cwd_flig = ',ED_val_cwd_flig
write(fates_log(),fmt0) 'ED_val_bbopt_c3 = ',ED_val_bbopt_c3
write(fates_log(),fmt0) 'ED_val_bbopt_c4 = ',ED_val_bbopt_c4
write(fates_log(),fmt0) 'ED_val_base_mr_20 = ', ED_val_base_mr_20
write(fates_log(),fmt0) 'ED_val_phen_drought_threshold = ',ED_val_phen_drought_threshold
write(fates_log(),fmt0) 'ED_val_phen_doff_time = ',ED_val_phen_doff_time
Expand All @@ -551,7 +542,8 @@ subroutine FatesReportParams(is_master)
write(fates_log(),fmt0) 'ED_val_cohort_size_fusion_tol = ',ED_val_cohort_size_fusion_tol
write(fates_log(),fmt0) 'ED_val_cohort_age_fusion_tol = ',ED_val_cohort_age_fusion_tol
write(fates_log(),fmt0) 'ED_val_patch_fusion_tol = ',ED_val_patch_fusion_tol
write(fates_log(),fmt0) 'ED_val_canopy_closure_thresh = ',ED_val_canopy_closure_thresh
write(fates_log(),fmt0) 'ED_val_canopy_closure_thresh = ',ED_val_canopy_closure_thresh
write(fates_log(),fmt0) 'stomatal_model = ',stomatal_model
write(fates_log(),fmt0) 'hydr_kmax_rsurf1 = ',hydr_kmax_rsurf1
write(fates_log(),fmt0) 'hydr_kmax_rsurf2 = ',hydr_kmax_rsurf2
write(fates_log(),fmt0) 'hydr_psi0 = ',hydr_psi0
Expand Down
Loading