diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 3089b6df8d..570b14d506 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -47,7 +47,7 @@ module EDPhysiologyMod use EDCohortDynamicsMod , only : InitPRTObject use FatesAllometryMod , only : tree_lai_sai use FatesAllometryMod , only : leafc_from_treelai - use FatesAllometryMod , only : decay_coeff_vcmax + use LeafBiophysicsMod , only : DecayCoeffVcmax use FatesLitterMod , only : litter_type use EDTypesMod , only : site_massbal_type use EDTypesMod , only : numlevsoil_max @@ -772,7 +772,7 @@ subroutine trim_canopy( currentSite ) ! Calculate sla_levleaf following the sla profile with overlying leaf area ! Scale for leaf nitrogen profile - kn = decay_coeff_vcmax(currentCohort%vcmax25top, & + kn = DecayCoeffVcmax(currentCohort%vcmax25top, & prt_params%leafn_vert_scaler_coeff1(ipft), & prt_params%leafn_vert_scaler_coeff2(ipft)) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 32fee8e39f..28cd59dde3 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -100,7 +100,8 @@ module FatesAllometryMod use FatesGlobals , only : FatesWarn,N2S,A2S,I2S use EDParamsMod , only : nlevleaf,dinc_vai,dlower_vai use DamageMainMod , only : GetCrownReduction - + use LeafBiophysicsMod, only : DecayCoeffVcmax + implicit none private @@ -116,7 +117,6 @@ module FatesAllometryMod public :: bdead_allom ! Generic bdead wrapper public :: carea_allom ! Generic crown area wrapper public :: bstore_allom ! Generic maximum storage carbon wrapper - public :: decay_coeff_vcmax ! vertical canopy decay rate, scaled on vcmax public :: ForceDBH ! Method to set DBH to sync with structure ! or fineroot biomass public :: CheckIntegratedAllometries @@ -717,9 +717,9 @@ real(r8) function tree_lai( leaf_c, pft, c_area, nplant, cl, canopy_lai, vcmax25 end if ! Coefficient for exponential decay of 1/sla with canopy depth: - kn = decay_coeff_vcmax(vcmax25top, & - prt_params%leafn_vert_scaler_coeff1(pft), & - prt_params%leafn_vert_scaler_coeff2(pft)) + kn = DecayCoeffVcmax(vcmax25top, & + prt_params%leafn_vert_scaler_coeff1(pft), & + prt_params%leafn_vert_scaler_coeff2(pft)) ! take PFT-level maximum SLA value, even if under a thick canopy (which has units of m2/gC), ! and put into units of m2/kgC @@ -940,9 +940,9 @@ real(r8) function leafc_from_treelai( treelai, treesai, pft, c_area, nplant, cl, sla_max = g_per_kg * prt_params%slamax(pft) ! Coefficient for exponential decay of 1/sla with canopy depth: - kn = decay_coeff_vcmax(vcmax25top, & - prt_params%leafn_vert_scaler_coeff1(pft), & - prt_params%leafn_vert_scaler_coeff2(pft)) + kn = DecayCoeffVcmax(vcmax25top, & + prt_params%leafn_vert_scaler_coeff1(pft), & + prt_params%leafn_vert_scaler_coeff2(pft)) if(treelai > 0.0_r8)then ! Leafc_per_unitarea at which sla_max is reached due to exponential sla profile in canopy: @@ -2979,42 +2979,6 @@ end subroutine jackson_beta_root_profile ! ===================================================================================== - real(r8) function decay_coeff_vcmax(vcmax25top,slope_param,intercept_param) - - ! --------------------------------------------------------------------------------- - ! This function estimates the decay coefficient used to estimate vertical - ! attenuation of properties in the canopy. - ! - ! Decay coefficient (kn) is a function of vcmax25top for each pft. - ! - ! Currently, this decay is applied to vcmax attenuation, SLA (optionally) - ! and leaf respiration (optionally w/ Atkin) - ! - ! --------------------------------------------------------------------------------- - - !ARGUMENTS - - real(r8),intent(in) :: vcmax25top - real(r8),intent(in) :: slope_param ! multiplies vcmax25top - real(r8),intent(in) :: intercept_param ! adds to vcmax25top - - - !LOCAL VARIABLES - ! ----------------------------------------------------------------------------------- - - ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used - ! kn = 0.11. Here, we derive kn from vcmax25 as in Lloyd et al - ! (2010) Biogeosciences, 7, 1833-1859 - ! This function is also used to vertically scale leaf maintenance - ! respiration. - - decay_coeff_vcmax = exp(slope_param * vcmax25top - intercept_param) - - return - end function decay_coeff_vcmax - - ! ===================================================================================== - subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d, h, bdead, bl ) ! ========================================================================= diff --git a/biogeophys/FatesLeafBiophysParamsMod.F90 b/biogeophys/FatesLeafBiophysParamsMod.F90 new file mode 100644 index 0000000000..57ac9365b8 --- /dev/null +++ b/biogeophys/FatesLeafBiophysParamsMod.F90 @@ -0,0 +1,289 @@ +module FatesLeafBiophysParamsMod + + use FatesConstantsMod , only: r8 => fates_r8 + use FatesConstantsMod , only: fates_check_param_set + use FatesParametersInterface, only : param_string_length + use FatesGlobals, only : fates_log + use FatesGlobals, only : endrun => fates_endrun + use shr_log_mod , only : errMsg => shr_log_errMsg + use LeafBiophysicsMod, only : lb_params,btran_on_gs_gs1,btran_on_ag_none + use FatesParametersInterface, only : fates_parameters_type + ! Register the parameters we want the host to provide, and + ! indicate whether they are fates parameters or host parameters + ! that need to be synced with host values. + use FatesParametersInterface, only : fates_parameters_type, param_string_length + use FatesParametersInterface, only : dimension_name_pft, dimension_shape_1d + use FatesParametersInterface, only : dimension_shape_scalar, dimension_name_scalar + use FatesUtilsMod, only : ArrayNint + + implicit none + private ! Modules are private by default + save + + public :: LeafBiophysRegisterParams + public :: LeafBiophysReceiveParams + public :: LeafBiophysReportParams + + character(len=*), parameter :: sourcefile = & + __FILE__ + + integer, parameter :: lower_bound_pft = 1 + +contains + + ! ===================================================================================== + + subroutine LeafBiophysRegisterParams(fates_params) + + + class(fates_parameters_type), intent(inout) :: fates_params + + character(len=param_string_length), parameter :: dim_names(1) = (/dimension_name_pft/) + integer, parameter :: dim_lower_bound(1) = (/ lower_bound_pft /) + character(len=param_string_length) :: name + character(len=param_string_length), parameter :: dim_names_scalar(1) = (/dimension_name_scalar/) + + + ! Register scalars + + name = 'fates_daylength_factor_switch' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + + name = 'fates_leaf_stomatal_model' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + + name = 'fates_leaf_stomatal_assim_model' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + + name = 'fates_leaf_photo_tempsens_model' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + + ! Register PFT dimensioned + + name = 'fates_leaf_c3psn' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_stomatal_btran_model' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_agross_btran_model' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_stomatal_slope_ballberry' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_stomatal_slope_medlyn' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_stomatal_intercept' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_maintresp_reduction_curvature' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_maintresp_reduction_intercept' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_maintresp_reduction_upthresh' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_maintresp_leaf_atkin2017_baserate' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_maintresp_leaf_ryan1991_baserate' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_vcmaxha' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_jmaxha' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_vcmaxhd' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_jmaxhd' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_vcmaxse' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_leaf_jmaxse' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + return + end subroutine LeafBiophysRegisterParams + + ! ===================================================================================== + + subroutine LeafBiophysReceiveParams(fates_params) + + + class(fates_parameters_type), intent(inout) :: fates_params + real(r8), allocatable :: tmpreal(:) ! Temporary variable to hold floats + real(r8) :: tmpscalar + character(len=param_string_length) :: name + + name = 'fates_daylength_factor_switch' + call fates_params%RetrieveParameter(name=name, & + data=tmpscalar) + lb_params%dayl_switch = nint(tmpscalar) + + name = 'fates_leaf_stomatal_model' + call fates_params%RetrieveParameter(name=name, & + data=tmpscalar) + lb_params%stomatal_model = nint(tmpscalar) + + name = 'fates_leaf_stomatal_assim_model' + call fates_params%RetrieveParameter(name=name, & + data=tmpscalar) + lb_params%stomatal_assim_model = nint(tmpscalar) + + name = 'fates_leaf_photo_tempsens_model' + call fates_params%RetrieveParameter(name=name, & + data=tmpscalar) + lb_params%photo_tempsens_model = nint(tmpscalar) + + name = 'fates_leaf_c3psn' + call fates_params%RetrieveParameterAllocate(name=name, & + data=tmpreal) + allocate(lb_params%c3psn(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,lb_params%c3psn) + deallocate(tmpreal) + + name = 'fates_leaf_stomatal_btran_model' + call fates_params%RetrieveParameterAllocate(name=name, & + data=tmpreal) + allocate(lb_params%stomatal_btran_model(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,lb_params%stomatal_btran_model) + deallocate(tmpreal) + + name = 'fates_leaf_agross_btran_model' + call fates_params%RetrieveParameterAllocate(name=name, & + data=tmpreal) + allocate(lb_params%agross_btran_model(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,lb_params%agross_btran_model) + deallocate(tmpreal) + + name = 'fates_leaf_stomatal_slope_medlyn' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%medlyn_slope) + + name = 'fates_leaf_stomatal_slope_ballberry' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%bb_slope) + + name = 'fates_leaf_stomatal_intercept' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%stomatal_intercept) + + name = 'fates_maintresp_leaf_ryan1991_baserate' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%maintresp_leaf_ryan1991_baserate) + + name = 'fates_maintresp_leaf_atkin2017_baserate' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%maintresp_leaf_atkin2017_baserate) + + name = 'fates_maintresp_reduction_curvature' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%maintresp_reduction_curvature) + + name = 'fates_maintresp_reduction_intercept' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%maintresp_reduction_intercept) + + name = 'fates_maintresp_reduction_upthresh' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%maintresp_reduction_upthresh) + + name = 'fates_leaf_vcmaxha' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%vcmaxha) + + name = 'fates_leaf_jmaxha' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%jmaxha) + + name = 'fates_leaf_vcmaxhd' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%vcmaxhd) + + name = 'fates_leaf_jmaxhd' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%jmaxhd) + + name = 'fates_leaf_vcmaxse' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%vcmaxse) + + name = 'fates_leaf_jmaxse' + call fates_params%RetrieveParameterAllocate(name=name, & + data=lb_params%jmaxse) + + + return + end subroutine LeafBiophysReceiveParams + + ! ==================================================================================== + + subroutine LeafBiophysReportParams(is_master) + + ! Argument + logical, intent(in) :: is_master ! Only log if this is the master proc + + logical, parameter :: debug_report = .false. + character(len=32),parameter :: fmt_rout = '(a,F16.8)' + character(len=32),parameter :: fmt_iout = '(a,I8)' + + integer :: npft,ipft + + if(debug_report .and. is_master) then + write(fates_log(),fmt_iout) 'fates_leaf_c3psn = ',lb_params%c3psn + write(fates_log(),fmt_iout) 'fates_leaf_stomatal_btran_model = ',lb_params%stomatal_btran_model + write(fates_log(),fmt_iout) 'fates_leaf_agross_btran_model = ',lb_params%agross_btran_model + write(fates_log(),fmt_rout) 'fates_leaf_vcmaxha = ',lb_params%vcmaxha + write(fates_log(),fmt_rout) 'fates_leaf_jmaxha = ',lb_params%jmaxha + write(fates_log(),fmt_rout) 'fates_leaf_vcmaxhd = ',lb_params%vcmaxhd + write(fates_log(),fmt_rout) 'fates_leaf_jmaxhd = ',lb_params%jmaxhd + write(fates_log(),fmt_rout) 'fates_leaf_vcmaxse = ',lb_params%vcmaxse + write(fates_log(),fmt_rout) 'fates_leaf_jmaxse = ',lb_params%jmaxse + write(fates_log(),fmt_iout) 'fates_daylength_factor_switch = ',lb_params%dayl_switch + write(fates_log(),fmt_iout) 'fates_leaf_stomatal_model = ',lb_params%stomatal_model + write(fates_log(),fmt_iout) 'fates_leaf_stomatal_assim_model = ',lb_params%stomatal_assim_model + write(fates_log(),fmt_iout) 'fates_leaf_photo_tempsens_model = ',lb_params%photo_tempsens_model + write(fates_log(),fmt_rout) 'fates_leaf_stomatal_slope_medlyn = ',lb_params%medlyn_slope + write(fates_log(),fmt_rout) 'fates_leaf_stomatal_slope_ballberry = ',lb_params%bb_slope + write(fates_log(),fmt_rout) 'fates_leaf_stomatal_intercept = ',lb_params%stomatal_intercept + write(fates_log(),fmt_rout) 'fates_maintresp_leaf_ryan1991_baserate = ',lb_params%maintresp_leaf_ryan1991_baserate + write(fates_log(),fmt_rout) 'fates_maintresp_leaf_atkin2017_baserate = ',lb_params%maintresp_leaf_atkin2017_baserate + write(fates_log(),fmt_rout) 'fates_maintresp_reduction_curvature = ',lb_params%maintresp_reduction_curvature + write(fates_log(),fmt_rout) 'fates_maintresp_reduction_intercept = ',lb_params%maintresp_reduction_intercept + write(fates_log(),fmt_rout) 'fates_maintresp_reduction_upthresh = ',lb_params%maintresp_reduction_upthresh + end if + + + end subroutine LeafBiophysReportParams + +end module FatesLeafBiophysParamsMod diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index d1413f1b15..b898c93e0a 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -28,12 +28,12 @@ module FATESPlantRespPhotosynthMod use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : molar_mass_ratio_vapdry use FatesConstantsMod, only : molar_mass_water - use FatesConstantsMod, only : rgas_J_K_mol + use FatesConstantsMod, only : rgas => rgas_J_K_kmol use FatesConstantsMod, only : fates_unset_r8 use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm + use FatesConstantsMod, only : umol_per_kmol + use FatesConstantsMod, only : wm2_to_umolm2s use FatesConstantsMod, only : nocomp_bareground - use FatesConstantsMod, only : photosynth_acclim_model_none - use FatesConstantsMod, only : photosynth_acclim_model_kumarathunge_etal_2019 use FatesInterfaceTypesMod, only : hlm_use_planthydro use FatesInterfaceTypesMod, only : hlm_parteh_mode use FatesInterfaceTypesMod, only : numpft @@ -61,19 +61,25 @@ module FATESPlantRespPhotosynthMod use PRTGenericMod, only : repro_organ use PRTGenericMod, only : struct_organ use EDParamsMod, only : maintresp_nonleaf_baserate - use EDParamsMod, only : stomatal_model - use EDParamsMod, only : stomatal_assim_model - use EDParamsMod, only : dayl_switch - use EDParamsMod, only : photo_tempsens_model use PRTParametersMod, only : prt_params use EDPftvarcon , only : EDPftvarcon_inst - use TemperatureType, only : temperature_type use FatesRadiationMemMod, only : norman_solver,twostr_solver use EDParamsMod, only : radiation_model use FatesRadiationMemMod, only : ipar use FatesTwoStreamUtilsMod, only : FatesGetCohortAbsRad use FatesAllometryMod , only : VegAreaLayer - use FatesAllometryMod, only : decay_coeff_vcmax + + use LeafBiophysicsMod, only : LeafLayerPhotosynthesis + use LeafBiophysicsMod, only : LeafHumidityStomaResis + use LeafBiophysicsMod, only : GetCanopyGasParameters + use LeafBiophysicsMod, only : LeafLayerMaintenanceRespiration_Ryan_1991 + use LeafBiophysicsMod, only : LeafLayerMaintenanceRespiration_Atkin_etal_2017 + use LeafBiophysicsMod, only : LeafLayerBiophysicalRates + use LeafBiophysicsMod, only : LowstorageMainRespReduction + use LeafBiophysicsMod, only : rsmax0 + use LeafBiophysicsMod, only : DecayCoeffVcmax + use LeafBiophysicsMod, only : VeloToMolarCF + use FatesRadiationMemMod, only : idirect ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -89,35 +95,9 @@ module FATESPlantRespPhotosynthMod character(len=1024) :: warn_msg ! for defining a warning message - !------------------------------------------------------------------------------------- - - ! maximum stomatal resistance [s/m] (used across several procedures) - 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 - - ! Constants used to define C3 versus C4 photosynth pathways - integer, parameter :: c3_path_index = 1 - integer, parameter :: c4_path_index = 0 - - - ! Constants used to define conductance models - integer, parameter :: medlyn_model = 2 - integer, parameter :: ballberry_model = 1 - - ! Alternatively, Gross Assimilation can be used to estimate - ! leaf co2 partial pressure and therefore conductance. The default - ! is to use anet - integer, parameter :: net_assim_model = 1 - integer, parameter :: gross_assim_model = 2 - logical, parameter :: preserve_b4b = .true. contains @@ -142,16 +122,12 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) use EDCanopyStructureMod, only : calc_areaindex use FatesConstantsMod, only : umolC_to_kgC use FatesConstantsMod, only : umol_per_mmol - use FatesConstantsMod, only : rgas => rgas_J_K_kmol + use FatesParameterDerivedMod, only : param_derived - use FatesAllometryMod, only : bleaf, bstore_allom use FatesAllometryMod, only : storage_fraction_of_target use FatesAllometryMod, only : set_root_fraction - - use DamageMainMod, only : GetCrownReduction - use FatesInterfaceTypesMod, only : hlm_use_tree_damage ! ARGUMENTS: @@ -186,47 +162,46 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! leaf maintenance (dark) respiration [umol CO2/m**2/s] real(r8) :: lmr_z(nlevleaf,maxpft,nclmax) - + ! stomatal resistance [s/m] real(r8) :: rs_z(nlevleaf,maxpft,nclmax) ! net leaf photosynthesis averaged over sun and shade leaves. [umol CO2/m**2/s] real(r8) :: anet_av_z(nlevleaf,maxpft,nclmax) - ! Photosynthesis [umol /m2 /s] - real(r8) :: psn_z(nlevleaf,maxpft,nclmax) + ! photsynthesis + real(r8) :: psn_z(nlevleaf,maxpft,nclmax) + + ! carbon 13 in newly assimilated carbon at leaf level + real(r8) :: c13disc_z(nlevleaf,maxpft,nclmax) ! Mask used to determine which leaf-layer biophysical rates have been ! used already logical :: rate_mask_z(nlevleaf,maxpft,nclmax) real(r8) :: vcmax_z ! leaf layer maximum rate of carboxylation - ! (umol co2/m**2/s) + ! (umol co2/m**2/s) real(r8) :: jmax_z ! leaf layer maximum electron transport rate - ! (umol electrons/m**2/s) + ! (umol electrons/m**2/s) real(r8) :: kp_z ! leaf layer initial slope of CO2 response - ! curve (C4 plants) - real(r8) :: c13disc_z(nclmax,maxpft,nlevleaf) ! carbon 13 in newly assimilated carbon at leaf level - + ! curve (C4 plants) 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) :: 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) :: 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]) - real(r8) :: ceair ! vapor pressure of air, constrained (Pa) real(r8) :: nscaler ! leaf nitrogen scaling coefficient + real(r8) :: rdark_scaler ! scaling coefficient for Atkin dark respiration real(r8) :: leaf_frac ! ratio of to leaf biomass to total alive biomass real(r8) :: tcsoi ! Temperature response function for root respiration. real(r8) :: tcwood ! Temperature response function for wood - real(r8) :: patch_la ! exposed leaf area (patch scale) real(r8) :: live_stem_n ! Live stem (above-ground sapwood) - ! nitrogen content (kgN/plant) + ! nitrogen content (kgN/plant) real(r8) :: live_croot_n ! Live coarse root (below-ground sapwood) - ! nitrogen content (kgN/plant) + ! nitrogen content (kgN/plant) real(r8) :: sapw_c ! Sapwood carbon (kgC/plant) real(r8) :: store_c_target ! Target storage carbon (kgC/plant) real(r8) :: fnrt_c ! Fine root carbon (kgC/plant) @@ -260,7 +235,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: cumulative_lai ! the cumulative LAI, top down, to the leaf layer of interest real(r8) :: leaf_psi ! leaf xylem matric potential [MPa] (only meaningful/used w/ hydro) real(r8) :: fnrt_mr_layer ! fine root maintenance respiation per layer [kgC/plant/s] - + integer :: isunsha ! Index for differentiating sun and shade real(r8) :: fnrt_mr_nfix_layer ! fineroot maintenance respiration specifically for symbiotic fixation [kgC/plant/layer/s] real(r8) :: nfix_layer ! Nitrogen fixed in each layer this timestep [kgN/plant/layer/timestep] real(r8), allocatable :: rootfr_ft(:,:) ! Root fractions per depth and PFT @@ -278,6 +253,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: rd_abs_leaf, rb_abs_leaf, r_abs_stem, r_abs_snow, rb_abs, rd_abs real(r8) :: fsun real(r8) :: par_per_sunla, par_per_shala ! PAR per sunlit and shaded leaf area [W/m2 leaf] + real(r8) :: ac_utest, aj_utest, ap_utest, co2_inter_c_utest real(r8),dimension(75) :: cohort_vaitop real(r8),dimension(75) :: cohort_vaibot real(r8),dimension(75) :: cohort_layer_elai @@ -287,8 +263,25 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: cohort_elai real(r8) :: cohort_esai real(r8) :: laisun,laisha + real(r8) :: leaf_area ! m2 leaf per m2 footprint, for either sunlit or shaded leaf + real(r8) :: par_abs ! absorbed PAR [umol photons/m2leaf/s] + real(r8) :: par_sun, par_sha real(r8) :: canopy_area real(r8) :: elai + real(r8) :: test_out + real(r8) :: area_frac ! this is either f_sun , or 1-f_sun, its for area weighting + real(r8) :: psn_ll ! Leaf level photosyntheis + real(r8) :: gstoma_ll ! leaf level stomatal conductance (separate for sun shade) [m/s] + real(r8) :: gstoma ! stomatal conductance at leaf bin (sun/shade combined) [m/s] + real(r8) :: anet_ll ! leaf level net assimilation [umol CO2/m**2/s] + real(r8) :: c13disc_ll ! leaf level c13 assimilation + real(r8) :: hydr_k_lwp ! inner leaf humidity scaling coefficient [-] + real(r8) :: gs0 ! stomatal intercept, possibly scaled by btran depending on hypothesis + real(r8) :: gs1 ! stomatal slope, possibly scaled by btran depending on hypothesis + real(r8) :: gs2 ! optional btran scaling factor for Medlyn conductance only, instead + ! of applying to gs1, this would scale the whole non-intercept + ! portion of the conductance equation + ! ----------------------------------------------------------------------------------- ! Keeping these two definitions in case they need to be added later ! @@ -300,9 +293,14 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) integer :: nv ! number of leaf layers integer :: NCL_p ! number of canopy layers in patch integer :: iage ! loop counter for leaf age classes - + integer :: solve_iter + integer :: ico + ! Parameters - ! + ! Absolute convergence tolerance on solving intracellular CO2 concentration [Pa] + + real(r8), parameter :: ci_tol = 0.5_r8 + ! Base rate is at 20C. Adjust to 25C using the CN Q10 = 1.5 ! (gC/gN/s) ! ------------------------------------------------------------------------ @@ -316,12 +314,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) associate( & - c3psn => EDPftvarcon_inst%c3psn , & slatop => prt_params%slatop , & ! specific leaf area at top of canopy, - ! projected area basis [m^2/gC] - woody => prt_params%woody, & ! Is vegetation woody or not? - stomatal_intercept => EDPftvarcon_inst%stomatal_intercept ) !Unstressed minimum stomatal conductance - + ! projected area basis [m^2/gC] + woody => prt_params%woody) ! Is vegetation woody or not? do s = 1,nsites @@ -356,8 +351,6 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) bc_out(s)%rssun_pa(ifp) = 0._r8 bc_out(s)%rssha_pa(ifp) = 0._r8 - psn_z(:,:,:) = 0._r8 - g_sb_leaves = 0._r8 patch_la = 0._r8 @@ -374,7 +367,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! Part III. Calculate the number of sublayers for each pft and layer. ! And then identify which layer/pft combinations have things in them. ! Output: - ! currentPatch%nleaf(:,:) + ! currentPatch%ncan(:,:) ! currentPatch%canopy_mask(:,:) call UpdateCanopyNCanNRadPresent(currentPatch) @@ -390,16 +383,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) call GetCanopyGasParameters(bc_in(s)%forc_pbot, & ! in bc_in(s)%oair_pa(ifp), & ! in bc_in(s)%t_veg_pa(ifp), & ! in - bc_in(s)%tgcm_pa(ifp), & ! in - bc_in(s)%eair_pa(ifp), & ! in - bc_in(s)%esat_tv_pa(ifp), & ! in - bc_in(s)%rb_pa(ifp), & ! in mm_kco2, & ! out mm_ko2, & ! out - co2_cpoint, & ! out - cf, & ! out - gb_mol, & ! out - ceair) ! out + co2_cpoint) ! ------------------------------------------------------------------------ ! Part VI: Loop over all leaf layers. @@ -421,16 +407,24 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! scratch space. ! ------------------------------------------------------------------------ rate_mask_z(:,1:numpft,:) = .false. - + psn_z(:,:,:) = 0._r8 + anet_av_z(:,:,:) = 0._r8 + c13disc_z(:,:,:) = 0._r8 + rs_z(:,:,:) = 0._r8 + lmr_z(:,:,:) = 0._r8 + if_any_cohorts: if(currentPatch%countcohorts > 0.0)then + currentCohort => currentPatch%tallest + ico = 0 do_cohort_drive: do while (associated(currentCohort)) ! Cohort loop ! Identify the canopy layer (cl), functional type (ft) ! and the leaf layer (IV) for this cohort ft = currentCohort%pft cl = currentCohort%canopy_layer - + ico = ico + 1 + ! Calculate the cohort specific elai profile ! And the top and bottom edges of the veg area index ! of each layer bin are. Note, if the layers @@ -476,7 +470,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) call storage_fraction_of_target(store_c_target, & currentCohort%prt%GetState(store_organ, carbon12_element), & frac) - call lowstorage_maintresp_reduction(frac,currentCohort%pft, & + call LowstorageMainRespReduction(frac,currentCohort%pft, & maintresp_reduction_factor) ! are there any leaves of this pft in this layer? @@ -504,9 +498,16 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) (nleafage > 1) .or. & (hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then + + ! These values are incremented, therefore since + ! sometimes we re-do layers, we need to re-zero them as well + ! since it is not an over-write + psn_z(iv,ft,cl) = 0._r8 + anet_av_z(iv,ft,cl) = 0._r8 + c13disc_z(iv,ft,cl) = 0._r8 + 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_vai(:) is the total vegetation area index of each "leaf" layer @@ -527,8 +528,6 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) else - 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 @@ -546,12 +545,11 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) btran_eff = btran_eff*currentPatch%bstress_sal_ft(ft) endif - ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al ! (2010) Biogeosciences, 7, 1833-1859 - kn = decay_coeff_vcmax(currentCohort%vcmax25top, & + kn = DecayCoeffVcmax(currentCohort%vcmax25top, & prt_params%leafn_vert_scaler_coeff1(ft), & prt_params%leafn_vert_scaler_coeff2(ft)) @@ -587,6 +585,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) end select + ! Part VII: Calculate dark respiration (leaf maintenance) for this layer select case (maintresp_leaf_model) @@ -601,9 +600,22 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) case (lmrmodel_atkin_etal_2017) - call LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, & ! in - cumulative_lai, & ! in - currentCohort%vcmax25top, & ! in + ! This uses the relationship between leaf N and respiration from Atkin et al + ! for the top of the canopy, but then scales through the canopy based on a rdark_scaler. + ! To assume proportionality with N through the canopy following Lloyd et al. 2010, use the + ! default parameter value of 2.43, which results in the scaling of photosynthesis and respiration + ! being proportional through the canopy. To have a steeper decrease in respiration than photosynthesis + ! this number can be smaller. There is some observational evidence for this being the case + ! in Lamour et al. 2023. + + kn = DecayCoeffVcmax(currentCohort%vcmax25top, & + EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff1(ft), & + EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff2(ft)) + + rdark_scaler = exp(-kn * cumulative_lai) + + call LeafLayerMaintenanceRespiration_Atkin_etal_2017( lnc_top, & ! in + rdark_scaler, & ! in ft, & ! in bc_in(s)%t_veg_pa(ifp), & ! in currentPatch%tveg_lpa%GetMean(), & ! in @@ -615,7 +627,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) call endrun(msg=errMsg(sourcefile, __LINE__)) end select - + ! Pre-process PAR absorbed per unit leaf area for different schemes ! par_per_sunla = [W absorbed beam+diffuse radiation / m2 of sunlit leaves] ! par_per_shala = [W absorbed diffuse radiation / m2 of shaded leaves] @@ -674,70 +686,130 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) end if if_radsolver - ! Part VII: Calculate (1) maximum rate of carboxylation (vcmax), - ! (2) maximum electron transport rate, (3) triose phosphate - ! utilization rate and (4) the initial slope of CO2 response curve - ! (C4 plants). Earlier we calculated their base rates as dictated - ! by their plant functional type and some simple scaling rules for - ! nitrogen limitation baesd on canopy position (not prognostic). - ! These rates are the specific rates used in the actual photosynthesis - ! calculations that take localized environmental effects (temperature) - ! into consideration. - - call LeafLayerBiophysicalRates(par_per_sunla, & ! in - ft, & ! in - currentCohort%vcmax25top, & ! in - currentCohort%jmax25top, & ! in - currentCohort%kp25top, & ! in - nscaler, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - bc_in(s)%dayl_factor_pa(ifp), & ! in - currentPatch%tveg_lpa%GetMean(), & ! in - currentPatch%tveg_longterm%GetMean(),& ! in - btran_eff, & ! in - vcmax_z, & ! out - jmax_z, & ! out - kp_z ) ! out - - ! Part IX: This call calculates the actual photosynthesis for the - ! leaf layer, as well as the stomatal resistance and the net assimilated carbon. - - call LeafLayerPhotosynthesis(fsun, & ! in - par_per_sunla, & ! in - par_per_shala, & ! in - laisun, & ! in - laisha, & ! in - canopy_area, & ! in - ft, & ! in - vcmax_z, & ! in - jmax_z, & ! in - kp_z, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - bc_in(s)%esat_tv_pa(ifp), & ! in - bc_in(s)%forc_pbot, & ! in - bc_in(s)%cair_pa(ifp), & ! in - bc_in(s)%oair_pa(ifp), & ! in - btran_eff, & ! in - stomatal_intercept_btran, & ! in - cf, & ! in - gb_mol, & ! in - ceair, & ! in - mm_kco2, & ! in - mm_ko2, & ! in - co2_cpoint, & ! in - lmr_z(iv,ft,cl), & ! in - leaf_psi, & ! in - bc_in(s)%rb_pa(ifp), & ! in - psn_z(iv,ft,cl), & ! out - rs_z(iv,ft,cl), & ! out - anet_av_z(iv,ft,cl), & ! out - c13disc_z(cl,ft,iv)) ! out - rate_mask_z(iv,ft,cl) = .true. + + + ! Perform photosynthesis calculations on sunlit and shaded leaves + ! --------------------------------------------------------------- + + ! Calculate leaf boundary layer conductance in molar form [umol/m2/s] + gb_mol = (1._r8/bc_in(s)%rb_pa(ifp)) * & + VeloToMolarCF(bc_in(s)%forc_pbot,bc_in(s)%t_veg_pa(ifp)) + + gstoma = 0._r8 + do_sunsha: do isunsha = 1,2 + + ! Part VII: Calculate (1) maximum rate of carboxylation (vcmax), + ! (2) maximum electron transport rate, (3) triose phosphate + ! utilization rate and (4) the initial slope of CO2 response curve + ! (C4 plants). Earlier we calculated their base rates as dictated + ! by their plant functional type and some simple scaling rules for + ! nitrogen limitation baesd on canopy position (not prognostic). + ! These rates are the specific rates used in the actual photosynthesis + ! calculations that take localized environmental effects (temperature) + ! into consideration. + + call LeafLayerBiophysicalRates(ft, & ! in + currentCohort%vcmax25top, & ! in + currentCohort%jmax25top, & ! in + currentCohort%kp25top, & ! in + nscaler, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + bc_in(s)%dayl_factor_pa(ifp), & ! in + currentPatch%tveg_lpa%GetMean(), & ! in + currentPatch%tveg_longterm%GetMean(),& ! in + btran_eff, & ! in + vcmax_z, & ! out + jmax_z, & ! out + kp_z, & ! out + gs0, & ! out + gs1, & ! out + gs2 ) ! out + + ! Part IX: This call calculates the actual photosynthesis for the + ! leaf layer, as well as the stomatal resistance and the net assimilated carbon. + + if(isunsha == idirect) then + leaf_area = laisun*canopy_area + par_abs = ConvertPar(leaf_area, par_per_sunla) + area_frac = fsun + else + leaf_area = laisha*canopy_area + par_abs = ConvertPar(leaf_area, par_per_shala) + area_frac = 1._r8 - fsun + end if + + if ( (hlm_use_planthydro.eq.itrue .and. EDPftvarcon_inst%hydr_k_lwp(ft)>nearzero) ) then + hydr_k_lwp = EDPftvarcon_inst%hydr_k_lwp(ft) + else + hydr_k_lwp = 1._r8 + end if + call LeafLayerPhotosynthesis(par_per_sunla, & ! + par_abs, & ! in + leaf_area, & ! in + ft, & ! in + vcmax_z, & ! in + jmax_z, & ! in + kp_z, & ! in + gs0, & ! in + gs1, & ! in + gs2, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + bc_in(s)%forc_pbot, & ! in + bc_in(s)%cair_pa(ifp), & ! in + bc_in(s)%oair_pa(ifp), & ! in + gb_mol, & ! in + bc_in(s)%eair_pa(ifp), & ! in + mm_kco2, & ! in + mm_ko2, & ! in + co2_cpoint, & ! in + lmr_z(iv,ft,cl), & ! in + ci_tol, & ! in + psn_ll, & ! out + gstoma_ll, & ! out + anet_ll, & ! out + c13disc_ll, & ! out + ac_utest, & ! out (unit tests) + aj_utest, & ! out (unit tests) + ap_utest, & ! out (unit tests) + co2_inter_c_utest, & ! out (unit tests) + solve_iter) ! out performance tracking + + ! Average output quantities across sunlit and shaded leaves + ! Convert from molar to velocity (umol /m**2/s) to (m/s) + gstoma = gstoma + area_frac* (gstoma_ll / VeloToMolarCF(bc_in(s)%forc_pbot,bc_in(s)%t_veg_pa(ifp))) + + psn_z(iv,ft,cl) = psn_z(iv,ft,cl) + area_frac * psn_ll + anet_av_z(iv,ft,cl) = anet_av_z(iv,ft,cl) + area_frac * anet_ll + c13disc_z(iv,ft,cl) = c13disc_z(iv,ft,cl) + area_frac * c13disc_ll + + end do do_sunsha + + + + ! Stomatal resistance of the leaf-layer + if ( (hlm_use_planthydro.eq.itrue .and. EDPftvarcon_inst%hydr_k_lwp(ft)>nearzero) ) then + + gstoma = max(gstoma,1._r8/rsmax0) + rs_z(iv,ft,cl) = LeafHumidityStomaResis(leaf_psi, EDPftvarcon_inst%hydr_k_lwp(ft), bc_in(s)%t_veg_pa(ifp), & + bc_in(s)%cair_pa(ifp),bc_in(s)%forc_pbot, bc_in(s)%rb_pa(ifp), gstoma, ft) + + else + ! if gstoma is truly zero, this will be weighted by zero + rs_z(iv,ft,cl)= 1._r8/max(gstoma,1._r8/rsmax0) + + end if + + rate_mask_z(iv,ft,cl) = .true. + end if rate_mask_if end do leaf_layer_loop + !if(maxval(psn_z(:,ft,cl))>nearzero .and. ico==2) then + ! print*,psn_z(1:8,ft,cl) + !end if + ! Zero cohort flux accumulators. currentCohort%resp_m_tstep = 0.0_r8 @@ -759,11 +831,11 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if(radiation_model.eq.norman_solver) then call ScaleLeafLayerFluxToCohort(nv, & !in - psn_z(1:nv,ft,cl), & !in + psn_z(1:nv,ft,cl), & !in lmr_z(1:nv,ft,cl), & !in rs_z(1:nv,ft,cl), & !in currentPatch%elai_profile(cl,ft,1:nv), & !in - c13disc_z(cl, ft, 1:nv), & !in + c13disc_z(1:nv,ft,cl), & !in currentCohort%c_area, & !in currentCohort%n, & !in bc_in(s)%rb_pa(ifp), & !in @@ -776,12 +848,14 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) else + + call ScaleLeafLayerFluxToCohort(nv, & !in - psn_z(1:nv,ft,cl), & !in + psn_z(1:nv,ft,cl), & !in lmr_z(1:nv,ft,cl), & !in rs_z(1:nv,ft,cl), & !in cohort_layer_elai(1:nv), & !in - c13disc_z(cl, ft, 1:nv), & !in + c13disc_z(1:nv,ft,cl), & !in currentCohort%c_area, & !in currentCohort%n, & !in bc_in(s)%rb_pa(ifp), & !in @@ -793,7 +867,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) cohort_eleaf_area) !out end if - + ! Net Uptake does not need to be scaled, just transfer directly currentCohort%ts_net_uptake(1:nv) = anet_av_z(1:nv,ft,cl) * umolC_to_kgC @@ -998,6 +1072,15 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) currentCohort => currentCohort%shorter enddo do_cohort_drive + !if(maxval(psn_z(:,1,1))>nearzero)then + ! currentCohort => currentPatch%tallest + ! do while (associated(currentCohort)) ! Cohort loop + ! print*,currentCohort%gpp_tstep,currentCohort%pft,currentCohort%canopy_layer + ! currentCohort => currentCohort%shorter + ! enddo + ! stop + !end if + end if if_any_cohorts ! Normalize canopy total conductance by the effective LAI @@ -1059,7 +1142,10 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! This value is used for diagnostics, the molar form of conductance ! is what is used in the field usually, so we track that form - currentPatch%c_stomata = cf / r_stomata + ! cf = s/ umol/m2 -> s/m + !real(r8) :: cf ! s m**2/umol -> s/m (ideal gas conversion) [umol/m3] + + currentPatch%c_stomata = VeloToMolarCF(bc_in(s)%forc_pbot,bc_in(s)%t_veg_pa(ifp)) / r_stomata else !if_any_lai @@ -1069,13 +1155,13 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! This value is used for diagnostics, the molar form of conductance ! is what is used in the field usually, so we track that form - currentPatch%c_stomata = cf / rsmax0 + currentPatch%c_stomata = VeloToMolarCF(bc_in(s)%forc_pbot,bc_in(s)%t_veg_pa(ifp)) / rsmax0 end if if_any_lai ! This value is used for diagnostics, the molar form of conductance ! is what is used in the field usually, so we track that form - currentPatch%c_lblayer = cf / bc_in(s)%rb_pa(ifp) + currentPatch%c_lblayer = VeloToMolarCF(bc_in(s)%forc_pbot,bc_in(s)%t_veg_pa(ifp)) / bc_in(s)%rb_pa(ifp) end if if_filter2 @@ -1151,605 +1237,21 @@ end subroutine RootLayerNFixation ! ======================================================================================= -subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in - parsun_lsl, & ! in - parsha_lsl, & ! in - laisun_lsl, & ! in - laisha_lsl, & ! in - canopy_area_lsl, & ! in - ft, & ! in - vcmax, & ! in - jmax, & ! in - co2_rcurve_islope, & ! in - veg_tempk, & ! in - veg_esat, & ! in - can_press, & ! in - can_co2_ppress, & ! in - can_o2_ppress, & ! in - btran, & ! in - stomatal_intercept_btran, & ! in - cf, & ! in - gb_mol, & ! in - ceair, & ! in - mm_kco2, & ! in - mm_ko2, & ! in - co2_cpoint, & ! in - lmr, & ! in - leaf_psi, & ! in - rb, & ! in - psn_out, & ! out - rstoma_out, & ! out - anet_av_out, & ! out - c13disc_z) ! out - - - ! ------------------------------------------------------------------------------------ - ! This subroutine calculates photosynthesis and stomatal conductance within each leaf - ! sublayer. - ! A note on naming conventions: As this subroutine is called for every - ! leaf-sublayer, many of the arguments are specific to that "leaf sub layer" - ! (LSL), those variables are given a dimension tag "_lsl" - ! Other arguments or variables may be indicative of scales broader than the LSL. - ! ------------------------------------------------------------------------------------ - - use EDParamsMod , only : theta_cj_c3, theta_cj_c4 - - - ! Arguments - ! ------------------------------------------------------------------------------------ - real(r8), intent(in) :: f_sun_lsl ! - real(r8), intent(in) :: parsun_lsl ! Absorbed PAR in sunlist leaves - real(r8), intent(in) :: parsha_lsl ! Absorved PAR in shaded leaves - real(r8), intent(in) :: laisun_lsl ! LAI in sunlit leaves - real(r8), intent(in) :: laisha_lsl ! LAI in shaded leaves - real(r8), intent(in) :: canopy_area_lsl ! - integer, intent(in) :: ft ! (plant) Functional Type Index - real(r8), intent(in) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) - real(r8), intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) - real(r8), intent(in) :: co2_rcurve_islope ! initial slope of CO2 response curve (C4 plants) - real(r8), intent(in) :: veg_tempk ! vegetation temperature - real(r8), intent(in) :: veg_esat ! saturation vapor pressure at veg_tempk (Pa) - - ! Important Note on the following gas pressures. This photosynthesis scheme will iteratively - ! solve for the co2 partial pressure at the leaf surface (ie in the stomata). The reference - ! point for these input values are NOT within that boundary layer that separates the stomata from - ! the canopy air space. The reference point for these is on the outside of that boundary - ! layer. This routine, which operates at the leaf scale, makes no assumptions about what the - ! scale of the refernce is, it could be lower atmosphere, it could be within the canopy - ! but most likely it is the closest value one can get to the edge of the leaf's boundary - ! layer. We use the convention "can_" because a reference point of within the canopy - ! ia a best reasonable scenario of where we can get that information from. - - real(r8), intent(in) :: can_press ! Air pressure NEAR the surface of the leaf (Pa) - 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) :: 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) - real(r8), intent(in) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) - real(r8), intent(in) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) - real(r8), intent(in) :: co2_cpoint ! CO2 compensation point (Pa) - real(r8), intent(in) :: lmr ! Leaf Maintenance Respiration (umol CO2/m**2/s) - real(r8), intent(in) :: leaf_psi ! Leaf water potential [MPa] - real(r8), intent(in) :: rb ! Boundary Layer resistance of leaf [s/m] - - real(r8), intent(out) :: psn_out ! carbon assimilated in this leaf layer umolC/m2/s - real(r8), intent(out) :: rstoma_out ! stomatal resistance (1/gs_lsl) (s/m) - real(r8), intent(out) :: anet_av_out ! net leaf photosynthesis (umol CO2/m**2/s) - ! averaged over sun and shade leaves. - real(r8), intent(out) :: c13disc_z ! carbon 13 in newly assimilated carbon - - - - - ! Locals - ! ------------------------------------------------------------------------ - integer :: c3c4_path_index ! Index for which photosynthetic pathway - ! is active. C4 = 0, C3 = 1 - integer :: sunsha ! Index for differentiating sun and shade - real(r8) :: gstoma ! Stomatal Conductance of this leaf layer (m/s) - real(r8) :: agross ! co-limited gross leaf photosynthesis (umol CO2/m**2/s) - real(r8) :: anet ! net leaf photosynthesis (umol CO2/m**2/s) - real(r8) :: a_gs ! The assimilation (a) for calculating conductance (gs) - ! is either = to anet or agross - real(r8) :: je ! electron transport rate (umol electrons/m**2/s) - real(r8) :: qabs ! PAR absorbed by PS II (umol photons/m**2/s) - real(r8) :: aquad,bquad,cquad ! terms for quadratic equations - real(r8) :: r1,r2 ! roots of quadratic equation - real(r8) :: co2_inter_c ! intercellular leaf CO2 (Pa) - real(r8) :: co2_inter_c_old ! intercellular leaf CO2 (Pa) (previous iteration) - logical :: loop_continue ! Loop control variable - integer :: niter ! iteration loop index - real(r8) :: gs_mol ! leaf stomatal conductance (umol H2O/m**2/s) - real(r8) :: gs ! leaf stomatal conductance (m/s) - real(r8) :: hs ! fractional humidity at leaf surface (dimensionless) - real(r8) :: gs_mol_err ! gs_mol for error check - real(r8) :: ac ! Rubisco-limited gross photosynthesis (umol CO2/m**2/s) - real(r8) :: aj ! RuBP-limited gross photosynthesis (umol CO2/m**2/s) - real(r8) :: ap ! product-limited (C3) or CO2-limited - ! (C4) gross photosynthesis (umol CO2/m**2/s) - real(r8) :: ai ! intermediate co-limited photosynthesis (umol CO2/m**2/s) - 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) :: term ! intermediate variable in Medlyn stomatal conductance model - real(r8) :: vpd ! water vapor deficit in Medlyn stomatal model (KPa) - - - ! Parameters - ! ------------------------------------------------------------------------ - ! Fraction of light absorbed by non-photosynthetic pigments - real(r8),parameter :: fnps = 0.15_r8 - - ! term accounting that two photons are needed to fully transport a single - ! electron in photosystem 2 - real(r8), parameter :: photon_to_e = 0.5_r8 - - ! Unit conversion of w/m2 to umol photons m-2 s-1 - real(r8), parameter :: wm2_to_umolm2s = 4.6_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 - - ! First guess on ratio between intercellular co2 and the atmosphere - ! an iterator converges on actual - real(r8),parameter :: init_a2l_co2_c3 = 0.7_r8 - real(r8),parameter :: init_a2l_co2_c4 = 0.4_r8 - - ! quantum efficiency, used only for C4 (mol CO2 / mol photons) (index 0) - real(r8),parameter,dimension(0:1) :: quant_eff = [0.05_r8,0.0_r8] - - ! empirical curvature parameter for ap photosynthesis co-limitation - real(r8),parameter :: theta_ip = 0.999_r8 - - ! minimum Leaf area to solve, too little has shown instability - real(r8), parameter :: min_la_to_solve = 0.0000000001_r8 - - 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)) - - if (c3c4_path_index == c3_path_index) then - init_co2_inter_c = init_a2l_co2_c3 * can_co2_ppress - else - init_co2_inter_c = init_a2l_co2_c4 * can_co2_ppress - end if - - ! Part III: Photosynthesis and Conductance - ! ---------------------------------------------------------------------------------- - - if_daytime: if ( parsun_lsl <= 0._r8 ) then ! night time - - anet_av_out = -lmr - psn_out = 0._r8 - - ! The cuticular conductance already factored in maximum resistance as a bound - ! no need to re-bound it - - 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 - - else ! day time (a little bit more complicated ...) - - ! Is there leaf area? - (NV can be larger than 0 with only stem area if deciduous) - - if_leafarea: if ( laisun_lsl + laisha_lsl > 0._r8 ) then - - !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. - anet_av_out = 0._r8 - gstoma = 0._r8 - - do sunsha = 1,2 - ! Electron transport rate for C3 plants. - ! Convert par from W/m2 to umol photons/m**2/s - ! Convert from units of par absorbed per unit ground area to par - ! absorbed per unit leaf area. - - if(sunsha == 1)then !sunlit - if(( laisun_lsl * canopy_area_lsl) > min_la_to_solve)then - - qabs = parsun_lsl / (laisun_lsl * canopy_area_lsl ) - qabs = qabs * photon_to_e * (1._r8 - fnps) * wm2_to_umolm2s - - else - qabs = 0.0_r8 - end if - else - - if( (parsha_lsl>nearzero) .and. (laisha_lsl * canopy_area_lsl) > min_la_to_solve ) then - - qabs = parsha_lsl / (laisha_lsl * canopy_area_lsl) - qabs = qabs * photon_to_e * (1._r8 - fnps) * wm2_to_umolm2s - else - ! The radiative transfer schemes are imperfect - ! they can sometimes generate negative values here - qabs = 0._r8 - end if - - end if - - !convert the absorbed par into absorbed par per m2 of leaf, - ! so it is consistant with the vcmax and lmr numbers. - aquad = theta_psii - bquad = -(qabs + jmax) - cquad = qabs * jmax - call QuadraticRoots(aquad, bquad, cquad, r1, r2) - je = min(r1,r2) - - ! Initialize intercellular co2 - co2_inter_c = init_co2_inter_c - - niter = 0 - loop_continue = .true. - iter_loop: do while(loop_continue) - ! Increment iteration counter. Stop if too many iterations - niter = niter + 1 - - ! Save old co2_inter_c - co2_inter_c_old = co2_inter_c - - ! Photosynthesis limitation rate calculations - if (c3c4_path_index == c3_path_index)then - - ! C3: Rubisco-limited photosynthesis - ac = vcmax * max(co2_inter_c-co2_cpoint, 0._r8) / & - (co2_inter_c+mm_kco2 * (1._r8+can_o2_ppress / mm_ko2 )) - - ! C3: RuBP-limited photosynthesis - aj = je * max(co2_inter_c-co2_cpoint, 0._r8) / & - (4._r8*co2_inter_c+8._r8*co2_cpoint) - - ! Gross photosynthesis smoothing calculations. Co-limit ac and aj. - aquad = theta_cj_c3 - bquad = -(ac + aj) - cquad = ac * aj - call QuadraticRoots(aquad, bquad, cquad, r1, r2) - agross = min(r1,r2) - - else - - ! C4: Rubisco-limited photosynthesis - ac = vcmax - - ! C4: RuBP-limited photosynthesis - if(sunsha == 1)then !sunlit - !guard against /0's in the night. - if((laisun_lsl * canopy_area_lsl) > 0.0000000001_r8) then - aj = quant_eff(c3c4_path_index) * parsun_lsl * wm2_to_umolm2s - !convert from per cohort to per m2 of leaf) - aj = aj / (laisun_lsl * canopy_area_lsl) - else - aj = 0._r8 - end if - else - aj = quant_eff(c3c4_path_index) * parsha_lsl * wm2_to_umolm2s - aj = aj / (laisha_lsl * canopy_area_lsl) - end if - - ! C4: PEP carboxylase-limited (CO2-limited) - ap = co2_rcurve_islope * max(co2_inter_c, 0._r8) / can_press - - ! Gross photosynthesis smoothing calculations. First co-limit ac and aj. Then co-limit ap - - aquad = theta_cj_c4 - bquad = -(ac + aj) - cquad = ac * aj - call QuadraticRoots(aquad, bquad, cquad, r1, r2) - ai = min(r1,r2) - - aquad = theta_ip - bquad = -(ai + ap) - cquad = ai * ap - call QuadraticRoots(aquad, bquad, cquad, r1, r2) - agross = min(r1,r2) - - end if - - ! Calculate anet, only exit iteration with negative anet when - ! using anet in calculating gs this is version B - anet = agross - lmr - - if ( stomatal_assim_model == gross_assim_model ) then - if ( stomatal_model == medlyn_model ) then - write (fates_log(),*) 'Gross Assimilation conductance is incompatible with the Medlyn model' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - a_gs = agross - else - if (anet < 0._r8) then - loop_continue = .false. - end if - a_gs = anet - end if - - ! With an <= 0, then gs_mol = stomatal_intercept_btran - leaf_co2_ppress = can_co2_ppress- h2o_co2_bl_diffuse_ratio/gb_mol * a_gs * can_press - leaf_co2_ppress = max(leaf_co2_ppress,1.e-06_r8) - - ! A note about the use of the quadratic equations for calculating stomatal conductance - ! ------------------------------------------------------------------------------------ - ! These two following models calculate the conductance between the intercellular leaf - ! space and the leaf surface, not the canopy air space. Transport between the leaf - ! surface and the canopy air space is governed by the leaf boundary layer conductance. - ! However, we need to estimate the properties at the surface of the leaf to solve for - ! the stomatal conductance. We do this by using Fick's law (gradient resistance - ! approximation of diffusion) to estimate the flux of water vapor across the - ! leaf boundary layer, and balancing that with the flux across the stomata. It - ! results in the following equation for leaf surface humidity: - ! - ! e_s = (e_i g_s + e_c g_b)/(g_b + g_s) - ! - ! The leaf surface humidity (e_s) becomes an expression of canopy humidity (e_c), - ! intercellular humidity (e_i, which is the saturation humidity at leaf temperature), - ! boundary layer conductance (g_b) (these are all known) and stomatal conductance - ! (g_s) (this is still unknown). This expression is substituted into the stomatal - ! conductance equation. The resulting form of these equations becomes a quadratic. - ! - ! For a detailed explanation, see the FATES technical note, section - ! "1.11 Stomatal Conductance" - ! - ! ------------------------------------------------------------------------------------ - - - if ( stomatal_model == medlyn_model ) 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 QuadraticRoots(aquad, bquad, cquad, r1, r2) - gs_mol = max(r1,r2) - - else if ( stomatal_model == ballberry_model ) 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) * a_gs * can_press - cquad = -gb_mol*(leaf_co2_ppress*stomatal_intercept_btran + & - bb_slope(ft)*anet*can_press * ceair/ veg_esat ) - - call QuadraticRoots(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 * & - (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 - ! convergence criteria of +/- 1 x 10**-6 ppm is met OR if at least ten - ! iterations (niter=10) are completed - - if ((abs(co2_inter_c-co2_inter_c_old)/can_press*1.e06_r8 <= 2.e-06_r8) & - .or. niter == 5) then - loop_continue = .false. - end if - end do iter_loop - - ! End of co2_inter_c iteration. Check for an < 0, in which case gs_mol = bbb - ! And Final estimates for leaf_co2_ppress and co2_inter_c - ! (needed for early exit of co2_inter_c iteration when an < 0) - if (anet < 0._r8) then - gs_mol = stomatal_intercept_btran - end if - - ! Final estimates for leaf_co2_ppress and co2_inter_c - 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 * & - (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 - - ! 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 - - ! Accumulate total photosynthesis umol/m2 ground/s-1. - ! weight per unit sun and sha leaves. - if(sunsha == 1)then !sunlit - psn_out = psn_out + agross * f_sun_lsl - anet_av_out = anet_av_out + anet * f_sun_lsl - gstoma = gstoma + 1._r8/(min(1._r8/gs, rsmax0)) * f_sun_lsl - else - psn_out = psn_out + agross * (1.0_r8-f_sun_lsl) - anet_av_out = anet_av_out + anet * (1.0_r8-f_sun_lsl) - gstoma = gstoma + & - 1._r8/(min(1._r8/gs, rsmax0)) * (1.0_r8-f_sun_lsl) - end if - - ! Make sure iterative solution is correct - if (gs_mol < 0._r8) then - write (fates_log(),*)'Negative stomatal conductance:' - write (fates_log(),*)'gs_mol= ',gs_mol - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - ! 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 - warn_msg = 'Stomatal conductance error check - weak convergence: '//trim(N2S(gs_mol))//' '//trim(N2S(gs_mol_err)) - call FatesWarn(warn_msg,index=1) - end if - - enddo !sunsha loop - - ! Stomatal resistance of the leaf-layer - if ( (hlm_use_planthydro.eq.itrue .and. EDPftvarcon_inst%hydr_k_lwp(ft)>nearzero) ) then - rstoma_out = LeafHumidityStomaResis(leaf_psi, veg_tempk, ceair, can_press, veg_esat, & - rb, gstoma, ft) - else - rstoma_out = 1._r8/gstoma - end if - - - else - - ! 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 - anet_av_out = 0._r8 - - rstoma_out = min(rsmax0,cf/(stem_cuticle_loss_frac*stomatal_intercept(ft))) - c13disc_z = 0.0_r8 - - end if if_leafarea !is there leaf area? - - - end if if_daytime ! night or day - - - end associate - return - end subroutine LeafLayerPhotosynthesis - - ! ======================================================================================= - - function LeafHumidityStomaResis(leaf_psi, veg_tempk, ceair, can_press, veg_esat, & - rb, gstoma, ft) result(rstoma_out) - - ! ------------------------------------------------------------------------------------- - ! This calculates inner leaf humidity as a function of mesophyll water potential - ! Adopted from Vesala et al., 2017 https://www.frontiersin.org/articles/10.3389/fpls.2017.00054/full - ! - ! Equation 1 in Vesala et al: - ! lwp_star = wi/w0 = exp( k_lwp*leaf_psi*molar_mass_water/(rgas_J_k_mol * veg_tempk) ) - ! - ! Terms: - ! leaf_psi: leaf water potential [MPa] - ! k_lwp: inner leaf humidity scaling coefficient [-] - ! rgas_J_K_mol: universal gas constant, [J/K/mol], 8.3144598 - ! molar_mass_water, molar mass of water, [g/mol]: 18.0 - ! - ! Unit conversions: - ! 1 Pa = 1 N/m2 = 1 J/m3 - ! density of liquid water [kg/m3] = 1000 - ! - ! units of equation 1: exp( [MPa]*[g/mol]/( [J/K/mol] * [K] ) ) - ! [MJ/m3]*[g/mol]*[m3/kg]*[kg/g]*[J/MJ] / ([J/mol]) - ! dimensionless: [J/g]*[g/mol]/([J/mol]) - ! - ! Note: unit conversions drop out b/c [m3/kg]*[kg/g]*[J/MJ] = 1e-3*1.e-3*1e6 = 1.0 - ! - ! Junyan Ding 2021 - ! ------------------------------------------------------------------------------------- - - ! Arguments - real(r8) :: leaf_psi ! Leaf water potential [MPa] - real(r8) :: veg_tempk ! Leaf temperature [K] - real(r8) :: ceair ! vapor pressure of air, constrained [Pa] - real(r8) :: can_press ! Atmospheric pressure of canopy [Pa] - real(r8) :: veg_esat ! Saturated vapor pressure at veg surf [Pa] - real(r8) :: rb ! Leaf Boundary layer resistance [s/m] - real(r8) :: gstoma ! Stomatal Conductance of this leaf layer [m/s] - integer :: ft ! Plant Functional Type - real(r8) :: rstoma_out ! Total Stomatal resistance (stoma and BL) [s/m] - - ! Locals - real(r8) :: k_lwp ! Scaling coefficient for the ratio of leaf xylem - ! water potential to mesophyll water potential - real(r8) :: qs ! Specific humidity [g/kg] - real(r8) :: qsat ! Saturation specific humidity [g/kg] - real(r8) :: qsat_adj ! Adjusted saturation specific humidity [g/kg] - real(r8) :: lwp_star ! leaf water potential scaling coefficient - ! for inner leaf humidity, 0 means total dehydroted - ! leaf, 1 means total saturated leaf - - ! Note: to disable this control, set k_lwp to zero, LWP_star will be 1 - k_lwp = EDPftvarcon_inst%hydr_k_lwp(ft) - if (leaf_psi<0._r8) then - lwp_star = exp(k_lwp*leaf_psi*molar_mass_water/(rgas_J_K_mol *veg_tempk)) - else - lwp_star = 1._r8 - end if - - ! compute specific humidity from vapor pressure - ! q = molar_mass_ratio_vapdry*e/(can_press - (1-molar_mass_ratio_vapdry)*e) - ! source https://cran.r-project.org/web/packages/humidity/vignettes/humidity-measures.html - ! now adjust inner leaf humidity by LWP_star - - qs = molar_mass_ratio_vapdry * ceair / (can_press - (1._r8-molar_mass_ratio_vapdry) * ceair) - qsat = molar_mass_ratio_vapdry * veg_esat / (can_press - (1._r8-molar_mass_ratio_vapdry) * veg_esat) - qsat_adj = qsat*lwp_star - - ! Adjusting gs (compute a virtual gs) that will be passed to host model - - if ( qsat_adj < qs ) then - - ! if inner leaf vapor pressure is less then or equal to that at leaf surface - ! then set stomata resistance to be very large to stop the transpiration or back flow of vapor - rstoma_out = rsmax0 - - else - - rstoma_out = (qsat-qs)*( 1/gstoma + rb)/(qsat_adj - qs)-rb - - end if - - if (rstoma_out < nearzero ) then - write (fates_log(),*) 'qsat:', qsat, 'qs:', qs - write (fates_log(),*) 'LWP :', leaf_psi - write (fates_log(),*) 'ceair:', ceair, 'veg_esat:', veg_esat - write (fates_log(),*) 'rstoma_out:', rstoma_out, 'rb:', rb - write (fates_log(),*) 'LWP_star', lwp_star - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - end function LeafHumidityStomaResis - - - ! ===================================================================================== - - subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv - psn_llz, & ! in %psn_z(1:currentCohort%nv,ft,cl) - lmr_llz, & ! in lmr_z(1:currentCohort%nv,ft,cl) - rs_llz, & ! in rs_z(1:currentCohort%nv,ft,cl) - elai_llz, & ! in %elai_profile(cl,ft,1:currentCohort%nv) - c13disc_llz, & ! in c13disc_z(cl, ft, 1:currentCohort%nv) - c_area, & ! in currentCohort%c_area - nplant, & ! in currentCohort%n - rb, & ! in bc_in(s)%rb_pa(ifp) - maintresp_reduction_factor, & ! in - g_sb_laweight, & ! out currentCohort%g_sb_laweight [m/s] [m2-leaf] - gpp, & ! out currentCohort%gpp_tstep - rdark, & ! out currentCohort%rdark - c13disc_clm, & ! out currentCohort%c13disc_clm - cohort_eleaf_area ) ! out [m2] + subroutine ScaleLeafLayerFluxToCohort(nv, & ! in + psn_llz, & ! in + lmr_llz, & ! in + rs_llz, & ! in + elai_llz, & ! in + c13disc_llz, & ! in + c_area, & ! in + nplant, & ! in + rb, & ! in + maintresp_reduction_factor, & ! in + g_sb_laweight, & ! out + gpp, & ! out + rdark, & ! out + c13disc_clm, & ! out + cohort_eleaf_area ) ! out ! ------------------------------------------------------------------------------------ ! This subroutine effectively integrates leaf carbon fluxes over the @@ -1779,7 +1281,7 @@ subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv real(r8), intent(out) :: cohort_eleaf_area ! Effective leaf area of the cohort [m2] real(r8), intent(out) :: c13disc_clm ! unpacked Cohort level c13 discrimination real(r8) :: sum_weight ! sum of weight for unpacking d13c flux (c13disc_z) from - ! (canopy_layer, pft, leaf_layer) matrix to cohort (c13disc_clm) + ! (canopy_layer, pft, leaf_layer) matrix to cohort (c13disc_clm) ! GPP IN THIS SUBROUTINE IS A RATE. THE CALLING ARGUMENT IS GPP_TSTEP. AFTER THIS ! CALL THE RATE WILL BE MULTIPLIED BY THE INTERVAL TO GIVE THE INTEGRATED QUANT. @@ -1866,106 +1368,20 @@ end subroutine ScaleLeafLayerFluxToCohort ! ===================================================================================== - function ft1_f(tl, ha) result(ans) - ! - !!DESCRIPTION: - ! photosynthesis temperature response - ! - ! !REVISION HISTORY - ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 - ! 7/23/16: Copied over from CLM by Ryan Knox - ! - !!USES - use FatesConstantsMod, only : rgas => rgas_J_K_kmol - - ! - ! !ARGUMENTS: - real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temperature function (K) - real(r8), intent(in) :: ha ! activation energy in photosynthesis temperature function (J/mol) - ! - ! !LOCAL VARIABLES: - real(r8) :: ans - !------------------------------------------------------------------------------- - - ans = exp( ha / (rgas*1.e-3_r8*(tfrz+25._r8)) * (1._r8 - (tfrz+25._r8)/tl) ) - - return - end function ft1_f - - ! ===================================================================================== - - function fth_f(tl,hd,se,scaleFactor) result(ans) - ! - !!DESCRIPTION: - !photosynthesis temperature inhibition - ! - ! !REVISION HISTORY - ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 - ! 7/23/16: Copied over from CLM by Ryan Knox - ! - use FatesConstantsMod, only : rgas => rgas_J_K_kmol - - ! - ! !ARGUMENTS: - real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temp function (K) - real(r8), intent(in) :: hd ! deactivation energy in photosynthesis temp function (J/mol) - real(r8), intent(in) :: se ! entropy term in photosynthesis temp function (J/mol/K) - real(r8), intent(in) :: scaleFactor ! scaling factor for high temp inhibition (25 C = 1.0) - ! - ! !LOCAL VARIABLES: - real(r8) :: ans - !------------------------------------------------------------------------------- - - ans = scaleFactor / ( 1._r8 + exp( (-hd+se*tl) / (rgas*1.e-3_r8*tl) ) ) - - return - end function fth_f - - ! ===================================================================================== - - function fth25_f(hd,se)result(ans) - ! - !!DESCRIPTION: - ! scaling factor for photosynthesis temperature inhibition - ! - ! !REVISION HISTORY: - ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 - ! 7/23/16: Copied over from CLM by Ryan Knox - ! - !!USES - - use FatesConstantsMod, only : rgas => rgas_J_K_kmol - - ! - ! !ARGUMENTS: - real(r8), intent(in) :: hd ! deactivation energy in photosynthesis temp function (J/mol) - real(r8), intent(in) :: se ! entropy term in photosynthesis temp function (J/mol/K) - ! - ! !LOCAL VARIABLES: - real(r8) :: ans - !------------------------------------------------------------------------------- - - ans = 1._r8 + exp( (-hd+se*(tfrz+25._r8)) / (rgas*1.e-3_r8*(tfrz+25._r8)) ) - - return - end function fth25_f - - ! ===================================================================================== - subroutine UpdateCanopyNCanNRadPresent(currentPatch) ! --------------------------------------------------------------------------------- ! This subroutine calculates two patch level quanities: - ! currentPatch%nleaf and + ! currentPatch%ncan and ! currentPatch%canopy_mask ! - ! currentPatch%nleaf(:,:) is a two dimensional array that indicates + ! currentPatch%ncan(:,:) is a two dimensional array that indicates ! the total number of leaf layers (including those that are not exposed to light) ! in each canopy layer and for each functional type. ! ! currentPatch%nrad(:,:) is a two dimensional array that indicates ! the total number of EXPOSED leaf layers, but for all intents and purposes - ! in the photosynthesis routine, this appears to be the same as %nleaf... + ! in the photosynthesis routine, this appears to be the same as %ncan... ! ! currentPatch%canopy_mask(:,:) has the same dimensions, is binary, and ! indicates whether or not leaf layers are present (by evaluating the canopy area @@ -2004,7 +1420,7 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) currentPatch%nrad = currentPatch%nleaf ! Now loop through and identify which layer and pft combo has scattering elements - do cl = 1,currentPatch%ncl_p + do cl = 1,nclmax do ft = 1,numpft currentPatch%canopy_mask(cl,ft) = 0 do iv = 1, currentPatch%nrad(cl,ft); @@ -2017,478 +1433,32 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) return end subroutine UpdateCanopyNCanNRadPresent + + ! ===================================================================================== + - ! ==================================================================================== - - subroutine GetCanopyGasParameters(can_press, & - can_o2_partialpress, & - veg_tempk, & - air_tempk, & - air_vpress, & - veg_esat, & - rb, & - mm_kco2, & - mm_ko2, & - co2_cpoint, & - cf, & - gb_mol, & - ceair) - - ! --------------------------------------------------------------------------------- - ! This subroutine calculates the specific Michaelis Menten Parameters (pa) for CO2 - ! and O2, as well as the CO2 compentation point. - ! --------------------------------------------------------------------------------- - - use FatesConstantsMod, only: umol_per_mol - use FatesConstantsMod, only: mmol_per_mol - use FatesConstantsMod, only: umol_per_kmol - use FatesConstantsMod, only : rgas => rgas_J_K_kmol - - ! Arguments - real(r8), intent(in) :: can_press ! Air pressure within the canopy (Pa) - real(r8), intent(in) :: can_o2_partialpress ! Partial press of o2 in the canopy (Pa) - real(r8), intent(in) :: veg_tempk ! The temperature of the vegetation (K) - real(r8), intent(in) :: air_tempk ! Temperature of canopy air (K) - real(r8), intent(in) :: air_vpress ! Vapor pressure of canopy air (Pa) - real(r8), intent(in) :: veg_esat ! Saturated vapor pressure at veg surf (Pa) - real(r8), intent(in) :: rb ! Leaf Boundary layer resistance (s/m) - - real(r8), intent(out) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) - real(r8), intent(out) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) - real(r8), intent(out) :: co2_cpoint ! CO2 compensation point (Pa) - real(r8), intent(out) :: cf ! conversion factor between molar form and velocity form - ! of conductance and resistance: [umol/m3] - real(r8), intent(out) :: gb_mol ! leaf boundary layer conductance (umol H2O/m**2/s) - real(r8), intent(out) :: ceair ! vapor pressure of air, constrained (Pa) - - ! Locals - real(r8) :: kc25 ! Michaelis-Menten constant for CO2 at 25C (Pa) - real(r8) :: ko25 ! Michaelis-Menten constant for O2 at 25C (Pa) - real(r8) :: sco ! relative specificity of rubisco - real(r8) :: cp25 ! CO2 compensation point at 25C (Pa) - - ! --------------------------------------------------------------------------------- - ! Intensive values (per mol of air) - ! kc, ko, currentPatch, from: Bernacchi et al (2001) - ! Plant, Cell and Environment 24:253-259 - ! --------------------------------------------------------------------------------- - - real(r8), parameter :: mm_kc25_umol_per_mol = 404.9_r8 - real(r8), parameter :: mm_ko25_mmol_per_mol = 278.4_r8 - real(r8), parameter :: co2_cpoint_umol_per_mol = 42.75_r8 - - ! Activation energy, from: - ! Bernacchi et al (2001) Plant, Cell and Environment 24:253-259 - ! Bernacchi et al (2003) Plant, Cell and Environment 26:1419-1430 - - real(r8), parameter :: kcha = 79430._r8 ! activation energy for kc (J/mol) - real(r8), parameter :: koha = 36380._r8 ! activation energy for ko (J/mol) - real(r8), parameter :: cpha = 37830._r8 ! activation energy for cp (J/mol) - - - ! Derive sco from currentPatch and O2 using present-day O2 (0.209 mol/mol) and re-calculate - ! currentPatch to account for variation in O2 using currentPatch = 0.5 O2 / sco - - ! FIXME (RGK 11-30-2016 THere are more constants here, but I don't have enough information - ! about what they are or do, so I can't give them more descriptive names. Someone please - ! fill this in when possible) - - kc25 = ( mm_kc25_umol_per_mol / umol_per_mol ) * can_press - ko25 = ( mm_ko25_mmol_per_mol / mmol_per_mol ) * can_press - sco = 0.5_r8 * 0.209_r8 / (co2_cpoint_umol_per_mol / umol_per_mol ) - cp25 = 0.5_r8 * can_o2_partialpress / sco - - if( veg_tempk.gt.150_r8 .and. veg_tempk.lt.350_r8 )then - mm_kco2 = kc25 * ft1_f(veg_tempk, kcha) - mm_ko2 = ko25 * ft1_f(veg_tempk, koha) - co2_cpoint = cp25 * ft1_f(veg_tempk, cpha) - else - mm_kco2 = 1.0_r8 - mm_ko2 = 1.0_r8 - co2_cpoint = 1.0_r8 - end if - ! --------------------------------------------------------------------------------- - ! - ! cf is the conversion factor between molar form and velocity form - ! of conductance and resistance: [umol/m3] + real(r8) function ConvertPar(leaf_area, par_wm2) result(par_umolm2s) ! - ! i.e. - ! [m/s] * [umol/m3] -> [umol/m2/s] + ! DESCRIPTION: + ! Convert par from W/m2 to umol photons/m2leaf/s ! - ! Breakdown of the conversion factor: [ umol / m3 ] - ! - ! Rgas [J /K /kmol] - ! Air Potential Temperature [ K ] - ! Canopy Pressure [ Pa ] - ! conversion: umol/kmol = 1e9 - ! - ! [ Pa * K * kmol umol/kmol / J K ] = [ Pa * umol / J ] - ! since: 1 Pa = 1 N / m2 - ! [ Pa * umol / J ] = [ N * umol / J m2 ] - ! since: 1 J = 1 N * m - ! [ N * umol / J m2 ] = [ N * umol / N m3 ] - ! [ umol / m3 ] - ! - ! -------------------------------------------------------------------------------- - - cf = can_press/(rgas * air_tempk )*umol_per_kmol - gb_mol = (1._r8/ rb) * cf - - ! Constrain eair >= 0.05*esat_tv so that solution does not blow up. This ensures - ! that hs does not go to zero. Also eair <= veg_esat so that hs <= 1 - ceair = min( max(air_vpress, 0.05_r8*veg_esat ),veg_esat ) - - - - return - end subroutine GetCanopyGasParameters - - ! ==================================================================================== - - subroutine LeafLayerMaintenanceRespiration_Ryan_1991(lnc_top, & - nscaler, & - ft, & - veg_tempk, & - lmr) - use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - use FatesConstantsMod, only : umolC_to_kgC - use FatesConstantsMod, only : g_per_kg - use EDPftvarcon , only : EDPftvarcon_inst - - ! ----------------------------------------------------------------------- - ! Base maintenance respiration rate for plant tissues maintresp_leaf_ryan1991_baserate - ! M. Ryan, 1991. Effects of climate change on plant respiration. - ! Ecological Applications, 1(2), 157-167. - ! Original expression is br = 0.0106 molC/(molN h) - ! Conversion by molecular weights of C and N gives 2.525e-6 gC/(gN s) - ! Which is the default value of maintresp_nonleaf_baserate - - ! Arguments - real(r8), intent(in) :: lnc_top ! Leaf nitrogen content per unit area at canopy top [gN/m2] - real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile - integer, intent(in) :: ft ! (plant) Functional Type Index - real(r8), intent(in) :: veg_tempk ! vegetation temperature - real(r8), intent(out) :: lmr ! Leaf Maintenance Respiration (umol CO2/m**2/s) - - ! Locals - real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) - real(r8) :: lmr25top ! canopy top leaf maint resp rate at 25C for this pft (umol CO2/m**2/s) - integer :: c3c4_path_index ! Index for which photosynthetic pathway - - ! Parameter - real(r8), parameter :: lmrha = 46390._r8 ! activation energy for lmr (J/mol) - real(r8), parameter :: lmrhd = 150650._r8 ! deactivation energy for lmr (J/mol) - real(r8), parameter :: lmrse = 490._r8 ! entropy term for lmr (J/mol/K) - real(r8), parameter :: lmrc = 1.15912391_r8 ! scaling factor for high - ! temperature inhibition (25 C = 1.0) - - lmr25top = EDPftvarcon_inst%maintresp_leaf_ryan1991_baserate(ft) * (1.5_r8 ** ((25._r8 - 20._r8)/10._r8)) - lmr25top = lmr25top * lnc_top / (umolC_to_kgC * g_per_kg) - - - ! Part I: Leaf Maintenance respiration: umol CO2 / m**2 [leaf] / s - ! ---------------------------------------------------------------------------------- - lmr25 = lmr25top * nscaler - - ! photosynthetic pathway: 0. = c4, 1. = c3 - c3c4_path_index = nint(EDPftvarcon_inst%c3psn(ft)) - - if (c3c4_path_index == c3_path_index) then - ! temperature sensitivity of C3 plants - lmr = lmr25 * ft1_f(veg_tempk, lmrha) * & - fth_f(veg_tempk, lmrhd, lmrse, lmrc) - else - ! temperature sensitivity of C4 plants - lmr = lmr25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) - lmr = lmr / (1._r8 + exp( 1.3_r8*(veg_tempk-(tfrz+55._r8)) )) - endif - - ! Any hydrodynamic limitations could go here, currently none - ! lmr = lmr * (nothing) - - end subroutine LeafLayerMaintenanceRespiration_Ryan_1991 - - ! ==================================================================================== - - subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, & - cumulative_lai, & - vcmax25top, & - ft, & - veg_tempk, & - tgrowth, & - lmr) - - use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - use FatesConstantsMod, only : umolC_to_kgC - use FatesConstantsMod, only : g_per_kg - use FatesConstantsMod, only : lmr_b - use FatesConstantsMod, only : lmr_c - use FatesConstantsMod, only : lmr_TrefC - use FatesConstantsMod, only : lmr_r_1 - use FatesConstantsMod, only : lmr_r_2 - use EDPftvarcon , only : EDPftvarcon_inst - - ! Arguments - real(r8), intent(in) :: lnc_top ! Leaf nitrogen content per unit area at canopy top [gN/m2] - integer, intent(in) :: ft ! (plant) Functional Type Index - real(r8), intent(in) :: vcmax25top ! top of canopy vcmax - real(r8), intent(in) :: cumulative_lai ! cumulative lai above the current leaf layer - real(r8), intent(in) :: veg_tempk ! vegetation temperature (degrees K) - real(r8), intent(in) :: tgrowth ! lagged vegetation temperature averaged over acclimation timescale (degrees K) - real(r8), intent(out) :: lmr ! Leaf Maintenance Respiration (umol CO2/m**2/s) - - ! Locals - real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) - real(r8) :: r_0 ! base respiration rate, PFT-dependent (umol CO2/m**2/s) - real(r8) :: r_t_ref ! acclimated ref respiration rate (umol CO2/m**2/s) - real(r8) :: lmr25top ! canopy top leaf maint resp rate at 25C for this pft (umol CO2/m**2/s) + ! ARGUMENTS: + real(r8), intent(in) :: leaf_area ! leaf area index [m2 leaf / m2 ground] + real(r8), intent(in) :: par_wm2 ! absorbed PAR [W/m2 ground] - real(r8) :: rdark_scaler ! negative exponential scaling of rdark - real(r8) :: kn ! decay coefficient - - ! parameter values of r_0 as listed in Atkin et al 2017: (umol CO2/m**2/s) - ! Broad-leaved trees 1.7560 - ! Needle-leaf trees 1.4995 - ! Shrubs 2.0749 - ! C3 herbs/grasses 2.1956 - ! In the absence of better information, we use the same value for C4 grasses as C3 grasses. - - ! r_0 currently put into the EDPftvarcon_inst%dev_arbitrary_pft - ! all figs in Atkin et al 2017 stop at zero Celsius so we will assume acclimation is fixed below that - r_0 = EDPftvarcon_inst%maintresp_leaf_atkin2017_baserate(ft) - - ! This code uses the relationship between leaf N and respiration from Atkin et al - ! for the top of the canopy, but then scales through the canopy based on a rdark_scaler. - ! To assume proportionality with N through the canopy following Lloyd et al. 2010, use the - ! default parameter value of 2.43, which results in the scaling of photosynthesis and respiration - ! being proportional through the canopy. To have a steeper decrease in respiration than photosynthesis - ! this number can be smaller. There is some observational evidence for this being the case - ! in Lamour et al. 2023. - - kn = decay_coeff_vcmax(vcmax25top, & - EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff1(ft), & - EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff2(ft)) - - rdark_scaler = exp(-kn * cumulative_lai) + ! minimum Leaf area to solve, too little has shown instability + real(r8), parameter :: min_la_to_solve = 0.0000000001_r8 - r_t_ref = max(0._r8, rdark_scaler * (r_0 + lmr_r_1 * lnc_top + lmr_r_2 * max(0._r8, (tgrowth - tfrz) )) ) - - if (r_t_ref .eq. 0._r8) then - warn_msg = 'Rdark is negative at this temperature and is capped at 0. tgrowth (C): '//trim(N2S(tgrowth-tfrz))//' pft: '//trim(I2S(ft)) - call FatesWarn(warn_msg,index=4) - end if - - lmr = r_t_ref * exp(lmr_b * (veg_tempk - tfrz - lmr_TrefC) + lmr_c * & - ((veg_tempk-tfrz)**2 - lmr_TrefC**2)) - - end subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017 - - ! ==================================================================================== - - subroutine LeafLayerBiophysicalRates( parsun_per_la, & - ft, & - vcmax25top_ft, & - jmax25top_ft, & - co2_rcurve_islope25top_ft, & - nscaler, & - veg_tempk, & - dayl_factor, & - t_growth, & - t_home, & - btran, & - vcmax, & - jmax, & - co2_rcurve_islope ) - - ! --------------------------------------------------------------------------------- - ! This subroutine calculates the localized rates of several key photosynthesis - ! rates. By localized, we mean specific to the plant type and leaf layer, - ! which factors in leaf physiology, as well as environmental effects. - ! This procedure should be called prior to iterative solvers, and should - ! have pre-calculated the reference rates for the pfts before this. - ! - ! The output biophysical rates are: - ! vcmax: maximum rate of carboxilation, - ! jmax: maximum electron transport rate, - ! co2_rcurve_islope: initial slope of CO2 response curve (C4 plants) - ! --------------------------------------------------------------------------------- - - use EDPftvarcon , only : EDPftvarcon_inst - - ! Arguments - ! ------------------------------------------------------------------------------ - - real(r8), intent(in) :: parsun_per_la ! PAR absorbed per sunlit leaves for this layer - integer, intent(in) :: ft ! (plant) Functional Type Index - real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile - real(r8), intent(in) :: vcmax25top_ft ! canopy top maximum rate of carboxylation at 25C - ! for this pft (umol CO2/m**2/s) - real(r8), intent(in) :: jmax25top_ft ! canopy top maximum electron transport rate at 25C - ! for this pft (umol electrons/m**2/s) - real(r8), intent(in) :: co2_rcurve_islope25top_ft ! initial slope of CO2 response curve - ! (C4 plants) at 25C, canopy top, this pft - real(r8), intent(in) :: veg_tempk ! vegetation temperature - real(r8), intent(in) :: dayl_factor ! daylength scaling factor (0-1) - real(r8), intent(in) :: t_growth ! T_growth (short-term running mean temperature) (K) - real(r8), intent(in) :: t_home ! T_home (long-term running mean temperature) (K) - real(r8), intent(in) :: btran ! transpiration wetness factor (0 to 1) - - real(r8), intent(out) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) - real(r8), intent(out) :: jmax ! maximum electron transport rate - ! (umol electrons/m**2/s) - real(r8), intent(out) :: co2_rcurve_islope ! initial slope of CO2 response curve (C4 plants) - - ! Locals - ! ------------------------------------------------------------------------------- - real(r8) :: vcmax25 ! leaf layer: maximum rate of carboxylation at 25C - ! (umol CO2/m**2/s) - real(r8) :: jmax25 ! leaf layer: maximum electron transport rate at 25C - ! (umol electrons/m**2/s) - real(r8) :: co2_rcurve_islope25 ! leaf layer: Initial slope of CO2 response curve - ! (C4 plants) at 25C - integer :: c3c4_path_index ! Index for which photosynthetic pathway - real(r8) :: dayl_factor_local ! Local version of daylength factor - - ! Parameters - ! --------------------------------------------------------------------------------- - real(r8) :: vcmaxha ! activation energy for vcmax (J/mol) - real(r8) :: jmaxha ! activation energy for jmax (J/mol) - real(r8) :: vcmaxhd ! deactivation energy for vcmax (J/mol) - real(r8) :: jmaxhd ! deactivation energy for jmax (J/mol) - real(r8) :: vcmaxse ! entropy term for vcmax (J/mol/K) - real(r8) :: jmaxse ! entropy term for jmax (J/mol/K) - real(r8) :: t_growth_celsius ! average growing temperature - real(r8) :: t_home_celsius ! average home temperature - real(r8) :: jvr ! ratio of Jmax25 / Vcmax25 - real(r8) :: vcmaxc ! scaling factor for high temperature inhibition (25 C = 1.0) - real(r8) :: jmaxc ! scaling factor for high temperature inhibition (25 C = 1.0) - - select case(photo_tempsens_model) - case (photosynth_acclim_model_none) !No temperature acclimation - vcmaxha = EDPftvarcon_inst%vcmaxha(FT) - jmaxha = EDPftvarcon_inst%jmaxha(FT) - vcmaxhd = EDPftvarcon_inst%vcmaxhd(FT) - jmaxhd = EDPftvarcon_inst%jmaxhd(FT) - vcmaxse = EDPftvarcon_inst%vcmaxse(FT) - jmaxse = EDPftvarcon_inst%jmaxse(FT) - case (photosynth_acclim_model_kumarathunge_etal_2019) !Kumarathunge et al. temperature acclimation, Thome=30-year running mean - t_growth_celsius = t_growth-tfrz - t_home_celsius = t_home-tfrz - vcmaxha = (42.6_r8 + (1.14_r8*t_growth_celsius))*1e3_r8 !J/mol - jmaxha = 40.71_r8*1e3_r8 !J/mol - vcmaxhd = 200._r8*1e3_r8 !J/mol - jmaxhd = 200._r8*1e3_r8 !J/mol - vcmaxse = (645.13_r8 - (0.38_r8*t_growth_celsius)) - jmaxse = 658.77_r8 - (0.84_r8*t_home_celsius) - 0.52_r8*(t_growth_celsius-t_home_celsius) - jvr = 2.56_r8 - (0.0375_r8*t_home_celsius)-(0.0202_r8*(t_growth_celsius-t_home_celsius)) - case default - write (fates_log(),*)'error, incorrect leaf photosynthesis temperature acclimation model specified' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end select - - vcmaxc = fth25_f(vcmaxhd, vcmaxse) - jmaxc = fth25_f(jmaxhd, jmaxse) - - if(parsun_per_la <= 0._r8) then - vcmax = 0._r8 - jmax = 0._r8 - co2_rcurve_islope = 0._r8 - else ! day time - - ! update the daylength factor local variable if the switch is on - if ( dayl_switch == itrue ) then - dayl_factor_local = dayl_factor - else - dayl_factor_local = 1.0_r8 - endif - - ! Vcmax25top was already calculated to derive the nscaler function - vcmax25 = vcmax25top_ft * nscaler * dayl_factor_local - select case(photo_tempsens_model) - case (photosynth_acclim_model_none) - jmax25 = jmax25top_ft * nscaler * dayl_factor_local - case (photosynth_acclim_model_kumarathunge_etal_2019) - jmax25 = vcmax25*jvr - case default - write (fates_log(),*)'error, incorrect leaf photosynthesis temperature acclimation model specified' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end select - - co2_rcurve_islope25 = co2_rcurve_islope25top_ft * nscaler - - ! Adjust for temperature - ! photosynthetic pathway: 0. = c4, 1. = c3 - c3c4_path_index = nint(EDPftvarcon_inst%c3psn(ft)) - - if (c3c4_path_index == c3_path_index) then - vcmax = vcmax25 * ft1_f(veg_tempk, vcmaxha) * fth_f(veg_tempk, vcmaxhd, vcmaxse, vcmaxc) - else - vcmax = vcmax25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) - vcmax = vcmax / (1._r8 + exp( 0.2_r8*((tfrz+15._r8)-veg_tempk ) )) - vcmax = vcmax / (1._r8 + exp( 0.3_r8*(veg_tempk-(tfrz+40._r8)) )) - end if - - jmax = jmax25 * ft1_f(veg_tempk, jmaxha) * fth_f(veg_tempk, jmaxhd, jmaxse, jmaxc) - - !q10 response of product limited psn. - co2_rcurve_islope = co2_rcurve_islope25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) + if (par_wm2 > nearzero .and. leaf_area > min_la_to_solve) then + par_umolm2s = par_wm2/leaf_area*wm2_to_umolm2s + else + ! The radiative transfer schemes are imperfect + ! they can sometimes generate negative values here if par or leaf area is 0.0 + par_umolm2s = 0.0_r8 end if - ! Adjust for water limitations - vcmax = vcmax * btran - - return - - end subroutine LeafLayerBiophysicalRates - - subroutine lowstorage_maintresp_reduction(frac, pft, maintresp_reduction_factor) - - ! This subroutine reduces maintenance respiration rates when storage pool is low. The premise - ! of this is that mortality of plants increases when storage is low because they are not able - ! to repair tissues, generate defense compounds, etc. This reduction is reflected in a reduced - ! maintenance demand. The output of this function takes the form of a curve between 0 and 1, - ! and the curvature of the function is determined by a parameter. - - ! Uses - use EDPftvarcon , only : EDPftvarcon_inst - - ! Arguments - ! ------------------------------------------------------------------------------ - real(r8), intent(in) :: frac ! ratio of storage to target leaf biomass - integer, intent(in) :: pft ! what pft is this cohort? - real(r8), intent(out) :: maintresp_reduction_factor ! the factor by which to reduce maintenance respiration - - ! -------------------------------------------------------------------------------- - ! Parameters are at the PFT level: - ! fates_maintresp_reduction_curvature controls the curvature of this. - ! If this parameter is zero, then there is no reduction until the plant dies at storage = 0. - ! If this parameter is one, then there is a linear reduction in respiration below the storage point. - ! Intermediate values will give some (concave-downwards) curvature. - ! - ! maintresp_reduction_intercept controls the maximum amount of throttling. - ! zero means no throttling at any point, so it turns this mechanism off completely and so - ! allows an entire cohort to die via negative carbon-induced termination mortality. - ! one means complete throttling, so no maintenance respiration at all, when out of carbon. - ! --------------------------------------------------------------------------------- - - if( frac .lt. 1._r8 )then - if ( abs(EDPftvarcon_inst%maintresp_reduction_curvature(pft)-1._r8) > nearzero ) then - maintresp_reduction_factor = (1._r8 - EDPftvarcon_inst%maintresp_reduction_intercept(pft)) + & - EDPftvarcon_inst%maintresp_reduction_intercept(pft) * & - (1._r8 - EDPftvarcon_inst%maintresp_reduction_curvature(pft)**frac) & - / (1._r8-EDPftvarcon_inst%maintresp_reduction_curvature(pft)) - else ! avoid nan answer for linear case - maintresp_reduction_factor = (1._r8 - EDPftvarcon_inst%maintresp_reduction_intercept(pft)) + & - EDPftvarcon_inst%maintresp_reduction_intercept(pft) * frac - endif - - else - maintresp_reduction_factor = 1._r8 - endif - - - end subroutine lowstorage_maintresp_reduction + end function ConvertPar end module FATESPlantRespPhotosynthMod diff --git a/biogeophys/LeafBiophysicsMod.F90 b/biogeophys/LeafBiophysicsMod.F90 new file mode 100644 index 0000000000..cce0e0dbf2 --- /dev/null +++ b/biogeophys/LeafBiophysicsMod.F90 @@ -0,0 +1,2275 @@ +module LeafBiophysicsMod + + !------------------------------------------------------------------------------------- + ! !DESCRIPTION: + ! + ! This module contains routines for leaf-level biophyics. These + ! routines are all called with primitive arguments to facilitate + ! use accross models, with the exception of an internally defined + ! set of constants associated with plant functional type. + ! + ! Assumptions: + ! + ! Units are always in micro-moles (umol), square meters (m2) and seconds + ! + ! Calculations are per-square-meter of leaf tissue. + ! + ! + ! + ! ------------------------------------------------------------------------------------ + + use shr_log_mod, only : errMsg => shr_log_errMsg + use shr_sys_mod, only : shr_sys_abort + use FatesConstantsMod, only : r8 => fates_r8 + use FatesGlobals, only : endrun => fates_endrun + use FatesGlobals, only : fates_log + use FatesGlobals, only : FatesWarn,N2S,A2S,I2S + use FatesConstantsMod, only : itrue + use FatesConstantsMod, only : nearzero + use FatesConstantsMod, only : molar_mass_ratio_vapdry + use FatesConstantsMod, only : molar_mass_water + use FatesConstantsMod, only : fates_unset_r8 + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm + use FatesConstantsMod, only : nocomp_bareground + use FatesConstantsMod, only : lmrmodel_ryan_1991 + use FatesConstantsMod, only : lmrmodel_atkin_etal_2017 + use FatesConstantsMod, only : kpa_per_pa + use FatesConstantsMod, only : umol_per_kmol + use FatesUtilsMod, only : QuadraticRoots => QuadraticRootsSridharachary + use FatesConstantsMod, only : rgas_J_K_kmol + use FatesConstantsMod, only : rgas_J_K_mol + use FatesConstantsMod, only : g_per_kg + use FatesConstantsMod, only : umolC_to_kgC + + implicit none + private + + public :: LeafLayerPhotosynthesis + public :: LeafHumidityStomaResis + public :: GetCanopyGasParameters + public :: LeafLayerMaintenanceRespiration_Ryan_1991 + public :: LeafLayerMaintenanceRespiration_Atkin_etal_2017 + public :: LeafLayerBiophysicalRates + public :: LowstorageMainRespReduction + public :: GetConstrainedVPress + public :: DecayCoeffVcmax + public :: QSat + public :: AgrossRubiscoC3 + public :: AgrossRuBPC3 + public :: AgrossRuBPC4 + public :: AgrossPEPC4 + public :: StomatalCondMedlyn + public :: StomatalCondBallBerry + public :: VeloToMolarCF + public :: CiMinMax + public :: CiFunc + public :: CiBisection + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + + + character(len=1024) :: warn_msg ! for defining a warning message + + !------------------------------------------------------------------------------------- + + ! maximum stomatal resistance [s/m] (used across several procedures) + real(r8),public, parameter :: rsmax0 = 2.e8_r8 + + ! for non cuticular (ie through branches) + ! minimum allowable stomatal conductance for numerics purposes [m/s] + real(r8),parameter :: gsmin0 = 1._r8/rsmax0 + + ! minimum allowable conductance for numerics purposes at 20C (293.15K) + ! and 1 standard atmosphere (101325 Pa) [umol/m2/s] + ! this follows the same math as FatesPlantRespPhotosynthMod:FetMolarVeloCF() + real(r8),parameter :: gsmin0_20C1A_mol = gsmin0 * 101325.0_r8/(rgas_J_K_kmol * 293.15 )*umol_per_kmol + + ! Set this to true to perform debugging + logical,parameter :: 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 + + ! Constants used to define C3 versus C4 photosynth pathways + integer, public, parameter :: c3_path_index = 1 + integer, public, parameter :: c4_path_index = 0 + + + ! Constants used to define conductance models + integer, parameter :: medlyn_model = 2 + integer, parameter :: ballberry_model = 1 + + ! Alternatively, Gross Assimilation can be used to estimate + ! leaf co2 partial pressure and therefore conductance. The default + ! is to use anet + integer, parameter :: net_assim_model = 1 + integer, parameter :: gross_assim_model = 2 + + ! Constants defining the photosynthesis temperature acclimation model + integer, parameter :: photosynth_acclim_model_none = 1 + integer, parameter :: photosynth_acclim_model_kumarathunge_etal_2019 = 2 + + ! Rdark constants from Atkin et al., 2017 https://doi.org/10.1007/978-3-319-68703-2_6 + ! and Heskel et al., 2016 https://doi.org/10.1073/pnas.1520282113 + real(r8), parameter :: lmr_b = 0.1012_r8 ! (degrees C**-1) + real(r8), parameter :: lmr_c = -0.0005_r8 ! (degrees C**-2) + real(r8), parameter :: lmr_TrefC = 25._r8 ! (degrees C) + + ! These two are public for error checking during parameter read-in + real(r8), parameter, public :: lmr_r_1 = 0.2061_r8 ! (umol CO2/m**2/s / (gN/(m2 leaf))) + real(r8), parameter, public :: lmr_r_2 = -0.0402_r8 ! (umol CO2/m**2/s/degree C) + + ! Fraction of light absorbed by non-photosynthetic pigments + real(r8),parameter :: fnps = 0.15_r8 + + ! term accounting that two photons are needed to fully transport a single + ! electron in photosystem 2 + real(r8), parameter :: photon_to_e = 0.5_r8 + + ! empirical curvature parameter for electron transport rate + real(r8),parameter :: theta_psii = 0.7_r8 + + + ! curvature parameter for quadratic smoothing of C4 gross assimilation + real(r8),parameter :: theta_ip_c4 = 0.95_r8 !0.95 is from Collatz 91, 0.999 was api 36 + real(r8),parameter :: theta_cj_c4 = 0.98_r8 !0.98 from Collatz 91, 0.099 was api 36 + + + ! There is a minimum stomatal conductance, below which we just don't + ! allow, this is well below reasonable ranges of cuticular conductance + real(r8),parameter :: gs0_min = 100.0_r8 + + + ! Set this to true to match results with main + logical, parameter :: base_compare_revert = .false. + + + ! 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 + + + ! There seems little evidence that jmax and vcmax ever truly reach zero + ! Due to leaf nitrogen concentration gradients, and water stress, + ! model estimated values of vcmax and jmax could reach zero. Until + ! a better, perhaps asymtotic method of preventing zero values is generated + ! + real(r8),parameter :: min_vcmax_frac = 0.10_r8 + real(r8),parameter :: min_jmax_frac = 0.10_r8 + + + ! The stomatal slope can either be scaled by btran or not. FATES had + ! a precedent of using this into 2024, but discussions here: https://github.com/NGEET/fates/issues/719 + ! suggest we should try other hypotheses + + integer, parameter,public :: btran_on_gs_none = 0 ! do not apply btran to conductance terms + integer, parameter,public :: btran_on_gs_gs0 = 1 ! apply btran to stomatal intercept (API 36.1) + integer, parameter,public :: btran_on_gs_gs1 = 2 ! apply btran to stomatal slope + integer, parameter,public :: btran_on_gs_gs01 = 3 ! apply btran to both stomatal slope and intercept + integer, parameter,public :: btran_on_gs_gs2 = 4 ! apply btran to the whole non-intercept portion + ! of conductance equation. (NOTE! for Medlyn, + ! this is different than btran_on_gs_gs1, + ! for Ball-Berry, this is the SAME as + ! btran_on_gs_gs1 + integer, parameter,public :: btran_on_gs_gs02 = 5 ! same as btran_on_gs_gs2, but also apply to the intercept + + integer, parameter,public :: btran_on_ag_none = 0 ! do not apply btran to vcmax or jmax + integer, parameter,public :: btran_on_ag_vcmax = 1 ! apply btran to vcmax (API 36.1) + integer, parameter,public :: btran_on_ag_vcmax_jmax = 2 ! apply btran to vcmax and jmax + + ! These are parameter constants read in externally + ! some are differentiated by PFT, others are not + ! ------------------------------------------------------------------------------------- + + type, public :: leafbiophys_params_type + + integer, allocatable :: c3psn(:) ! Photosynthetic pathway index (C3=1,C4=0,etc) + real(r8),allocatable :: medlyn_slope(:) ! Stomatal Slope, Medlyn, e.g. g1 [-] + real(r8),allocatable :: bb_slope(:) ! Stomatal Slope, Ball-Berry, e.g. g1 [-] + real(r8),allocatable :: stomatal_intercept(:) ! Stomatal int, BB or Medlyn, e.g. g0, [-] + real(r8),allocatable :: maintresp_leaf_ryan1991_baserate(:) ! Base maintenance resp rate M.Ryan 1991 [gC gN-1 s-1] + real(r8),allocatable :: maintresp_leaf_atkin2017_baserate(:) ! Base maintenance resp rate Atkin 2017 [umol CO2 m-2 s-1] + real(r8),allocatable :: maintresp_reduction_curvature(:) ! curvature of MR reduction as f(carbon storage), + ! 1=linear, 0=very curved + real(r8),allocatable :: maintresp_reduction_intercept(:) ! intercept of MR reduction as f(carbon storage), + ! 0=no throttling, 1=max throttling + real(r8),allocatable :: maintresp_reduction_upthresh (:) ! Upper threshold for storage biomass (relative + ! to leaf biomass) above which MR is not reduced + real(r8),allocatable :: vcmaxha(:) ! activation energy for vcmax (J/mol) + real(r8),allocatable :: jmaxha(:) ! activation energy for jmax (J/mol) + real(r8),allocatable :: vcmaxhd(:) ! deactivation energy for vcmax (J/mol) + real(r8),allocatable :: jmaxhd(:) ! deactivation energy for jmax (J/mol) + real(r8),allocatable :: vcmaxse(:) ! entropy term for vcmax (J/mol/K) + real(r8),allocatable :: jmaxse(:) ! entropy term for jmax (J/mol/K) + integer :: dayl_switch ! switch for turning on or off day length factor + ! scaling for photosynthetic parameters + integer :: photo_tempsens_model ! switch for choosing the model that defines the temperature + ! sensitivity of photosynthetic parameters (vcmax, jmax). + ! 1=non-acclimating, 2=Kumarathunge et al., 2019 + integer :: stomatal_assim_model ! Switch designating whether to use net or + ! gross assimilation in the stomata model + integer :: stomatal_model ! switch for choosing between stomatal conductance models, + ! for Ball-Berry, 2 for Medlyn + integer,allocatable :: stomatal_btran_model(:) ! index for how btran effects conductance + ! 0: btran does not scale the stomatal slope or intercept + ! 1: btran scales the stomatal intercept only + ! 2: btran scales the stomatal slope only + ! 3: btran scales both stomatal slope and intercept + integer,allocatable :: agross_btran_model(:) ! index for how btran scales gross assimilation processes + ! 0: btran does not scale vcmax or jmax + ! 1: btran scales only vcmax + ! 2: btran scales both vcmax and jmax + + ! ------------------------------------------------------------------------------------- + ! Note the omission of several parameter constants: + ! + ! Vcmax at 25C, canopy top + ! Jmax at 25C, at canopy + ! Kp (initial co2 response slope) at 25C, canopy top + ! + ! These parameters are omitted because some models (like FATES) have these as functions + ! of leaf age. By setting these as arguments in our routines, we add a degree of + ! flexibility. So these routines can be run on mixed-age leaves if desired. + !-------------------------------------------------------------------------------------- + + end type leafbiophys_params_type + + + type(leafbiophys_params_type),public :: lb_params + + + + ! A possible sequence of calls for leaf biophysics is as follows: + ! 1) determine any gas quantities or parameters that are derived and + ! are applicable to a super-set of leaf-layers (like MM, and compensation points) + ! 2) loop over discrete portions of leaves, perhaps differentiated by spatial position and pft + ! 3) determine if this particular leaf is present and has been solved + ! 4) solve for leaf maintenance respiration (e.g. dark respiration) + ! 5) update leaf-level rates, such as vcmax, jmax + ! 6) solve for photosynthesis + ! 7) solve for maintenance respiration of other tissues (other library?) + +contains + + subroutine StomatalCondMedlyn(anet,veg_esat,can_vpress,gs0,gs1,gs2,leaf_co2_ppress,can_press,gb,gs) + + ! ----------------------------------------------------------------------------------- + ! Calculate stomatal conductance (gs [umol/m2/s]) via the Medlyn approach, described in: + ! Medlyn et al. Reconciling the optimal and empirical approaches to modelling stomatal + ! conductance. Global Change Biology (2011) 17, 2134–2144, doi: 10.1111/j.1365-2486.2010.02375.x + ! + ! Original implementation in FATES is described in: + ! Li, Q. and Serbin, S. P. and Lamour, J. and Davidson, K. J. and Ely, K. S. and Rogers, A. + ! Implementation and evaluation of the unified stomatal optimization approach in the Functionally + ! Assembled Terrestrial Ecosystem Simulator (FATES). Geoscientific Model Development, 15(11), 2022. + ! doi: 10.5194/gmd-15-4313-2022. + ! + ! The implementation, and adaptation to include a leaf boundary layer in serial + ! resitance was taken from CLM5.0 + ! ----------------------------------------------------------------------------------- + + + ! Input + real(r8), intent(in) :: anet ! net leaf photosynthesis (umol CO2/m**2/s) + real(r8), intent(in) :: veg_esat ! saturation vapor pressure at veg_tempk (Pa) + real(r8), intent(in) :: can_press ! Air pressure NEAR the surface of the leaf (Pa) + real(r8), intent(in) :: gb ! leaf boundary layer conductance (umol /m**2/s) + real(r8), intent(in) :: can_vpress ! vapor pressure of canopy air (Pa) + real(r8), intent(in) :: gs0,gs1,gs2 ! stomatal intercept, and two different + ! mutually exclusive slopes + ! gs1 is a slope that is possibly multiplied by btran + ! gs2 is an alternate location for btran scaling + ! that applies to the non-intercept portion of + ! the equation + real(r8), intent(in) :: leaf_co2_ppress ! CO2 partial pressure at the leaf surface (Pa) + ! at the boundary layer interface with the leaf + + ! Output + real(r8),intent(out) :: gs ! leaf stomatal conductance (umol H2O/m**2/s) + + ! locals + real(r8) :: vpd ! water vapor deficit in Medlyn stomatal model [KPa] + real(r8) :: term ! intermediate term used to simplify equations + real(r8) :: aquad,bquad,cquad ! quadradic solve terms + real(r8) :: r1,r2 ! quadradic solve roots + real(r8),parameter :: min_vpd_pa = 50._r8 ! Minimum allowable vpd [Pa] + real(r8),parameter :: anet_min = 0.001_r8 + logical :: err + + ! Evaluate trival solution, if there is no positive net assimiolation + ! the stomatal conductance is the intercept conductance + if (anet <= anet_min) then + gs = gs0 + return + end if + + ! Trivial case (gs2 near 0) + if(gs2<0.01_r8) then + gs = gs0 + return + end if + + ! Trivial case (gs1 near 0) + if(gs1<0.01_r8)then + gs = gs0 + h2o_co2_stoma_diffuse_ratio * anet/(leaf_co2_ppress/can_press) + return + end if + + vpd = max((veg_esat - can_vpress), min_vpd_pa) * kpa_per_pa + term = gs2 * h2o_co2_stoma_diffuse_ratio * anet / (leaf_co2_ppress / can_press) + aquad = 1.0_r8 + bquad = -(2.0 * (gs0 + term) + (gs1 * term)**2 / (gb * vpd )) + cquad = gs0*gs0 + (2.0*gs0 + term * & + (1.0 - gs1*gs1 / vpd)) * term + + call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + gs = max(r1,r2) + + if(err)then + write(fates_log(),*) "medquadfail:",anet,veg_esat,can_vpress,gs0,gs1,gs2,leaf_co2_ppress,can_press + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! RGK: re-derived solution to check units. + ! For 50*10*10*5*50 operations in unit tests the original method is + ! about 2% faster. Alternative derivation below. + + ! vpd = max((veg_esat - can_vpress), min_vpd_pa) + ! term = h2o_co2_stoma_diffuse_ratio * anet * can_press / leaf_co2_ppress + ! term2 = pa_per_kpa*(gs1*term)**2._r8 + ! aquad = gb*vpd + ! bquad = -2*gb*vpd*(gs0 + term) - term2 + ! cquad = gb*vpd*(gs0 + term)**2._r8 - gb*term2 + + ! call QuadraticRoots(aquad, bquad, cquad, r3, r4) + + ! if( abs(max(r1,r2)-max(r3,r4))>1.e-5 ) then + ! print*,"Math check failed",r1,r2,r3,r4,anet,veg_esat,can_vpress,gs0,gs1,leaf_co2_ppress,can_press,gb + ! stop + ! end if + + return + end subroutine StomatalCondMedlyn + + ! ======================================================================================= + + subroutine StomatalCondBallBerry(a_gs,a_net,veg_esat,can_vpress,gs0,gs1,leaf_co2_ppress,can_press, gb, gs) + + + ! Input + real(r8), intent(in) :: veg_esat ! saturation vapor pressure at veg_tempk (Pa) + real(r8), intent(in) :: can_press ! Air pressure near the surface of the leaf (Pa) + real(r8), intent(in) :: gb ! leaf boundary layer conductance (umol /m**2/s) + real(r8), intent(in) :: can_vpress ! vapor pressure of the canopy air (Pa) + real(r8), intent(in) :: gs0,gs1 ! Stomatal intercept and slope (umol H2O/m**2/s) + real(r8), intent(in) :: leaf_co2_ppress ! CO2 partial pressure at leaf surface (Pa) + real(r8), intent(in) :: a_gs ! The assimilation (a) for calculating conductance (gs) + ! is either = to anet or agross + real(r8), intent(in) :: a_net ! Net assimilation rate of co2 (umol/m2/s) + ! Output + real(r8), intent(out) :: gs ! leaf stomatal conductance (umol H2O/m**2/s) + + ! locals + real(r8) :: ceair ! constrained vapor pressure (Pa) + real(r8) :: aquad,bquad,cquad ! quadradic solve terms + real(r8) :: r1,r2 ! quadradic solve roots + logical :: err + + if (a_gs <= nearzero) then + gs = gs0 + return + end if + + ! Trivial case (gs1 near 0) + if(gs1<0.01_r8)then + gs = gs0 + return + end if + + ! Apply a constraint to the vapor pressure + ceair = GetConstrainedVPress(can_vpress,veg_esat) + !ceair = can_vpress + + aquad = leaf_co2_ppress + bquad = leaf_co2_ppress*(gb - gs0) - gs1 * a_gs * can_press + + if(.not.base_compare_revert) then + cquad = -gb*(leaf_co2_ppress * gs0 + gs1 * a_gs * can_press * ceair/ veg_esat ) + else + cquad = -gb*(leaf_co2_ppress * gs0 + gs1 * a_net* can_press * ceair/ veg_esat ) + end if + + call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + if(err)then + write(fates_log(),*) "bbquadfail:",a_net,a_gs,veg_esat,can_vpress,gs0,gs1,leaf_co2_ppress,can_press + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + + gs = max(r1,r2) + + return + end subroutine StomatalCondBallBerry + + ! ===================================================================================== + + function AgrossRubiscoC3(vcmax,ci,can_o2_ppress,co2_cpoint,mm_kco2,mm_ko2) result(ac) + + ! Input + real(r8) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) + real(r8) :: ci ! intracellular leaf CO2 (Pa) + real(r8) :: co2_cpoint ! CO2 compensation point (Pa) + real(r8) :: can_o2_ppress ! Partial pressure of O2 near the leaf surface (Pa) + real(r8) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) + real(r8) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) + + ! Output + real(r8) :: ac ! Rubisco-limited gross photosynthesis (umol CO2/m**2/s) + + ac = vcmax * max(ci-co2_cpoint, 0._r8) / & + (ci+mm_kco2 * (1._r8+can_o2_ppress / mm_ko2 )) + + end function AgrossRubiscoC3 + + ! ===================================================================================== + + function GetJe(par_abs,jmax) result(je) + + ! Input + real(r8) :: par_abs ! Absorbed PAR per leaf area [umol photons/m**2/s] + real(r8) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8) :: je ! electron transport rate (umol electrons/m**2/s) + real(r8) :: aquad,bquad,cquad ! terms for quadratic equations + real(r8) :: r1,r2 ! roots of quadratic equation + real(r8) :: jpar ! absorbed photons in photocenters as an electron + ! rate (umol electrons/m**2/s) + logical :: err + + ! Electron transport rate for C3 plants. + ! Convert absorbed photon density [umol/m2 leaf /s] to + ! that absorbed only by the photocenters (fnps) and also + ! convert from photon energy into electron transport rate (photon_to_e) + jpar = par_abs*photon_to_e*(1.0_r8 - fnps) + + ! convert the absorbed par into absorbed par per m2 of leaf, + ! so it is consistant with the vcmax and lmr numbers. + aquad = theta_psii + bquad = -(jpar + jmax) + cquad = jpar * jmax + call QuadraticRoots(aquad, bquad, cquad, r1, r2, err) + if(err)then + write(fates_log(),*) "jequadfail:",par_abs,jpar,jmax + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + je = min(r1,r2) + + end function GetJe + + ! ===================================================================================== + + function AgrossRuBPC3(par_abs,jmax,ci,co2_cpoint) result(aj) + + ! Input + real(r8) :: par_abs ! Absorbed PAR per leaf area [umol photons/m2leaf/s ] + real(r8) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8) :: ci ! intracellular leaf CO2 (Pa) + real(r8) :: co2_cpoint ! CO2 compensation point (Pa) + + ! Output + real(r8) :: aj ! RuBP-limited gross photosynthesis (umol CO2/m**2/s) + + ! locals + real(r8) :: je ! actual electron transport rate (umol electrons/m**2/s) + + ! Get the smoothed (quadratic between J and Jmax) electron transport rate + je = GetJe(par_abs,jmax) + + + aj = je * max(ci-co2_cpoint, 0._r8) / & + (4._r8*ci+8._r8*co2_cpoint) + + + end function AgrossRuBPC3 + + ! ======================================================================================= + + function AgrossRuBPC4(par_abs) result(aj) + + real(r8) :: par_abs ! Absorbed PAR per leaf area [umol photons/m2leaf/s ] + real(r8) :: aj ! RuBP-limited gross photosynthesis (umol CO2/m**2/s) + + ! quantum efficiency, used only for C4 (mol CO2 / mol photons) + real(r8),parameter :: c4_quant_eff = 0.05_r8 + + aj = c4_quant_eff*par_abs*photon_to_e*(1.0_r8 - fnps) + + end function AgrossRuBPC4 + + ! ======================================================================================= + + function AgrossPEPC4(ci,kp,can_press) result(ap) + + real(r8) :: ci ! intracellular leaf CO2 (Pa) + real(r8) :: kp ! initial co2 response slope + real(r8) :: can_press ! Air pressure near the surface of the leaf (Pa) + real(r8) :: ap ! PEP limited gross assimilation rate (umol co2/m2/s) + + ap = kp * max(ci, 0._r8) / can_press + + end function AgrossPEPC4 + + ! ======================================================================================= + + + subroutine CiMinMax(ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, & + can_co2_ppress,can_o2_ppress,can_press,lmr,par_abs, & + gb,gs0,ci_min,ci_max) + + ! This routine is used to find the first values of Ci that are the bounds + ! for the bisection algorithm. It finds the values associated with minimum + ! and maximum conductance when equating the source and sink limitations + ! on net photosynthesis. + + ! input + integer, intent(in) :: ft ! plant functional type index + real(r8), intent(in) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) + real(r8), intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8), intent(in) :: kp ! initial slope of CO2 response curve (C4 plants) + real(r8), intent(in) :: co2_cpoint ! CO2 compensation point (Pa) + real(r8), intent(in) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) + real(r8), intent(in) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) + 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) :: can_press ! Air pressure near the surface of the leaf (Pa) + real(r8), intent(in) :: lmr ! leaf maintenance respiration rate (umol CO2/m**2/s) + real(r8), intent(in) :: par_abs ! par absorbed per unit leaf area [umol photons/m2leaf/s ] + real(r8), intent(in) :: gb ! leaf boundary layer conductance (umol H2O/m**2/s) + real(r8), intent(in) :: gs0 ! stomatal intercept + + ! output + real(r8), intent(out) :: ci_max ! Intracellular Co2 at maximum conductance [Pa] + real(r8), intent(out) :: ci_min ! Intracellular CO2 at minimum conductance [Pa] + + ! Intermediate terms + real(r8), dimension(3) :: ci ! Ci for Rubisco,RuBP and PEP respectively + real(r8), dimension(3) :: ag ! Agross for Rubisco, Rubp and PEP + + ! These are compound terms used to solve the equation that balances + ! net assimilation with the gradient flux equation + real(r8) :: a,b,c,d,e,f,g + real(r8) :: ap,ac,aj,ai + real(r8) :: Je ! Electron tranport rate (umol e/m2/s) + real(r8) :: rmin ,rmax ! Maximum and minimum resistance [s/umol h2o/m2] + real(r8) :: aquad,bquad,cquad,r1,r2 + logical :: err + + ! Minimum possible resistance (with a little buffer) + rmin = 0.75_r8*h2o_co2_bl_diffuse_ratio/gb + ! Maximum possible resistance (with a littel buffer) + rmax = 1.25_r8*(h2o_co2_bl_diffuse_ratio/gb + h2o_co2_stoma_diffuse_ratio/gs0) + + if (lb_params%c3psn(ft) == c4_path_index)then + + ! Maximum conductance (minimum resistance) + + ag(1) = vcmax + ag(2) = AgrossRuBPC4(par_abs) + + ! C4: Rubisco-limited photosynthesis + ac = vcmax + + ! C4: RuBP-limited photosynthesis + aj = AgrossRuBPC4(par_abs) + + + + aquad = theta_cj_c4 + bquad = -(ac + aj) + cquad = ac * aj + call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + if(err)then + write(fates_log(),*) "c41quadfail_minmax1:",par_abs,can_press,vcmax + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + ai = min(r1,r2) + + a = rmin*can_press + + aquad = theta_ip_c4/a**2.0_r8 + kp/(can_press*a) + bquad = - theta_ip_c4*2.0_r8*can_co2_ppress/(a*a) & + - theta_ip_c4*2._r8*lmr/a & + - kp*lmr/can_press & + + ai/a & + - kp*can_co2_ppress/(can_press*a) & + + ai*kp/can_press + cquad = theta_ip_c4*can_co2_ppress*can_co2_ppress/(a*a) & + + 2._r8*lmr*can_co2_ppress*theta_ip_c4/a & + + theta_ip_c4*lmr*lmr & + - ai*can_co2_ppress/a & + - ai*lmr + + call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + if(err)then + write(fates_log(),*) "c41quadfail_minmax2:",par_abs,can_press,vcmax + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ci(3) = max(r1,r2) + + !! C4: PEP carboxylase-limited (CO2-limited) + ap = AgrossPEPC4(ci(3),kp,can_press) + + aquad = theta_ip_c4 + bquad = -(ai + ap) + cquad = ai * ap + call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + if(err)then + write(fates_log(),*) "c42quadfail:",par_abs,ci,kp,can_press,vcmax + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + ag(:) = min(r1,r2) + + if(debug) then + if( abs(ci(3)-(can_co2_ppress - (ag(3)-lmr)*can_press*rmin))>0.001_r8) then + write(fates_log(),*) "c4 ci conv check fail" + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + + ci(1) = can_co2_ppress - (ag(1)-lmr)*can_press*rmin + ci(2) = can_co2_ppress - (ag(1)-lmr)*can_press*rmin + + ci_max = ci(minloc(ag,DIM=1)) + + ! Minimum conductance (maximum resistance) + + a = rmax*can_press + + aquad = theta_ip_c4/a**2.0_r8 + kp/(can_press*a) + bquad = - theta_ip_c4*2.0_r8*can_co2_ppress/(a*a) & + - theta_ip_c4*2._r8*lmr/a & + - kp*lmr/can_press & + + ai/a & + - kp*can_co2_ppress/(can_press*a) & + + ai*kp/can_press + cquad = theta_ip_c4*can_co2_ppress*can_co2_ppress/(a*a) & + + 2._r8*lmr*can_co2_ppress*theta_ip_c4/a & + + theta_ip_c4*lmr*lmr & + - ai*can_co2_ppress/a & + - ai*lmr + + call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + if(err)then + write(fates_log(),*) "c41quadfail_minmax2:",par_abs,can_press,vcmax + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + ci(3) = max(r1,r2) + + ap = AgrossPEPC4(ci(3),kp,can_press) + + aquad = theta_ip_c4 + bquad = -(ai + ap) + cquad = ai * ap + call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + if(err)then + write(fates_log(),*) "c42quadfail:",par_abs,ci,kp,can_press,vcmax + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + ag(:) = min(r1,r2) + + ci(1) = can_co2_ppress - (ag(1)-lmr)*can_press*rmax + ci(2) = can_co2_ppress - (ag(1)-lmr)*can_press*rmax + + if(debug) then + if( abs(ci(3)-(can_co2_ppress - (ag(3)-lmr)*can_press*rmax))>0.001_r8) then + write(fates_log(),*) "c4 ci conv check fail" + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + ci_min = ci(minloc(ag,DIM=1)) + + return + end if + + ! Get the maximum e tranport rate for when we solve for RuBP (twice) + Je = GetJe(par_abs,jmax) + + ! Find ci at maximum conductance (1/inf = 0) + + a = can_co2_ppress + b = can_press*rmin + c = vcmax + d = vcmax*co2_cpoint + e = 1._r8 + f = mm_kco2*(1._r8+can_o2_ppress / mm_ko2 ) + g = lmr + ci(1) = CiFromAnetDiffGrad(a,b,c,d,e,f,g) + ag(1) = AgrossRubiscoC3(vcmax,ci(1),can_o2_ppress,co2_cpoint,mm_kco2,mm_ko2) + + + ! check the math + if(debug)then + if ( abs((can_co2_ppress-ci(1))/b - (ag(1)-lmr) ) > 1.e-3_r8 ) then + write(fates_log(),*) 'incorrect ci, max cond:',ci(1),(can_co2_ppress-ci(1))/b,ag(1)-lmr + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + c = Je + d = Je*co2_cpoint + e = 4._r8 + f = 8._r8*co2_cpoint + g = lmr + ci(2) = CiFromAnetDiffGrad(a,b,c,d,e,f,g) + ag(2) = AgrossRuBPC3(par_abs,jmax,ci(2),co2_cpoint) + + if(debug)then + if ( abs((can_co2_ppress-ci(2))/b - (ag(2)-lmr)) > 1.e-3_r8 ) then + write(fates_log(),*)'incorrect ci, max cond:',ci(2),(can_co2_ppress-ci(2))/b,ag(2)-lmr + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + ! Take the Ci that matches the minimizing gross assimilation + ci_max = ci(minloc(ag(1:2),DIM=1)) + + ! Find ci at minimum conductance (1/g0) (max conductance) + + a = can_co2_ppress + b = can_press*rmax + c = vcmax + d = vcmax*co2_cpoint + e = 1._r8 + f = mm_kco2*(1._r8+can_o2_ppress / mm_ko2 ) + g = lmr + ci(1) = CiFromAnetDiffGrad(a,b,c,d,e,f,g) + ag(1) = AgrossRubiscoC3(vcmax,ci(1),can_o2_ppress,co2_cpoint,mm_kco2,mm_ko2) + + ! check the math + if(debug)then + if ( abs((can_co2_ppress-ci(1))/b - (ag(1)-lmr) ) > 1.e-3_r8 ) then + write(fates_log(),*)'incorrect ci, min cond:',ci(1),(can_co2_ppress-ci(1))/b,ag(1)-lmr + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + c = Je + d = Je*co2_cpoint + e = 4._r8 + f = 8._r8*co2_cpoint + g = lmr + ci(2) = CiFromAnetDiffGrad(a,b,c,d,e,f,g) + ag(2) = AgrossRuBPC3(par_abs,jmax,ci(2),co2_cpoint) + + if(debug)then + if ( abs((can_co2_ppress-ci(2))/b -(ag(2)-lmr)) > 1.e-3_r8 ) then + write(fates_log(),*)'incorrect ci, min cond:',ci(2),(can_co2_ppress-ci(2))/b,ag(2)-lmr + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + ci_min = ci(minloc(ag(1:2),DIM=1)) + + return + end subroutine CiMinMax + + ! ===================================================================================== + + function CiFromAnetDiffGrad(a,b,c,d,e,f,g) result(ci) + + ! When the equation for net photosynthesis is coupled with the diffusion + ! flux equation, for both Rubisco and RuBP limited assimilation + ! the form is like so: + ! + ! (a-ci)/b = (c*ci - d)/(e*ci + f) - g + ! + ! The expression below simply solves for ci + ! This function is called to find endpoints, where conductance is maximized + ! and minimized, to perform a binary search + + real(r8) :: a,b,c,d,e,f,g ! compound terms to solve the coupled Anet + ! and diffusive flux gradient equations + real(r8) :: ci ! intracellular co2 [Pa] + real(r8) :: r1,r2 ! roots for quadratic + real(r8) :: aquad,bquad,cquad + logical :: err + + aquad = -(e/b) + bquad = (e*a-f)/b + e*g - c + cquad = (f*a)/b + f*g + d + call QuadraticRoots(aquad, bquad, cquad, r1, r2, err) + + ci = max(r1,r2) + + end function CiFromAnetDiffGrad + + ! ======================================================================================= + + subroutine CiFunc(ci, & + ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, & + can_co2_ppress,can_o2_ppress,can_press,can_vpress,lmr,par_abs,gb,veg_tempk, & + gs0, gs1, gs2, & + anet,agross,gs,fval) + + ! ----------------------------------------------------------------------------------- + ! DESCRIPTION: + ! The photosynthesis solver must solve 3 equations simultaneously, + ! Anet, Stomatal Conductance and the update to Ci. To solve this, + ! we search values of Ci that satisfy the equations, and create an evaluation + ! function "fval" that is the difference between the input value and the predicted + ! value of ci. + ! + ! fval = ci_input - ci_updated + ! + ! The updated ci: + ! ci = ca - (1.37rb+1.65rs))*patm*anet + ! + ! Thus: + ! fval = ci_input - [ca-(1.37rb+1.65rs))*patm*anet] + ! + ! This method borrows from Jinyun Tang's 2011 code for CLM + ! ----------------------------------------------------------------------------------- + + + real(r8), intent(in) :: ci ! Input (trial) intracellular leaf CO2 (Pa) + integer, intent(in) :: ft ! plant functional type index + real(r8), intent(in) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) + real(r8), intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8), intent(in) :: kp ! initial slope of CO2 response curve (C4 plants) + real(r8), intent(in) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) + real(r8), intent(in) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) + real(r8), intent(in) :: co2_cpoint ! CO2 compensation point (Pa) + real(r8), intent(in) :: can_press ! Air pressure near the surface of the leaf (Pa) + 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) :: can_vpress ! vapor pressure of the canopy air (Pa) + real(r8), intent(in) :: lmr ! leaf maintenance respiration rate (umol CO2/m**2/s) + real(r8), intent(in) :: par_abs ! par absorbed per unit lai [umol photons/m2leaf/s ] + real(r8), intent(in) :: gb ! leaf boundary layer conductance (umol H2O/m**2/s) + real(r8), intent(in) :: veg_tempk ! vegetation temperature + real(r8), intent(in) :: gs0 ! conductance intercept (umol H20/m2/s) + real(r8), intent(in) :: gs1 ! conductance slope (could be multiplied by btran) + real(r8), intent(in) :: gs2 ! for Medlyn only: either 1 or btran + real(r8), intent(out) :: anet ! net leaf photosynthesis (umol CO2/m**2/s) + real(r8), intent(out) :: agross ! gross leaf photosynthesis (umol CO2/m**2/s) + real(r8), intent(out) :: gs ! stomatal conductance (umol h2o/m2/s) + real(r8), intent(out) :: fval ! ci_input - ci_updated (Pa) + + real(r8) :: veg_esat ! Saturation vapor pressure at leaf-surface [Pa] + real(r8) :: veg_qs ! DUMMY, specific humidity at leaf-surface [kg/kg] + real(r8) :: a_gs ! The assimilation (a) for calculating conductance (gs) + ! is either = to anet or agross + real(r8) :: ac ! Rubisco-limited gross photosynthesis (umol CO2/m**2/s) + real(r8) :: aj ! RuBP-limited gross photosynthesis (umol CO2/m**2/s) + real(r8) :: ap ! product-limited (C3) or CO2-limited + ! (C4) gross photosynthesis (umol CO2/m**2/s) + real(r8) :: leaf_co2_ppress ! CO2 partial pressure at leaf surface (Pa) + real(r8) :: ai ! For C4, smoothed co-limited assimilation of Rubisco/RuBP + real(r8) :: aquad,bquad,cquad ! a,b,c terms in the quadratic equation + real(r8) :: r1,r2 ! The roots from the quadratic + logical :: err + !------------------------------------------------------------------------------ + + ! Photosynthesis limitation rate calculations + if (lb_params%c3psn(ft) == c3_path_index)then + + ! C3: Rubisco-limited photosynthesis + ac = AgrossRubiscoC3(vcmax,ci,can_o2_ppress,co2_cpoint,mm_kco2,mm_ko2) + + ! C3: RuBP-limited photosynthesis + aj = AgrossRuBPC3(par_abs,jmax,ci,co2_cpoint ) + + if(base_compare_revert) then + ! Gross photosynthesis smoothing calculations. Co-limit ac and aj. + aquad = 0.999_r8 + bquad = -(ac + aj) + cquad = ac * aj + call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + agross = min(r1,r2) + else + ! Take the minimum, no smoothing + agross = min(ac,aj) + end if + + else + + ! C4: Rubisco-limited photosynthesis + ac = vcmax + + ! C4: RuBP-limited photosynthesis + aj = AgrossRuBPC4(par_abs) + + ! C4: PEP carboxylase-limited (CO2-limited) + ap = AgrossPEPC4(ci,kp,can_press) + + aquad = theta_cj_c4 + bquad = -(ac + aj) + cquad = ac * aj + call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + if(err)then + write(fates_log(),*) "c41quadfail:",par_abs,ci,kp,can_press,vcmax + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + ai = min(r1,r2) + + aquad = theta_ip_c4 + bquad = -(ai + ap) + cquad = ai * ap + call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + if(err)then + write(fates_log(),*) "c42quadfail:",par_abs,ci,kp,can_press,vcmax + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + agross = min(r1,r2) + + + end if + + ! Calculate anet + anet = agross - lmr + + if ( lb_params%stomatal_assim_model == gross_assim_model ) then + if ( lb_params%stomatal_model == medlyn_model ) then + write (fates_log(),*) 'Gross Assimilation conductance is incompatible with the Medlyn model' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + a_gs = agross + else + a_gs = anet + end if + + ! A note about the use of the quadratic equations for calculating stomatal conductance + ! ------------------------------------------------------------------------------------ + ! These two following models calculate the conductance between the intracellular leaf + ! space and the leaf surface, not the canopy air space. Transport between the leaf + ! surface and the canopy air space is governed by the leaf boundary layer conductance. + ! However, we need to estimate the properties at the surface of the leaf to solve for + ! the stomatal conductance. We do this by using Fick's law (gradient resistance + ! approximation of diffusion) to estimate the flux of water vapor across the + ! leaf boundary layer, and balancing that with the flux across the stomata. It + ! results in the following equation for leaf surface humidity: + ! + ! e_s = (e_i g_s + e_c g_b)/(g_b + g_s) + ! + ! The leaf surface humidity (e_s) becomes an expression of canopy humidity (e_c), + ! intracellular humidity (e_i, which is the saturation humidity at leaf temperature), + ! boundary layer conductance (g_b) (these are all known) and stomatal conductance + ! (g_s) (this is still unknown). This expression is substituted into the stomatal + ! conductance equation. The resulting form of these equations becomes a quadratic. + ! + ! For a detailed explanation, see the FATES technical note, section + ! "1.11 Stomatal Conductance" + ! + ! ------------------------------------------------------------------------------------ + + leaf_co2_ppress = can_co2_ppress - h2o_co2_bl_diffuse_ratio/gb * a_gs * can_press + leaf_co2_ppress = max(leaf_co2_ppress,1.e-06_r8) + + ! Determine saturation vapor pressure at the leaf surface, from temp and atm-pressure + call QSat(veg_tempk, can_press, veg_qs, veg_esat) + + if ( lb_params%stomatal_model == medlyn_model ) then + call StomatalCondMedlyn(anet,veg_esat,can_vpress,gs0,gs1,gs2, & + leaf_co2_ppress,can_press,gb,gs) + else + call StomatalCondBallBerry(a_gs,anet,veg_esat,can_vpress,gs0,gs1, & + leaf_co2_ppress,can_press,gb,gs) + end if + + ! Derive new estimate for ci + ! ci = can_co2_ppress - anet * can_press * & + ! (h2o_co2_bl_diffuse_ratio/gb + h2o_co2_stoma_diffuse_ratio/gs_out) + + ! fval = ci_input - ci_predicted + fval = ci - (can_co2_ppress - anet * can_press * & + (h2o_co2_bl_diffuse_ratio*gs + h2o_co2_stoma_diffuse_ratio*gb)/(gb*gs)) + + if(base_compare_revert) then + if (anet < 0._r8) then + fval = 0 + end if + end if + + end subroutine CiFunc + + ! ==================================================================================== + + subroutine CiBisection(ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, & + can_co2_ppress,can_o2_ppress,can_press,can_vpress,lmr,par_abs,gb,veg_tempk, & + gs0,gs1,gs2,ci_tol, & + anet,agross,gs,ci,solve_iter) + + ! ----------------------------------------------------------------------------------- + ! + ! Co-solve for photosynthesis and stomatal conductance via a bisectional search for + ! intracellular CO2. This is an inefficient yet robust search method that we resort + ! to if the recursive solver cannot find an answer within its convergence criteria. + ! ----------------------------------------------------------------------------------- + + integer, intent(in) :: ft ! plant functional type index + real(r8), intent(in) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) + real(r8), intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8), intent(in) :: kp ! initial slope of CO2 response curve (C4 plants) + real(r8), intent(in) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) + real(r8), intent(in) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) + real(r8), intent(in) :: co2_cpoint ! CO2 compensation point (Pa) + real(r8), intent(in) :: can_press ! Air pressure near the surface of the leaf (Pa) + 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) :: can_vpress ! vapor pressure of the canopy air (Pa) + real(r8), intent(in) :: lmr ! leaf maintenance respiration rate (umol CO2/m**2/s) + real(r8), intent(in) :: par_abs ! par absorbed per unit lai [umol photons/m2leaf/s ] + real(r8), intent(in) :: gb ! leaf boundary layer conductance (umol H2O/m**2/s) + real(r8), intent(in) :: veg_tempk ! vegetation temperature [Kelvin] + real(r8), intent(in) :: gs0,gs1,gs2 ! stomatal intercept and slope and alternative btran + real(r8), intent(in) :: ci_tol ! Convergence tolerance for solutions for intracellular CO2 (Pa) + real(r8), intent(out) :: anet ! net leaf photosynthesis (umol CO2/m**2/s) + real(r8), intent(out) :: agross ! gross leaf photosynthesis (umol CO2/m**2/s) + real(r8), intent(out) :: gs ! stomatal conductance (umol h2o/m2/s) + real(r8), intent(out) :: ci ! Input (trial) intracellular leaf CO2 (Pa) + integer, intent(inout) :: solve_iter ! number of bisections required + + ! With bisection, we need to keep track of three different ci values at any given time + ! The high, the low and the bisection. + + real(r8) :: ci_h, fval_h + real(r8) :: ci_l, fval_l + real(r8) :: ci_b, fval_b + + logical :: loop_continue ! Continue bisecting until tolerance reached + + ! Maximum number of iterations on intracelluar co2 solver until is quits + integer, parameter :: max_iters = 200 + + ! Find the starting points (end-points) for bisection + ! We dont need stomatal slope because we just want the two extremes + ! which is the intercept and infinite conductance + call CiMinMax(ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, & + can_co2_ppress,can_o2_ppress,can_press,lmr,par_abs, & + gb,gs0,ci_l,ci_h) + + call CiFunc(ci_h, & + ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, & + can_co2_ppress,can_o2_ppress,can_press,can_vpress,lmr,par_abs,gb,veg_tempk, & + gs0,gs1,gs2, & + anet,agross,gs,fval_h) + + call CiFunc(ci_l, & + ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, & + can_co2_ppress,can_o2_ppress,can_press,can_vpress,lmr,par_abs,gb,veg_tempk, & + gs0,gs1,gs2, & + anet,agross,gs,fval_l) + + ! It is necessary that our starting points are on opposite sides of the root + if( nint(fval_h/abs(fval_h)) .eq. nint(fval_l/abs(fval_l)) ) then + + ! Try an exteremly large bisection range, if this doesn't work, then + ! fail the run + ci_h = 0.1_r8 + call CiFunc(ci_h, & + ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, & + can_co2_ppress,can_o2_ppress,can_press,can_vpress,lmr,par_abs,gb,veg_tempk, & + gs0,gs1,gs2, & + anet,agross,gs,fval_h) + + ci_l = 2000._r8 + call CiFunc(ci_l, & + ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, & + can_co2_ppress,can_o2_ppress,can_press,can_vpress,lmr,par_abs,gb,veg_tempk, & + gs0,gs1,gs2, & + anet,agross,gs,fval_l) + + ! It is necessary that our starting points are on opposite sides of the root + if( nint(fval_h/abs(fval_h)) .eq. nint(fval_l/abs(fval_l)) ) then + write(fates_log(),*)'While attempting bisection for Ci calculations,' + write(fates_log(),*)'the two starting values for Ci were on the same' + write(fates_log(),*)'side of the root. Try increasing and decreasing' + write(fates_log(),*)'init_ci_high and init_ci_low respectively' + write(fates_log(),*) "ci_h=",ci_h,"fval_h=",fval_h,"ci_l=",ci_l,"fval_l=",fval_l + write(fates_log(),*) "ft= ",ft,"is c3psn:",lb_params%c3psn(ft) == c3_path_index + write(fates_log(),*) "vcmax=",vcmax,"jmax=",jmax,"kp=",kp + write(fates_log(),*) "co2_cpoint=",co2_cpoint,"mm_kco2=",mm_kco2,"mm_ko2=",mm_ko2 + write(fates_log(),*) "can_co2_ppress=",can_co2_ppress,"can_o2_ppress=",can_o2_ppress,"can_press=",can_press + write(fates_log(),*) "can_vpress=",can_vpress,"lmr=",lmr,"par_abs=",par_abs,"gb=",gb + write(fates_log(),*) "veg_tempk=",veg_tempk,"gs0=",gs0,"gs1=",gs1,"gs2=",gs2,"ci_tol=",ci_tol + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + + loop_continue = .true. + ci_b = 0.5*(ci_l+ci_h) + bi_iter_loop: do while(loop_continue) + + solve_iter = solve_iter + 1 + + call CiFunc(ci_b, & + ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, & + can_co2_ppress,can_o2_ppress,can_press,can_vpress,lmr,par_abs,gb,veg_tempk, & + gs0,gs1,gs2, & + anet,agross,gs,fval_b) + + if(abs(fval_b) 1.e-2)) then + write (fates_log(),*) 'qsat from QSat():', qsat_alt + write (fates_log(),*) 'qsat localy :', qsat_loc + write (fates_log(),*) 'values of qsat are too different' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Adjusting gs (compute a virtual gs) that will be passed to host model + + if ( qsat_adj < qs ) then + + ! if inner leaf vapor pressure is less then or equal to that at leaf surface + ! then set stomata resistance to be very large to stop the transpiration or back flow of vapor + rstoma_out = rsmax0 + + else + + rstoma_out = (qsat_loc-qs)*( 1/gstoma + rb)/(qsat_adj - qs)-rb + + end if + + if (rstoma_out < nearzero ) then + write (fates_log(),*) 'qsat:', qsat_loc, 'qs:', qs + write (fates_log(),*) 'LWP :', leaf_psi + write (fates_log(),*) 'pft :', ft + write (fates_log(),*) 'ceair:', ceair, 'veg_esat:', veg_esat + write (fates_log(),*) 'rstoma_out:', rstoma_out, 'rb:', rb + write (fates_log(),*) 'LWP_star', lwp_star + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + end function LeafHumidityStomaResis + + ! ===================================================================================== + + function ft1_f(tl, ha) result(ans) + ! + !!DESCRIPTION: + ! photosynthesis temperature response + ! + ! This function, along with the fth1_f() function, scales the Vcmax25 and Jmax25 + ! constants so they are adjusted for leaf temperature. + ! This equation does not have a name, but is equation 9.10 in the CLM5.0 technical + ! manual. + ! + ! !REVISION HISTORY + ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 + ! 7/23/16: Copied over from CLM by Ryan Knox + ! + !!USES + + + ! + ! !ARGUMENTS: + real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temperature function (K) + real(r8), intent(in) :: ha ! activation energy in photosynthesis temperature function (J/mol) + ! + ! !LOCAL VARIABLES: + real(r8) :: ans + !------------------------------------------------------------------------------- + + ans = exp( ha / (rgas_J_K_mol*(tfrz+25._r8)) * (1._r8 - (tfrz+25._r8)/tl) ) + + return + end function ft1_f + + ! ===================================================================================== + + function fth_f(tl,hd,se,scaleFactor) result(ans) + ! + !!DESCRIPTION: + !photosynthesis temperature inhibition + ! + ! This function, along with the ft1_f() function, scales the Vcmax25 and Jmax25 + ! constants so they are adjusted for leaf temperature. + ! This equation does not have a name, but is equation 9.10 in the CLM5.0 technical + ! manual. + ! !REVISION HISTORY + ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 + ! 7/23/16: Copied over from CLM by Ryan Knox + ! + + ! + ! !ARGUMENTS: + real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temp function (K) + real(r8), intent(in) :: hd ! deactivation energy in photosynthesis temp function (J/mol) + real(r8), intent(in) :: se ! entropy term in photosynthesis temp function (J/mol/K) + real(r8), intent(in) :: scaleFactor ! scaling factor for high temp inhibition (25 C = 1.0) + ! + ! !LOCAL VARIABLES: + real(r8) :: ans + !------------------------------------------------------------------------------- + + ans = scaleFactor / ( 1._r8 + exp( (-hd+se*tl) / (rgas_J_K_mol*tl) ) ) + + + + return + end function fth_f + + ! ===================================================================================== + + function fth25_f(hd,se)result(ans) + ! + !!DESCRIPTION: + ! scaling factor for photosynthesis temperature inhibition + ! + ! !REVISION HISTORY: + ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 + ! 7/23/16: Copied over from CLM by Ryan Knox + ! + !!USES + + + ! + ! !ARGUMENTS: + real(r8), intent(in) :: hd ! deactivation energy in photosynthesis temp function (J/mol) + real(r8), intent(in) :: se ! entropy term in photosynthesis temp function (J/mol/K) + ! + ! !LOCAL VARIABLES: + real(r8) :: ans + !------------------------------------------------------------------------------- + + ans = 1._r8 + exp( (-hd+se*(tfrz+25._r8)) / (rgas_J_K_mol*(tfrz+25._r8)) ) + + return + end function fth25_f + + ! ===================================================================================== + + subroutine GetCanopyGasParameters(can_press, & + can_o2_partialpress, & + veg_tempk, & + mm_kco2, & + mm_ko2, & + co2_cpoint) + + ! --------------------------------------------------------------------------------- + ! This subroutine calculates the specific Michaelis Menten Parameters (pa) for CO2 + ! and O2, as well as the CO2 compentation point. + ! --------------------------------------------------------------------------------- + + use FatesConstantsMod, only: umol_per_mol + use FatesConstantsMod, only: mmol_per_mol + use FatesConstantsMod, only: umol_per_kmol + + ! Arguments + real(r8), intent(in) :: can_press ! Air pressure within the canopy (Pa) + real(r8), intent(in) :: can_o2_partialpress ! Partial press of o2 in the canopy (Pa) + real(r8), intent(in) :: veg_tempk ! The temperature of the vegetation (K) + + real(r8), intent(out) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) + real(r8), intent(out) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) + real(r8), intent(out) :: co2_cpoint ! CO2 compensation point (Pa) + ! of conductance and resistance: [umol/m3] + + ! Locals + real(r8) :: kc25 ! Michaelis-Menten constant for CO2 at 25C (Pa) + real(r8) :: ko25 ! Michaelis-Menten constant for O2 at 25C (Pa) + real(r8) :: sco ! relative specificity of rubisco + real(r8) :: cp25 ! CO2 compensation point at 25C (Pa) + + ! --------------------------------------------------------------------------------- + ! Intensive values (per mol of air) + ! kc, ko, currentPatch, from: Bernacchi et al (2001) + ! Plant, Cell and Environment 24:253-259 + ! --------------------------------------------------------------------------------- + + real(r8), parameter :: mm_kc25_umol_per_mol = 404.9_r8 + real(r8), parameter :: mm_ko25_mmol_per_mol = 278.4_r8 + real(r8), parameter :: co2_cpoint_umol_per_mol = 42.75_r8 + + ! Activation energy, from: + ! Bernacchi et al (2001) Plant, Cell and Environment 24:253-259 + ! Bernacchi et al (2003) Plant, Cell and Environment 26:1419-1430 + + real(r8), parameter :: kcha = 79430._r8 ! activation energy for kc (J/mol) + real(r8), parameter :: koha = 36380._r8 ! activation energy for ko (J/mol) + real(r8), parameter :: cpha = 37830._r8 ! activation energy for cp (J/mol) + + + ! Derive sco from currentPatch and O2 using present-day O2 (0.209 mol/mol) and re-calculate + ! currentPatch to account for variation in O2 using currentPatch = 0.5 O2 / sco + + ! FIXME (RGK 11-30-2016 THere are more constants here, but I don't have enough information + ! about what they are or do, so I can't give them more descriptive names. Someone please + ! fill this in when possible) + + kc25 = ( mm_kc25_umol_per_mol / umol_per_mol ) * can_press + ko25 = ( mm_ko25_mmol_per_mol / mmol_per_mol ) * can_press + sco = 0.5_r8 * 0.209_r8 / (co2_cpoint_umol_per_mol / umol_per_mol ) + cp25 = 0.5_r8 * can_o2_partialpress / sco + + if( veg_tempk.gt.150_r8 .and. veg_tempk.lt.350_r8 )then + mm_kco2 = kc25 * ft1_f(veg_tempk, kcha) + mm_ko2 = ko25 * ft1_f(veg_tempk, koha) + co2_cpoint = cp25 * ft1_f(veg_tempk, cpha) + else + mm_kco2 = 1.0_r8 + mm_ko2 = 1.0_r8 + co2_cpoint = 1.0_r8 + end if + + return + end subroutine GetCanopyGasParameters + + + ! ==================================================================================== + + subroutine LeafLayerMaintenanceRespiration_Ryan_1991(lnc_top, & + nscaler, & + ft, & + veg_tempk, & + lmr) + + + + + ! ----------------------------------------------------------------------- + ! Base maintenance respiration rate for plant tissues maintresp_leaf_ryan1991_baserate + ! M. Ryan, 1991. Effects of climate change on plant respiration. + ! Ecological Applications, 1(2), 157-167. + ! Original expression is br = 0.0106 molC/(molN h) + ! Conversion by molecular weights of C and N gives 2.525e-6 gC/(gN s) + ! Which is the default value of maintresp_nonleaf_baserate + + ! Arguments + real(r8), intent(in) :: lnc_top ! Leaf nitrogen content per unit area at canopy top [gN/m2] + real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile + integer, intent(in) :: ft ! (plant) Functional Type Index + real(r8), intent(in) :: veg_tempk ! vegetation temperature + real(r8), intent(out) :: lmr ! Leaf Maintenance Respiration (umol CO2/m**2/s) + + ! Locals + real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) + real(r8) :: lmr25top ! canopy top leaf maint resp rate at 25C for this pft (umol CO2/m**2/s) + + ! Parameter + real(r8), parameter :: lmrha = 46390._r8 ! activation energy for lmr (J/mol) + real(r8), parameter :: lmrhd = 150650._r8 ! deactivation energy for lmr (J/mol) + real(r8), parameter :: lmrse = 490._r8 ! entropy term for lmr (J/mol/K) + real(r8), parameter :: lmrc = 1.15912391_r8 ! scaling factor for high + + ! temperature inhibition (25 C = 1.0) + + lmr25top = lb_params%maintresp_leaf_ryan1991_baserate(ft) * (1.5_r8 ** ((25._r8 - 20._r8)/10._r8)) + lmr25top = lmr25top * lnc_top / (umolC_to_kgC * g_per_kg) + + + ! Part I: Leaf Maintenance respiration: umol CO2 / m**2 [leaf] / s + ! ---------------------------------------------------------------------------------- + lmr25 = lmr25top * nscaler + + + if (lb_params%c3psn(ft) == c3_path_index) then + ! temperature sensitivity of C3 plants + lmr = lmr25 * ft1_f(veg_tempk, lmrha) * & + fth_f(veg_tempk, lmrhd, lmrse, lmrc) + else + ! temperature sensitivity of C4 plants + lmr = lmr25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) + lmr = lmr / (1._r8 + exp( 1.3_r8*(veg_tempk-(tfrz+55._r8)) )) + endif + + ! Any hydrodynamic limitations could go here, currently none + ! lmr = lmr * (nothing) + + end subroutine LeafLayerMaintenanceRespiration_Ryan_1991 + + ! ==================================================================================== + + subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, & + rdark_scaler, & + ft, & + veg_tempk, & + tgrowth, & + lmr) + + + + ! Arguments + real(r8), intent(in) :: lnc_top ! Leaf nitrogen content per unit area at canopy top [gN/m2] + real(r8), intent(in) :: rdark_scaler ! Decay coefficient for vertical scaling of respiration. See note 1 below + integer, intent(in) :: ft ! (plant) Functional Type Index + real(r8), intent(in) :: veg_tempk ! vegetation temperature (degrees K) + real(r8), intent(in) :: tgrowth ! lagged vegetation temperature averaged over acclimation timescale (degrees K) + real(r8), intent(out) :: lmr ! Leaf Maintenance Respiration (umol CO2/m**2/s) + + + ! Locals + real(r8) :: r_0 ! base respiration rate, PFT-dependent (umol CO2/m**2/s) + real(r8) :: r_t_ref ! acclimated ref respiration rate (umol CO2/m**2/s) + + ! parameter values of r_0 as listed in Atkin et al 2017: (umol CO2/m**2/s) + ! Broad-leaved trees 1.7560 + ! Needle-leaf trees 1.4995 + ! Shrubs 2.0749 + ! C3 herbs/grasses 2.1956 + ! In the absence of better information, we use the same value for C4 grasses as C3 grasses. + + ! r_0 currently put into the EDPftvarcon_inst%dev_arbitrary_pft + ! all figs in Atkin et al 2017 stop at zero Celsius so we will assume acclimation is fixed below that + r_0 = lb_params%maintresp_leaf_atkin2017_baserate(ft) + + ! Note 1: + ! This code uses the relationship between leaf N and respiration from Atkin et al + ! for the top of the canopy, but then scales through the canopy based on a rdark_scaler. + ! To assume proportionality with N through the canopy following Lloyd et al. 2010, use the + ! default parameter value of 2.43, which results in the scaling of photosynthesis and respiration + ! being proportional through the canopy. To have a steeper decrease in respiration than photosynthesis + ! this number can be smaller. There is some observational evidence for this being the case + ! in Lamour et al. 2023. + + r_t_ref = max(0._r8, rdark_scaler * (r_0 + lmr_r_1 * lnc_top + lmr_r_2 * max(0._r8, (tgrowth - tfrz) )) ) + + if (r_t_ref .eq. 0._r8) then + warn_msg = 'Rdark is negative at this temperature and is capped at 0. tgrowth (C): ' & + //trim(N2S(tgrowth-tfrz))//' pft: '//trim(I2S(ft)) + call FatesWarn(warn_msg,index=4) + end if + + lmr = r_t_ref * exp(lmr_b * (veg_tempk - tfrz - lmr_TrefC) + lmr_c * & + ((veg_tempk-tfrz)**2 - lmr_TrefC**2)) + + end subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017 + + ! ==================================================================================== + + subroutine LeafLayerBiophysicalRates( ft, & + vcmax25top_ft, & + jmax25top_ft, & + kp25_ft, & + nscaler, & + veg_tempk, & + dayl_factor, & + t_growth, & + t_home, & + btran, & + vcmax, & + jmax, & + kp, & + gs0, & + gs1, & + gs2) + + ! --------------------------------------------------------------------------------- + ! This subroutine calculates the localized values of several key photosynthesis + ! rates. By localized, we mean specific to the plant type and leaf layer, + ! which factors in leaf physiology, as well as environmental effects. + ! This procedure should be called prior to iterative solvers, and should + ! have pre-calculated the reference rates for the pfts before this. + ! + ! The output biophysical rates are: + ! vcmax: maximum rate of carboxilation, + ! jmax: maximum electron transport rate, + ! kp: initial slope of CO2 response curve (C4 plants) + ! gs0,gs1,gs2: Stomatal conductance slopes and intercepts (multiplied by btran) + ! --------------------------------------------------------------------------------- + + ! Arguments + ! ------------------------------------------------------------------------------ + + integer, intent(in) :: ft ! (plant) Functional Type Index + real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile + real(r8), intent(in) :: vcmax25top_ft ! canopy top maximum rate of carboxylation at 25C + ! for this pft (umol CO2/m**2/s) + real(r8), intent(in) :: jmax25top_ft ! canopy top maximum electron transport rate at 25C + ! for this pft (umol electrons/m**2/s) + real(r8), intent(in) :: kp25_ft ! initial slope of CO2 response curve + ! (C4 plants) at 25C, canopy top, this pft + real(r8), intent(in) :: veg_tempk ! vegetation temperature + real(r8), intent(in) :: dayl_factor ! daylength scaling factor (0-1) + real(r8), intent(in) :: t_growth ! T_growth (short-term running mean temperature) (K) + real(r8), intent(in) :: t_home ! T_home (long-term running mean temperature) (K) + real(r8), intent(in) :: btran ! transpiration wetness factor (0 to 1) + + real(r8), intent(out) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) + real(r8), intent(out) :: jmax ! maximum electron transport rate + ! (umol electrons/m**2/s) + real(r8), intent(out) :: kp ! initial slope of CO2 response curve (C4 plants) + real(r8), intent(out) :: gs0 ! effective stomatal intercept + real(r8), intent(out) :: gs1 ! effective stomatal slope + real(r8), intent(out) :: gs2 ! alternative btran term applied to Medlyn + ! conductance on whole non-intercept side of equation + + ! Locals + ! ------------------------------------------------------------------------------- + real(r8) :: vcmax25 ! leaf layer: maximum rate of carboxylation at 25C + ! (umol CO2/m**2/s) + real(r8) :: jmax25 ! leaf layer: maximum electron transport rate at 25C + ! (umol electrons/m**2/s) + real(r8) :: dayl_factor_local ! Local version of daylength factor + + ! Parameters + ! --------------------------------------------------------------------------------- + real(r8) :: vcmaxha ! activation energy for vcmax (J/mol) + real(r8) :: jmaxha ! activation energy for jmax (J/mol) + real(r8) :: vcmaxhd ! deactivation energy for vcmax (J/mol) + real(r8) :: jmaxhd ! deactivation energy for jmax (J/mol) + real(r8) :: vcmaxse ! entropy term for vcmax (J/mol/K) + real(r8) :: jmaxse ! entropy term for jmax (J/mol/K) + real(r8) :: t_growth_celsius ! average growing temperature + real(r8) :: t_home_celsius ! average home temperature + real(r8) :: jvr ! ratio of Jmax25 / Vcmax25 + real(r8) :: vcmaxc ! scaling factor for high temperature inhibition (25 C = 1.0) + real(r8) :: jmaxc ! scaling factor for high temperature inhibition (25 C = 1.0) + + + ! update the daylength factor local variable if the switch is on + if ( lb_params%dayl_switch == itrue ) then + dayl_factor_local = dayl_factor + else + dayl_factor_local = 1.0_r8 + endif + + select case(lb_params%photo_tempsens_model) + case (photosynth_acclim_model_none) !No temperature acclimation + vcmaxha = lb_params%vcmaxha(FT) + jmaxha = lb_params%jmaxha(FT) + vcmaxhd = lb_params%vcmaxhd(FT) + jmaxhd = lb_params%jmaxhd(FT) + vcmaxse = lb_params%vcmaxse(FT) + jmaxse = lb_params%jmaxse(FT) + + case (photosynth_acclim_model_kumarathunge_etal_2019) !Kumarathunge et al. temperature acclimation, Thome=30-year running mean + t_growth_celsius = t_growth-tfrz + t_home_celsius = t_home-tfrz + vcmaxha = (42.6_r8 + (1.14_r8*t_growth_celsius))*1e3_r8 !J/mol + jmaxha = 40.71_r8*1e3_r8 !J/mol + vcmaxhd = 200._r8*1e3_r8 !J/mol + jmaxhd = 200._r8*1e3_r8 !J/mol + vcmaxse = (645.13_r8 - (0.38_r8*t_growth_celsius)) + jmaxse = 658.77_r8 - (0.84_r8*t_home_celsius) - 0.52_r8*(t_growth_celsius-t_home_celsius) + jvr = 2.56_r8 - (0.0375_r8*t_home_celsius)-(0.0202_r8*(t_growth_celsius-t_home_celsius)) + + + case default + write (fates_log(),*)'error, incorrect leaf photosynthesis temperature acclimation model specified' + write (fates_log(),*)'lb_params%photo_tempsens_model: ',lb_params%photo_tempsens_model + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + + vcmaxc = fth25_f(vcmaxhd, vcmaxse) + jmaxc = fth25_f(jmaxhd, jmaxse) + + ! Vcmax25top was already calculated to derive the nscaler function + vcmax25 = vcmax25top_ft * nscaler * dayl_factor_local + + select case( lb_params%photo_tempsens_model) + case (photosynth_acclim_model_none) + jmax25 = jmax25top_ft * nscaler * dayl_factor_local + case (photosynth_acclim_model_kumarathunge_etal_2019) + jmax25 = vcmax25*jvr + case default + write (fates_log(),*)'error, incorrect leaf photosynthesis temperature acclimation model specified' + write (fates_log(),*)'lb_params%photo_tempsens_model:',lb_params%photo_tempsens_model + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + + ! Adjust for temperature + ! photosynthetic pathway: 0. = c4, 1. = c3 + + if (lb_params%c3psn(ft) == c3_path_index) then + vcmax = vcmax25 * ft1_f(veg_tempk, vcmaxha) * fth_f(veg_tempk, vcmaxhd, vcmaxse, vcmaxc) + kp = -9999._r8 + else + vcmax = vcmax25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) + vcmax = vcmax / (1._r8 + exp( 0.2_r8*((tfrz+15._r8)-veg_tempk ) )) + vcmax = vcmax / (1._r8 + exp( 0.3_r8*(veg_tempk-(tfrz+40._r8)) )) + kp = kp25_ft * nscaler * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) + end if + + jmax = jmax25 * ft1_f(veg_tempk, jmaxha) * fth_f(veg_tempk, jmaxhd, jmaxse, jmaxc) + + ! Adjust various rates for water limitations + ! ----------------------------------------------------------------------------------- + + ! Apply water limitations to Vcmax + if(lb_params%agross_btran_model(ft) .ne. btran_on_ag_none) then + vcmax = vcmax * btran + end if + + ! Apply water limitations to Jmax + if(lb_params%agross_btran_model(ft) .eq. btran_on_ag_vcmax_jmax) then + jmax = jmax * btran + end if + + ! Make sure that vcmax and jmax do not drop below a lower + ! threshold, see where the constants are defined for an explanation. + ! ----------------------------------------------------------------------------------- + + vcmax = max(min_vcmax_frac*vcmax25top_ft,vcmax) + + jmax = max(min_jmax_frac*jmax25top_ft,jmax) + + + ! Apply water limitations to stomatal intercept (hypothesis dependent) + + if(lb_params%stomatal_btran_model(ft)==btran_on_gs_gs0 .or. & + lb_params%stomatal_btran_model(ft)==btran_on_gs_gs01 .or. & + lb_params%stomatal_btran_model(ft)==btran_on_gs_gs02 )then + gs0 = max(gs0_min,lb_params%stomatal_intercept(ft)*btran) + else + gs0 = max(gs0_min,lb_params%stomatal_intercept(ft)) + end if + + ! Apply water limitations to stomatal slope (hypothesis dependent) + gs2 = 1._r8 + if(lb_params%stomatal_btran_model(ft)==btran_on_gs_gs1 .or. & + lb_params%stomatal_btran_model(ft)==btran_on_gs_gs01)then + if(lb_params%stomatal_model.eq.medlyn_model ) then + gs1 = lb_params%medlyn_slope(ft)*btran + else + gs1 = lb_params%bb_slope(ft)*btran + end if + elseif(lb_params%stomatal_btran_model(ft)==btran_on_gs_gs2 .or. & + lb_params%stomatal_btran_model(ft)==btran_on_gs_gs02)then + if(lb_params%stomatal_model.eq.medlyn_model ) then + gs2 = btran + gs1 = lb_params%medlyn_slope(ft) + else + gs1 = lb_params%bb_slope(ft) + end if + else + if(lb_params%stomatal_model.eq.medlyn_model ) then + gs1 = lb_params%medlyn_slope(ft) + else + gs1 = lb_params%bb_slope(ft) + end if + end if + + + + + return + end subroutine LeafLayerBiophysicalRates + + ! ===================================================================================== + + function DecayCoeffVcmax(vcmax25top,slope_param,intercept_param) result(decay_coeff_vcmax) + + ! --------------------------------------------------------------------------------- + ! This function estimates the decay coefficient used to estimate vertical + ! attenuation of properties in the canopy. + ! + ! Decay coefficient (kn) is a function of vcmax25top for each pft. + ! + ! Currently, this decay is applied to vcmax attenuation, SLA (optionally) + ! and leaf respiration (optionally w/ Atkin) + ! + ! --------------------------------------------------------------------------------- + + !ARGUMENTS + + real(r8),intent(in) :: vcmax25top + real(r8),intent(in) :: slope_param ! multiplies vcmax25top + real(r8),intent(in) :: intercept_param ! adds to vcmax25top + + real(r8) :: decay_coeff_vcmax + + !LOCAL VARIABLES + ! ----------------------------------------------------------------------------------- + + ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used + ! kn = 0.11. Here, we derive kn from vcmax25 as in Lloyd et al + ! (2010) Biogeosciences, 7, 1833-1859 + ! This function is also used to vertically scale leaf maintenance + ! respiration. + + decay_coeff_vcmax = exp(slope_param * vcmax25top - intercept_param) + + return + end function DecayCoeffVcmax + + ! ===================================================================================== + + subroutine LowstorageMainRespReduction(frac, ft, maintresp_reduction_factor) + + ! This subroutine reduces maintenance respiration rates when storage pool is low. The premise + ! of this is that mortality of plants increases when storage is low because they are not able + ! to repair tissues, generate defense compounds, etc. This reduction is reflected in a reduced + ! maintenance demand. The output of this function takes the form of a curve between 0 and 1, + ! and the curvature of the function is determined by a parameter. + + ! Arguments + ! ------------------------------------------------------------------------------ + real(r8), intent(in) :: frac ! ratio of storage to target leaf biomass + integer, intent(in) :: ft ! what pft is this cohort? + real(r8), intent(out) :: maintresp_reduction_factor ! the factor by which to reduce maintenance respiration + + ! -------------------------------------------------------------------------------- + ! Parameters are at the PFT level: + ! fates_maintresp_reduction_curvature controls the curvature of this. + ! If this parameter is zero, then there is no reduction until the plant dies at storage = 0. + ! If this parameter is one, then there is a linear reduction in respiration below the storage point. + ! Intermediate values will give some (concave-downwards) curvature. + ! + ! maintresp_reduction_intercept controls the maximum amount of throttling. + ! zero means no throttling at any point, so it turns this mechanism off completely and so + ! allows an entire cohort to die via negative carbon-induced termination mortality. + ! one means complete throttling, so no maintenance respiration at all, when out of carbon. + ! --------------------------------------------------------------------------------- + + if( frac .lt. 1._r8 )then + if ( abs(lb_params%maintresp_reduction_curvature(ft)-1._r8) > nearzero ) then + maintresp_reduction_factor = (1._r8 - lb_params%maintresp_reduction_intercept(ft)) + & + lb_params%maintresp_reduction_intercept(ft) * & + (1._r8 - lb_params%maintresp_reduction_curvature(ft)**frac) & + / (1._r8-lb_params%maintresp_reduction_curvature(ft)) + else ! avoid nan answer for linear case + maintresp_reduction_factor = (1._r8 - lb_params%maintresp_reduction_intercept(ft)) + & + lb_params%maintresp_reduction_intercept(ft) * frac + endif + + else + maintresp_reduction_factor = 1._r8 + endif + + + end subroutine LowstorageMainRespReduction + + ! ===================================================================================== + + real(r8) function GetConstrainedVPress(air_vpress,veg_esat) result(ceair) + + ! ----------------------------------------------------------------------------------- + ! Return a constrained vapor pressure [Pa] + ! + ! Vapor pressure may not be greater than saturation vapor pressure + ! and may not be less than 1% of saturation vapor pressure + ! + ! This function is used to make sure that Ball-Berry conductance is well + ! behaved. + ! ----------------------------------------------------------------------------------- + + real(r8) :: air_vpress ! vapor pressure of the air (unconstrained) [Pa] + real(r8) :: veg_esat ! saturated vapor pressure [Pa] + real(r8) :: min_frac_esat = 0.01_r8 ! We don't allow vapor pressures + ! below this fraction amount of saturation vapor pressure + + ceair = min( max(air_vpress, min_frac_esat*veg_esat ),veg_esat ) + + end function GetConstrainedVPress + + ! ===================================================================================== + + subroutine QSat (T, p, qs, es, qsdT, esdT) + ! + ! !DESCRIPTION: + ! + ! THIS IS AN EXACT CLONE OF QSat routine in QSatMod.F90 from CTSM, tag ctsm5.2.028 + ! https://github.com/ESCOMP/CTSM/blob/ctsm5.2.028/src/biogeophys/QSatMod.F90 + ! + ! Computes saturation mixing ratio and (optionally) the change in saturation mixing + ! ratio with respect to temperature. Mixing ratio and specific humidity are + ! approximately equal and can be treated as the same. + ! Reference: Polynomial approximations from: + ! Piotr J. Flatau, et al.,1992: Polynomial fits to saturation + ! vapor pressure. Journal of Applied Meteorology, 31, 1507-1513. + + ! !ARGUMENTS: + implicit none + real(r8), intent(in) :: T ! temperature (K) + real(r8), intent(in) :: p ! surface atmospheric pressure (pa) + real(r8), intent(out) :: qs ! humidity (kg/kg) + real(r8), intent(out), optional :: es ! vapor pressure (pa) + real(r8), intent(out), optional :: qsdT ! d(qs)/d(T) + real(r8), intent(out), optional :: esdT ! d(es)/d(T) + ! + ! !LOCAL VARIABLES: + real(r8) :: es_local ! local version of es (in case es is not present) + real(r8) :: esdT_local ! local version of esdT (in case esdT is not present) + real(r8) :: td,vp,vp1,vp2 + + ! For water vapor (temperature range 0C-100C) + real(r8), parameter :: a0 = 6.11213476_r8 + real(r8), parameter :: a1 = 0.444007856_r8 + real(r8), parameter :: a2 = 0.143064234e-01_r8 + real(r8), parameter :: a3 = 0.264461437e-03_r8 + real(r8), parameter :: a4 = 0.305903558e-05_r8 + real(r8), parameter :: a5 = 0.196237241e-07_r8 + real(r8), parameter :: a6 = 0.892344772e-10_r8 + real(r8), parameter :: a7 = -0.373208410e-12_r8 + real(r8), parameter :: a8 = 0.209339997e-15_r8 + ! For derivative:water vapor + real(r8), parameter :: b0 = 0.444017302_r8 + real(r8), parameter :: b1 = 0.286064092e-01_r8 + real(r8), parameter :: b2 = 0.794683137e-03_r8 + real(r8), parameter :: b3 = 0.121211669e-04_r8 + real(r8), parameter :: b4 = 0.103354611e-06_r8 + real(r8), parameter :: b5 = 0.404125005e-09_r8 + real(r8), parameter :: b6 = -0.788037859e-12_r8 + real(r8), parameter :: b7 = -0.114596802e-13_r8 + real(r8), parameter :: b8 = 0.381294516e-16_r8 + ! For ice (temperature range -75C-0C) + real(r8), parameter :: c0 = 6.11123516_r8 + real(r8), parameter :: c1 = 0.503109514_r8 + real(r8), parameter :: c2 = 0.188369801e-01_r8 + real(r8), parameter :: c3 = 0.420547422e-03_r8 + real(r8), parameter :: c4 = 0.614396778e-05_r8 + real(r8), parameter :: c5 = 0.602780717e-07_r8 + real(r8), parameter :: c6 = 0.387940929e-09_r8 + real(r8), parameter :: c7 = 0.149436277e-11_r8 + real(r8), parameter :: c8 = 0.262655803e-14_r8 + ! For derivative:ice + real(r8), parameter :: d0 = 0.503277922_r8 + real(r8), parameter :: d1 = 0.377289173e-01_r8 + real(r8), parameter :: d2 = 0.126801703e-02_r8 + real(r8), parameter :: d3 = 0.249468427e-04_r8 + real(r8), parameter :: d4 = 0.313703411e-06_r8 + real(r8), parameter :: d5 = 0.257180651e-08_r8 + real(r8), parameter :: d6 = 0.133268878e-10_r8 + real(r8), parameter :: d7 = 0.394116744e-13_r8 + real(r8), parameter :: d8 = 0.498070196e-16_r8 + + + !----------------------------------------------------------------------- + + td = min(100.0_r8, max(-75.0_r8, T - tfrz )) + + if (td >= 0.0_r8) then + es_local = a0 + td*(a1 + td*(a2 + td*(a3 + td*(a4 & + + td*(a5 + td*(a6 + td*(a7 + td*a8))))))) + else + es_local = c0 + td*(c1 + td*(c2 + td*(c3 + td*(c4 & + + td*(c5 + td*(c6 + td*(c7 + td*c8))))))) + endif + + es_local = es_local * 100._r8 ! pa + vp = 1.0_r8 / (p - 0.378_r8*es_local) + vp1 = 0.622_r8 * vp + qs = es_local * vp1 ! kg/kg + if (present(es)) then + es = es_local + end if + + if (present(qsdT) .or. present(esdT)) then + if (td >= 0.0_r8) then + esdT_local = b0 + td*(b1 + td*(b2 + td*(b3 + td*(b4 & + + td*(b5 + td*(b6 + td*(b7 + td*b8))))))) + else + esdT_local = d0 + td*(d1 + td*(d2 + td*(d3 + td*(d4 & + + td*(d5 + td*(d6 + td*(d7 + td*d8))))))) + end if + + esdT_local = esdT_local * 100._r8 ! pa/K + vp2 = vp1 * vp + if (present(qsdT)) then + qsdT = esdT_local * vp2 * p ! 1 / K + end if + if (present(esdT)) then + esdT = esdT_local + end if + end if + + end subroutine QSat + + ! ===================================================================================== + + real(r8) function VeloToMolarCF(press,tempk) result(cf) + + ! --------------------------------------------------------------------------------- + ! + ! "cf" is the conversion factor between molar form and velocity form + ! of conductance and resistance. The units on this factor are: [umol/m3] + ! This uses the ideal gas law. This routine is necessary because the + ! photosynthesis module uses units of micromoles, while the land-energy + ! balance solvers in the host models use units of velocity. + ! + ! i.e. + ! [m/s] * [umol/m3] -> [umol/m2/s] + ! + ! Breakdown of the conversion factor: [ umol / m3 ] + ! + ! Rgas [J /K /kmol] + ! Air Potential Temperature [ K ] + ! Air Pressure [ Pa ] + ! conversion: umol/kmol = 1e9 + ! + ! [ Pa * K * kmol umol/kmol / J K ] = [ Pa * umol / J ] + ! since: 1 Pa = 1 N / m2 + ! [ Pa * umol / J ] = [ N * umol / J m2 ] + ! since: 1 J = 1 N * m + ! [ N * umol / J m2 ] = [ N * umol / N m3 ] + ! [ umol / m3 ] + ! + ! -------------------------------------------------------------------------------- + + ! Arguments + real(r8) :: press ! air pressure at point of interest [Pa] + real(r8) :: tempk ! temperature at point of interest [K] + + cf = press/(rgas_J_K_kmol * tempk )*umol_per_kmol + + end function VeloToMolarCF + + ! ===================================================================================== + +end module LeafBiophysicsMod diff --git a/functional_unit_testing/leaf_biophys/CtypesLeafBiophys.py b/functional_unit_testing/leaf_biophys/CtypesLeafBiophys.py new file mode 100644 index 0000000000..40b16457e5 --- /dev/null +++ b/functional_unit_testing/leaf_biophys/CtypesLeafBiophys.py @@ -0,0 +1,89 @@ +import re +import subprocess +import sys +import ctypes +from ctypes import * +from operator import add +sys.path.append('../shared/py_src') +from PyF90Utils import c8, ci, cchar, c8_arr, ci_arr, ccharnb + +# Subroutines +# ======================================================================================= + +def GetModSymbol(mod_path,symbol): + + # This routine uses some minor string magic and a system call using + # "nm" to list the symbols inside a fortran object. Symbols are simply + # the text strings the define the routines and data structures inside + # the compiled objects. When symbols are created inside a compiled object, + # different compilers do different things, like force to lower or upper + # case, or attach trailing or preceding underscores. This routine + # helps identify the exact symbol name, from the name of the fortran + # routine (as it appears in code) that we want to match with. + + subprocess.run(["nm", mod_path],capture_output=True) + + list = str(subprocess.run(["nm",mod_path],capture_output=True)).split() + p = re.compile(symbol, re.IGNORECASE) + modlist = [ s for s in list if p.search(s)] + + # If nothing came up, thats bad + if (len(modlist)==0): + print('Failed to find the right module symbol for:{} in module {}'.format(symbol,mod_path)) + print(list) + print('Exiting') + exit(2) + + # Its possible a routine name could also be part of another routine name, + # such as alloc_such_and_such vs dealloc_such_and_such, we want the shorter + mod_symbols = [] + slen = 10000 + for item in modlist: + ms_str = item.split('\\')[0] + if len(ms_str)10cm/s? + # Convert to molar form using 1 standard atm at 25C + # Units: umol/m2/s + # 1 mol/m2/s was roughly the middle of the distribution + # when generating conductances at BCI using a static canopy + # and local driver data. + + bl_cond_const = 1.e6 + + # Lets look at canopy top + # kn = DecayCoeffVcmax(currentCohort%vcmax25top, & + # EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff1(ft), & + # EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff2(ft)) + # rdark_scaler = exp(-kn * cumulative_lai) + + rdark_scaler = 1.0 + + # kn = DecayCoeffVcmax(currentCohort%vcmax25top, & + # prt_params%leafn_vert_scaler_coeff1(ft), & + # prt_params%leafn_vert_scaler_coeff2(ft)) + + nscaler = 1.0 + + # Set convergence tolerance + #ci_tol = 2.*can_press_1atm + ci_tol = 0.1 + + # Initialize fortran return values + # these are all temps and doubles + vcmax_f = c_double(-9.0) + jmax_f = c_double(-9.0) + kp_f = c_double(-9.0) + agross_f = c_double(-9.0) + gstoma_f = c_double(-9.0) + anet_f = c_double(-9.0) + lmr_f = c_double(-9.0) + c13_f = c_double(-9.0) + ac_f = c_double(-9.0) + aj_f = c_double(-9.0) + ap_f = c_double(-9.0) + co2_interc_f = c_double(-9.0) + veg_qs_f = c_double(-9.0) + veg_es_f = c_double(-9.0) + mm_kco2_f = c_double(-9.0) + mm_ko2_f = c_double(-9.0) + co2_cpoint_f = c_double(-9.0) + qsdt_dummy_f = c_double(-9.0) + esdt_dummy_f = c_double(-9.0) + solve_iter_f = c_int(-9) + gs0_f = c_double(-9.0) + gs1_f = c_double(-9.0) + gs2_f = c_double(-9.0) + + print('Prepping Canopy Gas Parameters') + + for it, leaf_tempc in enumerate(leaf_tempc_vec): + + leaf_tempk = leaf_tempc + tfrz_1atm + + iret = f90_qsat_sub(c8(leaf_tempk),c8(can_press_1atm), \ + byref(veg_qs_f),byref(veg_es_f), \ + byref(qsdt_dummy_f),byref(esdt_dummy_f)) + + + veg_qsat_vec[it] = veg_qs_f.value + veg_esat_vec[it] = veg_es_f.value + + iret = f90_cangas_sub(c8(can_press_1atm), \ + c8(o2_ppress_209kppm), \ + c8(leaf_tempk), \ + byref(mm_kco2_f), \ + byref(mm_ko2_f), \ + byref(co2_cpoint_f)) + + mm_kco2_vec[it] = mm_kco2_f.value + mm_ko2_vec[it] = mm_ko2_f.value + co2_cpoint_vec[it] = co2_cpoint_f.value + + + + + print('\n') + print('Experiment 1: Evaluating Photosynthesis Equations by pft/Tl/RH/PR') + + elapsed_time = 0.0 + for pft in range(numpft): + + if(do_evalvjkbytemp): + print('\n') + print('Experiment 1: Evaluating Vcmax,Jmax,Kp by Temperature') + EvalVJKByTemp(pft,fates_leaf_vcmax25top[pft],leaf_c3psn[pft],pdf) + + + print('Evaluating PFT {}'.format(pft+1)) + + jmax25_top,kp25_top = GetJmaxKp25Top(fates_leaf_vcmax25top[pft]) + + vcmax = np.zeros([leaf_tempc_n]) + jmax = np.zeros([leaf_tempc_n]) + kp = np.zeros([leaf_tempc_n]) + lmr = np.zeros([leaf_tempc_n]) + agross = np.zeros([leaf_tempc_n,rh_n,par_abs_n,gb_n,btran_n]) + gstoma = np.zeros([leaf_tempc_n,rh_n,par_abs_n,gb_n,btran_n]) + co2_interc = np.zeros([leaf_tempc_n,rh_n,par_abs_n,gb_n,btran_n]) + anet = np.zeros([leaf_tempc_n,rh_n,par_abs_n,gb_n]) + ac = np.zeros([leaf_tempc_n,rh_n,par_abs_n,gb_n]) + aj = np.zeros([leaf_tempc_n,rh_n,par_abs_n,gb_n]) + ap = np.zeros([leaf_tempc_n,rh_n,par_abs_n,gb_n]) + iters = np.zeros([leaf_tempc_n,rh_n,par_abs_n,gb_n]) + + # When calling component limitations exclusively + # using an approximation of interstitial co2 as + # 0.7*canopy_co2 + + ac2 = np.zeros([leaf_tempc_n]) + aj2 = np.zeros([leaf_tempc_n,par_abs_n]) + ap2 = np.zeros([leaf_tempc_n]) + + # Leaf Nitrogen Concentration at the top + lnc_top = fates_stoich_nitr[pft]/fates_leaf_slatop[pft] + + for it, leaf_tempc in enumerate(leaf_tempc_vec): + + leaf_tempk = leaf_tempc + tfrz_1atm + + iret = f90_biophysrate_sub(ci(pft+1), c8(fates_leaf_vcmax25top[pft]), \ + c8(jmax25_top), c8(kp25_top), \ + c8(nscaler), c8(leaf_tempk), c8(dayl_factor_full), \ + c8(t_growth_kum),c8(t_home_kum),c8(btran_nolimit), \ + byref(vcmax_f), byref(jmax_f), byref(kp_f), byref(gs0_f), byref(gs1_f), byref(gs2_f)) + + vcmax[it] = vcmax_f.value + jmax[it] = jmax_f.value + kp[it] = kp_f.value + + if(leaf_c3psn[pft] == 1): + ap2[it] = 0.0 + else: + ap2[it] = f90_agross_pepc4(c8(co2_ppress_400ppm),c8(kp[it]),c8(can_press_1atm)) + + # Leaf Maintenance Respiration (temp and pft dependent) + if(fates_maintresp_leaf_model==1): + iret = f90_lmr_ryan_sub(c8(lnc_top),c8(nscaler), ci(pft+1), c8(leaf_tempk), byref(lmr_f)) + elif(fates_maintresp_leaf_model==2): + iret = f90_lmr_atkin_sub(c8(lnc_top),c8(rdark_scaler),c8(leaf_tempk),c8(atkin_mean_leaf_tempk),byref(lmr_f) ) + else: + print('unknown leaf respiration model') + exit(1) + + lmr[it] = lmr_f.value + + if(leaf_c3psn[pft] == 1): + + ac2[it] = f90_agross_rubiscoc3(c8(vcmax[it]),c8(ci_ppress_static),c8(o2_ppress_209kppm), \ + c8(co2_cpoint_vec[it]),c8(mm_kco2_vec[it]),c8(mm_ko2_vec[it])) + else: + ac2[it] = vcmax[it] + + + for ip, par_abs in enumerate(par_abs_vec): + + if(leaf_c3psn[pft] == 1): + aj2[it,ip] = f90_agross_rubpc3(c8(par_abs),c8(jmax[it]),c8(ci_ppress_static),c8(co2_cpoint_vec[it])) + else: + aj2[it,ip] = f90_agross_rubpc4(c8(par_abs)) + + for ir, rh in enumerate(rh_vec): + + vpress = rh * veg_esat_vec[it] + + for ig, gb in enumerate(gb_vec): + + for ib, btran in enumerate(btran_vec): + + iret = f90_biophysrate_sub(ci(pft+1), c8(fates_leaf_vcmax25top[pft]), \ + c8(jmax25_top), c8(kp25_top), \ + c8(nscaler), c8(leaf_tempk), c8(dayl_factor_full), \ + c8(t_growth_kum),c8(t_home_kum),c8(btran), \ + byref(vcmax_f), byref(jmax_f), byref(kp_f), byref(gs0_f), byref(gs1_f), byref(gs2_f)) + + iret = f90_leaflayerphoto_sub(c8(par_abs), \ + c8(par_abs), \ + c8(1.0), \ + ci(pft+1), \ + c8(vcmax_f.value), \ + c8(jmax_f.value), \ + c8(kp_f.value), \ + c8(gs0_f.value), \ + c8(gs1_f.value), \ + c8(gs2_f.value), \ + c8(leaf_tempk), \ + c8(can_press_1atm), \ + c8(co2_ppress_400ppm), \ + c8(o2_ppress_209kppm), \ + c8(gb), \ + c8(vpress), \ + c8(mm_kco2_vec[it]), \ + c8(mm_ko2_vec[it]), \ + c8(co2_cpoint_vec[it]), \ + c8(lmr[it]), \ + c8(ci_tol), \ + byref(agross_f), \ + byref(gstoma_f), \ + byref(anet_f), \ + byref(c13_f), \ + byref(ac_f), \ + byref(aj_f), \ + byref(ap_f), \ + byref(co2_interc_f), \ + byref(solve_iter_f) ) + + # Call the medlyn solve to get timing info + iret = f90_qsat_sub(c8(leaf_tempk),c8(can_press_1atm), \ + byref(veg_qs_f),byref(veg_es_f), \ + byref(qsdt_dummy_f),byref(esdt_dummy_f)) + + + + agross[it,ir,ip,ig,ib] = agross_f.value + gstoma[it,ir,ip,ig,ib] = gstoma_f.value + anet[it,ir,ip,ig] = anet_f.value + ac[it,ir,ip,ig] = ac_f.value + aj[it,ir,ip,ig] = aj_f.value + ap[it,ir,ip,ig] = ap_f.value + co2_interc[it,ir,ip,ig] = co2_interc_f.value + iters[it,ir,ip,ig] = solve_iter_f.value + time0 = time.process_time() + + iret = f90_gs_medlyn(c8(anet_f.value),c8(veg_es_f.value),c8(vpress),c8(gs0_f.value),c8(gs1_f.value), \ + c8(gs2_f.value),c8(co2_ppress_400ppm),c8(can_press_1atm),c8(gb),byref(gstoma_f)) + elapsed_time = elapsed_time + (time.process_time() - time0) + + + print("\nElapsed Time for Conductance: {}\n".format(elapsed_time)) + + # Plot out component gross assimilation rates + # by temperature with constant Ci, and by Ci with + # constant temperature + + fig2,ax1 = plt.subplots(1,1,figsize=(6.5,5.5)) + + ax1.plot(leaf_tempc_vec,ac2,label='Ac') + ix2_10 = int(par_abs_n*0.1) + ix2_50 = int(par_abs_n*0.5) + ix2_90 = int(par_abs_n*0.9) + ax1.plot(leaf_tempc_vec,aj2[:,ix2_10],color=[0.5,0.5,0.5],linestyle='dotted',label='Aj apar=%4.1f'%(par_abs_vec[ix2_10])) + ax1.plot(leaf_tempc_vec,aj2[:,ix2_50],color=[0.5,0.5,0.5],linestyle='dashed',label='Aj apar=%4.1f'%(par_abs_vec[ix2_50])) + ax1.plot(leaf_tempc_vec,aj2[:,ix2_90],color=[0.5,0.5,0.5],linestyle='solid',label='Aj apar=%4.1f'%(par_abs_vec[ix2_90])) + ax1.plot(leaf_tempc_vec,ap2[:],color='orange',linestyle='solid',label='Ap') + ax1.plot(leaf_tempc_vec, lmr, color=[0.7,0.5,0.3],linestyle='solid',label='Rdark') + ax1.set_ylabel('[umol/m2/s]') + ax1.set_xlabel('Leaf Temperature [C]') + ax1.set_title('PFT: %3i, Vcmax25: %4.1f, Jmax25: %4.1f, Ci: %4.1f'%(pft+1,fates_leaf_vcmax25top[pft],jmax25_top,ci_ppress_static)) + ax1.grid(True) + fig2.legend(loc='upper left') + if(pdf): + pdf.savefig(fig2) + plt.close(fig2) + else: + plt.savefig('images/rates_pft%2.2i.png'%(pft)) + + + # Lets plot metrics by temperature, using the + # 10th, and 90th percentiles of both RH and PAR + # Agross, Anet, Gstoma, Ac, Aj, Ap + + fig5, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2,figsize=(8.,7.)) + gb_id = 4 + bt_id = -1 + LinePlotY3dM1(ax1,leaf_tempc_vec,rh_vec,par_abs_vec, \ + agross[:,:,:,gb_id,bt_id].reshape([leaf_tempc_n,rh_n,par_abs_n]),'RH','APAR',True) + ax1.set_ylabel('Agross \n [umol/m2/s]') + ax1.set_xticklabels([]) + LinePlotY3dM1(ax2,leaf_tempc_vec,rh_vec,par_abs_vec,\ + gstoma[:,:,:,gb_id,bt_id].reshape([leaf_tempc_n,rh_n,par_abs_n])*1.e-6,'RH','APAR',False) + ax2.set_ylabel('Gs \n [mol/m2/s]') + ax2.yaxis.set_label_position("right") + ax2.yaxis.tick_right() + + LinePlotY3dM1(ax3,leaf_tempc_vec,rh_vec,par_abs_vec,\ + co2_interc[:,:,:,gb_id,bt_id].reshape([leaf_tempc_n,rh_n,par_abs_n]),'RH','APAR',False) + ax3.set_ylabel('Ci [Pa]') + ax3.axhline(y=co2_ppress_400ppm,color=[0.3,0.3,0.3]) + ax3.plot(leaf_tempc_vec,co2_cpoint_vec) + + ax3.set_xlabel('Leaf Temperature [C]') + ax1.set_xlabel('Leaf Temperature [C]') + + + ax4.text(0.25,0.6,'PFT: %2i \ng_b = %4.2f [mol/m2/s]'%(pft+1,gb_vec[gb_id]*1.e-6)) + ax4.axis('off') + + plt.tight_layout() + plt.subplots_adjust(wspace=0.02, hspace=0.03) + #fig3.legend(loc='lower right',labelspacing = 0.2, fontsize=12) + fig5.legend(loc='lower right', bbox_to_anchor=(0.87, 0.06),labelspacing = 0.2, fontsize=12, ncol=1) #, fancybox=True, shadow=True) + if(pdf): + pdf.savefig(fig5) + plt.close(fig5) + else: + plt.savefig('images/AgGsCI_temp_v2_pft%2.2i.png'%(pft)) + + + # Metrics across the boundary layer conductivity gradient + + # Metrics across the humidity Gradient (fixed conductance) + fig4, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2,figsize=(8.,7.)) + gb_id = 3 + bt_id = -1 + LinePlotY3dM1(ax1,rh_vec,leaf_tempc_vec,par_abs_vec, \ + np.transpose(agross,(1,0,2,3,4 ))[:,:,:,gb_id,bt_id].reshape([rh_n,leaf_tempc_n,par_abs_n]),'Tveg','APAR',True) + + ax1.set_ylabel('Agross \n [umol/m2/s]') + ax1.set_xticklabels([]) + LinePlotY3dM1(ax2,rh_vec,leaf_tempc_vec,par_abs_vec,\ + np.transpose(gstoma,(1,0,2,3,4 ))[:,:,:,gb_id,bt_id].reshape([rh_n,leaf_tempc_n,par_abs_n])*1.e-6,'Tveg','APAR',False) + ax2.set_ylabel('Gs \n [mol/m2/s]') + ax2.yaxis.set_label_position("right") + ax2.yaxis.tick_right() + + LinePlotY3dM1(ax3,rh_vec,leaf_tempc_vec,par_abs_vec,\ + np.transpose(co2_interc,(1,0,2,3,4 ))[:,:,:,gb_id,bt_id].reshape([rh_n,leaf_tempc_n,par_abs_n]),'Tveg','APAR',False) + ax3.set_ylabel('Ci [Pa]') + ax3.axhline(y=co2_ppress_400ppm,color=[0.3,0.3,0.3]) + + ax3.set_xlabel('Relative Humidity [%]') + ax1.set_xlabel('Relative Humidity [%]') + + + ax4.text(0.25,0.6,'PFT: %2i \ng_b = %4.2f [mol/m2/s]'%(pft+1,gb_vec[gb_id]*1.e-6)) + ax4.axis('off') + + plt.tight_layout() + plt.subplots_adjust(wspace=0.02, hspace=0.03) + fig4.legend(loc='lower right', bbox_to_anchor=(0.87, 0.06),labelspacing = 0.2, fontsize=12, ncol=1) #, fancybox=True, shadow=True) + if(pdf): + pdf.savefig(fig4) + plt.close(fig4) + else: + plt.savefig('images/AgGsCI_temp_v3_pft%2.2i.png'%(pft)) + + + + # Look at the btran gradient + + fig44, ((ax1,ax2)) = plt.subplots(1,2,figsize=(8.,5.)) + gb_id = 3 + rh_id = np.argmin(np.abs(rh_vec-1.0)) # ID at 100% humidity + par200_id = np.argmin(np.abs(par_abs_vec-200.)) + tv30_id = np.argmin(np.abs(leaf_tempc_vec-30.)) + par50_id = np.argmin(np.abs(par_abs_vec-50.)) + tv5_id = np.argmin(np.abs(leaf_tempc_vec-5.)) + + ax1.plot(btran_vec,agross[tv30_id,rh_id,par200_id,gb_id,:],label='PAR=%5.1f W Tv=%3.1f C'%(par_abs_vec[par200_id],leaf_tempc_vec[tv30_id]),color = 'r',linestyle = '-') + ax1.plot(btran_vec,agross[tv30_id,rh_id,par50_id,gb_id,:],label='PAR=%5.1f W Tv=%3.1f C'%(par_abs_vec[par50_id],leaf_tempc_vec[tv30_id]),color = 'r',linestyle = '--') + ax1.plot(btran_vec,agross[tv5_id,rh_id,par200_id,gb_id,:],label='PAR=%5.1f W Tv=%3.1f C'%(par_abs_vec[par200_id],leaf_tempc_vec[tv5_id]),color = 'b',linestyle = '-') + ax1.plot(btran_vec,agross[tv5_id,rh_id,par50_id,gb_id,:],label='PAR=%5.1f W Tv=%3.1f C'%(par_abs_vec[par50_id],leaf_tempc_vec[tv5_id]),color = 'b',linestyle = '--') + ax1.set_title('Ag [umol/m2/s]') + ax1.set_xlabel('BTRAN') + ax1.grid('on') + + ax2.plot(btran_vec,gstoma[tv30_id,rh_id,par200_id,gb_id,:]*1.e-6,color = 'r',linestyle = '-') + ax2.plot(btran_vec,gstoma[tv30_id,rh_id,par50_id,gb_id,:]*1.e-6,color = 'r',linestyle = '--') + ax2.plot(btran_vec,gstoma[tv5_id,rh_id,par200_id,gb_id,:]*1.e-6,color = 'b',linestyle = '-') + ax2.plot(btran_vec,gstoma[tv5_id,rh_id,par50_id,gb_id,:]*1.e-6,color = 'b',linestyle = '--') + ax2.set_title('gs [mol/m2/s]') + ax2.set_xlabel('BTRAN') + ax2.grid('on') + fig44.legend(loc='lower right', bbox_to_anchor=(0.87, 0.56),labelspacing = 0.2, fontsize=12, ncol=1) + + + + + + # Show for each pft + if(not pdf): + plt.show() + + + if(pdf): + pdf.close() + + print('Deallocating parameter space') + iret = f90_dealloc_leaf_param_sub() + + print('Functional Unit Testing Complete') + exit(0) + + +# ======================================================================================= +# This is the actual call to main + +if __name__ == "__main__": + main(sys.argv) diff --git a/functional_unit_testing/leaf_biophys/DriveSmokeTestLeafBiophys.py b/functional_unit_testing/leaf_biophys/DriveSmokeTestLeafBiophys.py new file mode 100644 index 0000000000..c1042af0b1 --- /dev/null +++ b/functional_unit_testing/leaf_biophys/DriveSmokeTestLeafBiophys.py @@ -0,0 +1,440 @@ +# ======================================================================================= +# +# For usage: $python LeafBiophysDriver.py --help +# +# This script runs unit tests on the leaf biophysics functions +# +# +# ======================================================================================= + +import matplotlib as mpl +#mpl.use('Agg') +import matplotlib.pyplot as plt +from datetime import datetime +import argparse +from matplotlib.backends.backend_pdf import PdfPages +import platform +import xml.etree.ElementTree as et +import numpy as np +import matplotlib +import os +import sys +import getopt +import code # For development: code.interact(local=dict(globals(), **locals())) +import time +import importlib +import csv +import subprocess +import re +import CtypesLeafBiophys +import ctypes +from ctypes import * +from operator import add +sys.path.append('../shared/py_src') +from PyF90Utils import c8, ci, cchar, c8_arr, ci_arr, ccharnb + +font = {'family' : 'sans-serif', + 'weight' : 'normal', + 'size' : 11} + +matplotlib.rc('font', **font) + + +# Global constants to use in all Leaf Biophysics unit testing +# ======================================================================================= + + +# For debugging +dump_parameters = False + +# Should we evaluate vcmax, jmax and kp actual? +do_evalvjkbytemp = False + +# Do an analysis on convergence +do_test_citol = False + +# Freezing point of water in Kelvin (at standard atmosphere) +tfrz_1atm = 273.15 + +# 25 degrees C in Kelvin (used because T25 functions) +leaf_tempk25 = tfrz_1atm + 25.0 + +# Daylight limitations can be imposed on Vcmax, a value of +# 1 means daylight length is at its maximum +dayl_factor_full = 1.0 + +# If Kumerathunga respiration is used, it requires moving averages +# of leaf temperature +t_growth_kum = -999 +t_home_kum = -999 + +# Simple conversion, number of micro-moles in a mole +umol_per_mol = 1.e6 +mol_per_umol = 1.e-6 + +# 1 standard atmosphere in [Pa] +can_press_1atm = 101325.0 + +# Atmospheric CO2 partial pressure [Pa] at 400 ppm +co2_ppress_400ppm = 0.0004*can_press_1atm + +# Atmospheric O2 partial pressure [Pa] %29.5 of atmosphere +o2_ppress_209kppm = 0.2095*can_press_1atm + +# 70% of atmospheric CO2 is a reasonablish guess for +# intercellular CO2 concentration during primary production +# We can use this to test the gross assimilation routines +# directly, without having to solve for the equilibrium +# intercellular CO2 +ci_ppress_static = 0.7*co2_ppress_400ppm + +# When there is hydrualic limitation on photosynthesis +# (via Vcmax reductions), then the btran factor is 1 +btran_nolimit = 1.0 + +# Respiration scaler at canopy top +rdark_scaler_top = 1.0 + +# Nitrogen scaler at canopy top +nscaler_top = 1.0 + + +# Create aliases for the ctype Fortran objects +# ======================================================================================= + +exec(open("CtypesLeafBiophys.py").read()) + + +# Subroutines +# ======================================================================================= + + +def GetJmaxKp25Top(vcmax25_top): + + # Calculate Jmax and Kp at the canopy top at 25C + # they scale off of vcmax + # + # jmax25_top: Canopy top maximum electron transport + # rate at 25C (umol electrons/m**2/s) + # + # kp25top Canopy top initial slope of CO2 response + # curve (C4 plants) at 25C + + jmax25_top = 1.67 * vcmax25_top + kp25_top = 20000. * vcmax25_top + + # q10 response of product limited psn. + # co2_rcurve_islope = co2_rcurve_islope25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) + + return jmax25_top, kp25_top + + + +def main(argv): + + parser = argparse.ArgumentParser(description='Parse command line arguments to this script.') + parser.add_argument('--fin', dest='xmlfile', type=str, help="XML control file, required.", required=True) + parser.add_argument('--smoketest', action='store_true') + args = parser.parse_args() + + # Load the xml control file + # ----------------------------------------------------------------------------------- + xmlroot = et.parse(args.xmlfile).getroot() + + # We will allocate 1 token pft to hold data, we will change the values as needed + numpft = 1 + + # Allocating parameters + print('Allocating parameter space for {} pfts'.format(numpft)) + iret = f90_alloc_leaf_param_sub(ci(numpft)) + + + + # Push scalar parameters + print('Pushing parameters from the xml file to the f90 lb_params datastructure') + scalar_root = xmlroot.find('f90_params').find('scalar_dim') + for param in scalar_root.iter('param'): + iret = f90_set_leaf_param_sub(c8(float(param.text.split(',')[0])),ci(0),*ccharnb(param.attrib['name'].strip())) + + # Push pft parameters to fortran instantiations + pft_root = xmlroot.find('f90_params').find('pft_dim') + leaf_c3psn = [] + leaf_stomatal_intercept = [] + for param in pft_root.iter('param'): + for pft in range(numpft): + iret = f90_set_leaf_param_sub(c8(float(param.text.split(',')[pft])),ci(pft+1),*ccharnb(param.attrib['name'].strip())) + if(param.attrib['name'].strip() == 'fates_leaf_c3psn'): + leaf_c3psn.append(int(param.text.split(',')[pft])) + if(param.attrib['name'].strip() == 'fates_leaf_stomatal_intercept'): + leaf_stomatal_intercept.append(int(param.text.split(',')[pft])) + + # Read in non-fortran parameters from the xml file + fates_leaf_vcmax25top = [] + fates_stoich_nitr = [] + fates_leaf_slatop = [] + + print('Reading non-fortran pft parameters') + py_pft_root = xmlroot.find('py_params').find('pft_dim') + for param in py_pft_root.iter('param'): + for pft in range(numpft): + if (param.attrib['name']=='fates_leaf_vcmax25top'): + fates_leaf_vcmax25top.append(np.float64(param.text.split(',')[pft])) + if (param.attrib['name']=='fates_leaf_slatop'): + fates_leaf_slatop.append(np.float64(param.text.split(',')[pft])) + if (param.attrib['name']=='fates_stoich_nitr'): + fates_stoich_nitr.append(np.float64(param.text.split(',')[pft])) + + print('Reading non-fortran scalar parameters') + py_scalar_root = xmlroot.find('py_params').find('scalar_dim') + for param in py_scalar_root.iter('param'): + if (param.attrib['name']=='fates_maintresp_leaf_model'): + fates_maintresp_leaf_model = int(param.text.split(',')[0]) + + + # Axes to test: + # Temperature + # PAR + # Humidity + # Boundary Layer conductance + # BTRAN + # vcmax (via btran) + # stomatal intercept (via btran) + # stomatal slope 1 (via btran) + # stomatal slope 2 (via btran/medlyn) + # Switches to test + # Stomtal Model (medlyn and BB) + # C3/C4 + # fates_leaf_stomatal_btran_model (0,1,2,3) + # fates_leaf_agross_btran_model + + # Axes implicitly tested + # daylength factors (scale vcmax, but covered by btran scaling) + # + + leaf_stomatal_btran_models = [0,1,2,3,4,5] #4,5 medlyn only + leaf_agross_btran_models = [0,1,2] + + + # Leaf temperature ranges [C] + leaf_tempc_min = -50.0 + leaf_tempc_max = 60.0 + leaf_tempc_n = 15 + leaf_tempc_vec = np.linspace(leaf_tempc_min,leaf_tempc_max,num=leaf_tempc_n) + + # Relative Humidity Ranges + rh_max = 1.00 + rh_min = 0.001 + rh_n = 15 + rh_vec = np.linspace(rh_min,rh_max,num=rh_n) + + # Absorbed PAR ranges [W/m2] + par_abs_min = 0.0 + par_abs_max = 800 + par_abs_n = 15 + par_abs_vec = np.linspace(par_abs_min,par_abs_max,num=par_abs_n) + + # Boundary Conductance ranges [umol/m2/s] + gb_min = 50000.0 # Lower limit imposed by CLM/ELM 0.5 mol/m2/s + gb_max = 8500000.0 # 50% larger than Roughly largestthe largest values seen at BCI (which are 2.5mol/m2/s) + gb_n = 10 + gb_vec = np.linspace(gb_min,gb_max,num=gb_n) + + # btran ranges + + btran_n = 15 + btran_min = -5 + btran_max = 0 + btran_vec = np.logspace(btran_min,btran_max,num=btran_n) + + # vcmax25top ranges + vcmax25t_n = 10 + vcmax25t_min = 1 + vcmax25t_max = 250 + vcmax25t_vec = np.linspace(vcmax25t_min,vcmax25t_max,num=vcmax25t_n) + + + # Set convergence tolerance + ci_tol = 0.1 + + # Initialize fortran return values + # these are all temps and doubles + vcmax_f = c_double(-9.0) + jmax_f = c_double(-9.0) + kp_f = c_double(-9.0) + agross_f = c_double(-9.0) + gstoma_f = c_double(-9.0) + anet_f = c_double(-9.0) + lmr_f = c_double(-9.0) + c13_f = c_double(-9.0) + ac_f = c_double(-9.0) + aj_f = c_double(-9.0) + ap_f = c_double(-9.0) + co2_interc_f = c_double(-9.0) + veg_qs_f = c_double(-9.0) + veg_es_f = c_double(-9.0) + mm_kco2_f = c_double(-9.0) + mm_ko2_f = c_double(-9.0) + co2_cpoint_f = c_double(-9.0) + qsdt_dummy_f = c_double(-9.0) + esdt_dummy_f = c_double(-9.0) + solve_iter_f = c_int(-9) + gs0_f = c_double(-9.0) + gs1_f = c_double(-9.0) + gs2_f = c_double(-9.0) + + pfails = 0 + ptests = 0 + printfail = True + + # unit conversion of W/m2 to umol photons/m^2/s + wm2_to_umolm2s = 4.6 + + print('\nStarting smoke tests with the following nested combinations:\n') + + print(' 2 conductance models, Medlyn and Ball-Berry') + print(' 2 pathway (C3/C4) models') + print(' {} stomatal btran options from {} to {}'.format(len(leaf_stomatal_btran_models),leaf_stomatal_btran_models[0],leaf_stomatal_btran_models[-1])) + print(' {} agross btran options from {} to {}'.format(len(leaf_agross_btran_models),leaf_agross_btran_models[0],leaf_agross_btran_models[-1])) + print(' {} leaf temperature values [C] from {} to {}'.format(leaf_tempc_n,leaf_tempc_min,leaf_tempc_max)) + print(' {} RH values [fraction] from {} to {}'.format(rh_n,rh_min,rh_max)) + print(' {} PAR Abs [W/m2] values from {} to {}'.format(par_abs_n,par_abs_min,par_abs_max)) + print(' {} BL conductance (gb) [umol/m2/s] values from {} to {}'.format(gb_n,gb_min,gb_max)) + print(' {} BTRAN values [fraction] from {} to {}'.format(btran_n,np.exp(btran_min),np.exp(btran_max))) + print(' {} Vcmax 25 top values [umol/m2/s] from {} to {}'.format(vcmax25t_n,vcmax25t_min,vcmax25t_max)) + + #print(' {} values from {} to {}'.format()) + + ntests = 2*len(leaf_agross_btran_models)*vcmax25t_n*leaf_tempc_n*btran_n*gb_n*par_abs_n*rh_n*(4 + 6) + + # Report every 5% + ntestmod = int(ntests/20) + + + print('\nRunning a total of {} tests: \n'.format(ntests)) + time0 = time.process_time() + # Switch between Medlyn and BB + for ism in [1,2]: + + # Push the conductance choice to the fortran code + iret = f90_set_leaf_param_sub(c8(float( ism )),ci(0),*ccharnb('fates_leaf_stomatal_model')) + + if ism==1: leaf_stomatal_btran_models = [0,1,2,3] + else: + leaf_stomatal_btran_models = [0,1,2,3,4,5] #4,5 medlyn only + + for isb in leaf_stomatal_btran_models: + + iret = f90_set_leaf_param_sub(c8(float( isb )),ci(1),*ccharnb('fates_leaf_stomatal_btran_model')) + + for iab in leaf_agross_btran_models: + + iret = f90_set_leaf_param_sub(c8(float( iab )),ci(1),*ccharnb('fates_leaf_agross_btran_model')) + + for ic3 in [0,1]: + + iret = f90_set_leaf_param_sub(c8(float( ic3 )),ci(1),*ccharnb('fates_leaf_c3psn')) + + for vcmax25_top in vcmax25t_vec: + + jmax25_top,kp25_top = GetJmaxKp25Top(vcmax25_top) + + for leaf_tempc in leaf_tempc_vec: + + leaf_tempk = leaf_tempc + tfrz_1atm + + iret = f90_qsat_sub(c8(leaf_tempk),c8(can_press_1atm), \ + byref(veg_qs_f),byref(veg_es_f), \ + byref(qsdt_dummy_f),byref(esdt_dummy_f)) + + iret = f90_cangas_sub(c8(can_press_1atm), \ + c8(o2_ppress_209kppm), \ + c8(leaf_tempk), \ + byref(mm_kco2_f), \ + byref(mm_ko2_f), \ + byref(co2_cpoint_f)) + + # Leaf Nitrogen Concentration at the top + lnc_top = fates_stoich_nitr[0]/fates_leaf_slatop[0] + + # Leaf Maintenance Respiration (temp and pft dependent) + if(fates_maintresp_leaf_model==1): + iret = f90_lmr_ryan_sub(c8(lnc_top),c8(nscaler_top), ci(1), c8(leaf_tempk), byref(lmr_f)) + elif(fates_maintresp_leaf_model==2): + iret = f90_lmr_atkin_sub(c8(lnc_top),c8(rdark_scaler_top),c8(leaf_tempk),c8(atkin_mean_leaf_tempk),byref(lmr_f) ) + else: + print('unknown leaf respiration model') + exit(1) + + for btran in btran_vec: + iret = f90_biophysrate_sub(ci(1), \ + c8(vcmax25_top), c8(jmax25_top), c8(kp25_top), \ + c8(nscaler_top), c8(leaf_tempk), c8(dayl_factor_full), \ + c8(t_growth_kum), c8(t_home_kum), c8(btran), \ + byref(vcmax_f), byref(jmax_f), byref(kp_f), byref(gs0_f), byref(gs1_f), byref(gs2_f)) + + for gb in gb_vec: + for par_abs in par_abs_vec: + par_abs_umol = par_abs*wm2_to_umolm2s + for rh in rh_vec: + vpress = rh * veg_es_f.value + ptests = ptests + 1 + try: + iret = f90_leaflayerphoto_sub(c8(par_abs_umol), \ + c8(par_abs_umol), \ + c8(1.0), \ + ci(1), \ + c8(vcmax_f.value), \ + c8(jmax_f.value), \ + c8(kp_f.value), \ + c8(gs0_f.value), \ + c8(gs1_f.value), \ + c8(gs2_f.value), \ + c8(leaf_tempk), \ + c8(can_press_1atm), \ + c8(co2_ppress_400ppm), \ + c8(o2_ppress_209kppm), \ + c8(gb), \ + c8(vpress), \ + c8(mm_kco2_f.value), \ + c8(mm_ko2_f.value), \ + c8(co2_cpoint_f.value), \ + c8(lmr_f.value), \ + c8(ci_tol), \ + byref(agross_f), \ + byref(gstoma_f), \ + byref(anet_f), \ + byref(c13_f), \ + byref(ac_f), \ + byref(aj_f), \ + byref(ap_f), \ + byref(co2_interc_f), \ + byref(solve_iter_f) ) + except: + pfails = pfails+1 + printfail=True + + if (np.mod(ptests,ntestmod)==0): + print('Completed {} tests -- {} percent complete'.format(ptests,100*float(ptests)/float(ntests))) + + if (pfails>0 and np.mod(pfails,100)==0 and printfail): + printfail=False + print('\n{} fails so far\n'.format(pfails)) + + + + print("\nCompleted Photosynthesis Smoke Test\n") + print("\nElapsed Time [s]: {}\n".format(time.process_time() - time0)) + print("{} Failures out of {} Encountered; {}% of Tests\n".format(pfails,ptests,float(pfails)/float(ptests))) + + print('Deallocating parameter space') + iret = f90_dealloc_leaf_param_sub() + + print('Functional Unit Testing Complete') + exit(0) + + +# ======================================================================================= +# This is the actual call to main + +if __name__ == "__main__": + main(sys.argv) diff --git a/functional_unit_testing/leaf_biophys/bld/README b/functional_unit_testing/leaf_biophys/bld/README new file mode 100644 index 0000000000..eba79cab57 --- /dev/null +++ b/functional_unit_testing/leaf_biophys/bld/README @@ -0,0 +1 @@ +Placeholder for directory diff --git a/functional_unit_testing/leaf_biophys/build_leafbiophys_objects.sh b/functional_unit_testing/leaf_biophys/build_leafbiophys_objects.sh new file mode 100755 index 0000000000..16e7813ddb --- /dev/null +++ b/functional_unit_testing/leaf_biophys/build_leafbiophys_objects.sh @@ -0,0 +1,29 @@ +#!/bin/bash + +# Path to FATES src + +FC='gfortran' + +F_OPTS="-shared -fPIC -O3" + +#F_OPTS="-shared -fPIC -O0 -g -ffpe-trap=zero,overflow,underflow -fbacktrace -fbounds-check -Wall" + +MOD_FLAG="-J" + +rm -f bld/*.o +rm -f bld/*.mod + +# Build the new file with constants + +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/FatesConstantsMod.o ../../main/FatesConstantsMod.F90 +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/WrapShrMod.o f90_src/WrapShrMod.F90 +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/FatesUtilsMod.o ../../main/FatesUtilsMod.F90 +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/LeafBiophysicsMod.o ../../biogeophys/LeafBiophysicsMod.F90 +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/LeafBiophysSuppMod.o f90_src/LeafBiophysSuppMod.F90 +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/FatesRadiationMemMod.o ../../radiation/FatesRadiationMemMod.F90 +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/TwoStreamMLPEMod.o ../../radiation/TwoStreamMLPEMod.F90 +${FC} ${F_OPTS} -I bld/ ${MOD_FLAG} bld/ -o bld/RadiationWrapMod.o ../radiation/f90_src/RadiationWrapMod.F90 + + + + diff --git a/functional_unit_testing/leaf_biophys/f90_src/LeafBiophysSuppMod.F90 b/functional_unit_testing/leaf_biophys/f90_src/LeafBiophysSuppMod.F90 new file mode 100644 index 0000000000..94c13ae531 --- /dev/null +++ b/functional_unit_testing/leaf_biophys/f90_src/LeafBiophysSuppMod.F90 @@ -0,0 +1,150 @@ +module LeafBiophysSuppMod + + use iso_c_binding, only : c_char + use iso_c_binding, only : c_int + use iso_c_binding, only : r8 => c_double + + use LeafBiophysicsMod, only : lb_params + + implicit none + public + save + + integer(kind=c_int), parameter :: param_string_length = 32 + logical, parameter :: debug = .true. + +contains + + subroutine AllocLeafParam(numpft) + + integer, intent(in) :: numpft + + allocate(lb_params%c3psn(numpft)) + allocate(lb_params%stomatal_btran_model(numpft)) + allocate(lb_params%agross_btran_model(numpft)) + allocate(lb_params%medlyn_slope(numpft)) + allocate(lb_params%bb_slope(numpft)) + allocate(lb_params%stomatal_intercept(numpft)) + allocate(lb_params%maintresp_leaf_ryan1991_baserate(numpft)) + allocate(lb_params%maintresp_leaf_atkin2017_baserate(numpft)) + allocate(lb_params%maintresp_reduction_curvature(numpft)) + allocate(lb_params%maintresp_reduction_intercept(numpft)) + allocate(lb_params%maintresp_reduction_upthresh(numpft)) + allocate(lb_params%vcmaxha(numpft)) + allocate(lb_params%jmaxha(numpft)) + allocate(lb_params%vcmaxhd(numpft)) + allocate(lb_params%jmaxhd(numpft)) + allocate(lb_params%vcmaxse(numpft)) + allocate(lb_params%jmaxse(numpft)) + + return + end subroutine AllocLeafParam + + subroutine DeallocLeafParam() + + deallocate(lb_params%c3psn) + deallocate(lb_params%stomatal_btran_model) + deallocate(lb_params%agross_btran_model) + deallocate(lb_params%medlyn_slope) + deallocate(lb_params%bb_slope) + deallocate(lb_params%stomatal_intercept) + deallocate(lb_params%maintresp_leaf_ryan1991_baserate) + deallocate(lb_params%maintresp_leaf_atkin2017_baserate) + deallocate(lb_params%maintresp_reduction_curvature) + deallocate(lb_params%maintresp_reduction_intercept) + deallocate(lb_params%maintresp_reduction_upthresh) + deallocate(lb_params%vcmaxha) + deallocate(lb_params%jmaxha) + deallocate(lb_params%vcmaxhd) + deallocate(lb_params%jmaxhd) + deallocate(lb_params%vcmaxse) + deallocate(lb_params%jmaxse) + + end subroutine DeallocLeafParam + + ! ===================================================================================== + + subroutine SetLeafParam(val,pft,pname) + + real(r8), intent(in) :: val + character(kind=c_char,len=*), intent(in) :: pname + integer(kind=c_int), intent(in) :: pft + + select case(trim(pname)) + case('fates_daylength_factor_switch') + lb_params%dayl_switch = nint(val) + case('fates_leaf_stomatal_model') + lb_params%stomatal_model = nint(val) + case('fates_leaf_stomatal_assim_model') + lb_params%stomatal_assim_model = nint(val) + case('fates_leaf_photo_tempsens_model') + lb_params%photo_tempsens_model = nint(val) + case('fates_leaf_c3psn') + lb_params%c3psn(pft) = nint(val) + case('fates_leaf_stomatal_btran_model') + lb_params%stomatal_btran_model(pft) = nint(val) + case('fates_leaf_agross_btran_model') + lb_params%agross_btran_model(pft) = nint(val) + case('fates_leaf_stomatal_slope_ballberry') + lb_params%bb_slope(pft) = val + case('fates_leaf_stomatal_slope_medlyn') + lb_params%medlyn_slope(pft) = val + case('fates_leaf_stomatal_intercept') + lb_params%stomatal_intercept(pft) = val + case('fates_maintresp_reduction_curvature') + lb_params%maintresp_reduction_curvature(pft) = val + case('fates_maintresp_reduction_intercept') + lb_params%maintresp_reduction_intercept(pft) = val + case('fates_maintresp_reduction_upthresh') + lb_params%maintresp_reduction_upthresh(pft) = val + case('fates_maintresp_leaf_atkin2017_baserate') + lb_params%maintresp_leaf_atkin2017_baserate(pft) = val + case('fates_maintresp_leaf_ryan1991_baserate') + lb_params%maintresp_leaf_ryan1991_baserate(pft) = val + case('fates_leaf_vcmaxha') + lb_params%vcmaxha(pft) = val + case('fates_leaf_jmaxha') + lb_params%jmaxha(pft) = val + case('fates_leaf_vcmaxhd') + lb_params%vcmaxhd(pft) = val + case('fates_leaf_jmaxhd') + lb_params%jmaxhd(pft) = val + case('fates_leaf_vcmaxse') + lb_params%vcmaxse(pft) = val + case('fates_leaf_jmaxse') + lb_params%jmaxse(pft) = val + case default + print*,"An unknown parameter name was sent to the parameter" + print*,"initialization function." + print*,"name:--",trim(pname),"--" + stop + end select + + end subroutine SetLeafParam + + subroutine DumpParams() + + print*,"Leaf Biophys Parameters In Structure lb_params: " + print*,'fates_daylength_factor_switch: ',lb_params%dayl_switch + print*,'fates_leaf_stomatal_model: ',lb_params%stomatal_model + print*,'fates_leaf_stomatal_assim_model: ',lb_params%stomatal_assim_model + print*,'fates_leaf_photo_tempsens_model: ',lb_params%photo_tempsens_model + print*,'fates_leaf_c3psn: ',lb_params%c3psn + print*,'fates_leaf_stomatal_slope_ballberry: ',lb_params%bb_slope + print*,'fates_leaf_stomatal_slope_medlyn: ',lb_params%medlyn_slope + print*,'fates_leaf_stomatal_intercept: ',lb_params%stomatal_intercept + print*,'fates_maintresp_reduction_curvature: ',lb_params%maintresp_reduction_curvature + print*,'fates_maintresp_reduction_intercept: ',lb_params%maintresp_reduction_intercept + print*,'fates_maintresp_reduction_upthresh: ',lb_params%maintresp_reduction_upthresh + print*,'fates_maintresp_leaf_atkin2017_baserate: ',lb_params%maintresp_leaf_atkin2017_baserate + print*,'fates_maintresp_leaf_ryan1991_baserate: ',lb_params%maintresp_leaf_ryan1991_baserate + print*,'fates_leaf_vcmaxha: ',lb_params%vcmaxha + print*,'fates_leaf_jmaxha: ',lb_params%jmaxha + print*,'fates_leaf_vcmaxhd: ',lb_params%vcmaxhd + print*,'fates_leaf_jmaxhd: ',lb_params%jmaxhd + print*,'fates_leaf_vcmaxse: ',lb_params%vcmaxse + print*,'fates_leaf_jmaxse: ',lb_params%jmaxse + + end subroutine DumpParams + +end module LeafBiophysSuppMod diff --git a/functional_unit_testing/leaf_biophys/f90_src/WrapShrMod.F90 b/functional_unit_testing/leaf_biophys/f90_src/WrapShrMod.F90 new file mode 100644 index 0000000000..80f3e41d0f --- /dev/null +++ b/functional_unit_testing/leaf_biophys/f90_src/WrapShrMod.F90 @@ -0,0 +1,108 @@ +module shr_log_mod + use iso_c_binding, only : c_char + use iso_c_binding, only : c_int + + public :: shr_log_errMsg + + contains + function shr_log_errMsg(source, line) result(ans) + character(kind=c_char,len=*), intent(in) :: source + integer(c_int), intent(in) :: line + character(kind=c_char,len=4) :: cline ! character version of int + character(kind=c_char,len=128) :: ans + + write(cline,'(I4)') line + ans = "source: " // trim(source) // " line: "// trim(cline) + + end function shr_log_errMsg + +end module shr_log_mod + +module shr_sys_mod + + public :: shr_sys_abort + +contains + + subroutine shr_sys_abort + call exit(0) + end subroutine shr_sys_abort + +end module shr_sys_mod + +module FatesGlobals + + use iso_c_binding, only : c_char + use iso_c_binding, only : c_int + use FatesConstantsMod, only : r8 => fates_r8 + + integer :: stdo_unit = 6 + +contains + + integer function fates_log() + fates_log = -1 + end function fates_log + + subroutine fates_endrun(msg) + + implicit none + character(len=*), intent(in) :: msg ! string to be printed + + write(stdo_unit,*) msg + + stop + + end subroutine fates_endrun + + subroutine FatesWarn(msg,index) + + character(len=*), intent(in) :: msg ! string to be printed + integer,optional,intent(in) :: index ! warning index + + + write(stdo_unit,*) index, msg + + stop + + end subroutine FatesWarn + + ! ===================================================================================== + + function N2S(real_in) result(str) + + real(r8) :: real_in + character(len=16) :: str + + !write(str,*) real_in + write(str,'(E12.6)') real_in + + end function N2S + + ! ===================================================================================== + + function I2S(int_in) result(str) + + integer :: int_in + character(len=16) :: str + + !write(str,*) real_in + write(str,'(I15)') int_in + + end function I2S + + ! ===================================================================================== + + function A2S(reals_in) result(str) + + real(r8) :: reals_in(:) + character(len=1024) :: str + integer :: i + + str = ', ' + do i = 1,ubound(reals_in,1) + str = trim(str)//', '//N2S(reals_in(i)) + end do + + end function A2S +end module FatesGlobals diff --git a/functional_unit_testing/leaf_biophys/images/README b/functional_unit_testing/leaf_biophys/images/README new file mode 100644 index 0000000000..3a61e9b89b --- /dev/null +++ b/functional_unit_testing/leaf_biophys/images/README @@ -0,0 +1 @@ +Placeholder for folder \ No newline at end of file diff --git a/functional_unit_testing/leaf_biophys/leaf_biophys_controls.xml b/functional_unit_testing/leaf_biophys/leaf_biophys_controls.xml new file mode 100644 index 0000000000..4fb11aa465 --- /dev/null +++ b/functional_unit_testing/leaf_biophys/leaf_biophys_controls.xml @@ -0,0 +1,46 @@ + + + + + base_pfts.pdf + 12 + + + 1 + 2 + 1 + 1 + + + 2,2,2,2,2,2,2,2,2,2,2,2 + 2,0,0,0,0,0,0,0,0,0,0,0 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0 + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8 + 4.1, 2.3, 2.3, 4.1, 4.4, 4.4, 4.7, 4.7, 4.7, 2.2, 5.3, 1.6 + 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 40000 + 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 + 1.756, 1.4995, 1.4995, 1.756, 1.756, 1.756, 2.0749, 2.0749, 2.0749, 2.1956, 2.1956, 2.1956 + 2.525e-06, 2.525e-06, 2.525e-06, 2.525e-06, 2.525e-06, 2.525e-06, 2.525e-06, 2.525e-06, 2.525e-06, 2.525e-06, 2.525e-06, 2.525e-06 + 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330 + 43540, 43540, 43540, 43540, 43540, 43540, 43540, 43540, 43540, 43540, 43540, 43540 + 149250, 149250, 149250, 149250, 149250, 149250, 149250, 149250, 149250, 149250, 149250, 149250 + 149250, 149250, 149250, 149250, 149250, 149250, 149250, 149250, 149250, 149250, 149250, 149250 + 485, 485, 485, 485, 485, 485, 485, 485, 485, 485, 485, 485 + 495, 495, 495, 495, 495, 495, 495, 495, 495, 495, 495, 495 + + + + + 1 + + + 0.00963, 0.00963, 0.00963, 0.00963, 0.00963, 0.00963, 0.00963, 0.00963, 0.00963, 0.00963, 0.00963, 0.00963 + 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43, 2.43 + 50, 62, 39, 61, 58, 58, 62, 54, 54, 78, 78, 78 + 0.033, 0.029, 0.04, 0.033, 0.04, 0.04, 0.033, 0.04, 0.04, 0.04, 0.04, 0.04 + 0.012, 0.005, 0.024, 0.009, 0.03, 0.03, 0.012, 0.03, 0.03, 0.03, 0.03, 0.03 + + + diff --git a/functional_unit_testing/radiation/RadiationUTestDriver.py b/functional_unit_testing/radiation/RadiationUTestDriver.py index b7c32d74d3..76aba42136 100644 --- a/functional_unit_testing/radiation/RadiationUTestDriver.py +++ b/functional_unit_testing/radiation/RadiationUTestDriver.py @@ -20,7 +20,7 @@ import os import sys import getopt -import code # For development: code.interact(local=locals()) code.interact(local=dict(globals(), **locals())) +import code import time import importlib import csv @@ -50,7 +50,8 @@ alloc_radparams_call = f90_twostr_obj.__twostreammlpemod_MOD_allocateradparams set_radparams_call = f90_wrap_obj.__radiationwrapmod_MOD_setradparam set_radparams_call.argtypes = [POINTER(c_double),POINTER(c_int),POINTER(c_int),c_char_p,c_long] -param_prep_call = f90_twostr_obj.__twostreammlpemod_MOD_radparamprep + +param_prep_call = f90_twostr_obj.__twostreammlpemod_MOD_paramprep setup_canopy_call = f90_wrap_obj.__radiationwrapmod_MOD_setupcanopy setup_canopy_call.argtypes = [POINTER(c_int),POINTER(c_int),POINTER(c_int), \ diff --git a/functional_unit_testing/radiation/f90_src/RadiationWrapMod.F90 b/functional_unit_testing/radiation/f90_src/RadiationWrapMod.F90 index 90c31d21a8..ca0f102c62 100644 --- a/functional_unit_testing/radiation/f90_src/RadiationWrapMod.F90 +++ b/functional_unit_testing/radiation/f90_src/RadiationWrapMod.F90 @@ -189,7 +189,6 @@ subroutine WrapSolve(ib,boundary_type,Rbeam_atm,Rdiff_atm, & ipiv, & albedo_beam, & albedo_diff, & - err_solve, & err_consv, & frac_abs_can_beam, & frac_abs_can_diff, & diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index cc906fecef..33015b7fcc 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -46,9 +46,7 @@ module EDParamsMod real(r8),protected, public :: sdlng2sap_par_timescale !Length of the window for the exponential !moving average of par at the seedling layer used to !calculate seedling to sapling transition rates - integer,protected, public :: photo_tempsens_model ! switch for choosing the model that defines the temperature - ! sensitivity of photosynthetic parameters (vcmax, jmax). - ! 1=non-acclimating, 2=Kumarathunge et al., 2019 + integer,protected, public :: radiation_model ! Switch betrween Norman (1) and Two-stream (2) radiation models @@ -76,8 +74,6 @@ module EDParamsMod real(r8),protected, public :: ED_val_cohort_age_fusion_tol ! minimum fraction in differece in cohort age between cohorts real(r8),protected, public :: ED_val_patch_fusion_tol ! minimum fraction in difference in profiles between patches 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 - integer,protected, public :: dayl_switch ! switch for turning on or off day length factor scaling for photosynthetic parameters integer,protected, public :: regeneration_model ! Switch for choosing between regeneration models: ! (1) for Fates default ! (2) for the Tree Recruitment Scheme (Hanbury-Brown et al., 2022) @@ -90,11 +86,7 @@ module EDParamsMod real(r8), protected, public :: cg_strikes ! fraction of cloud to ground lightning strikes (0-1) character(len=param_string_length),parameter :: fates_name_cg_strikes="fates_fire_cg_strikes" - ! empirical curvature parameters for ac, aj photosynthesis co-limitation, c3 and c4 plants respectively - real(r8),protected,public :: theta_cj_c3 ! Empirical curvature parameter for ac, aj photosynthesis co-limitation in c3 plants - real(r8),protected,public :: theta_cj_c4 ! Empirical curvature parameter for ac, aj photosynthesis co-limitation in c4 plants - - ! Global identifier of how nutrients interact with the host land model + ! Global identifier of how nutrients interact with the host land model ! either they are fully coupled, or they generate uptake rates synthetically ! in prescribed mode. In the latter, there is both NO mass removed from the HLM's soil ! BGC N and P pools, and there is also none removed. @@ -144,7 +136,6 @@ module EDParamsMod integer, protected,allocatable,public :: hydr_htftype_node(:) character(len=param_string_length),parameter,public :: ED_name_photo_temp_acclim_timescale = "fates_leaf_photo_temp_acclim_timescale" character(len=param_string_length),parameter,public :: ED_name_photo_temp_acclim_thome_time = "fates_leaf_photo_temp_acclim_thome_time" - character(len=param_string_length),parameter,public :: name_photo_tempsens_model = "fates_leaf_photo_tempsens_model" character(len=param_string_length),parameter,public :: name_maintresp_model = "fates_maintresp_leaf_model" character(len=param_string_length),parameter,public :: name_radiation_model = "fates_rad_model" character(len=param_string_length),parameter,public :: ED_name_hydr_htftype_node = "fates_hydro_htftype_node" @@ -169,13 +160,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" - character(len=param_string_length),parameter,public :: ED_name_dayl_switch= "fates_daylength_factor_switch" character(len=param_string_length),parameter,public :: ED_name_regeneration_model= "fates_regeneration_model" - - character(len=param_string_length),parameter,public :: name_theta_cj_c3 = "fates_leaf_theta_cj_c3" - character(len=param_string_length),parameter,public :: name_theta_cj_c4 = "fates_leaf_theta_cj_c4" - character(len=param_string_length),parameter :: fates_name_q10_mr="fates_q10_mr" character(len=param_string_length),parameter :: fates_name_q10_froz="fates_q10_froz" @@ -222,10 +207,6 @@ module EDParamsMod real(r8),protected,public :: bgc_soil_salinity ! site-level soil salinity for FATES when not coupled to dynamic soil BGC of salinity character(len=param_string_length),parameter,public :: bgc_name_soil_salinity= "fates_soil_salinity" - ! Switch designating whether to use net or gross assimilation in the stomata model - integer, protected, public :: stomatal_assim_model - character(len=param_string_length), parameter, public :: stomatal_assim_name = "fates_leaf_stomatal_assim_model" - ! Integer code that options how damage events are structured integer, protected, public :: damage_event_code character(len=param_string_length), parameter, public :: damage_name_event_code = "fates_damage_event_code" @@ -320,7 +301,6 @@ subroutine FatesParamsInit() sdlng_mdd_timescale = nan sdlng2sap_par_timescale = nan photo_temp_acclim_thome_time = nan - photo_tempsens_model = -9 maintresp_leaf_model = -9 radiation_model = -9 fates_mortality_disturbance_fraction = nan @@ -344,10 +324,7 @@ subroutine FatesParamsInit() ED_val_cohort_age_fusion_tol = nan ED_val_patch_fusion_tol = nan ED_val_canopy_closure_thresh = nan - stomatal_model = -9 - dayl_switch = -9 regeneration_model = -9 - stomatal_assim_model = -9 max_cohort_per_patch = -9 hydr_kmax_rsurf1 = nan hydr_kmax_rsurf2 = nan @@ -366,8 +343,6 @@ subroutine FatesParamsInit() eca_plant_escalar = nan q10_mr = nan q10_froz = nan - theta_cj_c3 = nan - theta_cj_c4 = nan dev_arbitrary = nan damage_event_code = -9 damage_canopy_layer_code = -9 @@ -419,21 +394,12 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=ED_name_photo_temp_acclim_thome_time, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=name_photo_tempsens_model,dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=name_radiation_model,dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) call fates_params%RegisterParameter(name=name_maintresp_model,dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=name_theta_cj_c3, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - - call fates_params%RegisterParameter(name=name_theta_cj_c4, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=ED_name_mort_disturb_frac, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -497,18 +463,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=ED_name_dayl_switch, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=ED_name_regeneration_model, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=stomatal_assim_name, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=maxcohort_name, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -647,10 +604,6 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetrieveParameter(name=ED_name_photo_temp_acclim_thome_time, & data=photo_temp_acclim_thome_time) - call fates_params%RetrieveParameter(name=name_photo_tempsens_model, & - data=tmpreal) - photo_tempsens_model = nint(tmpreal) - call fates_params%RetrieveParameter(name=name_radiation_model, & data=tmpreal) radiation_model = nint(tmpreal) @@ -723,22 +676,10 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetrieveParameter(name=ED_name_canopy_closure_thresh, & data=ED_val_canopy_closure_thresh) - call fates_params%RetrieveParameter(name=ED_name_stomatal_model, & - data=tmpreal) - stomatal_model = nint(tmpreal) - - call fates_params%RetrieveParameter(name=ED_name_dayl_switch, & - data=tmpreal) - dayl_switch = nint(tmpreal) - call fates_params%RetrieveParameter(name=ED_name_regeneration_model, & data=tmpreal) regeneration_model = nint(tmpreal) - call fates_params%RetrieveParameter(name=stomatal_assim_name, & - data=tmpreal) - stomatal_assim_model = nint(tmpreal) - call fates_params%RetrieveParameter(name=maxcohort_name, & data=tmpreal) max_cohort_per_patch = nint(tmpreal) @@ -791,12 +732,6 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetrieveParameter(name=eca_name_plant_escalar, & data=eca_plant_escalar) - - call fates_params%RetrieveParameter(name=name_theta_cj_c3, & - data=theta_cj_c3) - - call fates_params%RetrieveParameter(name=name_theta_cj_c4, & - data=theta_cj_c4) call fates_params%RetrieveParameter(name=fates_name_q10_mr, & data=q10_mr) @@ -909,9 +844,6 @@ subroutine FatesReportParams(is_master) 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) 'regeneration_model = ',regeneration_model - write(fates_log(),fmt0) 'dayl_switch = ',dayl_switch - write(fates_log(),fmt0) 'stomatal_model = ',stomatal_model - write(fates_log(),fmt0) 'stomatal_assim_model = ',stomatal_assim_model write(fates_log(),fmt0) 'hydro_kmax_rsurf1 = ',hydr_kmax_rsurf1 write(fates_log(),fmt0) 'hydro_kmax_rsurf2 = ',hydr_kmax_rsurf2 write(fates_log(),fmt0) 'hydro_psi0 = ',hydr_psi0 diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 3d91f2f1b9..07c6e46a40 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -13,6 +13,9 @@ module EDPftvarcon use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : itrue, ifalse use PRTParametersMod, only : prt_params + use LeafBiophysicsMod, only : lb_params + use LeafBiophysicsMod, only : btran_on_gs_gs02,btran_on_ag_vcmax_jmax + use LeafBiophysicsMod, only : lmr_r_1, lmr_r_2 use FatesGlobals, only : fates_log use FatesGlobals, only : endrun => fates_endrun use FatesLitterMod, only : ilabile,icellulose,ilignin @@ -56,10 +59,6 @@ module EDPftvarcon real(r8), allocatable :: initd(:) ! initial seedling density real(r8), allocatable :: seed_suppl(:) ! seeds that come from outside the gridbox. - real(r8), allocatable :: bb_slope(:) ! ball berry slope parameter - real(r8), allocatable :: medlyn_slope(:) ! Medlyn slope parameter KPa^0.5 - real(r8), allocatable :: stomatal_intercept(:) ! intercept of stomatal conductance model - real(r8), allocatable :: lf_flab(:) ! Leaf litter labile fraction [-] real(r8), allocatable :: lf_fcel(:) ! Leaf litter cellulose fraction [-] @@ -71,29 +70,11 @@ module EDPftvarcon real(r8), allocatable :: clumping_index(:) ! factor describing how much self-occlusion ! of leaf scattering elements ! decreases light interception - real(r8), allocatable :: c3psn(:) ! index defining the photosynthetic - ! pathway C4 = 0, C3 = 1 - real(r8), allocatable :: smpso(:) ! Soil water potential at full stomatal opening ! (non-HYDRO mode only) [mm] real(r8), allocatable :: smpsc(:) ! Soil water potential at full stomatal closure ! (non-HYDRO mode only) [mm] - - real(r8), allocatable :: maintresp_reduction_curvature(:) ! curvature of MR reduction as f(carbon storage), - ! 1=linear, 0=very curved - real(r8), allocatable :: maintresp_reduction_intercept(:) ! intercept of MR reduction as f(carbon storage), - ! 0=no throttling, 1=max throttling - real(r8), allocatable :: maintresp_reduction_upthresh (:) ! Upper threshold for storage biomass (relative - ! to leaf biomass) above which MR is not reduced - - real(r8), allocatable :: maintresp_leaf_atkin2017_baserate(:) ! leaf maintenance respiration base rate (r0) - ! per Atkin et al 2017 - - real(r8), allocatable :: maintresp_leaf_ryan1991_baserate(:) ! leaf maintenance respiration per Ryan et al 1991 - - - real(r8), allocatable :: maintresp_leaf_vert_scaler_coeff1(:) ! leaf maintenance respiration decrease through the canopy param 1 ! only with Atkin et al. 2017 respiration model real(r8), allocatable :: maintresp_leaf_vert_scaler_coeff2(:) ! leaf maintenance respiration decrease through the canopy param 2 @@ -109,12 +90,7 @@ module EDPftvarcon real(r8), allocatable :: mort_upthresh_cstarvation(:) ! threshold for storage biomass (relative to target leaf biomass) above which carbon starvation is zero real(r8), allocatable :: hf_sm_threshold(:) ! soil moisture (btran units) at which drought mortality begins for non-hydraulic model real(r8), allocatable :: hf_flc_threshold(:) ! plant fractional loss of conductivity at which drought mortality begins for hydraulic model - real(r8), allocatable :: vcmaxha(:) ! activation energy for vcmax - real(r8), allocatable :: jmaxha(:) ! activation energy for jmax - real(r8), allocatable :: vcmaxhd(:) ! deactivation energy for vcmax - real(r8), allocatable :: jmaxhd(:) ! deactivation energy for jmax - real(r8), allocatable :: vcmaxse(:) ! entropy term for vcmax - real(r8), allocatable :: jmaxse(:) ! entropy term for jmax + real(r8), allocatable :: germination_rate(:) ! Fraction of seed mass germinating per year (yr-1) real(r8), allocatable :: seed_decay_rate(:) ! Fraction of seed mass (both germinated and ! ungerminated), decaying per year (yr-1) @@ -406,18 +382,6 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_leaf_stomatal_slope_ballberry' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_leaf_stomatal_slope_medlyn' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_leaf_stomatal_intercept' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_frag_leaf_flab' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -450,10 +414,6 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_leaf_c3psn' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_nonhydro_smpso' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -462,28 +422,6 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_maintresp_reduction_curvature' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_maintresp_reduction_intercept' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_maintresp_reduction_upthresh' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_maintresp_leaf_atkin2017_baserate' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_maintresp_leaf_ryan1991_baserate' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - - name = 'fates_maintresp_leaf_vert_scaler_coeff1' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -608,30 +546,6 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_leaf_vcmaxha' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_leaf_jmaxha' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_leaf_vcmaxhd' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_leaf_jmaxhd' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_leaf_vcmaxse' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - - name = 'fates_leaf_jmaxse' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_recruit_seed_germination_rate' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -873,18 +787,6 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetrieveParameterAllocate(name=name, & data=this%seed_suppl) - name = 'fates_leaf_stomatal_slope_ballberry' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%bb_slope) - - name = 'fates_leaf_stomatal_slope_medlyn' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%medlyn_slope) - - name = 'fates_leaf_stomatal_intercept' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%stomatal_intercept) - name = 'fates_frag_leaf_flab' call fates_params%RetrieveParameterAllocate(name=name, & data=this%lf_flab) @@ -917,10 +819,6 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetrieveParameterAllocate(name=name, & data=this%clumping_index) - name = 'fates_leaf_c3psn' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%c3psn) - name = 'fates_nonhydro_smpso' call fates_params%RetrieveParameterAllocate(name=name, & data=this%smpso) @@ -929,26 +827,6 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetrieveParameterAllocate(name=name, & data=this%smpsc) - name = 'fates_maintresp_reduction_curvature' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%maintresp_reduction_curvature) - - name = 'fates_maintresp_reduction_intercept' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%maintresp_reduction_intercept) - - name = 'fates_maintresp_reduction_upthresh' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%maintresp_reduction_upthresh) - - name = 'fates_maintresp_leaf_atkin2017_baserate' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%maintresp_leaf_atkin2017_baserate) - - name = 'fates_maintresp_leaf_ryan1991_baserate' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%maintresp_leaf_ryan1991_baserate) - name = 'fates_maintresp_leaf_vert_scaler_coeff1' call fates_params%RetrieveParameterAllocate(name=name, & data=this%maintresp_leaf_vert_scaler_coeff1) @@ -1087,30 +965,6 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetrieveParameterAllocate(name=name, & data=this%hf_flc_threshold) - name = 'fates_leaf_vcmaxha' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%vcmaxha) - - name = 'fates_leaf_jmaxha' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%jmaxha) - - name = 'fates_leaf_vcmaxhd' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%vcmaxhd) - - name = 'fates_leaf_jmaxhd' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%jmaxhd) - - name = 'fates_leaf_vcmaxse' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%vcmaxse) - - name = 'fates_leaf_jmaxse' - call fates_params%RetrieveParameterAllocate(name=name, & - data=this%jmaxse) - name = 'fates_recruit_seed_germination_rate' call fates_params%RetrieveParameterAllocate(name=name, & data=this%germination_rate) @@ -1518,7 +1372,7 @@ subroutine Receive_PFT_leafage(this, fates_params) name = 'fates_leaf_vcmax25top' call fates_params%RetrieveParameterAllocate(name=name, & data=this%vcmax25top) - + return end subroutine Receive_PFT_leafage @@ -1716,9 +1570,6 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'crown_kill = ',EDPftvarcon_inst%crown_kill write(fates_log(),fmt0) 'initd = ',EDPftvarcon_inst%initd write(fates_log(),fmt0) 'seed_suppl = ',EDPftvarcon_inst%seed_suppl - write(fates_log(),fmt0) 'bb_slope = ',EDPftvarcon_inst%bb_slope - write(fates_log(),fmt0) 'medlyn_slope = ',EDPftvarcon_inst%medlyn_slope - write(fates_log(),fmt0) 'stomatal_intercept = ',EDPftvarcon_inst%stomatal_intercept write(fates_log(),fmt0) 'lf_flab = ',EDPftvarcon_inst%lf_flab write(fates_log(),fmt0) 'lf_fcel = ',EDPftvarcon_inst%lf_fcel write(fates_log(),fmt0) 'lf_flig = ',EDPftvarcon_inst%lf_flig @@ -1727,7 +1578,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'fr_flig = ',EDPftvarcon_inst%fr_flig write(fates_log(),fmt0) 'xl = ',EDPftvarcon_inst%xl write(fates_log(),fmt0) 'clumping_index = ',EDPftvarcon_inst%clumping_index - write(fates_log(),fmt0) 'c3psn = ',EDPftvarcon_inst%c3psn + write(fates_log(),fmt0) 'vcmax25top = ',EDPftvarcon_inst%vcmax25top write(fates_log(),fmt0) 'smpso = ',EDPftvarcon_inst%smpso write(fates_log(),fmt0) 'smpsc = ',EDPftvarcon_inst%smpsc @@ -1742,12 +1593,7 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'mort_upthresh_cstarvation = ',EDPftvarcon_inst%mort_upthresh_cstarvation write(fates_log(),fmt0) 'hf_sm_threshold = ',EDPftvarcon_inst%hf_sm_threshold write(fates_log(),fmt0) 'hf_flc_threshold = ',EDPftvarcon_inst%hf_flc_threshold - write(fates_log(),fmt0) 'vcmaxha = ',EDPftvarcon_inst%vcmaxha - write(fates_log(),fmt0) 'jmaxha = ',EDPftvarcon_inst%jmaxha - write(fates_log(),fmt0) 'vcmaxhd = ',EDPftvarcon_inst%vcmaxhd - write(fates_log(),fmt0) 'jmaxhd = ',EDPftvarcon_inst%jmaxhd - write(fates_log(),fmt0) 'vcmaxse = ',EDPftvarcon_inst%vcmaxse - write(fates_log(),fmt0) 'jmaxse = ',EDPftvarcon_inst%jmaxse + write(fates_log(),fmt0) 'germination_timescale = ',EDPftvarcon_inst%germination_rate write(fates_log(),fmt0) 'seed_decay_turnover = ',EDPftvarcon_inst%seed_decay_rate write(fates_log(),fmt0) 'seed_dispersal_pdf_scale = ',EDPftvarcon_inst%seed_dispersal_pdf_scale @@ -1823,11 +1669,10 @@ subroutine FatesCheckParams(is_master) use FatesConstantsMod , only : fates_check_param_set use FatesConstantsMod , only : itrue, ifalse use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - use FatesConstantsMod, only : lmr_r_1 - use FatesConstantsMod, only : lmr_r_2 + use EDParamsMod , only : logging_mechanical_frac, logging_collateral_frac use EDParamsMod , only : logging_direct_frac,logging_export_frac - use EDParamsMod , only : radiation_model, dayl_switch + use EDParamsMod , only : radiation_model use FatesInterfaceTypesMod, only : hlm_use_fixed_biogeog,hlm_use_sp, hlm_name use FatesInterfaceTypesMod, only : hlm_use_inventory_init use FatesInterfaceTypesMod, only : hlm_use_nocomp @@ -1877,10 +1722,10 @@ subroutine FatesCheckParams(is_master) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if(.not.any(dayl_switch == [itrue,ifalse])) then + if(.not.any(lb_params%dayl_switch == [itrue,ifalse])) then write(fates_log(),*) 'The only valid switch options for ' write(fates_log(),*) 'fates_daylength_factor_switch is 0 or 1 ...' - write(fates_log(),*) 'You specified fates_daylength_factor_switch = ',dayl_switch + write(fates_log(),*) 'You specified fates_daylength_factor_switch = ',lb_params%dayl_switch write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -2216,15 +2061,13 @@ subroutine FatesCheckParams(is_master) ! Check if photosynthetic pathway is neither C3/C4 ! ---------------------------------------------------------------------------------- - - if ( ( EDPftvarcon_inst%c3psn(ipft) < 0.0_r8 ) .or. & - ( EDPftvarcon_inst%c3psn(ipft) > 1.0_r8 ) ) then + if(.not.any(lb_params%c3psn(ipft) == [c3_path_index,c4_path_index])) then write(fates_log(),*) ' Two photosynthetic pathways are currently supported' - write(fates_log(),*) ' C4 plants have c3psn = 0' - write(fates_log(),*) ' C3 plants have c3psn = 1' + write(fates_log(),*) ' C4 plants have c3psn = ',c3_path_index + write(fates_log(),*) ' C3 plants have c3psn = ',c4_path_index write(fates_log(),*) ' PFT#: ',ipft - write(fates_log(),*) ' c3psn(pft): ',EDPftvarcon_inst%c3psn(ipft) + write(fates_log(),*) ' c3psn(pft): ',lb_params%c3psn(ipft) write(fates_log(),*) ' Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) @@ -2273,7 +2116,7 @@ subroutine FatesCheckParams(is_master) !------------------------------------------------------------------------------------ do ipft = 1,npft - r_0 = EDPftvarcon_inst%maintresp_leaf_atkin2017_baserate(ipft) + r_0 = lb_params%maintresp_leaf_atkin2017_baserate(ipft) lnc_top = prt_params%nitr_stoich_p1(ipft, prt_params%organ_param_id(leaf_organ)) @@ -2291,6 +2134,30 @@ subroutine FatesCheckParams(is_master) write(fates_log(),*) 'FatesPlantRespPhotosynthMod' end do ! ipft + + ! Check to make sure the btran limitation models are within expected ranges + + do ipft = 1,npft + + if( lb_params%stomatal_btran_model(ipft) < 0 .or. & + lb_params%stomatal_btran_model(ipft) > btran_on_gs_gs02 ) then + + write(fates_log(),*) 'PFT ', ipft + write(fates_log(),*) 'Undefined fates_leaf_stomatal_btran_model = ',lb_params%stomatal_btran_model + write(fates_log(),*) 'See biogeophys/LeafbiophysicsMod.F90 btran_on_gs_* for model types' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if( lb_params%agross_btran_model(ipft) < 0 .or. & + lb_params%agross_btran_model(ipft) > btran_on_ag_vcmax_jmax ) then + write(fates_log(),*) 'PFT ', ipft + write(fates_log(),*) 'Undefined fates_leaf_agross_btran_model = ',lb_params%agross_btran_model + write(fates_log(),*) 'See biogeophys/LeafbiophysicsMod.F90 btran_on_ag_* for model types' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + end do + diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 4aec41404d..8a1b94962b 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -135,10 +135,6 @@ module FatesConstantsMod ! based on average age of global ! secondary 1900s land in hurtt-2011 - ! integer labels for specifying harvest units - integer, parameter, public :: photosynth_acclim_model_none = 1 - integer, parameter, public :: photosynth_acclim_model_kumarathunge_etal_2019 = 2 - ! integer labels for specifying harvest units integer, parameter, public :: hlm_harvest_area_fraction = 1 ! Code for harvesting by area integer, parameter, public :: hlm_harvest_carbon = 2 ! Code for harvesting based on carbon extracted. @@ -317,6 +313,9 @@ module FatesConstantsMod ! Approximate molar mass of water vapor to dry air (-) real(fates_r8), parameter, public :: molar_mass_ratio_vapdry= 0.622_fates_r8 + ! unit conversion of W/m2 to umol photons/m^2/s + real(fates_r8), parameter, public :: wm2_to_umolm2s = 4.6_fates_r8 + ! Gravity constant on earth [m/s] real(fates_r8), parameter, public :: grav_earth = 9.8_fates_r8 @@ -326,10 +325,15 @@ module FatesConstantsMod ! Pascals to megapascals real(fates_r8), parameter, public :: mpa_per_pa = 1.e-6_fates_r8 + ! Pascals to kilopascals + real(fates_r8), parameter, public :: kpa_per_pa = 1.e-3_fates_r8 + ! Conversion: megapascals per mm H2O suction real(fates_r8), parameter, public :: mpa_per_mm_suction = dens_fresh_liquid_water * & grav_earth * 1.0E-9_fates_r8 + + ! For numerical inquiry real(fates_r8), parameter, public :: fates_huge = huge(g_per_kg) @@ -346,17 +350,7 @@ module FatesConstantsMod real(fates_r8), parameter, public :: pi_const = 3.14159265359_fates_r8 real(fates_r8), parameter, public :: rad_per_deg = pi_const/180.0_fates_r8 - ! Rdark constants from Atkin et al., 2017 https://doi.org/10.1007/978-3-319-68703-2_6 - ! and Heskel et al., 2016 https://doi.org/10.1073/pnas.1520282113 - real(fates_r8), parameter, public :: lmr_b = 0.1012_fates_r8 ! (degrees C**-1) - - real(fates_r8), parameter, public :: lmr_c = -0.0005_fates_r8 ! (degrees C**-2) - - real(fates_r8), parameter, public :: lmr_TrefC = 25._fates_r8 ! (degrees C) - - real(fates_r8), parameter, public :: lmr_r_1 = 0.2061_fates_r8 ! (umol CO2/m**2/s / (gN/(m2 leaf))) - real(fates_r8), parameter, public :: lmr_r_2 = -0.0402_fates_r8 ! (umol CO2/m**2/s/degree C) ! some integers related to termination mortality integer, parameter, public :: n_term_mort_types = 3 diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index b01cb442d1..6bb42be1b0 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -71,6 +71,7 @@ module FatesInterfaceMod use EDParamsMod , only : FatesRegisterParams, FatesReceiveParams use SFParamsMod , only : SpitFireRegisterParams, SpitFireReceiveParams use PRTInitParamsFATESMod , only : PRTRegisterParams, PRTReceiveParams + use FatesLeafBiophysParamsMod , only : LeafBiophysRegisterParams, LeafBiophysReceiveParams,LeafBiophysReportParams use FatesSynchronizedParamsMod, only : FatesSynchronizedParamsInst use EDParamsMod , only : p_uptake_mode use EDParamsMod , only : n_uptake_mode @@ -2085,6 +2086,7 @@ subroutine FatesReportParameters(masterproc) call FatesReportPFTParams(masterproc) call FatesReportParams(masterproc) + call LeafBiophysReportParams(masterproc) call PRTDerivedParams() ! Update PARTEH derived constants call FatesCheckParams(masterproc) ! Check general fates parameters call PRTCheckParams(masterproc) ! Check PARTEH parameters @@ -2524,6 +2526,7 @@ subroutine FatesReadParameters(param_reader) call FatesRegisterParams(fates_params) !EDParamsMod, only operates on fates_params class call SpitFireRegisterParams(fates_params) !SpitFire Mod, only operates of fates_params class call PRTRegisterParams(fates_params) ! PRT mod, only operates on fates_params class + call LeafBiophysRegisterParams(fates_params) call FatesSynchronizedParamsInst%RegisterParams(fates_params) !Synchronized params class in Synchronized params mod, only operates on fates_params class call param_reader%Read(fates_params) @@ -2531,6 +2534,7 @@ subroutine FatesReadParameters(param_reader) call FatesReceiveParams(fates_params) call SpitFireReceiveParams(fates_params) call PRTReceiveParams(fates_params) + call LeafBiophysReceiveParams(fates_params) call FatesSynchronizedParamsInst%ReceiveParams(fates_params) call fates_params%Destroy() diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 4b945741fd..e54ffe1b23 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -851,7 +851,7 @@ subroutine set_inventory_cohort_type1(csite,bc_in,css_file_unit,npatches, & read(css_file_unit,fmt=*,iostat=ios) c_time, p_name, c_dbh, & c_height, c_pft, c_nplant end if - + if( debug_inv) then write(*,fmt=wr_fmt) & c_time, p_name, c_dbh, c_height, c_pft, c_nplant diff --git a/main/FatesParametersInterface.F90 b/main/FatesParametersInterface.F90 index 0559ec3eb6..feebf503b7 100644 --- a/main/FatesParametersInterface.F90 +++ b/main/FatesParametersInterface.F90 @@ -9,10 +9,10 @@ module FatesParametersInterface implicit none private ! Modules are private by default - integer, parameter, public :: max_params = 250 + integer, parameter, public :: max_params = 500 integer, parameter, public :: max_dimensions = 2 integer, parameter, public :: max_used_dimensions = 25 - integer, parameter, public :: param_string_length = 40 + integer, parameter, public :: param_string_length = 80 ! NOTE(bja, 2017-02) these are the values returned from netcdf after ! inquiring about the number of dimensions integer, parameter, public :: dimension_shape_scalar = 0 diff --git a/main/FatesUtilsMod.F90 b/main/FatesUtilsMod.F90 index be955d5da0..720fc2ce9c 100644 --- a/main/FatesUtilsMod.F90 +++ b/main/FatesUtilsMod.F90 @@ -20,12 +20,27 @@ module FatesUtilsMod public :: FindIndex public :: QuadraticRootsNSWC public :: QuadraticRootsSridharachary + public :: ArrayNint character(len=*), parameter, private :: sourcefile = & __FILE__ contains + subroutine ArrayNint(realarr,intarr) + + real(r8),intent(in) :: realarr(:) + integer,intent(out) :: intarr(:) + integer :: i + + do i = 1,size(realarr,dim=1) + intarr(i) = nint(realarr(i)) + end do + + return + end subroutine ArrayNint + + ! ============================================================================================= function check_hlm_list(hlms,hlm_name) result(astatus) @@ -186,7 +201,7 @@ function FindIndex(input_string_array,string_to_match) result(array_index) end function FindIndex - subroutine QuadraticRootsNSWC(a,b,c,root1,root2) + subroutine QuadraticRootsNSWC(a,b,c,root1,root2,err) ! This code is based off of routines from the NSWC Mathematics Subroutine Library ! From the NSWC README (https://github.com/jacobwilliams/nswc) @@ -202,9 +217,10 @@ subroutine QuadraticRootsNSWC(a,b,c,root1,root2) real(r8),intent(in) :: a , b , c !! coefficients real(r8),intent(out) :: root1 ! sr !! real part of first root real(r8),intent(out) :: root2 ! lr !! real part of second root - + logical,intent(out) :: err real(r8) :: b1, d, e + err = .false. if ( abs(a)=0.0_r8 ) d = -d @@ -239,15 +256,17 @@ subroutine QuadraticRootsNSWC(a,b,c,root1,root2) end subroutine QuadraticRootsNSWC - subroutine QuadraticRootsSridharachary(a,b,c,root1,root2) + subroutine QuadraticRootsSridharachary(a,b,c,root1,root2,err) real(r8),intent(in) :: a , b , c !! coefficients real(r8),intent(out) :: root1 ! sr !! real part of first root real(r8),intent(out) :: root2 ! lr !! real part of second root + logical,intent(out) :: err real(r8) :: d ! discriminant real(r8) :: das ! sqrt(abs(d)) + err = .false. ! If a is 0, then equation is not quadratic, but linear if (abs(a) < nearzero ) then root2 = 0.0_r8 @@ -271,7 +290,8 @@ subroutine QuadraticRootsSridharachary(a,b,c,root1,root2) else write (fates_log(),*)'error, imaginary roots detected in quadratic solve' - call endrun(msg=errMsg(sourcefile, __LINE__)) + err = .true. + !call endrun(msg=errMsg(sourcefile, __LINE__)) end if diff --git a/parameter_files/archive/apiX_pr1262.xml b/parameter_files/archive/apiX_pr1262.xml new file mode 100644 index 0000000000..48e1eaeed5 --- /dev/null +++ b/parameter_files/archive/apiX_pr1262.xml @@ -0,0 +1,50 @@ + + + + + + + + + + + + + + + + + + + archive/base.cdl + fates_params_default.cdl + 1,2,3,4,5,6,7,8,9,10,11,12,13,14 + + + fates_leaf_theta_cj_c3 + + + fates_leaf_theta_cj_c4 + + + fates_leaf_stomatal_btran_model + fates_pft + index + model switch for how btran affects conductance. See LeafBiophysicsMod.F90, integer constants: btran_on_ + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 + + + fates_leaf_agross_btran_model + fates_pft + index + model switch for how gross assimilation affects conductance. See LeafBiophysicsMod.F90, integer constants: btran_on_ + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 + + + diff --git a/parameter_files/patch_default_bciopt224.xml b/parameter_files/patch_default_bciopt224.xml index 73ab594ff8..a86156122e 100644 --- a/parameter_files/patch_default_bciopt224.xml +++ b/parameter_files/patch_default_bciopt224.xml @@ -2,7 +2,7 @@ This parameter dataset was created by Ryan Knox rgknox@lbl.gov. Please contact if using in published work. The calibration uses the following datasets: [1] Ely et al. 2019. Leaf mass area, Panama. NGEE-Tropics data collection.http://dx.doi.org/10.15486/ngt/1411973 and [2] Condit et al. 2019. Complete data from the Barro Colorado 50-ha plot. https://doi.org/10.15146/5xcp-0d46. [3] Koven et al. 2019. Benchmarking and parameter sensitivity of physiological and vegetation dynamics using the functionally assembled terrestrial ecosystem simulator. Biogeosciences. The ECA nutrient aquisition parmeters are unconstrained, the file output naming convention vmn6phi is shorthand for vmax for nitrogen uptake is order e-6 and for phosphorus is excessively high. These parameters were calibrated with the special fates modification in main/EDTypesMod.F90: nclmax = 3 fates_params_default.cdl - fates_params_opt224_0424_api34.cdl + fates_params_opt224_1024_api35.cdl 1 diff --git a/radiation/TwoStreamMLPEMod.F90 b/radiation/TwoStreamMLPEMod.F90 index 25b8efdc1b..f4cdb05253 100644 --- a/radiation/TwoStreamMLPEMod.F90 +++ b/radiation/TwoStreamMLPEMod.F90 @@ -27,7 +27,6 @@ Module TwoStreamMLPEMod use shr_log_mod , only: errMsg => shr_log_errMsg use shr_sys_mod , only: shr_sys_abort use FatesConstantsMod, only : r8 => fates_r8 - use shr_infnan_mod, only : shr_infnan_isnan implicit none private @@ -372,8 +371,8 @@ function GetRdDn(this,ican,icol,ib,vai) result(r_diff_dn) if(debug)then ! if(isnan(r_diff_dn))then !RGK: NVHPC HAS A BUG IN THIS INTRINSIC (01-2024) - ! if(r_diff_dn /= r_diff_dn) then - if(shr_infnan_isnan(r_diff_dn)) then + if(r_diff_dn /= r_diff_dn) then + !if(shr_infnan_isnan(r_diff_dn)) then write(log_unit,*)"GETRDN" write(log_unit,*)scelg%Kb write(log_unit,*)scelb%a @@ -832,8 +831,8 @@ subroutine CanopyPrep(this,frac_snow) if(debug)then !if(isnan(scelb%betad))then !RGK: NVHPC HAS A BUG IN THIS INTRINSIC (01-2024) - !if(scelb%betad /= scelb%betad) then - if(shr_infnan_isnan(scelb%betad))then + if(scelb%betad /= scelb%betad) then + !if(shr_infnan_isnan(scelb%betad))then write(log_unit,*)"nans in canopy prep" write(log_unit,*) ib,ican,icol,ft write(log_unit,*) scelb%betad,scelb%om,lai,sai