diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 2792d02c7f..af20e8670b 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -40,8 +40,9 @@ module FatesPlantHydraulicsMod use FatesConstantsMod, only : cm2_per_m2 use FatesConstantsMod, only : g_per_kg - use EDParamsMod , only : hydr_kmax_rsurf - + use EDParamsMod , only : hydr_kmax_rsurf1 + ! use EDParamsMod , only : hydr_kmax_rsurf2 + use EDTypesMod , only : ed_site_type use EDTypesMod , only : ed_patch_type use EDTypesMod , only : ed_cohort_type @@ -141,6 +142,7 @@ module FatesPlantHydraulicsMod public :: CopyCohortHydraulics public :: FuseCohortHydraulics public :: updateSizeDepTreeHydProps + public :: updateWaterDepTreeHydProps public :: updateSizeDepTreeHydStates public :: initTreeHydStates public :: updateSizeDepRhizHydProps @@ -150,6 +152,7 @@ module FatesPlantHydraulicsMod public :: SavePreviousRhizVolumes public :: UpdateTreeHydrNodes public :: UpdateTreeHydrLenVolCond + public :: UpdateWaterDepTreeHydrCond public :: ConstrainRecruitNumber !------------------------------------------------------------------------------ @@ -325,7 +328,8 @@ subroutine initTreeHydStates(site_p, cc_p, bc_in) ! watsat(c,j), watres(c,j), alpha_VG(c,j), n_VG(c,j), m_VG(c,j), l_VG(c,j), & ! smp) !ccohort_hydr%psi_aroot(j) = smp - ccohort_hydr%psi_aroot(j) = csite%si_hydr%psisoi_liq_innershell(j) + !ccohort_hydr%psi_aroot(j) = csite%si_hydr%psisoi_liq_innershell(j) + ccohort_hydr%psi_aroot(j) = -0.2_r8 !do not assume the equalibrium between soil and root call th_from_psi(ft, 4, ccohort_hydr%psi_aroot(j), ccohort_hydr%th_aroot(j), csite%si_hydr, bc_in ) end do @@ -346,6 +350,7 @@ subroutine initTreeHydStates(site_p, cc_p, bc_in) !hydrostatic equilibrium with the water potential immediately below dz = ccohort_hydr%z_node_ag(n_hypool_ag) - ccohort_hydr%z_node_troot(1) ccohort_hydr%psi_ag(n_hypool_ag) = ccohort_hydr%psi_troot(1) - 1.e-6_r8*denh2o*grav*dz + if (ccohort_hydr%psi_ag(n_hypool_ag)>0.0_r8) ccohort_hydr%psi_ag(n_hypool_ag) = -0.01_r8 call th_from_psi(ft, 2, ccohort_hydr%psi_ag(n_hypool_ag), ccohort_hydr%th_ag(n_hypool_ag), csite%si_hydr, bc_in) do k=n_hypool_ag-1, 1, -1 dz = ccohort_hydr%z_node_ag(k) - ccohort_hydr%z_node_ag(k+1) @@ -543,9 +548,44 @@ subroutine updateSizeDepTreeHydProps(currentSite,ccohort,bc_in) ! volumes, and UpdateTreeHydrNodes is called prior to this. call UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) + + end subroutine updateSizeDepTreeHydProps + + ! ===================================================================================== + + subroutine updateWaterDepTreeHydProps(currentSite,ccohort,bc_in) + + + ! DESCRIPTION: Updates absorbing root length (total and its vertical distribution) + ! as well as the consequential change in the size of the 'representative' rhizosphere + ! shell radii, volumes, and compartment volumes of plant tissues + + ! !USES: + use shr_sys_mod , only : shr_sys_abort + + ! ARGUMENTS: + type(ed_site_type) , intent(in) :: currentSite ! Site stuff + type(ed_cohort_type) , intent(inout) :: ccohort ! current cohort pointer + type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions + + ! Locals + integer :: nlevsoi_hyd ! Number of total soil layers + type(ed_cohort_hydr_type), pointer :: ccohort_hydr + integer :: ft + + nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd + ccohort_hydr => ccohort%co_hydr + ft = ccohort%pft + ! This updates plant compartment volumes, lengths and + ! maximum conductances. Make sure for already + ! initialized vegetation, that SavePreviousCompartment + ! volumes, and UpdateTreeHydrNodes is called prior to this. - end subroutine updateSizeDepTreeHydProps + call UpdateWaterDepTreeHydrCond(currentSite,ccohort,nlevsoi_hyd,bc_in) + + + end subroutine updateWaterDepTreeHydProps ! ===================================================================================== @@ -616,7 +656,6 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) ! hydraulic conductance [kg s-1 MPa-1] real(r8) :: kmax_tot ! total tree (leaf to root tip) ! hydraulic conductance [kg s-1 MPa-1] - real(r8),parameter :: taper_exponent = 1._r8/3._r8 ! Savage et al. (2010) xylem taper exponent [-] ccohort_hydr => ccohort%co_hydr @@ -810,10 +849,84 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) ccohort_hydr%kmax_treebg_layer(j) = rootfr*ccohort_hydr%kmax_treebg_tot end do end if + end if !check for bleaf end subroutine UpdateTreeHydrLenVolCond + + + + !===================================================================================== + subroutine UpdateWaterDepTreeHydrCond(currentSite,ccohort,nlevsoi_hyd,bc_in) + + ! ----------------------------------------------------------------------------------- + ! This subroutine calculates update the conductivity for the soil-root interface, + ! depending on the plant water uptake/loss. + ! we assume that the conductivitity for water uptake is larger than + ! water loss due to composite regulation of resistance the roots + ! hydraulic vs osmostic with and without transpiration + ! Steudle, E. Water uptake by roots: effects of water deficit. + ! J Exp Bot 51, 1531-1542, doi:DOI 10.1093/jexbot/51.350.1531 (2000). + ! ----------------------------------------------------------------------------------- + + ! Arguments + type(ed_site_type) , intent(in) :: currentSite ! Site target + type(ed_cohort_type),intent(inout) :: ccohort ! cohort target + integer,intent(in) :: nlevsoi_hyd ! number of soil hydro layers + type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions + + type(ed_cohort_hydr_type),pointer :: ccohort_hydr ! Plant hydraulics structure + type(ed_site_hydr_type),pointer :: csite_hydr + + integer :: j,k + real(r8) :: hksat_s ! hksat converted to units of 10^6sec + real(r8) :: kmax_root_surf_total ! maximum conducitivity for total root surface(kg water/Mpa/s) + real(r8) :: kmax_soil_total ! maximum conducitivity for from root surface to soil shell(kg water/Mpa/s) + ! which is equiv to [kg m-1 s-1 MPa-1] + real(r8) :: kmax_root_surf ! maximum conducitivity for unit root surface (kg water/m2 root area/Mpa/s) + + ! (RGK 4-2019) THE FOLLOWING SHOULD BE ADDED TO THE PARAMETER FILE IN NEXT GROUPING + real(r8), parameter :: hydr_kmax_rsurf2 = 0.0001_r8 ! kg water/m2 root area/Mpa/s + + + ccohort_hydr => ccohort%co_hydr + csite_hydr => currentSite%si_hydr + k = 1 !only for the first soil shell + do j=1, nlevsoi_hyd + + hksat_s = bc_in%hksat_sisl(j) * 1.e-3_r8 * 1/grav * 1.e6_r8 + if(ccohort_hydr%psi_aroot(j) cPatch%older enddo !patch - kmax_root_surf = hydr_kmax_rsurf + csite_hydr%l_aroot_1D = sum( csite_hydr%l_aroot_layer(:)) ! update outer radii of column-level rhizosphere shells (same across patches and cohorts) @@ -1646,6 +1756,9 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in) enddo call shellGeom( csite_hydr%l_aroot_1D, csite_hydr%rs1(1), AREA, sum(bc_in%dz_sisl(1:nlevsoi_hyd)), & csite_hydr%r_out_shell_1D(:), csite_hydr%r_node_shell_1D(:), csite_hydr%v_shell_1D(:)) + + !update the conductitivity for first soil shell is done at subroutine UpdateWaterDepTreeHydrCond + !which is dependant on whether it is water uptake or loss for every 30 minutes do j = 1,csite_hydr%nlevsoi_hyd @@ -1654,50 +1767,8 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in) ! proceed only if the total absorbing root length (site-level) has changed in this layer if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then - do k = 1,nshell - if(k == 1) then - kmax_root_surf_total = kmax_root_surf*2._r8*pi_const *csite_hydr%rs1(j)* & - csite_hydr%l_aroot_layer(j) - if(csite_hydr%r_node_shell(j,k) <= csite_hydr%rs1(j)) then - !csite_hydr%kmax_upper_shell(j,k) = large_kmax_bound - !csite_hydr%kmax_bound_shell(j,k) = large_kmax_bound - !csite_hydr%kmax_lower_shell(j,k) = large_kmax_bound - csite_hydr%kmax_upper_shell(j,k) = kmax_root_surf_total - csite_hydr%kmax_bound_shell(j,k) = kmax_root_surf_total - csite_hydr%kmax_lower_shell(j,k) = kmax_root_surf_total - - else - - kmax_soil_total = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / & - log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s - - !csite_hydr%kmax_upper_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / & - ! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s - !csite_hydr%kmax_bound_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / & - ! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s - !csite_hydr%kmax_lower_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / & - ! log(csite_hydr%r_node_shell(j,k)/csite_hydr%rs1(j))*hksat_s - - csite_hydr%kmax_upper_shell(j,k) = (1._r8/kmax_root_surf_total + & - 1._r8/kmax_soil_total)**(-1._r8) - csite_hydr%kmax_bound_shell(j,k) = (1._r8/kmax_root_surf_total + & - 1._r8/kmax_soil_total)**(-1._r8) - csite_hydr%kmax_lower_shell(j,k) = (1._r8/kmax_root_surf_total + & - 1._r8/kmax_soil_total)**(-1._r8) - end if - if(j == 1) then - if(csite_hydr%r_node_shell(j,k) <= csite_hydr%rs1(j)) then - csite_hydr%kmax_upper_shell_1D(k) = csite_hydr%kmax_upper_shell(1,k) - csite_hydr%kmax_bound_shell_1D(k) = csite_hydr%kmax_bound_shell(1,k) - csite_hydr%kmax_lower_shell_1D(k) = csite_hydr%kmax_lower_shell(1,k) - else - csite_hydr%kmax_upper_shell_1D(k) = csite_hydr%kmax_upper_shell(1,k) - csite_hydr%kmax_bound_shell_1D(k) = csite_hydr%kmax_bound_shell(1,k) - csite_hydr%kmax_lower_shell_1D(k) = csite_hydr%kmax_lower_shell(1,k) - end if - end if - else - csite_hydr%kmax_upper_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / & + do k = 2,nshell + csite_hydr%kmax_upper_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / & log(csite_hydr%r_node_shell(j,k)/csite_hydr%r_out_shell(j,k-1))*hksat_s csite_hydr%kmax_bound_shell(j,k) = 2._r8*pi_const*csite_hydr%l_aroot_layer(j) / & log(csite_hydr%r_node_shell(j,k)/csite_hydr%r_node_shell(j,k-1))*hksat_s @@ -1711,16 +1782,14 @@ subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in) csite_hydr%kmax_lower_shell_1D(k) = 2._r8*pi_const*csite_hydr%l_aroot_1D / & log(csite_hydr%r_out_shell_1D( k)/csite_hydr%r_node_shell_1D(k ))*hksat_s end if - end if enddo ! loop over rhizosphere shells end if !has l_aroot_layer changed? enddo ! loop over soil layers - return end subroutine UpdateSizeDepRhizVolLenCon - + ! ===================================================================================== @@ -2445,7 +2514,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ! [mm H2O/cohort/s] = [mm H2O / patch / s] / [cohort/patch] !! qflx_tran_veg_patch_coh = qflx_trans_patch_vol * qflx_rel_tran_coh - + call updateWaterDepTreeHydProps(sites(s),ccohort,bc_in(s)) if(site_hydr%nlevsoi_hyd > 1) then ! BUCKET APPROXIMATION OF THE SOIL-ROOT HYDRAULIC GRADIENT (weighted average across layers) @@ -2516,6 +2585,9 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) kmax_lower( 1 : n_hypool_ag ) = ccohort_hydr%kmax_lower(:) kmax_upper(( n_hypool_ag+1) ) = ccohort_hydr%kmax_upper_troot if(site_hydr%nlevsoi_hyd == 1) then + site_hydr%kmax_upper_shell_1D(1) = ccohort_hydr%kmax_innershell(1) + site_hydr%kmax_lower_shell_1D(1) = ccohort_hydr%kmax_innershell(1) + site_hydr%kmax_bound_shell_1D(1) = ccohort_hydr%kmax_innershell(1) !! estimate troot-aroot and aroot-radial components as a residual: !! 25% each of total (surface of aroots to leaves) resistance kmax_bound(( n_hypool_ag+1):(n_hypool_ag+2 )) = 2._r8 * ccohort_hydr%kmax_treebg_tot @@ -2645,11 +2717,14 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) v_node_1l((n_hypool_ag+n_hypool_troot+1) ) = ccohort_hydr%v_aroot_layer(j) v_node_1l((n_hypool_tot-nshell+1):(n_hypool_tot)) = site_hydr%v_shell(j,:) * & ccohort_hydr%l_aroot_layer(j)/bc_in(s)%dz_sisl(j) + site_hydr%kmax_bound_shell(j,1)=ccohort_hydr%kmax_innershell(j) + site_hydr%kmax_upper_shell(j,1)=ccohort_hydr%kmax_innershell(j) + site_hydr%kmax_lower_shell(j,1)=ccohort_hydr%kmax_innershell(j) kmax_bound_1l(:) = 0._r8 kmax_bound_shell_1l(:) = site_hydr%kmax_bound_shell(j,:) * & ccohort_hydr%l_aroot_layer(j) / site_hydr%l_aroot_layer(j) - + ! transporting-to-absorbing root conductance: factor of 2 means one-half of the total ! belowground resistance in layer j kmax_bound_1l((n_hypool_ag+1)) = 2._r8 * ccohort_hydr%kmax_treebg_layer(j) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 09840da049..0baf150f4e 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -1,3 +1,4 @@ + module FATESPlantRespPhotosynthMod !------------------------------------------------------------------------------------- @@ -44,6 +45,7 @@ module FATESPlantRespPhotosynthMod use PRTGenericMod, only : store_organ use PRTGenericMod, only : repro_organ use PRTGenericMod, only : struct_organ + use EDParamsMod, only : ED_val_bbopt_c3, ED_val_bbopt_c4, ED_val_base_mr_20 ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -58,7 +60,7 @@ module FATESPlantRespPhotosynthMod !------------------------------------------------------------------------------------- ! maximum stomatal resistance [s/m] (used across several procedures) - real(r8),parameter :: rsmax0 = 2.e4_r8 + real(r8),parameter :: rsmax0 = 2.e8_r8 logical :: debug = .false. @@ -95,7 +97,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) use FatesConstantsMod, only : rgas => rgas_J_K_kmol use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use FatesParameterDerivedMod, only : param_derived - use EDParamsMod, only : ED_val_bbopt_c3, ED_val_bbopt_c4, ED_val_base_mr_20 + use FatesAllometryMod, only : bleaf use FatesAllometryMod, only : storage_fraction_of_target use FatesAllometryMod, only : set_root_fraction @@ -393,10 +395,10 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) (hlm_use_planthydro.eq.itrue) .or. & (nleafage > 1) .or. & (hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then - - if (hlm_use_planthydro.eq.itrue ) then - - bbb = max (bbbopt(nint(c3psn(ft)))*currentCohort%co_hydr%btran(1), 1._r8) + + if (hlm_use_planthydro.eq.itrue ) then + + bbb = max( cf/rsmax0, bbbopt(nint(c3psn(ft)))*currentCohort%co_hydr%btran(1) ) btran_eff = currentCohort%co_hydr%btran(1) ! dinc_ed is the total vegetation area index of each "leaf" layer @@ -413,7 +415,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) else - bbb = max (bbbopt(nint(c3psn(ft)))*currentPatch%btran_ft(ft), 1._r8) + bbb = max( cf/rsmax0, bbbopt(nint(c3psn(ft)))*currentPatch%btran_ft(ft) ) btran_eff = currentPatch%btran_ft(ft) ! For consistency sake, we use total LAI here, and not exposed ! if the plant is under-snow, it will be effectively dormant for @@ -831,7 +833,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in psn_out, & ! out rstoma_out, & ! out anet_av_out, & ! out - c13disc_z) ! out + c13disc_z) ! out ! ------------------------------------------------------------------------------------ ! This subroutine calculates photosynthesis and stomatal conductance within each leaf @@ -843,6 +845,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! ------------------------------------------------------------------------------------ use EDPftvarcon , only : EDPftvarcon_inst + ! Arguments ! ------------------------------------------------------------------------------------ @@ -916,11 +919,20 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in 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), dimension(0:1) :: bbbopt ! Cuticular conductance at full water potential (umol H2O /m2/s) + + + ! Parameters ! ------------------------------------------------------------------------ ! Fraction of light absorbed by non-photosynthetic pigments real(r8),parameter :: fnps = 0.15_r8 + ! For plants with no leaves, a miniscule amount of conductance + ! can happen through the stems, at a partial rate of cuticular conductance + real(r8),parameter :: stem_cuticle_loss_frac = 0.1_r8 + ! empirical curvature parameter for electron transport rate real(r8),parameter :: theta_psii = 0.7_r8 @@ -939,12 +951,16 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! empirical curvature parameter for ap photosynthesis co-limitation real(r8),parameter :: theta_ip = 0.999_r8 - + associate( bb_slope => EDPftvarcon_inst%BB_slope) ! slope of BB relationship + ! photosynthetic pathway: 0. = c4, 1. = c3 c3c4_path_index = nint(EDPftvarcon_inst%c3psn(ft)) + bbbopt(0) = ED_val_bbopt_c4 + bbbopt(1) = ED_val_bbopt_c3 + if (c3c4_path_index == 1) then init_co2_inter_c = init_a2l_co2_c3 * can_co2_ppress else @@ -955,22 +971,22 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! ---------------------------------------------------------------------------------- if ( parsun_lsl <= 0._r8 ) then ! night time - + anet_av_out = -lmr psn_out = 0._r8 - rstoma_out = min(rsmax0, 1._r8/bbb * cf) - c13disc_z = 0.0_r8 !carbon 13 discrimination in night time carbon flux, note value of 1.0 is used in CLM + + ! The cuticular conductance already factored in maximum resistance as a bound + ! no need to re-bound it + + rstoma_out = cf/bbb + + c13disc_z = 0.0_r8 !carbon 13 discrimination in night time carbon flux, note value of 1.0 is used in CLM else ! day time (a little bit more complicated ...) -! if ( debug ) write(fates_log(),*) 'EDphot 594 ',laisun_lsl -! if ( debug ) write(fates_log(),*) 'EDphot 595 ',laisha_lsl - - !is there leaf area? - (NV can be larger than 0 with only stem area if deciduous) - if ( laisun_lsl + laisha_lsl > 0._r8 ) then + !is there leaf area? - (NV can be larger than 0 with only stem area if deciduous) + if ( laisun_lsl + laisha_lsl > 0._r8 ) then -! if ( debug ) write(fates_log(),*) '600 in laisun, laisha loop ' - !Loop aroun shaded and unshaded leaves psn_out = 0._r8 ! psn is accumulated across sun and shaded leaves. rstoma_out = 0._r8 ! 1/rs is accumulated across sun and shaded leaves. @@ -1120,19 +1136,16 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! Convert gs_mol (umol /m**2/s) to gs (m/s) and then to rs (s/m) gs = gs_mol / cf - ! estimate carbon 13 discrimination in leaf level carbon flux Liang WEI and Hang ZHOU 2018, based on + ! estimate carbon 13 discrimination in leaf level carbon + ! flux Liang WEI and Hang ZHOU 2018, based on ! Ubierna and Farquhar, 2014 doi:10.1111/pce.12346, using the simplified model: ! $\Delta ^{13} C = \alpha_s + (b - \alpha_s) \cdot \frac{C_i}{C_a}$ ! just hard code b and \alpha_s for now, might move to parameter set in future ! b = 27.0 alpha_s = 4.4 ! TODO, not considering C4 or CAM right now, may need to address this - ! note co2_inter_c is intracelluar CO2, not intercelluar - c13disc_z = 4.4_r8 + (27.0_r8 - 4.4_r8) * min (can_co2_ppress, max (co2_inter_c, 0._r8)) / can_co2_ppress - - -! if ( debug ) write(fates_log(),*) 'EDPhoto 737 ', psn_out -! if ( debug ) write(fates_log(),*) 'EDPhoto 738 ', agross -! if ( debug ) write(fates_log(),*) 'EDPhoto 739 ', f_sun_lsl + ! note co2_inter_c is intracelluar CO2, not intercelluar + c13disc_z = 4.4_r8 + (27.0_r8 - 4.4_r8) * & + min (can_co2_ppress, max (co2_inter_c, 0._r8)) / can_co2_ppress ! Accumulate total photosynthesis umol/m2 ground/s-1. ! weight per unit sun and sha leaves. @@ -1146,10 +1159,6 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in gstoma = gstoma + & 1._r8/(min(1._r8/gs, rsmax0)) * (1.0_r8-f_sun_lsl) end if - -! if ( debug ) write(fates_log(),*) 'EDPhoto 758 ', psn_out -! if ( debug ) write(fates_log(),*) 'EDPhoto 759 ', agross -! if ( debug ) write(fates_log(),*) 'EDPhoto 760 ', f_sun_lsl ! Make sure iterative solution is correct if (gs_mol < 0._r8) then @@ -1171,19 +1180,25 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! This is the stomatal resistance of the leaf layer rstoma_out = 1._r8/gstoma - + else - !No leaf area. This layer is present only because of stems. + + ! No leaf area. This layer is present only because of stems. + ! Net assimilation is zero, not negative because there are + ! no leaves to even respire ! (leaves are off, or have reduced to 0) - psn_out = 0._r8 - rstoma_out = min(rsmax0, 1._r8/bbb * cf) - c13disc_z = 0.0_r8 + psn_out = 0._r8 + anet_av_out = 0._r8 + rstoma_out = min(rsmax0, cf/(stem_cuticle_loss_frac*bbbopt(c3c4_path_index))) + c13disc_z = 0.0_r8 - end if !is there leaf area? + end if !is there leaf area? end if ! night or day + + end associate return end subroutine LeafLayerPhotosynthesis diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 15103cdd29..badc5d457c 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -79,8 +79,11 @@ module EDParamsMod ! Hydraulics Control Parameters (ONLY RELEVANT WHEN USE_FATES_HYDR = TRUE) ! ---------------------------------------------------------------------------------------------- - real(r8),protected :: hydr_kmax_rsurf ! maximum conducitivity for unit root surface (kg water/m2 root area/Mpa/s) - character(len=param_string_length),parameter :: hydr_name_kmax_rsurf = "fates_hydr_kmax_rsurf" + real(r8),protected :: hydr_kmax_rsurf1 ! maximum conducitivity for unit root surface for water uptake(kg water/m2 root area/Mpa/s) +! real(r8),protected :: hydr_kmax_rsurf2 ! maximum conducitivity for unit root surface for water loss (kg water/m2 root area/Mpa/s) + + character(len=param_string_length),parameter :: hydr_name_kmax_rsurf1 = "fates_hydr_kmax_rsurf" +! character(len=param_string_length),parameter :: hydr_name_kmax_rsurf2 = "fates_hydr_kmax_rsurf2" real(r8),protected :: hydr_psi0 ! sapwood water potential at saturation (MPa) character(len=param_string_length),parameter :: hydr_name_psi0 = "fates_hydr_psi0" @@ -160,7 +163,8 @@ subroutine FatesParamsInit() ED_val_patch_fusion_tol = nan ED_val_canopy_closure_thresh = nan - hydr_kmax_rsurf = nan + hydr_kmax_rsurf1 = nan +! hydr_kmax_rsurf2 = nan hydr_psi0 = nan hydr_psicap = nan @@ -263,8 +267,11 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=ED_name_canopy_closure_thresh, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) - call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf, dimension_shape=dimension_shape_1d, & + call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf1, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) + +! call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf2, dimension_shape=dimension_shape_1d, & +! dimension_names=dim_names) call fates_params%RegisterParameter(name=hydr_name_psi0, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) @@ -387,8 +394,11 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=ED_name_canopy_closure_thresh, & data=ED_val_canopy_closure_thresh) - call fates_params%RetreiveParameter(name=hydr_name_kmax_rsurf, & - data=hydr_kmax_rsurf) + call fates_params%RetreiveParameter(name=hydr_name_kmax_rsurf1, & + data=hydr_kmax_rsurf1) + +! call fates_params%RetreiveParameter(name=hydr_name_kmax_rsurf2, & +! data=hydr_kmax_rsurf2) call fates_params%RetreiveParameter(name=hydr_name_psi0, & data=hydr_psi0) @@ -470,7 +480,8 @@ subroutine FatesReportParams(is_master) write(fates_log(),fmt0) 'ED_val_cohort_fusion_tol = ',ED_val_cohort_fusion_tol write(fates_log(),fmt0) 'ED_val_patch_fusion_tol = ',ED_val_patch_fusion_tol write(fates_log(),fmt0) 'ED_val_canopy_closure_thresh = ',ED_val_canopy_closure_thresh - write(fates_log(),fmt0) 'hydr_kmax_rsurf = ',hydr_kmax_rsurf + write(fates_log(),fmt0) 'hydr_kmax_rsurf1 = ',hydr_kmax_rsurf1 + ! write(fates_log(),fmt0) 'hydr_kmax_rsurf2 = ',hydr_kmax_rsurf2 write(fates_log(),fmt0) 'hydr_psi0 = ',hydr_psi0 write(fates_log(),fmt0) 'hydr_psicap = ',hydr_psicap write(fates_log(),fmt0) 'bgc_soil_salinity = ', bgc_soil_salinity diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90 index fad1790e03..93517abb44 100644 --- a/main/FatesHydraulicsMemMod.F90 +++ b/main/FatesHydraulicsMemMod.F90 @@ -130,6 +130,7 @@ module FatesHydraulicsMemMod real(r8),allocatable :: psisoi_liq_innershell(:) ! Matric potential of the inner rhizosphere shell (MPa) + real(r8),allocatable :: recruit_w_uptake(:) ! recruitment water uptake (kg H2o/m2/s) @@ -211,6 +212,8 @@ module FatesHydraulicsMemMod real(r8),allocatable :: v_aroot_layer_init(:) ! previous day's volume of absorbing roots by soil layer [m3] real(r8),allocatable :: v_aroot_layer(:) ! volume of absorbing roots by soil layer [m3] real(r8),allocatable :: l_aroot_layer(:) ! length of absorbing roots by soil layer [m] + + real(r8),allocatable :: kmax_innershell(:) ! Maximum hydraulic conductivity of the inner rhizosphere shell (kg s-1 MPa-1) ! BC PLANT HYDRAULICS - state variables real(r8) :: th_ag(n_hypool_ag) ! water in aboveground compartments [kgh2o/indiv] @@ -307,6 +310,7 @@ subroutine AllocateHydrCohortArrays(this,nlevsoil_hydr) allocate(this%flc_min_aroot(1:nlevsoil_hydr)) allocate(this%errh2o_growturn_aroot(1:nlevsoil_hydr)) allocate(this%errh2o_pheno_aroot(1:nlevsoil_hydr)) + allocate(this%kmax_innershell(1:nlevsoil_hydr)) return end subroutine AllocateHydrCohortArrays @@ -329,6 +333,7 @@ subroutine DeallocateHydrCohortArrays(this) deallocate(this%flc_min_aroot) deallocate(this%errh2o_growturn_aroot) deallocate(this%errh2o_pheno_aroot) + deallocate(this%kmax_innershell) return end subroutine DeallocateHydrCohortArrays diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 8d15d9e56a..1ab5311814 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -51,7 +51,7 @@ variables: fates_fire_nignitions:long_name = "number of daily ignitions (nfires = nignitions*FDI*area_scaling)" ; float fates_hydr_kmax_rsurf(fates_scalar) ; fates_hydr_kmax_rsurf:units = "kg water/m2 root area/Mpa/s" ; - fates_hydr_kmax_rsurf:long_name = "maximum conducitivity for unit root surface" ; + fates_hydr_kmax_rsurf:long_name = "maximum conducitivity for unit root surface for water uptake" ; float fates_hydr_psi0(fates_scalar) ; fates_hydr_psi0:units = "MPa" ; fates_hydr_psi0:long_name = "sapwood water potential at saturation" ;