From aa868207404553aef5a0417a49a0575d00d11237 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 12 Nov 2018 12:32:57 -0800 Subject: [PATCH 01/12] Starting to scratch out variable dependencies for hydro restarts --- main/FatesRestartInterfaceMod.F90 | 133 +++++++++++++++++++++++++++++- 1 file changed, 132 insertions(+), 1 deletion(-) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index a710fd922e..af9ec2a77b 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -846,7 +846,7 @@ subroutine define_restart_vars(this, initialize_variables) end subroutine define_restart_vars - ! ===================================================================================== + ! ===================================================================================== subroutine DefinePRTRestartVars(this,initialize_variables,ivar) @@ -978,6 +978,137 @@ end subroutine DefinePRTRestartVars ! ===================================================================================== + subroutine DefineHydroVecCohortRestartVars(this,initialize_variables,ivar) + + ! ---------------------------------------------------------------------------------- + ! Plant Hydraulics has many variables that are bound to the cohort + ! and are vectors. To prevent absolutely massive restart files + ! we will actually split the restart information of these different + ! vector indices into different saved variables. + ! ----------------------------------------------------------------------------------- + + + use FatesIOVariableKindMod, only : cohort_r8 + + class(fates_restart_interface_type) :: this + logical, intent(in) :: initialize_variables + integer,intent(inout) :: ivar ! global variable counter + + integer :: dummy_out ! dummy index for variable + ! position in global file + integer :: i_var ! loop counter for prt variables + integer :: i_pos ! loop counter for discrete position + + character(len=32) :: symbol_base ! Symbol name without position or flux type + character(len=128) :: name_base ! name without position or flux type + character(len=4) :: pos_symbol + character(len=128) :: symbol + character(len=256) :: long_name + + ! NEED FOR RESTARTS + 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(:) + + + ! CAN BE REMOVED ALTOGETHER + ccohort_hydr%z_node_aroot(1:nlevsoi_hyd) = -bc_in%z_sisl(1:nlevsoi_hyd) + + + do i_var = 1, prt_global%num_vars + + + + + + ! The base symbol name + symbol_base = prt_global%state_descriptor(i_var)%symbol + + ! The long name of the variable + name_base = prt_global%state_descriptor(i_var)%longname + + do i_pos = 1, prt_global%state_descriptor(i_var)%num_pos + + ! String describing the physical position of the variable + write(pos_symbol, '(I3.3)') i_pos + + ! Register the instantaneous state variable "val" + ! ---------------------------------------------------------------------------- + + ! The symbol that is written to file + symbol = trim(symbol_base)//'_val_'//trim(pos_symbol) + + ! The expanded long name of the variable + long_name = trim(name_base)//', state var, position:'//trim(pos_symbol) + + call this%set_restart_var(vname=trim(symbol), & + vtype=cohort_r8, & + long_name=trim(long_name), & + units='kg', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, & + ivar=ivar, index = dummy_out ) + + ! Register the turnover flux variables + ! ---------------------------------------------------------------------------- + + ! The symbol that is written to file + symbol = trim(symbol_base)//'_turn_'//trim(pos_symbol) + + ! The expanded long name of the variable + long_name = trim(name_base)//', turnover, position:'//trim(pos_symbol) + + call this%set_restart_var(vname=trim(symbol), & + vtype=cohort_r8, & + long_name=trim(long_name), & + units='kg', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, & + ivar=ivar, index = dummy_out ) + + + + ! Register the net allocation flux variable + ! ---------------------------------------------------------------------------- + + ! The symbol that is written to file + symbol = trim(symbol_base)//'_net_'//trim(pos_symbol) + + ! The expanded long name of the variable + long_name = trim(name_base)//', net allocation/transp, position:'//trim(pos_symbol) + + call this%set_restart_var(vname=trim(symbol), & + vtype=cohort_r8, & + long_name=trim(long_name), & + units='kg', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, & + ivar=ivar, index = dummy_out ) + + + + ! Register the burn flux variable + ! ---------------------------------------------------------------------------- + ! The symbol that is written to file + symbol = trim(symbol_base)//'_burned_'//trim(pos_symbol) + + ! The expanded long name of the variable + long_name = trim(name_base)//', burned mass:'//trim(pos_symbol) + + call this%set_restart_var(vname=symbol, & + vtype=cohort_r8, & + long_name=trim(long_name), & + units='kg', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, & + ivar=ivar, index = dummy_out ) + + end do + end do + + return + end subroutine DefineHydroRestartVars + + ! ===================================================================================== + + + subroutine set_restart_var(this,vname,vtype,long_name,units,flushval, & hlms,initialize,ivar,index) From 977c6c7e2744c974aad33ebbe81418278cecbbda Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 12 Nov 2018 14:41:31 -0800 Subject: [PATCH 02/12] Spacing out and listing hydro variables that will be added to restarts. --- main/FatesRestartInterfaceMod.F90 | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index af9ec2a77b..c1b2814b0c 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -1005,11 +1005,25 @@ subroutine DefineHydroVecCohortRestartVars(this,initialize_variables,ivar) character(len=128) :: symbol character(len=256) :: long_name - ! NEED FOR RESTARTS + ! From update props 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(:) + ! From update states + ccohort_hydr%th_ag(:) + ccohort_hydr%th_troot(k) + ccohort_hydr%th_aroot(j) + + ! From update rhiz props + 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(:,:) + + + ! From update rhiz states + csite_hydr%h2osoi_liqvol_shell(j,k) + ! CAN BE REMOVED ALTOGETHER ccohort_hydr%z_node_aroot(1:nlevsoi_hyd) = -bc_in%z_sisl(1:nlevsoi_hyd) From b84261ad2151cfdb979e31296f1c59ecbb3ac383 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 13 Nov 2018 12:26:04 -0800 Subject: [PATCH 03/12] First pass at incorporating hydraulics variables into restarts --- main/FatesHydraulicsMemMod.F90 | 6 +- main/FatesRestartInterfaceMod.F90 | 444 ++++++++++++++++++++---------- 2 files changed, 305 insertions(+), 145 deletions(-) diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90 index a1d2fbfd71..8bc0e2f8af 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 c1b2814b0c..19ad6f0ce8 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -14,6 +14,7 @@ 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 FatesSizeAgeTypeIndicesMod, only : get_sizeage_class_index use PRTGenericMod, only : prt_global @@ -134,6 +135,18 @@ module FatesRestartInterfaceMod integer, private :: ir_prt_base ! Base index for all PRT variables + ! Hydraulic indices + integer, private :: ir_hydro_v_ag_covec + integer, private :: ir_hydro_v_troot_covec + integer, private :: ir_hydro_v_aroot_layer_covec + 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_aroot_layer_si + integer, private :: ir_hydro_r_node_shell_si + integer, private :: ir_hydro_v_shell_si + integer, private :: ir_hydro_liqvol_shell_si + ! The number of variable dim/kind types we have defined (static) integer, parameter :: fates_restart_num_dimensions = 2 !(cohort,column) @@ -824,6 +837,65 @@ subroutine define_restart_vars(this, initialize_variables) hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_area_pa ) + if(hlm_use_planthydro==itrue) then + + call RegisterCohortVectors('fates_hydro_v_ag', vtype=cohort_r8, & + long_name='maximum storage volume of hydraulic compartments (above ground)', & + units='m3', n_hypool_ag, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_ag_covec) + + call RegisterCohortVectors('fates_hydro_v_troot', vtype=cohort_r8, & + long_name='maximum storage volume of transporting root compartments', & + units='m3', n_hypool_troot, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_troot_covec) + + call RegisterCohortVectors('fates_hydro_v_aroot_layer', vtype=cohort_r8, & + long_name='maximum storage volume of absorbing roots hydraulic compartments by soil layer', & + units='m3', nlevsoi_hyd_max, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_aroot_layer_covec) + + call RegisterCohortVectors('fates_hydro_th_ag', vtype=cohort_r8, & + long_name='water in aboveground compartments', & + units='kg/plant', n_hypool_ag, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_th_ag_covec) + + call RegisterCohortVectors('fates_hydro_th_troot', vtype=cohort_r8, & + long_name='water in transporting roots', & + units='kg/plant', n_hypool_troot, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_th_troot_covec) + + call RegisterCohortVectors('fates_hydro_th_aroot', vtype=cohort_r8, & + long_name='water in absorbing roots', & + units='kg/plant', nlevsoi_hyd_max, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_th_aroot_covec) + + ! Site-level absorbing root maximum volume + call this%set_restart_var(vname='fates_hydro_l_aroot_layer', vtype=cohort_r8, & + long_name='Total length (across cohorts) of absorbing roots by soil layer', & + units='m', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_l_aroot_layer_si ) + + ! Site-level Nodal Radius of rhizosphere compartments + call this%set_restart_var(vname='fates_hydro_r_node_shell', vtype=cohort_r8, & + long_name='Nodal Radius of rhizosphere compartments (layerxshell)', & + units='m', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_r_node_shell_si ) + + ! Site-level volume of rhizosphere compartments + call this%set_restart_var(vname='fates_hydro_v_shell', vtype=cohort_r8, & + long_name='Volume of rhizosphere compartments (layerxshell)', & + units='m3', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_shell_si ) + + ! 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 ) + + end if + + ! ! site x time level vars ! @@ -976,152 +1048,117 @@ subroutine DefinePRTRestartVars(this,initialize_variables,ivar) return end subroutine DefinePRTRestartVars - ! ===================================================================================== - - subroutine DefineHydroVecCohortRestartVars(this,initialize_variables,ivar) - - ! ---------------------------------------------------------------------------------- - ! Plant Hydraulics has many variables that are bound to the cohort - ! and are vectors. To prevent absolutely massive restart files - ! we will actually split the restart information of these different - ! vector indices into different saved variables. - ! ----------------------------------------------------------------------------------- - - use FatesIOVariableKindMod, only : cohort_r8 + subroutine RegisterCohortVector(this,symbol_base, vtype_in, long_name_base, & + units_in, veclength, flushval, hlms_in, & + initialize=initialize_variables, & + ivar=ivar_inout, index = index_out) - class(fates_restart_interface_type) :: this - logical, intent(in) :: initialize_variables - integer,intent(inout) :: ivar ! global variable counter - - integer :: dummy_out ! dummy index for variable - ! position in global file - integer :: i_var ! loop counter for prt variables - integer :: i_pos ! loop counter for discrete position - - character(len=32) :: symbol_base ! Symbol name without position or flux type - character(len=128) :: name_base ! name without position or flux type - character(len=4) :: pos_symbol - character(len=128) :: symbol - character(len=256) :: long_name - - ! From update props - 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(:) - - ! From update states - ccohort_hydr%th_ag(:) - ccohort_hydr%th_troot(k) - ccohort_hydr%th_aroot(j) - - ! From update rhiz props - 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(:,:) - - - ! From update rhiz states - csite_hydr%h2osoi_liqvol_shell(j,k) - - - ! CAN BE REMOVED ALTOGETHER - ccohort_hydr%z_node_aroot(1:nlevsoi_hyd) = -bc_in%z_sisl(1:nlevsoi_hyd) - - - do i_var = 1, prt_global%num_vars - - - - - - ! The base symbol name - symbol_base = prt_global%state_descriptor(i_var)%symbol - - ! The long name of the variable - name_base = prt_global%state_descriptor(i_var)%longname - - do i_pos = 1, prt_global%state_descriptor(i_var)%num_pos - - ! String describing the physical position of the variable - write(pos_symbol, '(I3.3)') i_pos - - ! Register the instantaneous state variable "val" - ! ---------------------------------------------------------------------------- - - ! The symbol that is written to file - symbol = trim(symbol_base)//'_val_'//trim(pos_symbol) - - ! The expanded long name of the variable - long_name = trim(name_base)//', state var, position:'//trim(pos_symbol) - - call this%set_restart_var(vname=trim(symbol), & - vtype=cohort_r8, & - long_name=trim(long_name), & - units='kg', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, & - ivar=ivar, index = dummy_out ) - - ! Register the turnover flux variables - ! ---------------------------------------------------------------------------- - - ! The symbol that is written to file - symbol = trim(symbol_base)//'_turn_'//trim(pos_symbol) - - ! The expanded long name of the variable - long_name = trim(name_base)//', turnover, position:'//trim(pos_symbol) - - call this%set_restart_var(vname=trim(symbol), & - vtype=cohort_r8, & - long_name=trim(long_name), & - units='kg', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, & - ivar=ivar, index = dummy_out ) - - - - ! Register the net allocation flux variable - ! ---------------------------------------------------------------------------- - - ! The symbol that is written to file - symbol = trim(symbol_base)//'_net_'//trim(pos_symbol) - - ! The expanded long name of the variable - long_name = trim(name_base)//', net allocation/transp, position:'//trim(pos_symbol) - - call this%set_restart_var(vname=trim(symbol), & - vtype=cohort_r8, & - long_name=trim(long_name), & - units='kg', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, & - ivar=ivar, index = dummy_out ) - - - - ! Register the burn flux variable - ! ---------------------------------------------------------------------------- - ! The symbol that is written to file - symbol = trim(symbol_base)//'_burned_'//trim(pos_symbol) - - ! The expanded long name of the variable - long_name = trim(name_base)//', burned mass:'//trim(pos_symbol) - - call this%set_restart_var(vname=symbol, & - vtype=cohort_r8, & - long_name=trim(long_name), & - units='kg', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, & - ivar=ivar, index = dummy_out ) - - end do - end do - - return - end subroutine DefineHydroRestartVars + + ! 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_in ! String defining variable type + character(*),intent(in) :: long_name_base ! name without position + character(*),intent(in) :: units_in ! units for this variable + integer,intent(in) :: veclength ! length of the vector + real(r8),intent(in) :: flushval_in ! Value to flush to + character(*),intent(in) :: hlms_in ! The HLMs this works in + logical, intent(in) :: init_var_in ! Is this registering or counting? + integer,intent(inout) :: ivar_inout ! global variable counter + integer,intent(out) :: index_out ! 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 + + ! We give each vector its own index + + index_out = ivar_inout + 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_in, & + long_name=trim(long_name), & + units=units_in, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, & + ivar=ivar, index = index_out ) + + end do + + end subroutine RegisterCohortVector ! ===================================================================================== + + subroutine RetrieveCohortRealVector(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 RetrieveCohortRealVector + + ! ===================================================================================== + + 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) @@ -1295,7 +1332,12 @@ 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_l_aroot_layer_si => this%rvars(ir_hydro_l_aroot_layer_si)%r81d, & + rio_hydro_r_node_shell_si => this%rvars(ir_hydro_r_node_shell_si)%r81d, & + rio_hydro_v_shell_si => this%rvars(ir_hydro_v_shell_si)%r81d, & + rio_hydro_liqvol_shell_si => this%rvars(ir_hydro_liqvol_shell_si)%r81d + ) totalCohorts = 0 @@ -1320,6 +1362,12 @@ 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 = io_idx_co_1st + io_idx_si_lyr_shell = io_idx_co_1st + + ! write seed_bank info(site-level, but PFT-resolved) do i = 1,numpft @@ -1385,6 +1433,19 @@ subroutine set_restart_vectors(this,nc,nsites,sites) end do end do + + if(hlm_use_planthydro==itrue)then + + ! Load the storage compartment volumes + call this%SetCohortRealVector(cohort%hydr_co%v_ag,n_hypool_ag,ir_hydro_v_ag_covec,io_idx_co) + call this%SetCohortRealVector(cohort%hydr_co%v_troot,n_hypool_troot,ir_hydro_v_troot_covec,io_idx_co) + call this%SetCohortRealVector(cohort%hydr_co%v_aroot_layer,site(s)%si_hydr%nlevsoi_hyd,ir_hydro_v_aroot_layer_covec,io_idx_co) + + ! Load the water contents + call this%SetCohortRealVector(cohort%hydr_co%th_ag,n_hypool_ag,ir_hydro_th_ag_covec,io_idx_co) + call this%SetCohortRealVector(cohort%hydr_co%th_troot,n_hypool_troot,ir_hydro_th_troot_covec,io_idx_co) + call this%SetCohortRealVector(cohort%hydr_co%th_aroot_layer,site(s)%si_hydr%nlevsoi_hyd,ir_hydro_th_aroot_layer_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 @@ -1536,6 +1597,35 @@ 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 + + ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell + do i = 1, site_hydr%nlevsoi_hyd + + ! Loop shells + do k = 1, nshell + + rio_hydro_r_node_shell_si(io_idx_si_lyr_shell) = & + sites(s)%si_hydr%r_node_shell(i,k) + + rio_hydro_v_shell_si(io_idx_si_lyr_shell) = & + sites(s)%si_hydr%v_shell(i,k) + + 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 + + rio_hydro_l_aroot_layer_si(io_idx_si_lyr) = sites(s)%si_hydr%l_aroot_layer(i) + io_idx_si_lyr = io_idx_si_lyr + 1 + end do + end if enddo @@ -1557,6 +1647,7 @@ 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 @@ -1695,6 +1786,9 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) b_fineroot = 0.0_r8 b_sapwood = 0.0_r8 site_spread = 0.5_r8 + + ! Hydraulics - if turned on, the hydraulics arrays are being allocated in create_cohort as well. + 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, & @@ -1871,7 +1965,12 @@ 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_l_aroot_layer_si => this%rvars(ir_hydro_l_aroot_layer_si)%r81d, & + rio_hydro_r_node_shell_si => this%rvars(ir_hydro_r_node_shell_si)%r81d, & + rio_hydro_v_shell_si => this%rvars(ir_hydro_v_shell_si)%r81d, & + rio_hydro_liqvol_shell_si => this%rvars(ir_hydro_liqvol_shell_si)%r81d + ) totalcohorts = 0 @@ -1895,6 +1994,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)) @@ -1942,6 +2043,28 @@ subroutine get_restart_vectors(this, nc, nsites, sites) this%rvars(ir_prt_var)%r81d(io_idx_co) end do end do + + + if(hlm_use_planthydro==itrue)then + + ! Load the storage compartment volumes + call this%RetrieveCohortRealVector(cohort%hydr_co%v_ag,n_hypool_ag, & + ir_hydro_v_ag_covec,io_idx_co) + call this%RetrieveCohortRealVector(cohort%hydr_co%v_troot,n_hypool_troot, & + ir_hydro_v_troot_covec,io_idx_co) + call this%RetrieveCohortRealVector(cohort%hydr_co%v_aroot_layer,site(s)%si_hydr%nlevsoi_hyd, & + ir_hydro_v_aroot_layer_covec,io_idx_co) + + ! Load the water contents + call this%RetrieveCohortRealVector(cohort%hydr_co%th_ag,n_hypool_ag, & + ir_hydro_th_ag_covec,io_idx_co) + call this%RetrieveCohortRealVector(cohort%hydr_co%th_troot,n_hypool_troot, & + ir_hydro_th_troot_covec,io_idx_co) + call this%RetrieveCohortRealVector(cohort%hydr_co%th_aroot_layer,site(s)%si_hydr%nlevsoi_hyd, & + ir_hydro_th_aroot_layer_covec,io_idx_co) + end if + + ccohort%canopy_layer = rio_canopy_layer_co(io_idx_co) ccohort%canopy_layer_yesterday = rio_canopy_layer_yesterday_co(io_idx_co) @@ -2069,6 +2192,43 @@ 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 + + ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell + do i = 1, site(s)%si_hydr%nlevsoi_hyd + + ! Loop shells + do k = 1, nshell + + sites(s)%si_hydr%r_node_shell(i,k) = & + rio_hydro_r_node_shell_si(io_idx_si_lyr_shell) + + sites(s)%si_hydr%v_shell(i,k) = & + rio_hydro_v_shell_si(io_idx_si_lyr_shell) + + 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 + + sites(s)%si_hydr%l_aroot_layer(i) = & + rio_hydro_l_aroot_layer_si(io_idx_si_lyr) + + io_idx_si_lyr = io_idx_si_lyr + 1 + 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) From 37291fc469da6f9ab5b28ccade507cbf6a83b4b4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 13 Nov 2018 14:40:21 -0800 Subject: [PATCH 04/12] syntax and bug fixes for hydro restarts, code compiles, un-tested --- main/FatesRestartInterfaceMod.F90 | 164 ++++++++++++++++++------------ 1 file changed, 101 insertions(+), 63 deletions(-) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 19ad6f0ce8..5e59e55137 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -15,8 +15,12 @@ module FatesRestartInterfaceMod 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 @@ -142,7 +146,7 @@ module FatesRestartInterfaceMod 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_aroot_layer_si + integer, private :: ir_hydro_l_aroot_layer_si integer, private :: ir_hydro_r_node_shell_si integer, private :: ir_hydro_v_shell_si integer, private :: ir_hydro_liqvol_shell_si @@ -221,6 +225,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 @@ -836,37 +843,52 @@ 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 - call RegisterCohortVectors('fates_hydro_v_ag', vtype=cohort_r8, & - long_name='maximum storage volume of hydraulic compartments (above ground)', & - units='m3', n_hypool_ag, flushval = flushzero, & + 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_v_ag', vtype=cohort_r8, & + long_name_base='maximum storage volume of hydraulic compartments (above ground)', & + units='m3', veclength=n_hypool_ag, flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_ag_covec) - call RegisterCohortVectors('fates_hydro_v_troot', vtype=cohort_r8, & - long_name='maximum storage volume of transporting root compartments', & - units='m3', n_hypool_troot, flushval = flushzero, & + call this%RegisterCohortVector(symbol_base='fates_hydro_v_troot', vtype=cohort_r8, & + long_name_base='maximum storage volume of transporting root compartments', & + units='m3', veclength=n_hypool_troot, flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_troot_covec) - call RegisterCohortVectors('fates_hydro_v_aroot_layer', vtype=cohort_r8, & - long_name='maximum storage volume of absorbing roots hydraulic compartments by soil layer', & - units='m3', nlevsoi_hyd_max, flushval = flushzero, & + call this%RegisterCohortVector(symbol_base='fates_hydro_v_aroot_layer', vtype=cohort_r8, & + long_name_base='maximum storage volume of absorbing roots hydraulic compartments by layer', & + units='m3', veclength=nlevsoi_hyd_max, flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_aroot_layer_covec) - call RegisterCohortVectors('fates_hydro_th_ag', vtype=cohort_r8, & - long_name='water in aboveground compartments', & - units='kg/plant', n_hypool_ag, flushval = flushzero, & + 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 RegisterCohortVectors('fates_hydro_th_troot', vtype=cohort_r8, & - long_name='water in transporting roots', & - units='kg/plant', n_hypool_troot, flushval = flushzero, & + 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 RegisterCohortVectors('fates_hydro_th_aroot', vtype=cohort_r8, & - long_name='water in absorbing roots', & - units='kg/plant', nlevsoi_hyd_max, flushval = flushzero, & + 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) ! Site-level absorbing root maximum volume @@ -1048,11 +1070,11 @@ subroutine DefinePRTRestartVars(this,initialize_variables,ivar) return end subroutine DefinePRTRestartVars + ! ===================================================================================== - subroutine RegisterCohortVector(this,symbol_base, vtype_in, long_name_base, & - units_in, veclength, flushval, hlms_in, & - initialize=initialize_variables, & - ivar=ivar_inout, index = index_out) + 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 @@ -1065,25 +1087,27 @@ subroutine RegisterCohortVector(this,symbol_base, vtype_in, long_name_base, & class(fates_restart_interface_type) :: this character(*),intent(in) :: symbol_base ! Symbol name without position - character(*),intent(in) :: vtype_in ! String defining variable type + character(*),intent(in) :: vtype ! String defining variable type character(*),intent(in) :: long_name_base ! name without position - character(*),intent(in) :: units_in ! units for this variable + character(*),intent(in) :: units ! units for this variable integer,intent(in) :: veclength ! length of the vector - real(r8),intent(in) :: flushval_in ! Value to flush to - character(*),intent(in) :: hlms_in ! The HLMs this works in - logical, intent(in) :: init_var_in ! Is this registering or counting? - integer,intent(inout) :: ivar_inout ! global variable counter - integer,intent(out) :: index_out ! The variable index for this variable + 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 + + ! We give each vector its own index that points to the first position - index_out = ivar_inout + 1 + index = ivar + 1 do i_pos = 1, veclength @@ -1097,11 +1121,11 @@ subroutine RegisterCohortVector(this,symbol_base, vtype_in, long_name_base, & long_name = trim(long_name_base)//', position:'//trim(pos_symbol) call this%set_restart_var(vname=trim(symbol), & - vtype=vtype_in, & + vtype=vtype, & long_name=trim(long_name), & - units=units_in, flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, & - ivar=ivar, index = index_out ) + units=units, flushval = flushval, & + hlms='CLM:ALM', initialize=initialize, & + ivar=ivar, index = dummy_index ) end do @@ -1109,8 +1133,8 @@ end subroutine RegisterCohortVector ! ===================================================================================== - subroutine RetrieveCohortRealVector(this, state_vector, len_state_vector, & - variable_index_base, co_global_index) + 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 @@ -1130,7 +1154,7 @@ subroutine RetrieveCohortRealVector(this, state_vector, len_state_vector, & ir_pos_var = ir_pos_var + 1 end do return - end subroutine RetrieveCohortRealVector + end subroutine GetCohortRealVector ! ===================================================================================== @@ -1248,7 +1272,9 @@ 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 ! site - layer x index + 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 @@ -1266,7 +1292,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, & @@ -1336,8 +1362,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_hydro_l_aroot_layer_si => this%rvars(ir_hydro_l_aroot_layer_si)%r81d, & rio_hydro_r_node_shell_si => this%rvars(ir_hydro_r_node_shell_si)%r81d, & rio_hydro_v_shell_si => this%rvars(ir_hydro_v_shell_si)%r81d, & - rio_hydro_liqvol_shell_si => this%rvars(ir_hydro_liqvol_shell_si)%r81d - ) + rio_hydro_liqvol_shell_si => this%rvars(ir_hydro_liqvol_shell_si)%r81d ) totalCohorts = 0 @@ -1437,14 +1462,20 @@ subroutine set_restart_vectors(this,nc,nsites,sites) if(hlm_use_planthydro==itrue)then ! Load the storage compartment volumes - call this%SetCohortRealVector(cohort%hydr_co%v_ag,n_hypool_ag,ir_hydro_v_ag_covec,io_idx_co) - call this%SetCohortRealVector(cohort%hydr_co%v_troot,n_hypool_troot,ir_hydro_v_troot_covec,io_idx_co) - call this%SetCohortRealVector(cohort%hydr_co%v_aroot_layer,site(s)%si_hydr%nlevsoi_hyd,ir_hydro_v_aroot_layer_covec,io_idx_co) + call this%SetCohortRealVector(ccohort%co_hydr%v_ag,n_hypool_ag, & + ir_hydro_v_ag_covec,io_idx_co) + call this%SetCohortRealVector(ccohort%co_hydr%v_troot,n_hypool_troot, & + ir_hydro_v_troot_covec,io_idx_co) + call this%SetCohortRealVector(ccohort%co_hydr%v_aroot_layer,sites(s)%si_hydr%nlevsoi_hyd, & + ir_hydro_v_aroot_layer_covec,io_idx_co) ! Load the water contents - call this%SetCohortRealVector(cohort%hydr_co%th_ag,n_hypool_ag,ir_hydro_th_ag_covec,io_idx_co) - call this%SetCohortRealVector(cohort%hydr_co%th_troot,n_hypool_troot,ir_hydro_th_troot_covec,io_idx_co) - call this%SetCohortRealVector(cohort%hydr_co%th_aroot_layer,site(s)%si_hydr%nlevsoi_hyd,ir_hydro_th_aroot_layer_covec,io_idx_co) + 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) end if rio_canopy_layer_co(io_idx_co) = ccohort%canopy_layer @@ -1605,7 +1636,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) if(hlm_use_planthydro==itrue)then ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell - do i = 1, site_hydr%nlevsoi_hyd + do i = 1, sites(s)%si_hydr%nlevsoi_hyd ! Loop shells do k = 1, nshell @@ -1654,6 +1685,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) 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 @@ -1888,6 +1920,8 @@ 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 ! site - layer x index + 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) @@ -1969,8 +2003,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_hydro_l_aroot_layer_si => this%rvars(ir_hydro_l_aroot_layer_si)%r81d, & rio_hydro_r_node_shell_si => this%rvars(ir_hydro_r_node_shell_si)%r81d, & rio_hydro_v_shell_si => this%rvars(ir_hydro_v_shell_si)%r81d, & - rio_hydro_liqvol_shell_si => this%rvars(ir_hydro_liqvol_shell_si)%r81d - ) + rio_hydro_liqvol_shell_si => this%rvars(ir_hydro_liqvol_shell_si)%r81d ) totalcohorts = 0 @@ -1984,7 +2017,12 @@ 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 = io_idx_co_1st + 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) @@ -2048,20 +2086,20 @@ subroutine get_restart_vectors(this, nc, nsites, sites) if(hlm_use_planthydro==itrue)then ! Load the storage compartment volumes - call this%RetrieveCohortRealVector(cohort%hydr_co%v_ag,n_hypool_ag, & + call this%GetCohortRealVector(ccohort%co_hydr%v_ag,n_hypool_ag, & ir_hydro_v_ag_covec,io_idx_co) - call this%RetrieveCohortRealVector(cohort%hydr_co%v_troot,n_hypool_troot, & + call this%GetCohortRealVector(ccohort%co_hydr%v_troot,n_hypool_troot, & ir_hydro_v_troot_covec,io_idx_co) - call this%RetrieveCohortRealVector(cohort%hydr_co%v_aroot_layer,site(s)%si_hydr%nlevsoi_hyd, & + call this%GetCohortRealVector(ccohort%co_hydr%v_aroot_layer,sites(s)%si_hydr%nlevsoi_hyd, & ir_hydro_v_aroot_layer_covec,io_idx_co) ! Load the water contents - call this%RetrieveCohortRealVector(cohort%hydr_co%th_ag,n_hypool_ag, & + call this%GetCohortRealVector(ccohort%co_hydr%th_ag,n_hypool_ag, & ir_hydro_th_ag_covec,io_idx_co) - call this%RetrieveCohortRealVector(cohort%hydr_co%th_troot,n_hypool_troot, & + call this%GetCohortRealVector(ccohort%co_hydr%th_troot,n_hypool_troot, & ir_hydro_th_troot_covec,io_idx_co) - call this%RetrieveCohortRealVector(cohort%hydr_co%th_aroot_layer,site(s)%si_hydr%nlevsoi_hyd, & - ir_hydro_th_aroot_layer_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) end if @@ -2204,7 +2242,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) if(hlm_use_planthydro==itrue)then ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell - do i = 1, site(s)%si_hydr%nlevsoi_hyd + do i = 1, sites(s)%si_hydr%nlevsoi_hyd ! Loop shells do k = 1, nshell From ffd5e4f4663413f30d9262211445b6212c51bdce Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 14 Nov 2018 12:35:58 -0800 Subject: [PATCH 05/12] SPlitting up hydraulic props routine --- biogeochem/EDCohortDynamicsMod.F90 | 4 +- biogeochem/FatesAllometryMod.F90 | 28 +++ biogeophys/FatesPlantHydraulicsMod.F90 | 274 ++++++++++++++++--------- 3 files changed, 203 insertions(+), 103 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 939360388a..f2fdf22d14 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -243,8 +243,8 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine 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) +! call updateSizeDepTreeHydProps(CurrentSite,new_cohort, bc_in) +! call initTreeHydStates(CurrentSite,new_cohort, bc_in) if(recruitstatus==1)then new_cohort%co_hydr%is_newly_recruited = .true. endif diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 81881a190a..d1a0eb8fa9 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(hite,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(cCohort%hite,0.1_r8) + + return + end subroutine CrownDepth + + + ! ============================================================================= ! Specific diameter to crown area allometries ! ============================================================================= diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 4a1404d292..17523e7947 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -173,15 +173,22 @@ end subroutine Hydraulics_Drive ! ===================================================================================== 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 @@ -257,14 +264,118 @@ subroutine initTreeHydStates(site_p, cc_p, bc_in) end subroutine initTreeHydStates + + ! ===================================================================================== + + + subroutine UpdateTreeHydrNodes(ccohort_hydr,pft,plant_height,nlevsoi_hyd) + + ! -------------------------------------------------------------------------------- + ! 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) :: pft ! plant functional type index + real(r8), intent(in) :: plant_height ! [m] + integer,intent(in) :: nlevsoi_hyd ! number of soil hydro layers + + ! 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) :: 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 + + 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 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 + ! 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 FatesConstantsMod , only : pi_const @@ -272,22 +383,19 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) ! ! !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 + type(ed_cohort_type) , intent(inout) :: ccohort ! current cohort pointer + type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions ! !LOCAL VARIABLES: - type(ed_cohort_type), pointer :: cCohort + type(ed_patch_type), pointer :: cPatch - integer :: i,j,k,FT ! indices + 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] @@ -326,142 +434,107 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) 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 + 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] integer :: nlevsoi_hyd ! Number of soil hydraulic layers integer :: nlevsoil ! Number of total soil layers type(ed_cohort_hydr_type), pointer :: ccohort_hydr !----------------------------------------------------------------------- + nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd nlevsoil = bc_in%nlevsoil - cCohort => cc_p + 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) + cPatch => ccohort%patchptr + ft = ccohort%pft + + ! This updates all of the z_node positions + call UpdateTreeHydrNodes(ccohort_hydr,ft,ccohort%hite,nlevsoi_hyd) - 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) + ! 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(:) + + + ! Pre-process biomass pools + + 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 + if(leaf_c > 0._r8) then + + 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_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) - sla = EDPftvarcon_inst%slatop(FT) * cm2_per_m2 ! m2/gC * cm2/m2 -> cm2/gC + sla = EDPftvarcon_inst%slatop(ft) * cm2_per_m2 ! m2/gC * cm2/m2 -> cm2/gC 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 + ccohort_hydr%v_ag(1:n_hypool_leaf) = v_canopy / real(n_hypool_leaf,r8) - ! 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 + 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 - call bsap_allom(cCohort%dbh,cCohort%pft,cCohort%canopy_trim,a_sapwood_target,bsw_target) + call bsap_allom(ccohort%dbh,ccohort%pft,ccohort%canopy_trim,a_sapwood_target,bsw_target) a_sapwood = a_sapwood_target ! 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 + ! 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 + !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 - ! 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 + !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) + 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 - ! ABSORBING ROOT DEPTH, LENGTH & VOLUME - ccohort_hydr%z_node_aroot(1:nlevsoi_hyd) = -bc_in%z_sisl(1:nlevsoi_hyd) + - 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) + 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 + ! 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)) @@ -473,9 +546,7 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) 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 + ! MAXIMUM (SIZE-DEPENDENT) HYDRAULIC CONDUCTANCES ! first estimate cumulative (petiole to node k) conductances without taper as well as the chi taper function @@ -486,12 +557,13 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) 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 + 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 + call CrownDepth(ccohort%hite,crown_depth) + 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 @@ -528,10 +600,10 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) 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 + 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,:) + ccohort_hydr%kmax_treebg_layer(:) = ccohort_hydr%kmax_treebg_tot * cPatch%rootfr_ft(ft,:) else do j=1,nlevsoi_hyd if(j == 1) then From e6679a850c03ae92849d08472fd60a5415520d53 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 14 Nov 2018 16:07:42 -0800 Subject: [PATCH 06/12] Re-arranging calls to initialize hydro cohorts --- biogeochem/EDCohortDynamicsMod.F90 | 33 +- biogeochem/EDPhysiologyMod.F90 | 5 +- biogeophys/FatesPlantHydraulicsMod.F90 | 488 ++++++++++++++----------- 3 files changed, 311 insertions(+), 215 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index f2fdf22d14..a24d75a773 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -98,6 +98,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: ! @@ -126,6 +133,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine 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_hydr_type), pointer :: ccohort_hydr integer :: tnull,snull ! are the tallest and shortest cohorts allocate !---------------------------------------------------------------------- @@ -242,12 +250,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) + + ccohort_hydr => ccohort%co_hydr + nlevsoi_hydr = currentSite%si_hydr%nlevsoi_hyd + + ! This allocates array spaces + call InitHydrCohort(currentSite,new_cohort) + + ! This calculates node heights + call UpdateTreeHydrNodes(ccohort_hydr,ft,new_cohort%hite,nlevsoi_hydr) + + ! This calculates volumes, lengths and max conductances + call UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hydr,bc_in) + + ! Since this is a newly initialized plant, we set the previous compartment-size + ! equal to the ones we just calculated. + call SavePreviousCompartmentVolumes(ccohort_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/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 17523e7947..41ec79b178 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -83,7 +83,7 @@ module FatesPlantHydraulicsMod ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg use shr_infnan_mod , only : isnan => shr_infnan_isnan - + implicit none @@ -115,6 +115,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 @@ -366,257 +372,319 @@ subroutine UpdateTreeHydrNodes(ccohort_hydr,pft,plant_height,nlevsoi_hyd) return end subroutine UpdateTreeHydrNodes - ! ===================================================================================== - + 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 - ! + ! !USES: use FatesConstantsMod , only : pi_const use shr_sys_mod , only : shr_sys_abort - ! - ! !ARGUMENTS: + + ! 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 - ! !LOCAL VARIABLES: - - - 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) :: 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 ! 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] - integer :: nlevsoi_hyd ! Number of soil hydraulic layers - integer :: nlevsoil ! Number of total soil layers + ! Locals + integer :: nlevsoil_hyd ! Number of total soil layers type(ed_cohort_hydr_type), pointer :: ccohort_hydr - !----------------------------------------------------------------------- - - - nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd - nlevsoil = bc_in%nlevsoil - ccohort => cc_p - ccohort_hydr => cc_p%co_hydr - cPatch => ccohort%patchptr + 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) + + ! 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_hydr,bc_in) + + + end subroutine updateSizeDepTreeHydProps - ! 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(:) + ! ===================================================================================== + subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hydr,bc_in) - ! Pre-process biomass pools + ! ----------------------------------------------------------------------------------- + ! 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 + ! + ! ----------------------------------------------------------------------------------- - 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) + ! 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 - !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 + type(ed_cohort_hydr_type), intent(inout) :: ccohort_hydr ! Plant hydraulics structure + integer :: pft ! Plant functional type index + 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) :: v_stem ! aboveground stem volume [m3/indiv] + 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) :: 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] + + 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._r8) then - 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) - sla = EDPftvarcon_inst%slatop(ft) * cm2_per_m2 ! m2/gC * cm2/m2 -> cm2/gC + ! ------------------------------------------------------------------------------ + ! Part 1. Set the volumes of the leaf, stem and root compartments + ! and lenghts of the roots + ! ------------------------------------------------------------------------------ - 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 / real(n_hypool_leaf,r8) + 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 + + ccohort_hydr%v_ag(1:n_hypool_leaf) = v_canopy / real(n_hypool_leaf,r8) - 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 + b_stem_carb = b_tot_carb - b_bg_carb - b_canopy_carb + b_stem_biom = b_stem_carb * C2B ! kg DM + + !BOC...may be needed for testing/comparison w/ v_sapwood + v_stem = b_stem_biom / (EDPftvarcon_inst%wood_density(ft)*1.e3_r8) - call bsap_allom(ccohort%dbh,ccohort%pft,ccohort%canopy_trim,a_sapwood_target,bsw_target) + ! m2 leaf = kg leaf DM * cm2/g * 1000g/1kg * 1m2/10000cm2 + a_leaf_tot = b_canopy_carb * sla * 1.e3_r8 / 1.e4_r8 + + call bsap_allom(ccohort%dbh,ccohort%pft,ccohort%canopy_trim,a_sapwood_target,bsw_target) - a_sapwood = a_sapwood_target + a_sapwood = a_sapwood_target + + ! 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 + + v_sapwood = a_sapwood * z_stem + ccohort_hydr%v_ag(n_hypool_leaf+1:n_hypool_ag) = v_sapwood / n_hypool_stem - ! 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 - - 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 + ! Estimate absorbing root total length (all layers) + ! ------------------------------------------------------------------------------ + ccohort_hydr%l_aroot_tot = fnrt_c*C2B*EDPftvarcon_inst%hydr_srl(ft) - !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 + ! 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 - 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 + + ! 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 + - + ! ------------------------------------------------------------------------------ + ! Part II. Set maximum (size-dependent) hydraulic conductances + ! ------------------------------------------------------------------------------ - - 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)) - & + ! 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 + call CrownDepth(ccohort%hite,crown_depth) + 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 * 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)) - & 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 - + 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 - ! 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 - call CrownDepth(ccohort%hite,crown_depth) - 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 * 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)) - & - 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 ! ===================================================================================== subroutine updateSizeDepTreeHydStates(currentSite,cc_p) @@ -953,6 +1021,7 @@ end subroutine FuseCohortHydraulics ! ===================================================================================== ! Initialization Routines ! ===================================================================================== + subroutine InitHydrCohort(currentSite,currentCohort) ! Arguments @@ -964,8 +1033,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 ! ===================================================================================== From 44d639db4deff42702efde7e3813987f4981c1b3 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 14 Nov 2018 18:51:08 -0800 Subject: [PATCH 07/12] Debugging new hydro restarts --- biogeochem/EDCohortDynamicsMod.F90 | 19 ++-- biogeochem/FatesAllometryMod.F90 | 4 +- biogeophys/FatesPlantHydraulicsMod.F90 | 151 ++++++++++++++++++++----- main/FatesRestartInterfaceMod.F90 | 147 +++++++++--------------- 4 files changed, 190 insertions(+), 131 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index a24d75a773..6845fddd88 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -132,8 +132,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_hydr_type), pointer :: ccohort_hydr + 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 !---------------------------------------------------------------------- @@ -251,22 +251,21 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine if( hlm_use_planthydro.eq.itrue ) then - ccohort_hydr => ccohort%co_hydr - nlevsoi_hydr = currentSite%si_hydr%nlevsoi_hyd - + nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd + ! This allocates array spaces call InitHydrCohort(currentSite,new_cohort) ! This calculates node heights - call UpdateTreeHydrNodes(ccohort_hydr,ft,new_cohort%hite,nlevsoi_hydr) + call UpdateTreeHydrNodes(new_cohort%co_hydr,new_cohort%pft,new_cohort%hite,nlevsoi_hyd) ! This calculates volumes, lengths and max conductances - call UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hydr,bc_in) - + 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(ccohort_hydr) - + 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) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index d1a0eb8fa9..45df34aebc 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -1847,7 +1847,7 @@ end subroutine h2d_martcano ! ===================================================================================== - subroutine CrownDepth(hite,crown_depth) + subroutine CrownDepth(height,crown_depth) ! ----------------------------------------------------------------------------------- ! This routine returns the depth of a plant's crown. Which is the length @@ -1864,7 +1864,7 @@ subroutine CrownDepth(hite,crown_depth) ! crown depth from Poorter, Bongers & Bongers ! crown_depth = exp(-1.169_r8)*cCohort%hite**1.098_r8 - crown_depth = min(cCohort%hite,0.1_r8) + crown_depth = min(height,0.1_r8) return end subroutine CrownDepth diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 41ec79b178..cc1d65d992 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -36,7 +36,8 @@ 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 EDParamsMod , only : hydr_kmax_rsurf use EDTypesMod , only : ed_site_type @@ -139,6 +140,10 @@ module FatesPlantHydraulicsMod public :: initTreeHydStates public :: updateSizeDepRhizHydProps public :: updateSizeDepRhizHydStates + public :: RestartHydrStates + public :: SavePreviousCompartmentVolumes + public :: UpdateTreeHydrNodes + public :: UpdateTreeHydrLenVolCond !------------------------------------------------------------------------------ ! 01/18/16: Created by Brad Christoffersen @@ -175,9 +180,79 @@ 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 + 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)) + + ! This calculates node heights + call UpdateTreeHydrNodes(ccohort%co_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%co_hydr) + + ccohort%co_hydr%errh2o_growturn_aroot(:) = 0.0_r8 + ccohort%co_hydr%errh2o_growturn_ag(:) = 0.0_r8 + ccohort%co_hydr%errh2o_growturn_troot(:) = 0.0_r8 + + ccohort => ccohort%taller + enddo + + cpatch => cpatch%younger + end do + + call updateSizeDepRhizHydProps(sites(s), bc_in(s) ) + + end do + + call UpdateH2OVeg(nsites,sites,bc_out(s)) + + return + end subroutine RestartHydrStates + + ! ==================================================================================== + subroutine initTreeHydStates(site_p, cc_p, bc_in) ! REQUIRED INPUTS: @@ -274,7 +349,7 @@ end subroutine initTreeHydStates ! ===================================================================================== - subroutine UpdateTreeHydrNodes(ccohort_hydr,pft,plant_height,nlevsoi_hyd) + subroutine UpdateTreeHydrNodes(ccohort_hydr,ft,plant_height,nlevsoi_hyd,bc_in) ! -------------------------------------------------------------------------------- ! This subroutine calculates the nodal heights critical to hydraulics in the plant @@ -294,9 +369,11 @@ subroutine UpdateTreeHydrNodes(ccohort_hydr,pft,plant_height,nlevsoi_hyd) ! Arguments type(ed_cohort_hydr_type), intent(inout) :: ccohort_hydr - integer,intent(in) :: pft ! plant functional type index + 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 @@ -305,6 +382,7 @@ subroutine UpdateTreeHydrNodes(ccohort_hydr,pft,plant_height,nlevsoi_hyd) 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] @@ -313,7 +391,10 @@ subroutine UpdateTreeHydrNodes(ccohort_hydr,pft,plant_height,nlevsoi_hyd) ! 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) @@ -401,7 +482,6 @@ subroutine updateSizeDepTreeHydProps(currentSite,ccohort,bc_in) ! shell radii, volumes, and compartment volumes of plant tissues ! !USES: - use FatesConstantsMod , only : pi_const use shr_sys_mod , only : shr_sys_abort ! ARGUMENTS: @@ -410,8 +490,9 @@ subroutine updateSizeDepTreeHydProps(currentSite,ccohort,bc_in) type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions ! Locals - integer :: nlevsoil_hyd ! Number of total soil layers + 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 @@ -423,21 +504,21 @@ subroutine updateSizeDepTreeHydProps(currentSite,ccohort,bc_in) call SavePreviousCompartmentVolumes(ccohort_hydr) ! This updates all of the z_node positions - call UpdateTreeHydrNodes(ccohort_hydr,ft,ccohort%hite,nlevsoi_hyd) + 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_hydr,bc_in) + call UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) end subroutine updateSizeDepTreeHydProps ! ===================================================================================== - subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hydr,bc_in) + subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) ! ----------------------------------------------------------------------------------- ! This subroutine calculates three attributes of a plant: @@ -456,11 +537,12 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hydr,bc_in) 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), intent(inout) :: ccohort_hydr ! Plant hydraulics structure - integer :: pft ! Plant functional type index + 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] @@ -471,11 +553,17 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hydr,bc_in) 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) :: a_leaf_tot ! total leaf area [m2/indiv] 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] @@ -499,11 +587,20 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hydr,bc_in) 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 + ft = ccohort%pft + 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 = EDPftvarcon_inst%roota_par(ft) + rootb = EDPftvarcon_inst%rootb_par(ft) + !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) @@ -556,7 +653,9 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hydr,bc_in) ! 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 @@ -619,10 +718,9 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hydr,bc_in) 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) + 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 - call CrownDepth(ccohort%hite,crown_depth) if(crown_depth == 0._r8) then write(fates_log(),*) 'do_kbound_upstream requires a nonzero canopy depth ' call endrun(msg=errMsg(sourcefile, __LINE__)) @@ -669,7 +767,7 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hydr,bc_in) 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,:) + 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 @@ -687,7 +785,7 @@ end subroutine UpdateTreeHydrLenVolCond ! ===================================================================================== - subroutine updateSizeDepTreeHydStates(currentSite,cc_p) + subroutine updateSizeDepTreeHydStates(currentSite,ccohort) ! ! !DESCRIPTION: ! @@ -696,11 +794,10 @@ subroutine updateSizeDepTreeHydStates(currentSite,cc_p) use EDTypesMod , only : dump_cohort ! !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 integer :: j,k,FT ! indices integer :: err_code = 0 @@ -708,11 +805,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 @@ -1360,7 +1456,6 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) ! the same. ! ! !USES: - use FatesConstantsMod , only : pi_const use EDTypesMod , only : AREA ! !ARGUMENTS: @@ -4879,7 +4974,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/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 5e59e55137..ed56123425 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -140,9 +140,6 @@ module FatesRestartInterfaceMod integer, private :: ir_prt_base ! Base index for all PRT variables ! Hydraulic indices - integer, private :: ir_hydro_v_ag_covec - integer, private :: ir_hydro_v_troot_covec - integer, private :: ir_hydro_v_aroot_layer_covec integer, private :: ir_hydro_th_ag_covec integer, private :: ir_hydro_th_troot_covec integer, private :: ir_hydro_th_aroot_covec @@ -860,22 +857,6 @@ subroutine define_restart_vars(this, initialize_variables) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - - call this%RegisterCohortVector(symbol_base='fates_hydro_v_ag', vtype=cohort_r8, & - long_name_base='maximum storage volume of hydraulic compartments (above ground)', & - units='m3', veclength=n_hypool_ag, flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_ag_covec) - - call this%RegisterCohortVector(symbol_base='fates_hydro_v_troot', vtype=cohort_r8, & - long_name_base='maximum storage volume of transporting root compartments', & - units='m3', veclength=n_hypool_troot, flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_troot_covec) - - call this%RegisterCohortVector(symbol_base='fates_hydro_v_aroot_layer', vtype=cohort_r8, & - long_name_base='maximum storage volume of absorbing roots hydraulic compartments by layer', & - units='m3', veclength=nlevsoi_hyd_max, flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_aroot_layer_covec) - 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, & @@ -1461,14 +1442,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) if(hlm_use_planthydro==itrue)then - ! Load the storage compartment volumes - call this%SetCohortRealVector(ccohort%co_hydr%v_ag,n_hypool_ag, & - ir_hydro_v_ag_covec,io_idx_co) - call this%SetCohortRealVector(ccohort%co_hydr%v_troot,n_hypool_troot, & - ir_hydro_v_troot_covec,io_idx_co) - call this%SetCohortRealVector(ccohort%co_hydr%v_aroot_layer,sites(s)%si_hydr%nlevsoi_hyd, & - ir_hydro_v_aroot_layer_covec,io_idx_co) - ! Load the water contents call this%SetCohortRealVector(ccohort%co_hydr%th_ag,n_hypool_ag, & ir_hydro_th_ag_covec,io_idx_co) @@ -1707,22 +1680,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 @@ -1753,7 +1722,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 @@ -1788,45 +1757,48 @@ 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) + 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 - - ! Hydraulics - if turned on, the hydraulics arrays are being allocated in create_cohort as well. - - 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 @@ -2083,27 +2055,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end do - if(hlm_use_planthydro==itrue)then - - ! Load the storage compartment volumes - call this%GetCohortRealVector(ccohort%co_hydr%v_ag,n_hypool_ag, & - ir_hydro_v_ag_covec,io_idx_co) - call this%GetCohortRealVector(ccohort%co_hydr%v_troot,n_hypool_troot, & - ir_hydro_v_troot_covec,io_idx_co) - call this%GetCohortRealVector(ccohort%co_hydr%v_aroot_layer,sites(s)%si_hydr%nlevsoi_hyd, & - ir_hydro_v_aroot_layer_covec,io_idx_co) - - ! 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) - end if - - - ccohort%canopy_layer = rio_canopy_layer_co(io_idx_co) ccohort%canopy_layer_yesterday = rio_canopy_layer_yesterday_co(io_idx_co) ccohort%canopy_trim = rio_canopy_trim_co(io_idx_co) @@ -2135,6 +2086,19 @@ 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) + end if io_idx_co = io_idx_co + 1 @@ -2295,7 +2259,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 @@ -2305,8 +2269,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 @@ -2320,7 +2286,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) From 015533ccb9f7624cfb5afd363c69e47307676dc3 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 15 Nov 2018 21:56:20 -0800 Subject: [PATCH 08/12] fates hydro restart refactors - syntax updates --- biogeochem/EDCohortDynamicsMod.F90 | 6 +++++- biogeophys/FatesPlantHydraulicsMod.F90 | 26 ++++++++++++++------------ main/FatesRestartInterfaceMod.F90 | 6 +++++- 3 files changed, 24 insertions(+), 14 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 6845fddd88..57ea9dca09 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 @@ -257,7 +260,8 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine call InitHydrCohort(currentSite,new_cohort) ! This calculates node heights - call UpdateTreeHydrNodes(new_cohort%co_hydr,new_cohort%pft,new_cohort%hite,nlevsoi_hyd) + 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) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index cc1d65d992..c4212432af 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -37,6 +37,8 @@ module FatesPlantHydraulicsMod 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 @@ -50,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 @@ -77,8 +80,6 @@ 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 @@ -448,8 +449,6 @@ subroutine UpdateTreeHydrNodes(ccohort_hydr,ft,plant_height,nlevsoi_hyd,bc_in) ccohort_hydr%z_node_troot(:) = ccohort_hydr%z_node_aroot(nlevsoi_hyd) end if - - return end subroutine UpdateTreeHydrNodes @@ -468,7 +467,6 @@ subroutine SavePreviousCompartmentVolumes(ccohort_hydr) ccohort_hydr%v_troot_init(:) = ccohort_hydr%v_troot(:) ccohort_hydr%v_aroot_layer_init(:) = ccohort_hydr%v_aroot_layer(:) - return end subroutine SavePreviousCompartmentVolumes @@ -636,11 +634,12 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) b_stem_carb = b_tot_carb - b_bg_carb - b_canopy_carb b_stem_biom = b_stem_carb * C2B ! kg DM - !BOC...may be needed for testing/comparison w/ v_sapwood - v_stem = b_stem_biom / (EDPftvarcon_inst%wood_density(ft)*1.e3_r8) + !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 ) ! m2 leaf = kg leaf DM * cm2/g * 1000g/1kg * 1m2/10000cm2 - a_leaf_tot = b_canopy_carb * sla * 1.e3_r8 / 1.e4_r8 + a_leaf_tot = b_canopy_carb * sla * g_per_kg / cm2_per_m2 call bsap_allom(ccohort%dbh,ccohort%pft,ccohort%canopy_trim,a_sapwood_target,bsw_target) @@ -660,8 +659,9 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) 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 + ! 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 @@ -677,7 +677,8 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) ! 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 + 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 @@ -767,7 +768,8 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) 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,:) + 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 diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index ed56123425..673c994122 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -22,7 +22,10 @@ module FatesRestartInterfaceMod 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 @@ -1772,6 +1775,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) allocate(new_cohort) call nan_cohort(new_cohort) + call zero_cohort(new_cohort) new_cohort%patchptr => newp ! If this is the first in the list, it is tallest From 6481b735adf70d733c291c9706b0031e358700cc Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 16 Nov 2018 11:47:36 -0800 Subject: [PATCH 09/12] Added some zero-ing initializations to hydraulics during restarts --- biogeophys/FatesPlantHydraulicsMod.F90 | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index f43e20d67a..6d4ad18eea 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -245,11 +245,18 @@ subroutine RestartHydrStates(sites,nsites,bc_in,bc_out) cpatch => cpatch%younger end do + sites(s)%si_hydr%h2oveg = 0.0_r8 + sites(s)%si_hydr%h2oveg_recruit = 0.0_r8 + sites(s)%si_hydr%h2oveg_dead = 0.0_r8 + sites(s)%si_hydr%h2oveg_growturn_err = 0.0_r8 + sites(s)%si_hydr%h2oveg_pheno_err = 0.0_r8 + sites(s)%si_hydr%h2oveg_hydro_err = 0.0_r8 + call updateSizeDepRhizHydProps(sites(s), bc_in(s) ) end do - call UpdateH2OVeg(nsites,sites,bc_out(s)) + call UpdateH2OVeg(nsites,sites,bc_out) return end subroutine RestartHydrStates @@ -563,7 +570,6 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) 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) :: a_leaf_tot ! total leaf area [m2/indiv] 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] @@ -639,12 +645,11 @@ subroutine UpdateTreeHydrLenVolCond(ccohort,nlevsoi_hyd,bc_in) ! kg / ( g cm-3 * cm3/m3 * kg/g ) -> m3 v_stem = b_stem_biom / (EDPftvarcon_inst%wood_density(ft)*1.e3_r8 ) - ! calculate the leaf area based on the leaf profile - a_leaf_tot = cCohort%treelai * cCohort%c_area/cCohort%n + ! 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 @@ -1476,7 +1481,7 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) 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 From 50e8141bbe2ce030f4d9b3c15ab8f3a8b9a2f2fe Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 19 Nov 2018 14:43:10 -0800 Subject: [PATCH 10/12] Corrected hydro-restart initialization procedure --- biogeophys/FatesPlantHydraulicsMod.F90 | 120 ++++++++++++++++++++----- main/FatesRestartInterfaceMod.F90 | 116 ++++++++++++++++++++++-- 2 files changed, 207 insertions(+), 29 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 6d4ad18eea..68e342752c 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -145,6 +145,7 @@ module FatesPlantHydraulicsMod public :: updateSizeDepRhizHydStates public :: RestartHydrStates public :: SavePreviousCompartmentVolumes + public :: SavePreviousRhizVolumes public :: UpdateTreeHydrNodes public :: UpdateTreeHydrLenVolCond @@ -212,9 +213,10 @@ subroutine RestartHydrStates(sites,nsites,bc_in,bc_out) ! locals ! ---------------------------------------------------------------------------------- ! LL pointers - type(ed_patch_type),pointer :: cpatch ! current patch - type(ed_cohort_type),pointer :: ccohort ! current cohort - integer :: s ! site loop counter + 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 @@ -224,8 +226,10 @@ subroutine RestartHydrStates(sites,nsites,bc_in,bc_out) ccohort => cpatch%shortest do while(associated(ccohort)) + ccohort_hydr => ccohort%co_hydr + ! This calculates node heights - call UpdateTreeHydrNodes(ccohort%co_hydr,ccohort%pft,ccohort%hite, & + call UpdateTreeHydrNodes(ccohort_hydr,ccohort%pft,ccohort%hite, & sites(s)%si_hydr%nlevsoi_hyd,bc_in(s)) ! This calculates volumes, lengths and max conductances @@ -233,26 +237,42 @@ subroutine RestartHydrStates(sites,nsites,bc_in,bc_out) ! Since this is a newly initialized plant, we set the previous compartment-size ! equal to the ones we just calculated. - call SavePreviousCompartmentVolumes(ccohort%co_hydr) - - ccohort%co_hydr%errh2o_growturn_aroot(:) = 0.0_r8 - ccohort%co_hydr%errh2o_growturn_ag(:) = 0.0_r8 - ccohort%co_hydr%errh2o_growturn_troot(:) = 0.0_r8 + 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(:) = -9.9e10_r8 + sites(s)%si_hydr%r_node_shell_init(:,:) = -9.9e10_r8 + sites(s)%si_hydr%v_shell_init(:,:) = -9.9e10_r8 + - sites(s)%si_hydr%h2oveg = 0.0_r8 - sites(s)%si_hydr%h2oveg_recruit = 0.0_r8 - sites(s)%si_hydr%h2oveg_dead = 0.0_r8 - sites(s)%si_hydr%h2oveg_growturn_err = 0.0_r8 - sites(s)%si_hydr%h2oveg_pheno_err = 0.0_r8 - sites(s)%si_hydr%h2oveg_hydro_err = 0.0_r8 - call updateSizeDepRhizHydProps(sites(s), bc_in(s) ) + ! 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 @@ -1435,7 +1455,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 @@ -1469,14 +1510,10 @@ 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 @@ -1523,8 +1560,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) / & @@ -1570,6 +1609,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 ! ================================================================================= diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 673c994122..3412e07ae5 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -150,7 +150,14 @@ module FatesRestartInterfaceMod integer, private :: ir_hydro_r_node_shell_si integer, private :: ir_hydro_v_shell_si 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) @@ -875,6 +882,22 @@ subroutine define_restart_vars(this, initialize_variables) 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 absorbing root maximum volume call this%set_restart_var(vname='fates_hydro_l_aroot_layer', vtype=cohort_r8, & long_name='Total length (across cohorts) of absorbing roots by soil layer', & @@ -898,6 +921,36 @@ subroutine define_restart_vars(this, initialize_variables) 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 @@ -1346,7 +1399,13 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_hydro_l_aroot_layer_si => this%rvars(ir_hydro_l_aroot_layer_si)%r81d, & rio_hydro_r_node_shell_si => this%rvars(ir_hydro_r_node_shell_si)%r81d, & rio_hydro_v_shell_si => this%rvars(ir_hydro_v_shell_si)%r81d, & - rio_hydro_liqvol_shell_si => this%rvars(ir_hydro_liqvol_shell_si)%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 @@ -1452,6 +1511,20 @@ subroutine set_restart_vectors(this,nc,nsites,sites) 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 @@ -1611,6 +1684,12 @@ subroutine set_restart_vectors(this,nc,nsites,sites) 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 @@ -1630,6 +1709,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) end do rio_hydro_l_aroot_layer_si(io_idx_si_lyr) = sites(s)%si_hydr%l_aroot_layer(i) + io_idx_si_lyr = io_idx_si_lyr + 1 end do end if @@ -1979,7 +2059,12 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_hydro_l_aroot_layer_si => this%rvars(ir_hydro_l_aroot_layer_si)%r81d, & rio_hydro_r_node_shell_si => this%rvars(ir_hydro_r_node_shell_si)%r81d, & rio_hydro_v_shell_si => this%rvars(ir_hydro_v_shell_si)%r81d, & - rio_hydro_liqvol_shell_si => this%rvars(ir_hydro_liqvol_shell_si)%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 @@ -2097,11 +2182,23 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! Load the water contents call this%GetCohortRealVector(ccohort%co_hydr%th_ag,n_hypool_ag, & - ir_hydro_th_ag_covec,io_idx_co) + 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) + 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) + 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 @@ -2209,6 +2306,12 @@ subroutine get_restart_vectors(this, nc, nsites, sites) 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 @@ -2234,7 +2337,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) 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) From 757d93b31d9fbbc2fff868942bc64c81f905b801 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 19 Nov 2018 15:25:26 -0800 Subject: [PATCH 11/12] Removed length, node and volumes from shell restarts as they are auto calculated --- biogeophys/FatesPlantHydraulicsMod.F90 | 7 ++- main/FatesRestartInterfaceMod.F90 | 61 -------------------------- 2 files changed, 3 insertions(+), 65 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 68e342752c..8db569c78e 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -259,10 +259,9 @@ subroutine RestartHydrStates(sites,nsites,bc_in,bc_out) cpatch => cpatch%younger end do - sites(s)%si_hydr%l_aroot_layer_init(:) = -9.9e10_r8 - sites(s)%si_hydr%r_node_shell_init(:,:) = -9.9e10_r8 - sites(s)%si_hydr%v_shell_init(:,:) = -9.9e10_r8 - + 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 diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 3412e07ae5..fa851b32fc 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -146,9 +146,6 @@ module FatesRestartInterfaceMod 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_l_aroot_layer_si - integer, private :: ir_hydro_r_node_shell_si - integer, private :: ir_hydro_v_shell_si integer, private :: ir_hydro_liqvol_shell_si integer, private :: ir_hydro_err_growturn_aroot_covec integer, private :: ir_hydro_err_growturn_ag_covec @@ -897,25 +894,6 @@ subroutine define_restart_vars(this, initialize_variables) 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 absorbing root maximum volume - call this%set_restart_var(vname='fates_hydro_l_aroot_layer', vtype=cohort_r8, & - long_name='Total length (across cohorts) of absorbing roots by soil layer', & - units='m', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_l_aroot_layer_si ) - - ! Site-level Nodal Radius of rhizosphere compartments - call this%set_restart_var(vname='fates_hydro_r_node_shell', vtype=cohort_r8, & - long_name='Nodal Radius of rhizosphere compartments (layerxshell)', & - units='m', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_r_node_shell_si ) - - ! Site-level volume of rhizosphere compartments - call this%set_restart_var(vname='fates_hydro_v_shell', vtype=cohort_r8, & - long_name='Volume of rhizosphere compartments (layerxshell)', & - units='m3', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_hydro_v_shell_si ) - ! 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)', & @@ -1309,7 +1287,6 @@ 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 ! site - layer x index integer :: io_idx_si_lyr_shell ! site - layer x shell index ! Some counters (for checking mostly) @@ -1396,9 +1373,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) 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_hydro_l_aroot_layer_si => this%rvars(ir_hydro_l_aroot_layer_si)%r81d, & - rio_hydro_r_node_shell_si => this%rvars(ir_hydro_r_node_shell_si)%r81d, & - rio_hydro_v_shell_si => this%rvars(ir_hydro_v_shell_si)%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, & @@ -1432,11 +1406,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_si_wmem = io_idx_co_1st ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell - io_idx_si_lyr = io_idx_co_1st io_idx_si_lyr_shell = io_idx_co_1st - - ! write seed_bank info(site-level, but PFT-resolved) do i = 1,numpft rio_seed_bank_sift(io_idx_co_1st+i-1) = sites(s)%seed_bank(i) @@ -1692,25 +1663,12 @@ subroutine set_restart_vectors(this,nc,nsites,sites) ! 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_r_node_shell_si(io_idx_si_lyr_shell) = & - sites(s)%si_hydr%r_node_shell(i,k) - - rio_hydro_v_shell_si(io_idx_si_lyr_shell) = & - sites(s)%si_hydr%v_shell(i,k) - 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 - - rio_hydro_l_aroot_layer_si(io_idx_si_lyr) = sites(s)%si_hydr%l_aroot_layer(i) - - io_idx_si_lyr = io_idx_si_lyr + 1 end do end if @@ -1976,7 +1934,6 @@ 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 ! site - layer x index integer :: io_idx_si_lyr_shell ! site - layer x shell index ! Some counters (for checking mostly) @@ -2056,9 +2013,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) 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_hydro_l_aroot_layer_si => this%rvars(ir_hydro_l_aroot_layer_si)%r81d, & - rio_hydro_r_node_shell_si => this%rvars(ir_hydro_r_node_shell_si)%r81d, & - rio_hydro_v_shell_si => this%rvars(ir_hydro_v_shell_si)%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, & @@ -2080,7 +2034,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_si_wmem = io_idx_co_1st ! Hydraulics counters lyr = hydraulic layer, shell = rhizosphere shell - io_idx_si_lyr = io_idx_co_1st io_idx_si_lyr_shell = io_idx_co_1st @@ -2314,26 +2267,12 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! 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%r_node_shell(i,k) = & - rio_hydro_r_node_shell_si(io_idx_si_lyr_shell) - - sites(s)%si_hydr%v_shell(i,k) = & - rio_hydro_v_shell_si(io_idx_si_lyr_shell) - 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 - - sites(s)%si_hydr%l_aroot_layer(i) = & - rio_hydro_l_aroot_layer_si(io_idx_si_lyr) - - io_idx_si_lyr = io_idx_si_lyr + 1 end do end if From 9d70ddd84e2b03d3c1cf69f721fe27bc8ea2fd03 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 20 Nov 2018 13:50:38 -0800 Subject: [PATCH 12/12] Added debug protection to an abundant water balance error write statement --- biogeophys/FatesPlantHydraulicsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 8db569c78e..78046c1aa8 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -3361,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