From 113a3cac59f199fe592b1d9c4f7e255b6144771d Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 25 Apr 2019 16:50:01 -0700 Subject: [PATCH 1/4] Recalculated bbb (cuticular conductance) so that the same value is used through photosynth code, and is easier to follow) --- biogeophys/FatesPlantRespPhotosynthMod.F90 | 78 +++++++++------------- 1 file changed, 30 insertions(+), 48 deletions(-) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 65fcb9e677..5225824bba 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -394,10 +394,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 (bbbopt(nint(c3psn(ft)))*currentCohort%co_hydr%btran(1), 1._r8) + + if (hlm_use_planthydro.eq.itrue ) then + + bbb = max( cf/rsmax0, bbbopt(nint(c3psn(ft)))*currentCohort%co_hydr%btran(1) ) btran_eff = currentCohort%co_hydr%btran(1) ! dinc_ed is the total vegetation area index of each "leaf" layer @@ -414,7 +414,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) else - bbb = max (bbbopt(nint(c3psn(ft)))*currentPatch%btran_ft(ft), 1._r8) + bbb = max( cf/rsmax0, bbbopt(nint(c3psn(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 @@ -943,15 +943,9 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! empirical curvature parameter for ap photosynthesis co-limitation real(r8),parameter :: theta_ip = 0.999_r8 - real(r8), dimension(0:1) :: bbbopt !cuticular conductance - associate( bb_slope => EDPftvarcon_inst%BB_slope) ! slope of BB relationship - - bbbopt(0) = ED_val_bbopt_c4 - bbbopt(1) = ED_val_bbopt_c3 - ! photosynthetic pathway: 0. = c4, 1. = c3 c3c4_path_index = nint(EDPftvarcon_inst%c3psn(ft)) @@ -968,25 +962,22 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in anet_av_out = -lmr psn_out = 0._r8 - if(btran>0._r8) then - rstoma_out = min(rsmax0, cf*1._r8/(bbbopt(c3c4_path_index)*btran)) - else - rstoma_out = rsmax0 - endif + + ! The cuticular conductance already factored in maximum resistance as a bound + ! no need to re-bound it + + rstoma_out = cf/bbb c13disc_z = 0.0_r8 !carbon 13 discrimination in night time carbon flux, note value of 1.0 is used in CLM else ! day time (a little bit more complicated ...) -! if ( debug ) write(fates_log(),*) 'EDphot 594 ',laisun_lsl -! if ( debug ) write(fates_log(),*) 'EDphot 595 ',laisha_lsl + !is there leaf area? - (NV can be larger than 0 with only stem area if deciduous) + if ( laisun_lsl + laisha_lsl > 0._r8 ) then - !is there leaf area? - (NV can be larger than 0 with only stem area if deciduous) - if ( laisun_lsl + laisha_lsl > 0._r8 ) then - if(bbb > 1.0_r8)then !only if stomata open larger than cuticular conductance + !only if stomata open larger than cuticular conductance + if(bbb > 1.0_r8)then !only if stomata open larger than cuticular conductance -! if ( debug ) write(fates_log(),*) '600 in laisun, laisha loop ' - !Loop aroun shaded and unshaded leaves psn_out = 0._r8 ! psn is accumulated across sun and shaded leaves. rstoma_out = 0._r8 ! 1/rs is accumulated across sun and shaded leaves. @@ -1136,19 +1127,16 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! Convert gs_mol (umol /m**2/s) to gs (m/s) and then to rs (s/m) gs = gs_mol / cf - ! estimate carbon 13 discrimination in leaf level carbon flux Liang WEI and Hang ZHOU 2018, based on + ! estimate carbon 13 discrimination in leaf level carbon + ! flux Liang WEI and Hang ZHOU 2018, based on ! Ubierna and Farquhar, 2014 doi:10.1111/pce.12346, using the simplified model: ! $\Delta ^{13} C = \alpha_s + (b - \alpha_s) \cdot \frac{C_i}{C_a}$ ! just hard code b and \alpha_s for now, might move to parameter set in future ! b = 27.0 alpha_s = 4.4 ! TODO, not considering C4 or CAM right now, may need to address this - ! note co2_inter_c is intracelluar CO2, not intercelluar - c13disc_z = 4.4_r8 + (27.0_r8 - 4.4_r8) * min (can_co2_ppress, max (co2_inter_c, 0._r8)) / can_co2_ppress - - -! if ( debug ) write(fates_log(),*) 'EDPhoto 737 ', psn_out -! if ( debug ) write(fates_log(),*) 'EDPhoto 738 ', agross -! if ( debug ) write(fates_log(),*) 'EDPhoto 739 ', f_sun_lsl + ! note co2_inter_c is intracelluar CO2, not intercelluar + c13disc_z = 4.4_r8 + (27.0_r8 - 4.4_r8) * & + min (can_co2_ppress, max (co2_inter_c, 0._r8)) / can_co2_ppress ! Accumulate total photosynthesis umol/m2 ground/s-1. ! weight per unit sun and sha leaves. @@ -1162,10 +1150,6 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in gstoma = gstoma + & 1._r8/(min(1._r8/gs, rsmax0)) * (1.0_r8-f_sun_lsl) end if - -! if ( debug ) write(fates_log(),*) 'EDPhoto 758 ', psn_out -! if ( debug ) write(fates_log(),*) 'EDPhoto 759 ', agross -! if ( debug ) write(fates_log(),*) 'EDPhoto 760 ', f_sun_lsl ! Make sure iterative solution is correct if (gs_mol < 0._r8) then @@ -1188,25 +1172,23 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! This is the stomatal resistance of the leaf layer rstoma_out = 1._r8/gstoma - else !!set the situations with only cuticular conductance and no photosynthesis - - psn_out = 0._r8 - anet_av_out = -lmr - if(btran>0._r8) then - rstoma_out = min(rsmax0, cf*1._r8/(bbbopt(c3c4_path_index)*btran)) - else - rstoma_out = rsmax0 - endif - endif + else !!set the situations with only cuticular conductance and no photosynthesis + + psn_out = 0._r8 + anet_av_out = -lmr + rstoma_out = cf/bbb + + endif ! If stomata are open more than cuticular conductance else + !No leaf area. This layer is present only because of stems. ! (leaves are off, or have reduced to 0) psn_out = 0._r8 - rstoma_out = min(rsmax0, cf*1._r8/(0.1_r8*bbbopt(c3c4_path_index))) !assume stem loss is only 10% of leaf culticular - c13disc_z = 0.0_r8 + rstoma_out = min(rsmax0, cf/(0.1_r8*bbbopt(c3c4_path_index))) !assume stem loss is only 10% of leaf culticular + c13disc_z = 0.0_r8 - end if !is there leaf area? + end if !is there leaf area? end if ! night or day From 3e71640a346e91f3618a51659babd03966a503f4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 25 Apr 2019 16:59:11 -0700 Subject: [PATCH 2/4] Removed condition on bbb per discussion with C. Xu. Added a parameter for stem leakage. --- biogeophys/FatesPlantRespPhotosynthMod.F90 | 25 ++++++++++------------ 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 5225824bba..fe722f854c 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -924,6 +924,10 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! Fraction of light absorbed by non-photosynthetic pigments real(r8),parameter :: fnps = 0.15_r8 + ! For plants with no leaves, a miniscule amount of conductance + ! can happen through the stems, at a partial rate of cuticular conductance + real(r8),parameter :: stem_cuticle_loss_frac = 0.1_r8 + ! empirical curvature parameter for electron transport rate real(r8),parameter :: theta_psii = 0.7_r8 @@ -975,9 +979,6 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in !is there leaf area? - (NV can be larger than 0 with only stem area if deciduous) if ( laisun_lsl + laisha_lsl > 0._r8 ) then - !only if stomata open larger than cuticular conductance - if(bbb > 1.0_r8)then !only if stomata open larger than cuticular conductance - !Loop aroun shaded and unshaded leaves psn_out = 0._r8 ! psn is accumulated across sun and shaded leaves. rstoma_out = 0._r8 ! 1/rs is accumulated across sun and shaded leaves. @@ -1172,20 +1173,16 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! This is the stomatal resistance of the leaf layer rstoma_out = 1._r8/gstoma - else !!set the situations with only cuticular conductance and no photosynthesis - - psn_out = 0._r8 - anet_av_out = -lmr - rstoma_out = cf/bbb - - endif ! If stomata are open more than cuticular conductance - else - !No leaf area. This layer is present only because of stems. + ! No leaf area. This layer is present only because of stems. + ! Net assimilation is zero, not negative because there are + ! no leaves to even respire ! (leaves are off, or have reduced to 0) - psn_out = 0._r8 - rstoma_out = min(rsmax0, cf/(0.1_r8*bbbopt(c3c4_path_index))) !assume stem loss is only 10% of leaf culticular + + psn_out = 0._r8 + anet_av_out = 0._r8 + rstoma_out = min(rsmax0, cf/(stem_cuticle_loss_frac*bbbopt(c3c4_path_index))) c13disc_z = 0.0_r8 end if !is there leaf area? From f311ef5faf133399aa4de2c4abacce4cd29720f3 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 25 Apr 2019 17:05:55 -0700 Subject: [PATCH 3/4] Temporarily setting hydr_kmax_rsurf2 to be a hard-wired parameter constant. Will add this to parameter file during next API change --- biogeophys/FatesPlantHydraulicsMod.F90 | 6 +++++- main/EDParamsMod.F90 | 20 ++++++++++---------- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index a45cc45c38..af20e8670b 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -41,7 +41,7 @@ module FatesPlantHydraulicsMod use FatesConstantsMod, only : g_per_kg use EDParamsMod , only : hydr_kmax_rsurf1 - use EDParamsMod , only : hydr_kmax_rsurf2 + ! use EDParamsMod , only : hydr_kmax_rsurf2 use EDTypesMod , only : ed_site_type use EDTypesMod , only : ed_patch_type @@ -886,6 +886,10 @@ subroutine UpdateWaterDepTreeHydrCond(currentSite,ccohort,nlevsoi_hyd,bc_in) ! which is equiv to [kg m-1 s-1 MPa-1] real(r8) :: kmax_root_surf ! maximum conducitivity for unit root surface (kg water/m2 root area/Mpa/s) + ! (RGK 4-2019) THE FOLLOWING SHOULD BE ADDED TO THE PARAMETER FILE IN NEXT GROUPING + real(r8), parameter :: hydr_kmax_rsurf2 = 0.0001_r8 ! kg water/m2 root area/Mpa/s + + ccohort_hydr => ccohort%co_hydr csite_hydr => currentSite%si_hydr k = 1 !only for the first soil shell diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 767f4c2b33..8e423f9e22 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -80,10 +80,10 @@ module EDParamsMod ! Hydraulics Control Parameters (ONLY RELEVANT WHEN USE_FATES_HYDR = TRUE) ! ---------------------------------------------------------------------------------------------- real(r8),protected :: hydr_kmax_rsurf1 ! maximum conducitivity for unit root surface for water uptake(kg water/m2 root area/Mpa/s) - real(r8),protected :: hydr_kmax_rsurf2 ! maximum conducitivity for unit root surface for water loss (kg water/m2 root area/Mpa/s) +! real(r8),protected :: hydr_kmax_rsurf2 ! maximum conducitivity for unit root surface for water loss (kg water/m2 root area/Mpa/s) - character(len=param_string_length),parameter :: hydr_name_kmax_rsurf1 = "fates_hydr_kmax_rsurf1" - character(len=param_string_length),parameter :: hydr_name_kmax_rsurf2 = "fates_hydr_kmax_rsurf2" + character(len=param_string_length),parameter :: hydr_name_kmax_rsurf1 = "fates_hydr_kmax_rsurf" +! character(len=param_string_length),parameter :: hydr_name_kmax_rsurf2 = "fates_hydr_kmax_rsurf2" real(r8),protected :: hydr_psi0 ! sapwood water potential at saturation (MPa) character(len=param_string_length),parameter :: hydr_name_psi0 = "fates_hydr_psi0" @@ -161,7 +161,7 @@ subroutine FatesParamsInit() ED_val_canopy_closure_thresh = nan hydr_kmax_rsurf1 = nan - hydr_kmax_rsurf2 = nan +! hydr_kmax_rsurf2 = nan hydr_psi0 = nan hydr_psicap = nan @@ -266,8 +266,8 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf1, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) - call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf2, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names) +! call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf2, dimension_shape=dimension_shape_1d, & +! dimension_names=dim_names) call fates_params%RegisterParameter(name=hydr_name_psi0, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) @@ -390,8 +390,8 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=hydr_name_kmax_rsurf1, & data=hydr_kmax_rsurf1) - call fates_params%RetreiveParameter(name=hydr_name_kmax_rsurf2, & - data=hydr_kmax_rsurf2) +! call fates_params%RetreiveParameter(name=hydr_name_kmax_rsurf2, & +! data=hydr_kmax_rsurf2) call fates_params%RetreiveParameter(name=hydr_name_psi0, & data=hydr_psi0) @@ -470,8 +470,8 @@ subroutine FatesReportParams(is_master) write(fates_log(),fmt0) 'ED_val_cohort_fusion_tol = ',ED_val_cohort_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) 'hydr_kmax_rsurf1 = ',hydr_kmax_rsurf1 - write(fates_log(),fmt0) 'hydr_kmax_rsurf2 = ',hydr_kmax_rsurf2 + 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 write(fates_log(),fmt0) 'hydr_psicap = ',hydr_psicap write(fates_log(),fmt0) 'bgc_soil_salinity = ', bgc_soil_salinity From b4d643e34478fdc1a30529f1400db1e58fa6c3cf Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 25 Apr 2019 17:32:01 -0700 Subject: [PATCH 4/4] Small fix on declaring bbbopt as a local --- biogeophys/FatesPlantRespPhotosynthMod.F90 | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index fe722f854c..0baf150f4e 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -45,6 +45,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 ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -96,7 +97,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) use FatesConstantsMod, only : rgas => rgas_J_K_kmol use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use FatesParameterDerivedMod, only : param_derived - use EDParamsMod, only : ED_val_bbopt_c3, ED_val_bbopt_c4, ED_val_base_mr_20 + use FatesAllometryMod, only : bleaf use FatesAllometryMod, only : storage_fraction_of_target use FatesAllometryMod, only : set_root_fraction @@ -844,7 +845,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! ------------------------------------------------------------------------------------ use EDPftvarcon , only : EDPftvarcon_inst - use EDParamsMod, only : ED_val_bbopt_c3, ED_val_bbopt_c4 + ! Arguments ! ------------------------------------------------------------------------------------ @@ -919,6 +920,10 @@ 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 ! ------------------------------------------------------------------------ ! Fraction of light absorbed by non-photosynthetic pigments @@ -953,6 +958,9 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! 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