diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index c20a91841b..e98549286e 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -31,6 +31,9 @@ module EDCohortDynamicsMod use FatesPlantHydraulicsMod, only : InitHydrCohort use FatesPlantHydraulicsMod, only : DeallocateHydrCohort use FatesPlantHydraulicsMod, only : AccumulateMortalityWaterStorage + use FatesPlantHydraulicsMod, only : UpdateTreeHydrNodes + use FatesPlantHydraulicsMod, only : UpdateTreeHydrLenVolCond + use FatesPlantHydraulicsMod, only : SavePreviousCompartmentVolumes use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index use FatesAllometryMod , only : bleaf use FatesAllometryMod , only : bfineroot @@ -98,6 +101,13 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine ! ! !DESCRIPTION: ! create new cohort + ! There are 4 places this is called + ! 1) Initializing new cohorts at the beginning of a cold-start simulation + ! 2) Initializing new recruits during dynamics + ! 3) Initializing new cohorts at the beginning of a inventory read + ! 4) Initializing new cohorts during restart + ! + ! It is assumed that in the first 3, this is called with a reasonable amount of starter information. ! ! !USES: ! @@ -125,7 +135,8 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine ! !LOCAL VARIABLES: type(ed_cohort_type), pointer :: new_cohort ! Pointer to New Cohort structure. type(ed_cohort_type), pointer :: storesmallcohort - type(ed_cohort_type), pointer :: storebigcohort + type(ed_cohort_type), pointer :: storebigcohort + integer :: nlevsoi_hyd ! number of hydraulically active soil layers integer :: tnull,snull ! are the tallest and shortest cohorts allocate !---------------------------------------------------------------------- @@ -242,12 +253,31 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine new_cohort%isnew = .true. if( hlm_use_planthydro.eq.itrue ) then - call InitHydrCohort(CurrentSite,new_cohort) - call updateSizeDepTreeHydProps(CurrentSite,new_cohort, bc_in) - call initTreeHydStates(CurrentSite,new_cohort, bc_in) + + nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd + + ! This allocates array spaces + call InitHydrCohort(currentSite,new_cohort) + + ! This calculates node heights + call UpdateTreeHydrNodes(new_cohort%co_hydr,new_cohort%pft, & + new_cohort%hite,nlevsoi_hyd,bc_in) + + ! This calculates volumes, lengths and max conductances + call UpdateTreeHydrLenVolCond(new_cohort,nlevsoi_hyd,bc_in) + + ! Since this is a newly initialized plant, we set the previous compartment-size + ! equal to the ones we just calculated. + call SavePreviousCompartmentVolumes(new_cohort%co_hydr) + + ! This comes up with starter suctions and then water contents + ! based on the soil values + call initTreeHydStates(currentSite,new_cohort, bc_in) + if(recruitstatus==1)then new_cohort%co_hydr%is_newly_recruited = .true. endif + endif call insert_cohort(new_cohort, patchptr%tallest, patchptr%shortest, tnull, snull, & diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 50b91c3513..89e3257537 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -1002,13 +1002,12 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) endif if (temp_cohort%n > 0.0_r8 )then - if ( debug ) write(fates_log(),*) 'EDPhysiologyMod.F90 call create_cohort ' + call create_cohort(currentSite,currentPatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & temp_cohort%laimemory, cohortstatus,recruitstatus, temp_cohort%canopy_trim, currentPatch%NCL_p, & currentSite%spread, bc_in) - - + ! keep track of how many individuals were recruited for passing to history currentSite%recruitment_rate(ft) = currentSite%recruitment_rate(ft) + temp_cohort%n diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 81881a190a..45df34aebc 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -117,6 +117,7 @@ module FatesAllometryMod public :: decay_coeff_kn public :: StructureResetOfDH ! Method to set DBH to sync with structure biomass public :: CheckIntegratedAllometries + public :: CrownDepth public :: set_root_fraction ! Generic wrapper to calculate normalized ! root profiles @@ -1843,6 +1844,33 @@ subroutine h2d_martcano(h,p1,p2,p3,d,dddh) return end subroutine h2d_martcano + ! ===================================================================================== + + + subroutine CrownDepth(height,crown_depth) + + ! ----------------------------------------------------------------------------------- + ! This routine returns the depth of a plant's crown. Which is the length + ! from the bottom of the crown to the top in the vertical dimension. + ! + ! This code may be used as a wrapper if different hypotheses are wished to be + ! optioned. + ! ----------------------------------------------------------------------------------- + + real(r8),intent(in) :: height ! The height of the plant [m] + real(r8),intent(out) :: crown_depth ! The depth of the crown [m] + + ! Alternative Hypothesis: + ! crown depth from Poorter, Bongers & Bongers + ! crown_depth = exp(-1.169_r8)*cCohort%hite**1.098_r8 + + crown_depth = min(height,0.1_r8) + + return + end subroutine CrownDepth + + + ! ============================================================================= ! Specific diameter to crown area allometries ! ============================================================================= diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index f8567bc785..78046c1aa8 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -36,7 +36,10 @@ module FatesPlantHydraulicsMod use FatesConstantsMod, only : denh2o => dens_fresh_liquid_water use FatesConstantsMod, only : grav => grav_earth use FatesConstantsMod, only : ifalse, itrue - + use FatesConstantsMod, only : pi_const + use FatesConstantsMod, only : cm2_per_m2 + use FatesConstantsMod, only : g_per_kg + use EDParamsMod , only : hydr_kmax_rsurf use EDTypesMod , only : ed_site_type @@ -49,6 +52,7 @@ module FatesPlantHydraulicsMod use FatesInterfaceMod , only : hlm_use_planthydro use FatesAllometryMod, only : bsap_allom + use FatesAllometryMod, only : CrownDepth use FatesHydraulicsMemMod, only: ed_site_hydr_type use FatesHydraulicsMemMod, only: ed_cohort_hydr_type @@ -76,14 +80,12 @@ module FatesPlantHydraulicsMod use clm_time_manager , only : get_step_size, get_nstep - use FatesConstantsMod, only: cm2_per_m2 - use EDPftvarcon, only : EDPftvarcon_inst ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : isnan => shr_infnan_isnan - + implicit none @@ -117,6 +119,12 @@ module FatesPlantHydraulicsMod character(len=*), parameter, private :: sourcefile = & __FILE__ + + + ! We use this parameter as the value for which we set un-initialized values + real(r8), parameter :: un_initialized = -9.9e32_r8 + + ! ! !PUBLIC MEMBER FUNCTIONS: public :: AccumulateMortalityWaterStorage @@ -135,6 +143,11 @@ module FatesPlantHydraulicsMod public :: initTreeHydStates public :: updateSizeDepRhizHydProps public :: updateSizeDepRhizHydStates + public :: RestartHydrStates + public :: SavePreviousCompartmentVolumes + public :: SavePreviousRhizVolumes + public :: UpdateTreeHydrNodes + public :: UpdateTreeHydrLenVolCond !------------------------------------------------------------------------------ ! 01/18/16: Created by Brad Christoffersen @@ -171,19 +184,121 @@ subroutine hydraulics_drive( nsites, sites, bc_in,bc_out,dtime ) end select end subroutine Hydraulics_Drive - + ! ===================================================================================== + subroutine RestartHydrStates(sites,nsites,bc_in,bc_out) + + ! It is assumed that the following state variables have been read in by + ! the restart machinery. + ! + ! co_hydr%th_ag + ! co_hydr%th_troot + ! co_hydr%th_aroot + ! si_hydr%r_node_shell + ! si_hydr%v_shell + ! si_hydr%h2osoi_liqvol_shell + ! si_hydr%l_aroot_layer + ! + ! The goal of this subroutine is to call + ! the correct sequence of hydraulics initializations to repopulate + ! information that relies on these key states, as well as other vegetation + ! states such as carbon pools and plant geometry. + + type(ed_site_type) , intent(inout), target :: sites(nsites) + integer , intent(in) :: nsites + type(bc_in_type) , intent(in) :: bc_in(nsites) + type(bc_out_type) , intent(inout) :: bc_out(nsites) + + ! locals + ! ---------------------------------------------------------------------------------- + ! LL pointers + type(ed_patch_type),pointer :: cpatch ! current patch + type(ed_cohort_type),pointer :: ccohort ! current cohort + type(ed_cohort_hydr_type),pointer :: ccohort_hydr + integer :: s ! site loop counter + + do s = 1,nsites + + cpatch => sites(s)%oldest_patch + do while(associated(cpatch)) + + ccohort => cpatch%shortest + do while(associated(ccohort)) + + ccohort_hydr => ccohort%co_hydr + + ! This calculates node heights + call UpdateTreeHydrNodes(ccohort_hydr,ccohort%pft,ccohort%hite, & + sites(s)%si_hydr%nlevsoi_hyd,bc_in(s)) + + ! This calculates volumes, lengths and max conductances + call UpdateTreeHydrLenVolCond(ccohort,sites(s)%si_hydr%nlevsoi_hyd,bc_in(s)) + + ! Since this is a newly initialized plant, we set the previous compartment-size + ! equal to the ones we just calculated. + call SavePreviousCompartmentVolumes(ccohort_hydr) + + ! Set some generic initial values + ccohort_hydr%refill_days = 3.0_r8 + ccohort_hydr%lwp_mem(:) = 0.0_r8 + ccohort_hydr%lwp_stable = 0.0_r8 + ccohort_hydr%lwp_is_unstable = .false. + ccohort_hydr%flc_ag(:) = 1.0_r8 + ccohort_hydr%flc_troot(:) = 1.0_r8 + ccohort_hydr%flc_aroot(:) = 1.0_r8 + ccohort_hydr%flc_min_ag(:) = 1.0_r8 + ccohort_hydr%flc_min_troot(:) = 1.0_r8 + ccohort_hydr%flc_min_aroot(:) = 1.0_r8 + ccohort_hydr%refill_thresh = -0.01_r8 + ccohort_hydr%refill_days = 3.0_r8 + + ccohort => ccohort%taller + enddo + + cpatch => cpatch%younger + end do + + sites(s)%si_hydr%l_aroot_layer_init(:) = un_initialized + sites(s)%si_hydr%r_node_shell_init(:,:) = un_initialized + sites(s)%si_hydr%v_shell_init(:,:) = un_initialized + + + ! Update static quantities related to the rhizosphere + call UpdateSizeDepRhizVolLenCon(sites(s), bc_in(s)) + + ! We update the "initial" values of the rhizosphere after + ! the previous call to make sure that the conductances are updated + ! Now we set the prevous to the current so that the water states + ! are not perturbed + call SavePreviousRhizVolumes(sites(s), bc_in(s)) + + end do + + call UpdateH2OVeg(nsites,sites,bc_out) + + return + end subroutine RestartHydrStates + + ! ==================================================================================== + subroutine initTreeHydStates(site_p, cc_p, bc_in) + + ! REQUIRED INPUTS: ! + ! csite%si_hydr%psisoi_liq_innershell(:) + ! ccohort_hydr%z_node_troot(:) + ! ccohort_hydr%z_node_aroot + ! ccohort_hydr%z_node_ag + ! ! !DESCRIPTION: ! ! !USES: ! !ARGUMENTS: - type(ed_site_type), intent(inout), target :: site_p ! current cohort pointer - type(ed_cohort_type), intent(inout), target :: cc_p ! current cohort pointer - type(bc_in_type) , intent(in) :: bc_in + type(ed_site_type), intent(inout), target :: site_p ! current cohort pointer + type(ed_cohort_type), intent(inout), target :: cc_p ! current cohort pointer + type(bc_in_type) , intent(in) :: bc_in ! ! !LOCAL VARIABLES: type(ed_cohort_type), pointer :: cCohort @@ -259,298 +374,444 @@ subroutine initTreeHydStates(site_p, cc_p, bc_in) end subroutine initTreeHydStates + + ! ===================================================================================== + + + subroutine UpdateTreeHydrNodes(ccohort_hydr,ft,plant_height,nlevsoi_hyd,bc_in) + + ! -------------------------------------------------------------------------------- + ! This subroutine calculates the nodal heights critical to hydraulics in the plant + ! + ! Inputs: Plant height + ! Plant functional type + ! Number of soil hydraulic layers + ! + ! Outputs: cohort_hydr%z_node_ag(:) + ! %z_lower_ag(:) + ! %z_upper_ag(:) + ! %z_node_troot(:) + ! %z_lower_troot(:) + ! %z_upper_troot(:) + ! %z_node_aroot(:) + ! -------------------------------------------------------------------------------- + + ! Arguments + type(ed_cohort_hydr_type), intent(inout) :: ccohort_hydr + integer,intent(in) :: ft ! plant functional type index + real(r8), intent(in) :: plant_height ! [m] + integer,intent(in) :: nlevsoi_hyd ! number of soil hydro layers + type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions + + + ! Locals + + real(r8) :: roota ! root profile parameter a zeng2001_crootfr + real(r8) :: rootb ! root profile parameter b zeng2001_crootfr + real(r8) :: crown_depth ! crown depth for the plant [m] + real(r8) :: dz_canopy ! discrete crown depth intervals [m] + real(r8) :: z_stem ! the height of the plants stem below crown [m] + real(r8) :: dz_stem ! vertical stem discretization [m] + real(r8) :: dcumul_rf ! cumulative root distribution discretization [-] + real(r8) :: cumul_rf ! cumulative root distribution where depth is determined [-] + real(r8) :: z_cumul_rf ! depth at which cumul_rf occurs [m] + integer :: k ! Loop counter for compartments + + ! Crown Nodes + ! in special case where n_hypool_leaf = 1, the node height of the canopy + ! water pool is 1/2 the distance from the bottom of the canopy to the top of the tree + + roota = EDPftvarcon_inst%roota_par(ft) + rootb = EDPftvarcon_inst%rootb_par(ft) + + call CrownDepth(plant_height,crown_depth) + + dz_canopy = crown_depth / real(n_hypool_leaf,r8) + do k=1,n_hypool_leaf + ccohort_hydr%z_lower_ag(k) = plant_height - dz_canopy*real(k,r8) + ccohort_hydr%z_node_ag(k) = ccohort_hydr%z_lower_ag(k) + 0.5_r8*dz_canopy + ccohort_hydr%z_upper_ag(k) = ccohort_hydr%z_lower_ag(k) + dz_canopy + enddo + + + ! Stem Nodes + ! in special case where n_hypool_stem = 1, the node height of the stem water pool is + ! 1/2 the height from the ground to the bottom of the canopy + z_stem = plant_height - crown_depth + dz_stem = z_stem / real(n_hypool_stem,r8) + do k=n_hypool_leaf+1,n_hypool_ag + ccohort_hydr%z_upper_ag(k) = real(n_hypool_stem - (k - 1 - n_hypool_leaf),r8)*dz_stem + ccohort_hydr%z_node_ag(k) = ccohort_hydr%z_upper_ag(k) - 0.5_r8*dz_stem + ccohort_hydr%z_lower_ag(k) = ccohort_hydr%z_upper_ag(k) - dz_stem + enddo + + ! Transporting Root Nodes + ! in special case where n_hypool_troot = 1, the node depth of the single troot pool + ! is the depth at which 50% total root distribution is attained + dcumul_rf = 1._r8/real(n_hypool_troot,r8) + + do k=1,n_hypool_troot + cumul_rf = dcumul_rf*real(k,r8) + call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, & + 0.001_r8, 0.001_r8, cumul_rf, z_cumul_rf) + z_cumul_rf = min(z_cumul_rf, abs(bc_in%zi_sisl(nlevsoi_hyd))) + ccohort_hydr%z_lower_troot(k) = -z_cumul_rf + call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, & + 0.001_r8, 0.001_r8, cumul_rf-0.5_r8*dcumul_rf, z_cumul_rf) + z_cumul_rf = min(z_cumul_rf, abs(bc_in%zi_sisl(nlevsoi_hyd))) + ccohort_hydr%z_node_troot(k) = -z_cumul_rf + call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, & + 0.001_r8, 0.001_r8, cumul_rf-1.0_r8*dcumul_rf+1.E-10_r8, z_cumul_rf) + z_cumul_rf = min(z_cumul_rf, abs(bc_in%zi_sisl(nlevsoi_hyd))) + ccohort_hydr%z_upper_troot(k) = -z_cumul_rf + enddo + + + ! Absorbing root depth + ccohort_hydr%z_node_aroot(1:nlevsoi_hyd) = -bc_in%z_sisl(1:nlevsoi_hyd) + + + ! Shouldn't this be updating the upper and lower values as well? + ! (RGK 12-2018) + if(nlevsoi_hyd == 1) then + ccohort_hydr%z_node_troot(:) = ccohort_hydr%z_node_aroot(nlevsoi_hyd) + end if + + return + end subroutine UpdateTreeHydrNodes + ! ===================================================================================== - subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) + subroutine SavePreviousCompartmentVolumes(ccohort_hydr) + + type(ed_cohort_hydr_type),intent(inout) :: ccohort_hydr + + + ! Saving the current compartment volumes into an "initial" save-space + ! allows us to see how the compartments change size when plants + ! change size and effect water contents + + ccohort_hydr%v_ag_init(:) = ccohort_hydr%v_ag(:) + ccohort_hydr%v_troot_init(:) = ccohort_hydr%v_troot(:) + ccohort_hydr%v_aroot_layer_init(:) = ccohort_hydr%v_aroot_layer(:) + + return + end subroutine SavePreviousCompartmentVolumes + + ! ===================================================================================== + + subroutine updateSizeDepTreeHydProps(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 - ! - ! !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 - ! ! !USES: - use FatesConstantsMod , only : pi_const use shr_sys_mod , only : shr_sys_abort - ! - ! !ARGUMENTS: - type(ed_site_type) , intent(in) :: currentSite ! Site stuff - type(ed_cohort_type) , intent(inout), target :: cc_p ! current cohort pointer - type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions - ! !LOCAL VARIABLES: + ! 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 - type(ed_cohort_type), pointer :: cCohort - type(ed_patch_type), pointer :: cPatch - integer :: i,j,k,FT ! indices - real(r8) :: b_tot_carb ! total individual biomass in carbon units [kgC/indiv] - real(r8) :: b_bg_carb ! belowground biomass (coarse + fine roots) in carbon units [kgC/indiv] - real(r8) :: roota, rootb ! parameters for root distribution [m-1] - real(r8) :: latosa ! leaf:sapwood area ratio [m2/cm2] - ! TRANSPORTING ROOT QUANTITIES - real(r8) :: dcumul_rf ! cumulative root distribution discretization [-] - real(r8) :: cumul_rf ! cumulative root distribution where depth is determined [-] - real(r8) :: z_cumul_rf ! depth at which cumul_rf occurs [m] - real(r8) :: b_troot_carb ! transporting root biomass in carbon units [kgC/indiv] - real(r8) :: b_troot_biom ! transporting root biomass in dry wt units [kg/indiv] - real(r8) :: v_troot ! transporting root volume [m3/indiv] - real(r8) :: rootfr - ! CANOPY or LEAF QUANTITIES - real(r8) :: sla ! specific leaf area [cm2/g] - real(r8) :: depth_canopy ! crown (canopy) depth [m] - real(r8) :: dz_canopy ! vertical canopy discretization [m] - real(r8) :: a_sapwood_target ! sapwood cross-section area at reference height, at target biomass [m2] - real(r8) :: bsw_target ! sapwood carbon, at target [kgC] - real(r8) :: a_leaf_tot ! total leaf area [m2/indiv] - real(r8) :: b_canopy_carb ! total leaf (canopy) biomass in carbon units [kgC/indiv] - real(r8) :: b_canopy_biom ! total leaf (canopy) biomass in dry wt units [kg/indiv] - real(r8) :: v_canopy ! total leaf (canopy) volume [m3/indiv] - real(r8) :: denleaf ! leaf dry mass per unit fresh leaf volume [kg/m3] - ! STEM OR SAPWOOD QUANTITIES - real(r8) :: a_sapwood ! sapwood area [m2] - real(r8) :: v_sapwood ! sapwood volume [m3] - real(r8) :: z_stem ! tree height, minus any crown depth [m] - real(r8) :: dz_stem ! vertical stem discretization [m] - real(r8) :: b_woody_carb ! total woody biomass in carbon units [kgC/indiv] - real(r8) :: b_woody_bg_carb ! belowground woody biomass in carbon units [kgC/indiv] - real(r8) :: b_stem_carb ! aboveground stem biomass in carbon units [kgC/indiv] - real(r8) :: b_stem_biom ! aboveground stem biomass in dry wt units [kg/indiv] - real(r8) :: v_stem ! aboveground stem volume [m3/indiv] - !HYDRAULIC MAXIMUM CONDUCTANCES and assoc vars - real(r8) :: p=1._r8/3._r8 ! Savage et al. (2010) xylem taper exponent [-] - real(r8) :: rint_jansenchoat=22._r8 ! conduit radius at branch location where kmax measured, tropical mean [um] - real(r8) :: rint_petiole=10._r8 ! petiole conduit radius (assumed invariant, sensu Savage et al. 2010) [um] - real(r8) :: kmax_node_petiole ! maximum hydraulic conductivity at petiole [kg m-1 s-1 MPa-1] - real(r8) :: kmax_node1_nodekplus1(n_hypool_ag) ! cumulative kmax, petiole to node k+1, conduit taper effects excluded [kg s-1 MPa-1] - real(r8) :: kmax_node1_lowerk(n_hypool_ag) ! cumulative kmax, petiole to upper boundary of node k, conduit taper effects excluded [kg s-1 MPa-1] - real(r8) :: chi_node1_nodekplus1(n_hypool_ag)! ratio of cumulative kmax with taper effects included to that without [-] - real(r8) :: chi_node1_lowerk(n_hypool_ag) ! ratio of cumulative kmax with taper effects included to that without [-] - real(r8) :: kmax_treeag_tot ! total stem (petiole to transporting root node) 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) :: dz_node1_nodekplus1 ! cumulative distance between canopy node and node k + 1 [m] - real(r8) :: dz_node1_lowerk ! cumulative distance between canopy node and upper boundary of node k [m] - real(r8) :: leaf_c - real(r8) :: fnrt_c - real(r8) :: sapw_c - real(r8) :: struct_c - integer :: nlevsoi_hyd ! Number of soil hydraulic layers - integer :: nlevsoil ! Number of total soil layers + ! 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 + + ! Save the current vegetation compartment volumes into + ! a save space so that it can be compared with the updated quantity. + + call SavePreviousCompartmentVolumes(ccohort_hydr) + + ! This updates all of the z_node positions + call UpdateTreeHydrNodes(ccohort_hydr,ft,ccohort%hite,nlevsoi_hyd,bc_in) + ! 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. + + call UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) + + + end subroutine updateSizeDepTreeHydProps - nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd - nlevsoil = bc_in%nlevsoil - cCohort => cc_p - ccohort_hydr => cc_p%co_hydr - cPatch => cCohort%patchptr - FT = cCohort%pft - roota = EDPftvarcon_inst%roota_par(FT) - rootb = EDPftvarcon_inst%rootb_par(FT) - - leaf_c = cCohort%prt%GetState(leaf_organ, all_carbon_elements) - sapw_c = cCohort%prt%GetState(sapw_organ, all_carbon_elements) - fnrt_c = cCohort%prt%GetState(fnrt_organ, all_carbon_elements) - struct_c = cCohort%prt%GetState(struct_organ, all_carbon_elements) - - !roota = 4.372_r8 ! TESTING: deep (see Zeng 2001 Table 1) - !rootb = 0.978_r8 ! TESTING: deep (see Zeng 2001 Table 1) - !roota = 8.992_r8 ! TESTING: shallow (see Zeng 2001 Table 1) - !rootb = 8.992_r8 ! TESTING: shallow (see Zeng 2001 Table 1) - if(leaf_c>0.0) then !only update when bleaf >0 - b_woody_carb = sapw_c + struct_c - b_woody_bg_carb = (1.0_r8-EDPftvarcon_inst%allom_agb_frac(FT)) * b_woody_carb - - b_tot_carb = sapw_c + struct_c + leaf_c + fnrt_c - b_canopy_carb = leaf_c - b_bg_carb = (1.0_r8-EDPftvarcon_inst%allom_agb_frac(FT)) * b_tot_carb - - ! SAVE INITIAL VOLUMES - ccohort_hydr%v_ag_init(:) = ccohort_hydr%v_ag(:) - ccohort_hydr%v_troot_init(:) = ccohort_hydr%v_troot(:) - ccohort_hydr%v_aroot_layer_init(:) = ccohort_hydr%v_aroot_layer(:) - - ! CANOPY HEIGHT & CANOPY LEAF VOLUME - !in special case where n_hypool_leaf = 1, the node height of the canopy water pool is - !1/2 the distance from the bottom of the canopy to the top of the tree - !depth_canopy = exp(-1.169_r8)*cCohort%hite**1.098_r8 !! crown depth from Poorter, Bongers & Bongers - depth_canopy = min(cCohort%hite,0.1_r8) ! 0.0_r8 was default, now changed 01/14/2017 (BOC) - dz_canopy = depth_canopy / n_hypool_leaf - do k=1,n_hypool_leaf - ccohort_hydr%z_lower_ag(k) = cCohort%hite - dz_canopy*k - ccohort_hydr%z_node_ag(k) = ccohort_hydr%z_lower_ag(k) + 0.5_r8*dz_canopy - ccohort_hydr%z_upper_ag(k) = ccohort_hydr%z_lower_ag(k) + dz_canopy - enddo - b_canopy_biom = b_canopy_carb * C2B + ! ===================================================================================== - ! NOTE: SLATOP currently does not use any vertical scaling functions - ! but that may not be so forever. ie sla = slatop (RGK-082017) - sla = EDPftvarcon_inst%slatop(FT) * cm2_per_m2 ! m2/gC * cm2/m2 -> cm2/gC + subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) + + ! ----------------------------------------------------------------------------------- + ! This subroutine calculates three attributes of a plant: + ! 1) the volumes of storage compartments in the plants + ! 2) the lenghts of the organs + ! 3) the conductances + ! These and are not dependent on the hydraulic state of the + ! plant, it is more about the structural characteristics and how much biomass + ! is present in the different tissues. + ! + ! Inputs, plant geometries, plant carbon pools, z_node values + ! + ! ----------------------------------------------------------------------------------- - denleaf = -2.3231_r8*sla/C2B + 781.899_r8 ! empirical regression data from leaves at Caxiuana (~ 8 spp) - v_canopy = b_canopy_biom / denleaf - ccohort_hydr%v_ag(1:n_hypool_leaf) = v_canopy / n_hypool_leaf + ! Arguments + type(ed_cohort_type),intent(inout) :: ccohort + 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 + integer :: j,k + integer :: ft ! Plant functional type index + real(r8) :: roota ! root profile parameter a zeng2001_crootfr + real(r8) :: rootb ! root profile parameter b zeng2001_crootfr + real(r8) :: leaf_c ! Current amount of leaf carbon in the plant [kg] + real(r8) :: fnrt_c ! Current amount of fine-root carbon in the plant [kg] + real(r8) :: sapw_c ! Current amount of sapwood carbon in the plant [kg] + real(r8) :: struct_c ! Current amount of structural carbon in the plant [kg] + real(r8) :: b_canopy_carb ! total leaf (canopy) biomass in carbon units [kgC/indiv] + real(r8) :: b_canopy_biom ! total leaf (canopy) biomass in dry wt units [kg/indiv] + real(r8) :: b_woody_carb ! total woody biomass in carbon units [kgC/indiv] + real(r8) :: b_woody_bg_carb ! belowground woody biomass in carbon units [kgC/indiv] + real(r8) :: b_stem_carb ! aboveground stem biomass in carbon units [kgC/indiv] + real(r8) :: b_stem_biom ! aboveground stem biomass in dry wt units [kg/indiv] + real(r8) :: b_bg_carb ! belowground biomass (coarse + fine roots) in carbon units [kgC/indiv] + real(r8) :: b_tot_carb ! total individual biomass in carbon units [kgC/indiv] + real(r8) :: v_stem ! aboveground stem volume [m3/indiv] + real(r8) :: z_stem ! the height of the plants stem below crown [m] + real(r8) :: sla ! specific leaf area [cm2/g] + real(r8) :: v_canopy ! total leaf (canopy) volume [m3/indiv] + real(r8) :: denleaf ! leaf dry mass per unit fresh leaf volume [kg/m3] + real(r8) :: a_sapwood ! sapwood area [m2] + real(r8) :: a_sapwood_target ! sapwood cross-section area at reference height, at target biomass [m2] + real(r8) :: bsw_target ! sapwood carbon, at target [kgC] + real(r8) :: v_sapwood ! sapwood volume [m3] + real(r8) :: b_troot_carb ! transporting root biomass in carbon units [kgC/indiv] + real(r8) :: b_troot_biom ! transporting root biomass in dry wt units [kg/indiv] + real(r8) :: v_troot ! transporting root volume [m3/indiv] + real(r8) :: rootfr ! mass fraction of roots in each layer [kg/kg] + real(r8) :: crown_depth ! Depth of the plant's crown [m] + real(r8) :: kmax_node1_nodekplus1(n_hypool_ag) ! cumulative kmax, petiole to node k+1, + ! conduit taper effects excluded [kg s-1 MPa-1] + real(r8) :: kmax_node1_lowerk(n_hypool_ag) ! cumulative kmax, petiole to upper boundary of node k, + ! conduit taper effects excluded [kg s-1 MPa-1] + real(r8) :: chi_node1_nodekplus1(n_hypool_ag) ! ratio of cumulative kmax with taper effects + ! included to that without [-] + real(r8) :: chi_node1_lowerk(n_hypool_ag) ! ratio of cumulative kmax with taper effects + ! included to that without [-] + real(r8) :: dz_node1_nodekplus1 ! cumulative distance between canopy + ! node and node k + 1 [m] + real(r8) :: dz_node1_lowerk ! cumulative distance between canopy + ! node and upper boundary of node k [m] + real(r8) :: kmax_treeag_tot ! total stem (petiole to transporting root node) + ! 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 [-] - ! STEM HEIGHT & VOLUME - !in special case where n_hypool_stem = 1, the node height of the stem water pool is - !1/2 the height from the ground to the bottom of the canopy - z_stem = cCohort%hite - depth_canopy - dz_stem = z_stem / n_hypool_stem - do k=n_hypool_leaf+1,n_hypool_ag - ccohort_hydr%z_upper_ag(k) = (n_hypool_stem - (k - 1 - n_hypool_leaf))*dz_stem - ccohort_hydr%z_node_ag(k) = ccohort_hydr%z_upper_ag(k) - 0.5_r8*dz_stem - ccohort_hydr%z_lower_ag(k) = ccohort_hydr%z_upper_ag(k) - dz_stem - enddo - b_stem_carb = b_tot_carb - b_bg_carb - b_canopy_carb - b_stem_biom = b_stem_carb * C2B ! kg DM - v_stem = b_stem_biom / (EDPftvarcon_inst%wood_density(FT)*1.e3_r8) !BOC...may be needed for testing/comparison w/ v_sapwood - !a_leaf_tot = b_canopy_carb * sla * 1.e3_r8 / 1.e4_r8 ! m2 leaf = kg leaf DM * cm2/g * 1000g/1kg * 1m2/10000cm2 - a_leaf_tot = cCohort%treelai * cCohort%c_area/cCohort%n ! calculate the leaf area based on the leaf profile + ccohort_hydr => ccohort%co_hydr + ft = ccohort%pft - call bsap_allom(cCohort%dbh,cCohort%pft,cCohort%canopy_trim,a_sapwood_target,bsw_target) - - a_sapwood = a_sapwood_target + leaf_c = ccohort%prt%GetState(leaf_organ, all_carbon_elements) + sapw_c = ccohort%prt%GetState(sapw_organ, all_carbon_elements) + fnrt_c = ccohort%prt%GetState(fnrt_organ, all_carbon_elements) + struct_c = ccohort%prt%GetState(struct_organ, all_carbon_elements) - ! or .... - ! a_sapwood = a_sapwood_target * ccohort%bsw / bsw_target + roota = EDPftvarcon_inst%roota_par(ft) + rootb = EDPftvarcon_inst%rootb_par(ft) - ! a_sapwood = a_leaf_tot / EDPftvarcon_inst%allom_latosa_int(FT)*1.e-4_r8 - ! m2 sapwood = m2 leaf * cm2 sapwood/m2 leaf *1.0e-4m2 - ! or ... - !a_sapwood = a_leaf_tot / ( 0.001_r8 + 0.025_r8 * cCohort%hite ) * 1.e-4_r8 - - v_sapwood = a_sapwood * z_stem - ccohort_hydr%v_ag(n_hypool_leaf+1:n_hypool_ag) = v_sapwood / n_hypool_stem + !roota = 4.372_r8 ! TESTING: deep (see Zeng 2001 Table 1) + !rootb = 0.978_r8 ! TESTING: deep (see Zeng 2001 Table 1) + !roota = 8.992_r8 ! TESTING: shallow (see Zeng 2001 Table 1) + !rootb = 8.992_r8 ! TESTING: shallow (see Zeng 2001 Table 1) + + if(leaf_c > 0._r8) then - ! TRANSPORTING ROOT DEPTH & VOLUME - !in special case where n_hypool_troot = 1, the node depth of the single troot pool - !is the depth at which 50% total root distribution is attained - dcumul_rf = 1._r8/real(n_hypool_troot,r8) - do k=1,n_hypool_troot - cumul_rf = dcumul_rf*k - call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, & - 0.001_r8, 0.001_r8, cumul_rf, z_cumul_rf) - z_cumul_rf = min(z_cumul_rf, abs(bc_in%zi_sisl(nlevsoi_hyd))) - ccohort_hydr%z_lower_troot(k) = -z_cumul_rf - call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, & - 0.001_r8, 0.001_r8, cumul_rf-0.5_r8*dcumul_rf, z_cumul_rf) - z_cumul_rf = min(z_cumul_rf, abs(bc_in%zi_sisl(nlevsoi_hyd))) - ccohort_hydr%z_node_troot(k) = -z_cumul_rf - call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, & - 0.001_r8, 0.001_r8, cumul_rf-1.0_r8*dcumul_rf+1.E-10_r8, z_cumul_rf) - z_cumul_rf = min(z_cumul_rf, abs(bc_in%zi_sisl(nlevsoi_hyd))) - ccohort_hydr%z_upper_troot(k) = -z_cumul_rf - enddo + ! ------------------------------------------------------------------------------ + ! Part 1. Set the volumes of the leaf, stem and root compartments + ! and lenghts of the roots + ! ------------------------------------------------------------------------------ + b_woody_carb = sapw_c + struct_c + b_woody_bg_carb = (1.0_r8-EDPftvarcon_inst%allom_agb_frac(ft)) * b_woody_carb + b_tot_carb = sapw_c + struct_c + leaf_c + fnrt_c + b_canopy_carb = leaf_c + b_bg_carb = (1.0_r8-EDPftvarcon_inst%allom_agb_frac(ft)) * b_tot_carb + b_canopy_biom = b_canopy_carb * C2B + + ! NOTE: SLATOP currently does not use any vertical scaling functions + ! but that may not be so forever. ie sla = slatop (RGK-082017) + ! m2/gC * cm2/m2 -> cm2/gC + sla = EDPftvarcon_inst%slatop(ft) * cm2_per_m2 + + ! empirical regression data from leaves at Caxiuana (~ 8 spp) + denleaf = -2.3231_r8*sla/C2B + 781.899_r8 + v_canopy = b_canopy_biom / denleaf - !Determine belowground biomass as a function of total (sapwood, heartwood, leaf, fine root) biomass - !then subtract out the fine root biomass to get coarse (transporting) root biomass + ccohort_hydr%v_ag(1:n_hypool_leaf) = v_canopy / real(n_hypool_leaf,r8) - b_troot_carb = b_woody_bg_carb - b_troot_biom = b_troot_carb * C2B - v_troot = b_troot_biom / (EDPftvarcon_inst%wood_density(FT)*1.e3_r8) - ccohort_hydr%v_troot(:) = v_troot / n_hypool_troot !! BOC not sure if/how we should multiply this by the sapwood fraction + + b_stem_carb = b_tot_carb - b_bg_carb - b_canopy_carb + b_stem_biom = b_stem_carb * C2B ! kg DM - ! ABSORBING ROOT DEPTH, LENGTH & VOLUME - ccohort_hydr%z_node_aroot(1:nlevsoi_hyd) = -bc_in%z_sisl(1:nlevsoi_hyd) + !BOC...may be needed for testing/comparison w/ v_sapwood + ! kg / ( g cm-3 * cm3/m3 * kg/g ) -> m3 + v_stem = b_stem_biom / (EDPftvarcon_inst%wood_density(ft)*1.e3_r8 ) + + ! calculate the sapwood cross-sectional area + call bsap_allom(ccohort%dbh,ccohort%pft,ccohort%canopy_trim,a_sapwood_target,bsw_target) + a_sapwood = a_sapwood_target + + ! Alternative ways to calculate sapwood cross section + ! or .... + ! a_sapwood = a_sapwood_target * ccohort%bsw / bsw_target + + ! a_sapwood = a_leaf_tot / EDPftvarcon_inst%allom_latosa_int(ft)*1.e-4_r8 + ! m2 sapwood = m2 leaf * cm2 sapwood/m2 leaf *1.0e-4m2 + ! or ... + !a_sapwood = a_leaf_tot / ( 0.001_r8 + 0.025_r8 * ccohort%hite ) * 1.e-4_r8 + + call CrownDepth(ccohort%hite,crown_depth) + z_stem = ccohort%hite - crown_depth + v_sapwood = a_sapwood * z_stem + ccohort_hydr%v_ag(n_hypool_leaf+1:n_hypool_ag) = v_sapwood / n_hypool_stem + + + ! Determine belowground biomass as a function of total (sapwood, heartwood, + ! leaf, fine root) biomass then subtract out the fine root biomass to get + ! coarse (transporting) root biomass + + b_troot_carb = b_woody_bg_carb + b_troot_biom = b_troot_carb * C2B + v_troot = b_troot_biom / (EDPftvarcon_inst%wood_density(ft)*1.e3_r8) + + !! BOC not sure if/how we should multiply this by the sapwood fraction + ccohort_hydr%v_troot(:) = v_troot / n_hypool_troot - ccohort_hydr%l_aroot_tot = fnrt_c*C2B*EDPftvarcon_inst%hydr_srl(FT) - !ccohort_hydr%v_aroot_tot = fnrt_c/EDecophyscon%ccontent(FT)/EDecophyscon%rootdens(FT) - ccohort_hydr%v_aroot_tot = pi_const*(EDPftvarcon_inst%hydr_rs2(FT)**2._r8)*ccohort_hydr%l_aroot_tot - !ccohort_hydr%l_aroot_tot = ccohort_hydr%v_aroot_tot/(pi_const*EDecophyscon%rs2(FT)**2) - if(nlevsoi_hyd == 1) then - ccohort_hydr%l_aroot_layer(nlevsoi_hyd) = ccohort_hydr%l_aroot_tot - ccohort_hydr%v_aroot_layer(nlevsoi_hyd) = ccohort_hydr%v_aroot_tot - else - ! ccohort_hydr%l_aroot_layer(:) = cPatch%rootfr_ft(FT,:)*ccohort_hydr%l_aroot_tot - ! ccohort_hydr%v_aroot_layer(:) = cPatch%rootfr_ft(FT,:)*ccohort_hydr%v_aroot_tot - do j=1,nlevsoi_hyd - if(j == 1) then - rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j)) - else - rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j)) - & + ! Estimate absorbing root total length (all layers) + ! ------------------------------------------------------------------------------ + ccohort_hydr%l_aroot_tot = fnrt_c*C2B*EDPftvarcon_inst%hydr_srl(ft) + + ! Estimate absorbing root volume (all layers) + ! ------------------------------------------------------------------------------ + ccohort_hydr%v_aroot_tot = pi_const * (EDPftvarcon_inst%hydr_rs2(ft)**2._r8) * & + ccohort_hydr%l_aroot_tot + + + ! Partition the total absorbing root lengths and volumes into the active soil layers + ! ------------------------------------------------------------------------------ + if(nlevsoi_hyd == 1) then + ccohort_hydr%l_aroot_layer(nlevsoi_hyd) = ccohort_hydr%l_aroot_tot + ccohort_hydr%v_aroot_layer(nlevsoi_hyd) = ccohort_hydr%v_aroot_tot + else + do j=1,nlevsoi_hyd + if(j == 1) then + rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j)) + else + rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j)) - & zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j-1)) - end if - ccohort_hydr%l_aroot_layer(j) = rootfr*ccohort_hydr%l_aroot_tot - ccohort_hydr%v_aroot_layer(j) = rootfr*ccohort_hydr%v_aroot_tot - end do - end if - if(nlevsoi_hyd == 1) then - ccohort_hydr%z_node_troot(:) = ccohort_hydr%z_node_aroot(nlevsoi_hyd) - end if + end if + ccohort_hydr%l_aroot_layer(j) = rootfr*ccohort_hydr%l_aroot_tot + ccohort_hydr%v_aroot_layer(j) = rootfr*ccohort_hydr%v_aroot_tot + end do + end if + - ! MAXIMUM (SIZE-DEPENDENT) HYDRAULIC CONDUCTANCES - ! first estimate cumulative (petiole to node k) conductances without taper as well as the chi taper function - do k=n_hypool_leaf,n_hypool_ag - dz_node1_lowerk = ccohort_hydr%z_node_ag(n_hypool_leaf) - ccohort_hydr%z_lower_ag(k) - if(k < n_hypool_ag) then - dz_node1_nodekplus1 = ccohort_hydr%z_node_ag(n_hypool_leaf) - ccohort_hydr%z_node_ag(k+1) - else - dz_node1_nodekplus1 = ccohort_hydr%z_node_ag(n_hypool_leaf) - ccohort_hydr%z_node_troot(1) - end if - kmax_node1_nodekplus1(k) = EDPftvarcon_inst%hydr_kmax_node(FT,2) * a_sapwood / dz_node1_nodekplus1 - kmax_node1_lowerk(k) = EDPftvarcon_inst%hydr_kmax_node(FT,2) * a_sapwood / dz_node1_lowerk - chi_node1_nodekplus1(k) = xylemtaper(p, dz_node1_nodekplus1) - chi_node1_lowerk(k) = xylemtaper(p, dz_node1_lowerk) - if(.not.do_kbound_upstream) then - if(depth_canopy == 0._r8) then - write(fates_log(),*) 'do_kbound_upstream requires a nonzero canopy depth ' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end if - enddo - ! then calculate the conductances at node boundaries as the difference of cumulative conductances - do k=n_hypool_leaf,n_hypool_ag - if(k == n_hypool_leaf) then - ccohort_hydr%kmax_bound(k) = kmax_node1_nodekplus1(k) * chi_node1_nodekplus1(k) - ccohort_hydr%kmax_lower(k) = kmax_node1_lowerk(k) * chi_node1_lowerk(k) - else - ccohort_hydr%kmax_bound(k) = ( 1._r8/(kmax_node1_nodekplus1(k) *chi_node1_nodekplus1(k) ) - & - 1._r8/(kmax_node1_nodekplus1(k-1)*chi_node1_nodekplus1(k-1)) ) ** (-1._r8) - ccohort_hydr%kmax_lower(k) = ( 1._r8/(kmax_node1_lowerk(k) *chi_node1_lowerk(k) ) - & - 1._r8/(kmax_node1_nodekplus1(k-1)*chi_node1_nodekplus1(k-1)) ) ** (-1._r8) - end if - if(k < n_hypool_ag) then - ccohort_hydr%kmax_upper(k+1) = ( 1._r8/(kmax_node1_nodekplus1(k) *chi_node1_nodekplus1(k) ) - & - 1._r8/(kmax_node1_lowerk(k) *chi_node1_lowerk(k) ) ) ** (-1._r8) - else if(k == n_hypool_ag) then - ccohort_hydr%kmax_upper_troot = ( 1._r8/(kmax_node1_nodekplus1(k) *chi_node1_nodekplus1(k) ) - & - 1._r8/(kmax_node1_lowerk(k) *chi_node1_lowerk(k) ) ) ** (-1._r8) - end if - !!!!!!!!!! FOR TESTING ONLY - !ccohort_hydr%kmax_bound(:) = 0.02_r8 ! Diurnal lwp variation in coldstart: -0.1 MPa - ! Diurnal lwp variation in large-tree (50cmDBH) coldstart: less than -0.01 MPa - !ccohort_hydr%kmax_bound(:) = 0.0016_r8 ! Diurnal lwp variation in coldstart: -0.8 - 1.0 MPa - ! Diurnal lwp variation in large-tree (50cmDBH) coldstart: -1.5 - 2.0 MPa [seemingly unstable] - !ccohort_hydr%kmax_bound(:) = 0.0008_r8 ! Diurnal lwp variation in coldstart: -1.5 - 2.0 MPa - ! Diurnal lwp variation in large-tree (50cmDBH) coldstart: -2.0 - 3.0 MPa [seemingly unstable] - !ccohort_hydr%kmax_bound(:) = 0.0005_r8 ! Diurnal lwp variation in coldstart: -2.0 - 3.0 MPa and one -5 MPa outlier - ! Diurnal lwp variation in large-tree (50cmDBH) coldstart: -3.0 - 4.0 MPa and one -10 MPa outlier [Unstable] - !!!!!!!!!! - enddo - ! finally, estimate the remaining tree conductance belowground as a residual - kmax_treeag_tot = sum(1._r8/ccohort_hydr%kmax_bound(n_hypool_leaf:n_hypool_ag))**(-1._r8) - kmax_tot = EDPftvarcon_inst%hydr_rfrac_stem(FT) * kmax_treeag_tot - ccohort_hydr%kmax_treebg_tot = ( 1._r8/kmax_tot - 1._r8/kmax_treeag_tot ) ** (-1._r8) - if(nlevsoi_hyd == 1) then - ccohort_hydr%kmax_treebg_layer(:) = ccohort_hydr%kmax_treebg_tot * cPatch%rootfr_ft(FT,:) - else - do j=1,nlevsoi_hyd - if(j == 1) then - rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j)) - else - rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j)) - & + ! ------------------------------------------------------------------------------ + ! Part II. Set maximum (size-dependent) hydraulic conductances + ! ------------------------------------------------------------------------------ + + ! first estimate cumulative (petiole to node k) conductances + ! without taper as well as the chi taper function + + do k=n_hypool_leaf,n_hypool_ag + dz_node1_lowerk = ccohort_hydr%z_node_ag(n_hypool_leaf) & + - ccohort_hydr%z_lower_ag(k) + if(k < n_hypool_ag) then + dz_node1_nodekplus1 = ccohort_hydr%z_node_ag(n_hypool_leaf) & + - ccohort_hydr%z_node_ag(k+1) + else + dz_node1_nodekplus1 = ccohort_hydr%z_node_ag(n_hypool_leaf) & + - ccohort_hydr%z_node_troot(1) + end if + kmax_node1_nodekplus1(k) = EDPftvarcon_inst%hydr_kmax_node(ft,2) * a_sapwood / dz_node1_nodekplus1 + kmax_node1_lowerk(k) = EDPftvarcon_inst%hydr_kmax_node(ft,2) * a_sapwood / dz_node1_lowerk + chi_node1_nodekplus1(k) = xylemtaper(taper_exponent, dz_node1_nodekplus1) + chi_node1_lowerk(k) = xylemtaper(taper_exponent, dz_node1_lowerk) + if(.not.do_kbound_upstream) then + if(crown_depth == 0._r8) then + write(fates_log(),*) 'do_kbound_upstream requires a nonzero canopy depth ' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + enddo + + + ! then calculate the conductances at node boundaries as the difference of cumulative conductances + do k=n_hypool_leaf,n_hypool_ag + if(k == n_hypool_leaf) then + ccohort_hydr%kmax_bound(k) = kmax_node1_nodekplus1(k) * chi_node1_nodekplus1(k) + ccohort_hydr%kmax_lower(k) = kmax_node1_lowerk(k) * chi_node1_lowerk(k) + else + ccohort_hydr%kmax_bound(k) = ( 1._r8/(kmax_node1_nodekplus1(k) *chi_node1_nodekplus1(k) ) - & + 1._r8/(kmax_node1_nodekplus1(k-1)*chi_node1_nodekplus1(k-1)) ) ** (-1._r8) + ccohort_hydr%kmax_lower(k) = ( 1._r8/(kmax_node1_lowerk(k) *chi_node1_lowerk(k) ) - & + 1._r8/(kmax_node1_nodekplus1(k-1)*chi_node1_nodekplus1(k-1)) ) ** (-1._r8) + end if + if(k < n_hypool_ag) then + ccohort_hydr%kmax_upper(k+1) = ( 1._r8/(kmax_node1_nodekplus1(k) *chi_node1_nodekplus1(k) ) - & + 1._r8/(kmax_node1_lowerk(k) *chi_node1_lowerk(k) ) ) ** (-1._r8) + else if(k == n_hypool_ag) then + ccohort_hydr%kmax_upper_troot = ( 1._r8/(kmax_node1_nodekplus1(k) *chi_node1_nodekplus1(k) ) - & + 1._r8/(kmax_node1_lowerk(k) *chi_node1_lowerk(k) ) ) ** (-1._r8) + end if + + !!!!!!!!!! FOR TESTING ONLY + !ccohort_hydr%kmax_bound(:) = 0.02_r8 ! Diurnal lwp variation in coldstart: -0.1 MPa + ! Diurnal lwp variation in large-tree (50cmDBH) coldstart: less than -0.01 MPa + !ccohort_hydr%kmax_bound(:) = 0.0016_r8 ! Diurnal lwp variation in coldstart: -0.8 - 1.0 MPa + ! Diurnal lwp variation in large-tree (50cmDBH) coldstart: -1.5 - 2.0 MPa [seemingly unstable] + !ccohort_hydr%kmax_bound(:) = 0.0008_r8 ! Diurnal lwp variation in coldstart: -1.5 - 2.0 MPa + ! Diurnal lwp variation in large-tree (50cmDBH) coldstart: -2.0 - 3.0 MPa [seemingly unstable] + !ccohort_hydr%kmax_bound(:) = 0.0005_r8 ! Diurnal lwp variation in coldstart: -2.0 - 3.0 MPa and one -5 MPa outlier + ! Diurnal lwp variation in large-tree (50cmDBH) coldstart: -3.0 - 4.0 MPa and one -10 MPa outlier [Unstable] + !!!!!!!!!! + + enddo + + ! finally, estimate the remaining tree conductance belowground as a residual + kmax_treeag_tot = sum(1._r8/ccohort_hydr%kmax_bound(n_hypool_leaf:n_hypool_ag))**(-1._r8) + kmax_tot = EDPftvarcon_inst%hydr_rfrac_stem(ft) * kmax_treeag_tot + ccohort_hydr%kmax_treebg_tot = ( 1._r8/kmax_tot - 1._r8/kmax_treeag_tot ) ** (-1._r8) + + if(nlevsoi_hyd == 1) then + ccohort_hydr%kmax_treebg_layer(:) = ccohort_hydr%kmax_treebg_tot * & + ccohort%patchptr%rootfr_ft(ft,:) + else + do j=1,nlevsoi_hyd + if(j == 1) then + rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j)) + else + rootfr = zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j)) - & zeng2001_crootfr(roota, rootb, bc_in%zi_sisl(j-1)) - end if - ccohort_hydr%kmax_treebg_layer(j) = rootfr*ccohort_hydr%kmax_treebg_tot - end do - end if - end if !check for bleaf - end subroutine updateSizeDepTreeHydProps + end if + 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 updateSizeDepTreeHydStates(currentSite,cc_p) + subroutine updateSizeDepTreeHydStates(currentSite,ccohort) ! ! !DESCRIPTION: ! @@ -560,11 +821,10 @@ subroutine updateSizeDepTreeHydStates(currentSite,cc_p) use EDTypesMod , only : AREA ! !ARGUMENTS: - type(ed_site_type) , intent(in) :: currentSite ! Site stuff - type(ed_cohort_type) , intent(inout), target :: cc_p ! current cohort pointer + type(ed_site_type) , intent(in) :: currentSite ! Site stuff + type(ed_cohort_type) , intent(inout) :: ccohort ! ! !LOCAL VARIABLES: - type(ed_cohort_type), pointer :: cCohort type(ed_cohort_hydr_type), pointer :: ccohort_hydr type(ed_site_hydr_type),pointer :: csite_hydr integer :: j,k,FT ! indices @@ -573,11 +833,10 @@ subroutine updateSizeDepTreeHydStates(currentSite,cc_p) real(r8) :: th_troot_uncorr(n_hypool_troot) ! uncorrected transporting root water content[m3 m-3] real(r8) :: th_aroot_uncorr(currentSite%si_hydr%nlevsoi_hyd) ! uncorrected absorbing root water content[m3 m-3] real(r8), parameter :: small_theta_num = 1.e-7_r8 ! avoids theta values equalling thr or ths [m3 m-3] - integer :: nstep !number of time steps + integer :: nstep !number of time steps !----------------------------------------------------------------------- - cCohort => cc_p - ccohort_hydr => cCohort%co_hydr + ccohort_hydr => ccohort%co_hydr FT = cCohort%pft ! MAYBE ADD A NAN CATCH? If updateSizeDepTreeHydProps() was not called twice prior to the first @@ -879,6 +1138,7 @@ end subroutine FuseCohortHydraulics ! ===================================================================================== ! Initialization Routines ! ===================================================================================== + subroutine InitHydrCohort(currentSite,currentCohort) ! Arguments @@ -890,8 +1150,9 @@ subroutine InitHydrCohort(currentSite,currentCohort) allocate(ccohort_hydr) currentCohort%co_hydr => ccohort_hydr call ccohort_hydr%AllocateHydrCohortArrays(currentSite%si_hydr%nlevsoi_hyd) - ccohort_hydr%is_newly_recruited = .false. + ccohort_hydr%is_newly_recruited = .false. + end subroutine InitHydrCohort ! ===================================================================================== @@ -1193,7 +1454,28 @@ end subroutine RecruitWUptake ! ===================================================================================== - subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) + subroutine SavePreviousRhizVolumes(currentSite, bc_in) + + ! !ARGUMENTS: + type(ed_site_type) , intent(inout), target :: currentSite + type(bc_in_type) , intent(in) :: bc_in + type(ed_site_hydr_type), pointer :: csite_hydr + integer :: nlevsoi_hyd + + csite_hydr => currentSite%si_hydr + nlevsoi_hyd = csite_hydr%nlevsoi_hyd + + csite_hydr%l_aroot_layer_init(:) = csite_hydr%l_aroot_layer(:) + csite_hydr%r_node_shell_init(:,:) = csite_hydr%r_node_shell(:,:) + csite_hydr%v_shell_init(:,:) = csite_hydr%v_shell(:,:) + + return + end subroutine SavePreviousRhizVolumes + + ! ====================================================================================== + + subroutine UpdateSizeDepRhizVolLenCon(currentSite, bc_in) + ! ! !DESCRIPTION: Updates size of 'representative' rhizosphere -- node radii, volumes. ! As fine root biomass (and thus absorbing root length) increases, this characteristic @@ -1201,7 +1483,6 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) ! the same. ! ! !USES: - use FatesConstantsMod , only : pi_const use EDTypesMod , only : AREA ! !ARGUMENTS: @@ -1228,19 +1509,15 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) integer :: nlevsoi_hyd !----------------------------------------------------------------------- - + csite_hydr => currentSite%si_hydr nlevsoi_hyd = csite_hydr%nlevsoi_hyd - csite_hydr%l_aroot_layer_init(:) = csite_hydr%l_aroot_layer(:) - csite_hydr%r_node_shell_init(:,:) = csite_hydr%r_node_shell(:,:) - csite_hydr%v_shell_init(:,:) = csite_hydr%v_shell(:,:) - ! update cohort-level root length density and accumulate it across cohorts and patches to the column level csite_hydr%l_aroot_layer(:) = 0._r8 cPatch => currentSite%youngest_patch do while(associated(cPatch)) - cCohort => cPatch%tallest + cCohort => cPatch%tallest do while(associated(cCohort)) ccohort_hydr => cCohort%co_hydr csite_hydr%l_aroot_layer(:) = csite_hydr%l_aroot_layer(:) + ccohort_hydr%l_aroot_layer(:)*cCohort%n @@ -1282,8 +1559,10 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) 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) / & @@ -1329,6 +1608,43 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) end if !has l_aroot_layer changed? enddo ! loop over soil layers + + + return + end subroutine UpdateSizeDepRhizVolLenCon + + + ! ===================================================================================== + + + subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) + ! + ! !DESCRIPTION: Updates size of 'representative' rhizosphere -- node radii, volumes. + ! As fine root biomass (and thus absorbing root length) increases, this characteristic + ! rhizosphere shrinks even though the total volume of soil tapped by fine roots remains + ! the same. + ! + ! !USES: + + use EDTypesMod , only : AREA + + ! !ARGUMENTS: + type(ed_site_type) , intent(inout), target :: currentSite + type(bc_in_type) , intent(in) :: bc_in + + + ! Save current volumes, lenghts and nodes to an "initial" + ! used to calculate effects in states later on. + + call SavePreviousRhizVolumes(currentSite, bc_in) + + ! Update the properties of the vegetation-soil hydraulic environment + ! these are independent on the water state + + call UpdateSizeDepRhizVolLenCon(currentSite, bc_in) + + + return end subroutine updateSizeDepRhizHydProps ! ================================================================================= @@ -3045,7 +3361,7 @@ subroutine Hydraulics_1DSolve(cc_p, ft, z_node, v_node, ths_node, thr_node, kmax ccohort_hydr%iterh2 = real(iterh2) ! WATER BALANCE ERROR-HANDLING - if (abs(we_local) > thresh) then + if ( (abs(we_local) > thresh) .and. debug) then write(fates_log(),*)'WARNING: plant hydraulics water balance error exceeds threshold of ',& thresh else if (abs(we_local) > thresh_break) then @@ -4727,7 +5043,7 @@ subroutine shellGeom(l_aroot, rs1, area, dz, r_out_shell, r_node_shell, v_shell) ! the same. ! ! !USES: - use FatesConstantsMod, only : pi_const + ! ! !ARGUMENTS: real(r8) , intent(in) :: l_aroot diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90 index eb233c66d1..06545b3f10 100644 --- a/main/FatesHydraulicsMemMod.F90 +++ b/main/FatesHydraulicsMemMod.F90 @@ -81,9 +81,9 @@ module FatesHydraulicsMemMod ! may or may not cross that with a simple or ! non-simple layering - real(r8),allocatable :: v_shell(:,:) ! Volume of rhizosphere compartment (m) - real(r8),allocatable :: v_shell_init(:,:) ! Previous volume of rhizosphere compartment (m) - real(r8),allocatable :: v_shell_1D(:) ! Volume of rhizosphere compartment (m) + real(r8),allocatable :: v_shell(:,:) ! Volume of rhizosphere compartment (m3) + real(r8),allocatable :: v_shell_init(:,:) ! Previous volume of rhizosphere compartment (m3) + real(r8),allocatable :: v_shell_1D(:) ! Volume of rhizosphere compartment (m3) real(r8),allocatable :: r_node_shell(:,:) ! Nodal radius of rhizosphere compartment (m) real(r8),allocatable :: r_node_shell_init(:,:) ! Previous Nodal radius of rhizosphere compartment (m) real(r8),allocatable :: l_aroot_layer(:) ! Total length (across cohorts) of absorbing diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index a710fd922e..fa851b32fc 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -14,10 +14,18 @@ module FatesRestartInterfaceMod use FatesRestartVariableMod, only : fates_restart_variable_type use FatesInterfaceMod, only : bc_in_type use FatesInterfaceMod, only : bc_out_type + use FatesInterfaceMod, only : hlm_use_planthydro + use FatesInterfaceMod, only : fates_maxElementsPerSite use FatesSizeAgeTypeIndicesMod, only : get_sizeage_class_index - + use FatesHydraulicsMemMod, only : nshell + use FatesHydraulicsMemMod, only : n_hypool_ag + use FatesHydraulicsMemMod, only : n_hypool_troot + use FatesHydraulicsMemMod, only : nlevsoi_hyd_max use PRTGenericMod, only : prt_global - + use EDCohortDynamicsMod, only : nan_cohort + use EDCohortDynamicsMod, only : zero_cohort + use EDCohortDynamicsMod, only : InitPRTCohort + use FatesPlantHydraulicsMod, only : InitHydrCohort ! CIME GLOBALS use shr_log_mod , only : errMsg => shr_log_errMsg @@ -134,6 +142,19 @@ module FatesRestartInterfaceMod integer, private :: ir_prt_base ! Base index for all PRT variables + ! Hydraulic indices + integer, private :: ir_hydro_th_ag_covec + integer, private :: ir_hydro_th_troot_covec + integer, private :: ir_hydro_th_aroot_covec + integer, private :: ir_hydro_liqvol_shell_si + integer, private :: ir_hydro_err_growturn_aroot_covec + integer, private :: ir_hydro_err_growturn_ag_covec + integer, private :: ir_hydro_err_growturn_troot_covec + integer, private :: ir_hydro_recruit_si + integer, private :: ir_hydro_dead_si + integer, private :: ir_hydro_growturn_err_si + integer, private :: ir_hydro_pheno_err_si + integer, private :: ir_hydro_hydro_err_si ! The number of variable dim/kind types we have defined (static) integer, parameter :: fates_restart_num_dimensions = 2 !(cohort,column) @@ -208,6 +229,9 @@ module FatesRestartInterfaceMod procedure, private :: define_restart_vars procedure, private :: set_restart_var procedure, private :: DefinePRTRestartVars + procedure, private :: GetCohortRealVector + procedure, private :: SetCohortRealVector + procedure, private :: RegisterCohortVector end type fates_restart_interface_type @@ -823,6 +847,91 @@ subroutine define_restart_vars(this, initialize_variables) long_name='are of the ED patch', units='m2', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_area_pa ) + + ! Only register hydraulics restart variables if it is turned on! + + if(hlm_use_planthydro==itrue) then + + if ( fates_maxElementsPerSite < (nshell * nlevsoi_hyd_max) ) then + write(fates_log(), *) ' Ftes plant hydraulics needs space to store site-level hydraulics info.' + write(fates_log(), *) ' It uses array spaces typically reserved for cohorts to hold this.' + write(fates_log(), *) ' However, that space defined by fates_maxElementsPerSite must be larger' + write(fates_log(), *) ' than the product of maximum soil layers x rhizosphere shells' + write(fates_log(), *) ' See FatesInterfaceMod.F90 for how this array is set' + write(fates_log(), *) ' fates_maxElementsPerSite = ',fates_maxElementsPerSite + write(fates_log(), *) ' nshell = ',nshell + write(fates_log(), *) ' nlevsoi_hyd_max = ',nlevsoi_hyd_max + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + call this%RegisterCohortVector(symbol_base='fates_hydro_th_ag', vtype=cohort_r8, & + long_name_base='water in aboveground compartments', & + units='kg/plant', veclength=n_hypool_ag, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_th_ag_covec) + + call this%RegisterCohortVector(symbol_base='fates_hydro_th_troot', vtype=cohort_r8, & + long_name_base='water in transporting roots', & + units='kg/plant', veclength=n_hypool_troot, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_th_troot_covec) + + call this%RegisterCohortVector(symbol_base='fates_hydro_th_aroot', vtype=cohort_r8, & + long_name_base='water in absorbing roots', & + units='kg/plant', veclength=nlevsoi_hyd_max, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_th_aroot_covec) + + call this%RegisterCohortVector(symbol_base='fates_hydro_err_aroot', vtype=cohort_r8, & + long_name_base='error in plant-hydro balance in absorbing roots', & + units='kg/plant', veclength=nlevsoi_hyd_max, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_err_growturn_aroot_covec) + + call this%RegisterCohortVector(symbol_base='fates_hydro_err_ag', vtype=cohort_r8, & + long_name_base='error in plant-hydro balance above ground', & + units='kg/plant', veclength=n_hypool_ag, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_err_growturn_ag_covec) + + call this%RegisterCohortVector(symbol_base='fates_hydro_err_troot', vtype=cohort_r8, & + long_name_base='error in plant-hydro balance above ground', & + units='kg/plant', veclength=n_hypool_troot, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_err_growturn_troot_covec) + + ! Site-level volumentric liquid water content (shell x layer) + call this%set_restart_var(vname='fates_hydro_liqvol_shell', vtype=cohort_r8, & + long_name='Volumetric water content of rhizosphere compartments (layerxshell)', & + units='m3/m3', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_liqvol_shell_si ) + + ! Site-level water bound in new recruits + call this%set_restart_var(vname='fates_hydro_recruit_h2o', vtype=site_r8, & + long_name='Site level water mass used for new recruits', & + units='kg', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_recruit_si ) + + ! Site-level water bound in dead plants + call this%set_restart_var(vname='fates_hydro_dead_h2o', vtype=site_r8, & + long_name='Site level water bound in dead plants', & + units='kg', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_dead_si ) + + ! Site-level water balance error due to growth/turnover + call this%set_restart_var(vname='fates_hydro_growturn_err', vtype=site_r8, & + long_name='Site level error for hydraulics due to growth/turnover', & + units='kg', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_growturn_err_si ) + + ! Site-level water balance error due to phenology? + call this%set_restart_var(vname='fates_hydro_pheno_err', vtype=site_r8, & + long_name='Site level error for hydraulics due to phenology', & + units='kg', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_pheno_err_si ) + + ! Site-level water balance error in vegetation + call this%set_restart_var(vname='fates_hydro_hydro_err', vtype=site_r8, & + long_name='Site level error for hydrodynamics', & + units='kg', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_hydro_err_si ) + + end if + ! ! site x time level vars @@ -846,7 +955,7 @@ subroutine define_restart_vars(this, initialize_variables) end subroutine define_restart_vars - ! ===================================================================================== + ! ===================================================================================== subroutine DefinePRTRestartVars(this,initialize_variables,ivar) @@ -977,7 +1086,119 @@ subroutine DefinePRTRestartVars(this,initialize_variables,ivar) end subroutine DefinePRTRestartVars ! ===================================================================================== - + + subroutine RegisterCohortVector(this,symbol_base, vtype, long_name_base, & + units, veclength, flushval, hlms, & + initialize, ivar, index) + + + ! The basic idea here is that instead of saving cohorts with vector data + ! as long arrays in the restart file, we give each index of the vector + ! its own variable. This helps reduce the size of the restart files + ! considerably. + + + use FatesIOVariableKindMod, only : cohort_r8 + + class(fates_restart_interface_type) :: this + character(*),intent(in) :: symbol_base ! Symbol name without position + character(*),intent(in) :: vtype ! String defining variable type + character(*),intent(in) :: long_name_base ! name without position + character(*),intent(in) :: units ! units for this variable + integer,intent(in) :: veclength ! length of the vector + real(r8),intent(in) :: flushval ! Value to flush to + character(*),intent(in) :: hlms ! The HLMs this works in + logical, intent(in) :: initialize ! Is this registering or counting? + integer,intent(inout) :: ivar ! global variable counter + integer,intent(out) :: index ! The variable index for this variable + + ! Local Variables + character(len=4) :: pos_symbol ! vectors need text strings for each position + character(len=128) :: symbol ! symbol name written to file + character(len=256) :: long_name ! long name written to file + integer :: i_pos ! loop counter for discrete position + integer :: dummy_index + + + ! We give each vector its own index that points to the first position + + index = ivar + 1 + + do i_pos = 1, veclength + + ! String describing the physical position of the variable + write(pos_symbol, '(I3.3)') i_pos + + ! The symbol that is written to file + symbol = trim(symbol_base)//'_vec_'//trim(pos_symbol) + + ! The expanded long name of the variable + long_name = trim(long_name_base)//', position:'//trim(pos_symbol) + + call this%set_restart_var(vname=trim(symbol), & + vtype=vtype, & + long_name=trim(long_name), & + units=units, flushval = flushval, & + hlms='CLM:ALM', initialize=initialize, & + ivar=ivar, index = dummy_index ) + + end do + + end subroutine RegisterCohortVector + + ! ===================================================================================== + + subroutine GetCohortRealVector(this, state_vector, len_state_vector, & + variable_index_base, co_global_index) + + ! This subroutine walks through global cohort vector indices + ! and pulls from the different associated restart variables + + class(fates_restart_interface_type) , intent(inout) :: this + real(r8),intent(inout) :: state_vector(len_state_vector) + integer,intent(in) :: len_state_vector + integer,intent(in) :: variable_index_base + integer,intent(in) :: co_global_index + + integer :: i_pos ! vector position loop index + integer :: ir_pos_var ! global variable index + + ir_pos_var = variable_index_base + do i_pos = 1, len_state_vector + state_vector(i_pos) = this%rvars(ir_pos_var)%r81d(co_global_index) + ir_pos_var = ir_pos_var + 1 + end do + return + end subroutine GetCohortRealVector + + ! ===================================================================================== + + subroutine SetCohortRealVector(this, state_vector, len_state_vector, & + variable_index_base, co_global_index) + + ! This subroutine walks through global cohort vector indices + ! and pushes into the restart arrays the different associated restart variables + + class(fates_restart_interface_type) , intent(inout) :: this + real(r8),intent(in) :: state_vector(len_state_vector) + integer,intent(in) :: len_state_vector + integer,intent(in) :: variable_index_base + integer,intent(in) :: co_global_index + + integer :: i_pos ! vector position loop index + integer :: ir_pos_var ! global variable index + + ir_pos_var = variable_index_base + do i_pos = 1, len_state_vector + this%rvars(ir_pos_var)%r81d(co_global_index) = state_vector(i_pos) + ir_pos_var = ir_pos_var + 1 + end do + return + end subroutine SetCohortRealVector + + + ! ===================================================================================== + subroutine set_restart_var(this,vname,vtype,long_name,units,flushval, & hlms,initialize,ivar,index) @@ -1066,7 +1287,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: io_idx_pa_cwd ! each cwd class within each patch (pa_cwd) integer :: io_idx_pa_ib ! each SW band (vis/ir) per patch (pa_ib) integer :: io_idx_si_wmem ! each water memory class within each site - + integer :: io_idx_si_lyr_shell ! site - layer x shell index + ! Some counters (for checking mostly) integer :: totalcohorts ! total cohort count on this thread (diagnostic) integer :: patchespersite ! number of patches per site @@ -1084,7 +1306,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) associate( rio_npatch_si => this%rvars(ir_npatch_si)%int1d, & - rio_old_stock_si => this%rvars(ir_oldstock_si)%r81d, & + rio_old_stock_si => this%rvars(ir_oldstock_si)%r81d, & rio_cd_status_si => this%rvars(ir_cd_status_si)%r81d, & rio_dd_status_si => this%rvars(ir_dd_status_si)%r81d, & rio_nchill_days_si => this%rvars(ir_nchill_days_si)%r81d, & @@ -1150,7 +1372,14 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_livegrass_pa => this%rvars(ir_livegrass_pa)%r81d, & rio_age_pa => this%rvars(ir_age_pa)%r81d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & - rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d ) + rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d, & + rio_hydro_liqvol_shell_si => this%rvars(ir_hydro_liqvol_shell_si)%r81d, & + rio_hydro_recruit_si => this%rvars(ir_hydro_recruit_si)%r81d, & + rio_hydro_dead_si => this%rvars(ir_hydro_dead_si)%r81d, & + rio_hydro_growturn_err_si => this%rvars(ir_hydro_growturn_err_si)%r81d, & + rio_hydro_pheno_err_si => this%rvars(ir_hydro_pheno_err_si)%r81d, & + rio_hydro_hydro_err_si => this%rvars(ir_hydro_hydro_err_si)%r81d) + totalCohorts = 0 @@ -1175,6 +1404,9 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_pa_cwd = io_idx_co_1st io_idx_pa_ib = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st + + ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell + io_idx_si_lyr_shell = io_idx_co_1st ! write seed_bank info(site-level, but PFT-resolved) do i = 1,numpft @@ -1240,6 +1472,31 @@ subroutine set_restart_vectors(this,nc,nsites,sites) end do end do + + if(hlm_use_planthydro==itrue)then + + ! Load the water contents + call this%SetCohortRealVector(ccohort%co_hydr%th_ag,n_hypool_ag, & + ir_hydro_th_ag_covec,io_idx_co) + call this%SetCohortRealVector(ccohort%co_hydr%th_troot,n_hypool_troot, & + ir_hydro_th_troot_covec,io_idx_co) + call this%SetCohortRealVector(ccohort%co_hydr%th_aroot,sites(s)%si_hydr%nlevsoi_hyd, & + ir_hydro_th_aroot_covec,io_idx_co) + + ! Load the error terms + call this%setCohortRealVector(ccohort%co_hydr%errh2o_growturn_aroot, & + sites(s)%si_hydr%nlevsoi_hyd, & + ir_hydro_err_growturn_aroot_covec,io_idx_co) + + call this%setCohortRealVector(ccohort%co_hydr%errh2o_growturn_troot, & + n_hypool_troot, & + ir_hydro_err_growturn_troot_covec,io_idx_co) + + call this%setCohortRealVector(ccohort%co_hydr%errh2o_growturn_ag, & + n_hypool_ag, & + ir_hydro_err_growturn_ag_covec,io_idx_co) + + end if rio_canopy_layer_co(io_idx_co) = ccohort%canopy_layer rio_canopy_layer_yesterday_co(io_idx_co) = ccohort%canopy_layer_yesterday @@ -1391,6 +1648,29 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_watermem_siwm( io_idx_si_wmem ) = sites(s)%water_memory(i) io_idx_si_wmem = io_idx_si_wmem + 1 end do + + ! ----------------------------------------------------------------------------- + ! Set site-level hydraulics arrays + ! ----------------------------------------------------------------------------- + + if(hlm_use_planthydro==itrue)then + + rio_hydro_recruit_si(io_idx_si) = sites(s)%si_hydr%h2oveg_recruit + rio_hydro_dead_si(io_idx_si) = sites(s)%si_hydr%h2oveg_dead + rio_hydro_growturn_err_si(io_idx_si) = sites(s)%si_hydr%h2oveg_growturn_err + rio_hydro_pheno_err_si(io_idx_si) = sites(s)%si_hydr%h2oveg_pheno_err + rio_hydro_hydro_err_si(io_idx_si) = sites(s)%si_hydr%h2oveg_hydro_err + + ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell + do i = 1, sites(s)%si_hydr%nlevsoi_hyd + ! Loop shells + do k = 1, nshell + rio_hydro_liqvol_shell_si(io_idx_si_lyr_shell) = & + sites(s)%si_hydr%h2osoi_liqvol_shell(i,k) + io_idx_si_lyr_shell = io_idx_si_lyr_shell + 1 + end do + end do + end if enddo @@ -1412,12 +1692,14 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) ! subroutine is called prior to the transfer of the restart vectors into the ! linked-list state structure. ! --------------------------------------------------------------------------------- + use EDTypesMod, only : ed_site_type use EDTypesMod, only : ed_cohort_type use EDTypesMod, only : ed_patch_type use EDTypesMod, only : ncwd use EDTypesMod, only : maxSWb use FatesInterfaceMod, only : fates_maxElementsPerPatch + use EDTypesMod, only : maxpft use EDTypesMod, only : area use EDPatchDynamicsMod, only : zero_patch @@ -1439,22 +1721,18 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) ! local variables type(ed_patch_type) , pointer :: newp - type(ed_cohort_type), allocatable :: temp_cohort + type(ed_cohort_type), pointer :: new_cohort + type(ed_cohort_type), pointer :: prev_cohort real(r8) :: cwd_ag_local(ncwd) real(r8) :: cwd_bg_local(ncwd) real(r8) :: leaf_litter_local(maxpft) real(r8) :: root_litter_local(maxpft) real(r8) :: patch_age integer :: cohortstatus - integer :: s ! site index + integer :: s ! site index integer :: idx_pa ! local patch index integer :: io_idx_si ! global site index in IO vector integer :: io_idx_co_1st ! global cohort index in IO vector - real(r8) :: b_dead ! dummy structural biomass (kgC) - real(r8) :: b_store ! dummy storage carbon (kgC) - real(r8) :: b_leaf ! leaf biomass dummy var (kgC) - real(r8) :: b_fineroot ! fineroot dummy var (kgC) - real(r8) :: b_sapwood ! sapwood dummy var (kgC) real(r8) :: site_spread ! site sprea dummy var (0-1) integer :: fto integer :: ft @@ -1485,7 +1763,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) call init_site_vars( sites(s) ) call zero_site( sites(s) ) - + ! ! set a few items that are necessary on restart for ED but not on the ! restart file @@ -1520,42 +1798,49 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) ! give this patch a unique patch number newp%patchno = idx_pa - - do fto = 1, rio_ncohort_pa( io_idx_co_1st ) - allocate(temp_cohort) - - temp_cohort%n = 700.0_r8 - - temp_cohort%laimemory = 0.0_r8 - temp_cohort%canopy_trim = 1.0_r8 - temp_cohort%canopy_layer = 1.0_r8 - temp_cohort%canopy_layer_yesterday = 1.0_r8 - temp_cohort%pft = 1 ! Give it a nominal PFT value for allocation + ! Iterate over the number of cohorts + ! the file says are associated with this patch + ! we are just allocating space here, so we do + ! a simple list filling routine + + newp%tallest => null() + newp%shortest => null() + prev_cohort => null() - cohortstatus = 2 ! status of 2 means leaves are out (dummy var) + do fto = 1, rio_ncohort_pa( io_idx_co_1st ) - temp_cohort%hite = 1.25_r8 - ! Solve for diameter from height - call h2d_allom(temp_cohort%hite,temp_cohort%pft,temp_cohort%dbh) + allocate(new_cohort) + call nan_cohort(new_cohort) + call zero_cohort(new_cohort) + new_cohort%patchptr => newp - if (debug) then - write(fates_log(),*) 'EDRestVectorMod.F90::createPatchCohortStructure call create_cohort ' + ! If this is the first in the list, it is tallest + if (.not.associated(newp%tallest)) then + newp%tallest => new_cohort + endif + + ! Every cohort's taller is the one that came before + ! (unless it is first) + if(associated(prev_cohort)) then + new_cohort%taller => prev_cohort + prev_cohort%shorter => new_cohort end if - - b_dead = 0.0_r8 - b_store = 0.0_r8 - b_leaf = 0.0_r8 - b_fineroot = 0.0_r8 - b_sapwood = 0.0_r8 - site_spread = 0.5_r8 - call create_cohort(sites(s),newp, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & - b_leaf, b_fineroot, b_sapwood, b_dead, b_store, & - temp_cohort%laimemory, cohortstatus,recruitstatus, temp_cohort%canopy_trim, newp%NCL_p, & - site_spread, bc_in(s)) - deallocate(temp_cohort) + ! Ever cohort added takes over as shortest + newp%shortest => new_cohort + + ! Initialize the PRT environment (allocate/choose hypothesis only) + call InitPRTCohort(new_cohort) + + ! Allocate hydraulics arrays + if( hlm_use_planthydro.eq.itrue ) then + call InitHydrCohort(sites(s),new_cohort) + end if + + ! Update the previous + prev_cohort => new_cohort enddo ! ends loop over fto @@ -1649,6 +1934,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: io_idx_pa_cwd ! each cwd class within each patch (pa_cwd) integer :: io_idx_pa_ib ! each SW radiation band per patch (pa_ib) integer :: io_idx_si_wmem ! each water memory class within each site + integer :: io_idx_si_lyr_shell ! site - layer x shell index ! Some counters (for checking mostly) integer :: totalcohorts ! total cohort count on this thread (diagnostic) @@ -1726,7 +2012,13 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_livegrass_pa => this%rvars(ir_livegrass_pa)%r81d, & rio_age_pa => this%rvars(ir_age_pa)%r81d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & - rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d ) + rio_watermem_siwm => this%rvars(ir_watermem_siwm)%r81d, & + rio_hydro_liqvol_shell_si => this%rvars(ir_hydro_liqvol_shell_si)%r81d, & + rio_hydro_recruit_si => this%rvars(ir_hydro_recruit_si)%r81d, & + rio_hydro_dead_si => this%rvars(ir_hydro_dead_si)%r81d, & + rio_hydro_growturn_err_si => this%rvars(ir_hydro_growturn_err_si)%r81d, & + rio_hydro_pheno_err_si => this%rvars(ir_hydro_pheno_err_si)%r81d, & + rio_hydro_hydro_err_si => this%rvars(ir_hydro_hydro_err_si)%r81d) totalcohorts = 0 @@ -1740,7 +2032,11 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_pa_cwd = io_idx_co_1st io_idx_pa_ib = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st - + + ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell + io_idx_si_lyr_shell = io_idx_co_1st + + ! read seed_bank info(site-level, but PFT-resolved) do i = 1,numpft sites(s)%seed_bank(i) = rio_seed_bank_sift(io_idx_co_1st+i-1) @@ -1750,6 +2046,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! Perform a check on the number of patches per site patchespersite = 0 + + cpatch => sites(s)%oldest_patch do while(associated(cpatch)) @@ -1797,6 +2095,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) this%rvars(ir_prt_var)%r81d(io_idx_co) end do end do + ccohort%canopy_layer = rio_canopy_layer_co(io_idx_co) ccohort%canopy_layer_yesterday = rio_canopy_layer_yesterday_co(io_idx_co) @@ -1829,6 +2128,31 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%pft = rio_pft_co(io_idx_co) ccohort%status_coh = rio_status_co(io_idx_co) ccohort%isnew = ( rio_isnew_co(io_idx_co) .eq. new_cohort ) + + ! Initialize Plant Hydraulics + + if(hlm_use_planthydro==itrue)then + + ! Load the water contents + call this%GetCohortRealVector(ccohort%co_hydr%th_ag,n_hypool_ag, & + ir_hydro_th_ag_covec,io_idx_co) + call this%GetCohortRealVector(ccohort%co_hydr%th_troot,n_hypool_troot, & + ir_hydro_th_troot_covec,io_idx_co) + call this%GetCohortRealVector(ccohort%co_hydr%th_aroot,sites(s)%si_hydr%nlevsoi_hyd, & + ir_hydro_th_aroot_covec,io_idx_co) + + call this%GetCohortRealVector(ccohort%co_hydr%errh2o_growturn_aroot, & + sites(s)%si_hydr%nlevsoi_hyd, & + ir_hydro_err_growturn_aroot_covec,io_idx_co) + + call this%GetCohortRealVector(ccohort%co_hydr%errh2o_growturn_troot, & + n_hypool_troot, & + ir_hydro_err_growturn_troot_covec,io_idx_co) + + call this%GetCohortRealVector(ccohort%co_hydr%errh2o_growturn_ag, & + n_hypool_ag, & + ir_hydro_err_growturn_ag_covec,io_idx_co) + end if io_idx_co = io_idx_co + 1 @@ -1924,6 +2248,34 @@ subroutine get_restart_vectors(this, nc, nsites, sites) sites(s)%water_memory(i) = rio_watermem_siwm( io_idx_si_wmem ) io_idx_si_wmem = io_idx_si_wmem + 1 end do + + ! ----------------------------------------------------------------------------- + ! Retrieve site-level hydraulics arrays + ! Note that Hydraulics structures, their allocations, and the length + ! declaration nlevsoi_hyd should be allocated early on when the code first + ! allocates sites (before restart info), and when the soils layer is + ! first known. + ! ----------------------------------------------------------------------------- + + if(hlm_use_planthydro==itrue)then + + sites(s)%si_hydr%h2oveg_recruit = rio_hydro_recruit_si(io_idx_si) + sites(s)%si_hydr%h2oveg_dead = rio_hydro_dead_si(io_idx_si) + sites(s)%si_hydr%h2oveg_growturn_err = rio_hydro_growturn_err_si(io_idx_si) + sites(s)%si_hydr%h2oveg_pheno_err = rio_hydro_pheno_err_si(io_idx_si) + sites(s)%si_hydr%h2oveg_hydro_err = rio_hydro_hydro_err_si(io_idx_si) + + ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell + do i = 1, sites(s)%si_hydr%nlevsoi_hyd + ! Loop shells + do k = 1, nshell + sites(s)%si_hydr%h2osoi_liqvol_shell(i,k) = & + rio_hydro_liqvol_shell_si(io_idx_si_lyr_shell) + io_idx_si_lyr_shell = io_idx_si_lyr_shell + 1 + end do + end do + + end if sites(s)%old_stock = rio_old_stock_si(io_idx_si) sites(s)%status = rio_cd_status_si(io_idx_si) @@ -1952,7 +2304,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) sites(s)%resources_management%trunk_product_site = rio_trunk_product_si(io_idx_si) end do - + if ( debug ) then write(fates_log(),*) 'CVTL total cohorts ',totalCohorts end if @@ -1962,8 +2314,10 @@ end subroutine get_restart_vectors ! ==================================================================================== + + - subroutine update_3dpatch_radiation(this, nc, nsites, sites, bc_out) + subroutine update_3dpatch_radiation(this, nsites, sites, bc_out) ! ------------------------------------------------------------------------- ! This subroutine populates output boundary conditions related to radiation @@ -1977,7 +2331,6 @@ subroutine update_3dpatch_radiation(this, nc, nsites, sites, bc_out) ! !ARGUMENTS: class(fates_restart_interface_type) , intent(inout) :: this - integer , intent(in) :: nc integer , intent(in) :: nsites type(ed_site_type) , intent(inout), target :: sites(nsites) type(bc_out_type) , intent(inout) :: bc_out(nsites)