From 01019caf8910109fbcd739997f491626a8734599 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 7 May 2018 14:53:11 -0700 Subject: [PATCH 01/12] removed hlm_numlevsoil and started converting most arrays to the site defined quantity bc_in(s)%nlevsoil. Part-way through converting decomposition stuff. --- biogeochem/EDPatchDynamicsMod.F90 | 70 +++++++++----- biogeochem/EDPhysiologyMod.F90 | 18 ++-- biogeophys/EDBtranMod.F90 | 45 +++++---- biogeophys/FatesPlantRespPhotosynthMod.F90 | 9 +- main/EDInitMod.F90 | 2 +- main/FatesInterfaceMod.F90 | 107 ++++++++++----------- main/FatesInventoryInitMod.F90 | 2 +- main/FatesRestartInterfaceMod.F90 | 2 +- 8 files changed, 137 insertions(+), 118 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 992e798ec7..2364407bf7 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -18,8 +18,6 @@ module EDPatchDynamicsMod use EDTypesMod , only : dtype_ilog use EDTypesMod , only : dtype_ifire use FatesInterfaceMod , only : hlm_use_planthydro - use FatesInterfaceMod , only : hlm_numlevgrnd - use FatesInterfaceMod , only : hlm_numlevsoil use FatesInterfaceMod , only : hlm_numSWb use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_days_per_year @@ -345,7 +343,7 @@ subroutine spawn_patches( currentSite, bc_in) allocate(new_patch) call create_patch(currentSite, new_patch, age, site_areadis, & cwd_ag_local, cwd_bg_local, leaf_litter_local, & - root_litter_local) + root_litter_local, bc_in%nlevsoil) new_patch%tallest => null() new_patch%shortest => null() @@ -1116,7 +1114,7 @@ end subroutine mortality_litter_fluxes ! ============================================================================ subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_local, & - leaf_litter_local,root_litter_local) + leaf_litter_local,root_litter_local,nlevsoil) ! ! !DESCRIPTION: ! Set default values for creating a new patch @@ -1126,12 +1124,13 @@ subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_l ! !ARGUMENTS: type(ed_site_type) , intent(inout), target :: currentSite type(ed_patch_type), intent(inout), target :: new_patch - real(r8), intent(in) :: age ! notional age of this patch in years - real(r8), intent(in) :: areap ! initial area of this patch in m2. - real(r8), intent(in) :: cwd_ag_local(:) ! initial value of above ground coarse woody debris. KgC/m2 - real(r8), intent(in) :: cwd_bg_local(:) ! initial value of below ground coarse woody debris. KgC/m2 - real(r8), intent(in) :: root_litter_local(:)! initial value of root litter. KgC/m2 - real(r8), intent(in) :: leaf_litter_local(:)! initial value of leaf litter. KgC/m2 + real(r8), intent(in) :: age ! notional age of this patch in years + real(r8), intent(in) :: areap ! initial area of this patch in m2. + real(r8), intent(in) :: cwd_ag_local(:) ! initial value of above ground coarse woody debris. KgC/m2 + real(r8), intent(in) :: cwd_bg_local(:) ! initial value of below ground coarse woody debris. KgC/m2 + real(r8), intent(in) :: root_litter_local(:) ! initial value of root litter. KgC/m2 + real(r8), intent(in) :: leaf_litter_local(:) ! initial value of leaf litter. KgC/m2 + integer, intent(in) :: nlevsoil ! number of soil layers ! ! !LOCAL VARIABLES: !--------------------------------------------------------------------- @@ -1144,8 +1143,8 @@ subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_l allocate(new_patch%fabi(hlm_numSWb)) allocate(new_patch%sabs_dir(hlm_numSWb)) allocate(new_patch%sabs_dif(hlm_numSWb)) - allocate(new_patch%rootfr_ft(numpft,hlm_numlevgrnd)) - allocate(new_patch%rootr_ft(numpft,hlm_numlevgrnd)) + allocate(new_patch%rootfr_ft(numpft,nlevsoil) + allocate(new_patch%rootr_ft(numpft,nlevsoil) call zero_patch(new_patch) !The nan value in here is not working?? @@ -1892,36 +1891,55 @@ end function countPatches ! ==================================================================================== - subroutine set_root_fraction( cpatch , zi ) + subroutine set_root_fraction( root_fraction , zi ) ! ! !DESCRIPTION: ! Calculates the fractions of the root biomass in each layer for each pft. + ! It assumes an exponential decay. If the soil depth is shallower than + ! then exponential attenuation function, then it will normalize + ! the profile and divide through. ! ! !USES: ! ! !ARGUMENTS - type(ed_patch_type),intent(inout), target :: cpatch - real(r8),intent(in) :: zi(0:hlm_numlevsoil) + real(r8),intent(out) :: root_fraction(:,:) + real(r8),intent(in) :: zi(:) ! ! !LOCAL VARIABLES: integer :: lev,p,c,ft + integer :: nlevsoil + real(r8) :: sum_rootfr !---------------------------------------------------------------------- + if(lbound(zi,1).ne.0) then + write(fates_log(),*) 'layer interface levels should have 0 index' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + nlevsoil = ubound(zi,1) + do ft = 1,numpft - do lev = 1, hlm_numlevgrnd - cpatch%rootfr_ft(ft,lev) = 0._r8 - enddo + + sum_rootfr = 0.0_r8 - do lev = 1, hlm_numlevsoil-1 - cpatch%rootfr_ft(ft,lev) = .5_r8*( & - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev-1)) & - + exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev-1)) & - - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev)) & - - exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev))) + do lev = 1, nlevsoil + + root_fraction(ft,lev) = .5_r8*( & + exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev-1)) & + + exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev-1)) & + - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev)) & + - exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev))) + + sum_rootfr = sum_rootfr + cpatch%rootfr_ft(ft,lev) + end do + + ! Normalize the root profile + root_fraction(ft,:) = cpatch%rootfr_ft(ft,:)/sum_rootfr + end do - - end subroutine set_root_fraction + + end subroutine set_root_fraction end module EDPatchDynamicsMod diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 8ae081d9c2..5bdc6c9001 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -454,7 +454,7 @@ subroutine phenology( currentSite, bc_in ) do i = 1,numWaterMem-1 !shift memory along one currentSite%water_memory(numWaterMem+1-i) = currentSite%water_memory(numWaterMem-i) enddo - currentSite%water_memory(1) = bc_in%h2o_liqvol_gl(1) !waterstate_inst%h2osoi_vol_col(coli,1) + currentSite%water_memory(1) = bc_in%h2o_liqvol_sl(1) !waterstate_inst%h2osoi_vol_col(coli,1) !In drought phenology, we often need to force the leaves to stay on or off as moisture fluctuates... timesincedleafoff = 0 @@ -2005,7 +2005,6 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) use EDTypesMod, only : AREA - use FatesInterfaceMod, only : hlm_numlevdecomp_full use FatesInterfaceMod, only : hlm_numlevdecomp use EDPftvarcon, only : EDPftvarcon_inst use FatesConstantsMod, only : sec_per_day @@ -2034,11 +2033,14 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) integer :: begp,endp integer :: begc,endc !bounds !------------------------------------------------------------------------ - real(r8) :: cinput_rootfr(1:maxpft, 1:hlm_numlevdecomp_full) ! column by pft root fraction used for calculating inputs - real(r8) :: croot_prof_perpatch(1:hlm_numlevdecomp_full) - real(r8) :: surface_prof(1:hlm_numlevdecomp_full) + ! The following scratch arrays are allocated for maximum possible + ! pft and layer usage + real(r8) :: cinput_rootfr(1:maxpft, 1:hlm_numlevgrnd) ! column by pft root fraction used for calculating inputs + real(r8) :: croot_prof_perpatch(1:hlm_numlevgrnd) + real(r8) :: surface_prof(1:hlm_numlevgrnd) integer :: ft - real(r8) :: rootfr_tot(1:maxpft), biomass_bg_ft(1:maxpft) + real(r8) :: rootfr_tot(1:maxpft) + real(r8) :: biomass_bg_ft(1:maxpft) real(r8) :: surface_prof_tot, leaf_prof_sum, stem_prof_sum, froot_prof_sum, biomass_bg_tot real(r8) :: delta @@ -2126,6 +2128,10 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) end do endif else + + + + do ft = 1,numpft do j = 1, hlm_numlevdecomp ! use standard CLM root fraction profiles; diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 28593ad955..7f046dfe68 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -12,7 +12,6 @@ module EDBtranMod ed_patch_type, & ed_cohort_type, & maxpft - use FatesInterfaceMod , only : hlm_numlevgrnd use shr_kind_mod , only : r8 => shr_kind_r8 use FatesInterfaceMod , only : bc_in_type, & bc_out_type, & @@ -66,11 +65,11 @@ subroutine get_active_suction_layers(nsites, sites, bc_in, bc_out) do s = 1,nsites if (bc_in(s)%filter_btran) then - do j = 1,hlm_numlevgrnd - bc_out(s)%active_suction_gl(j) = check_layer_water( bc_in(s)%h2o_liqvol_gl(j),bc_in(s)%tempk_gl(j) ) + do j = 1,bc_in(s)%nlevsoil + bc_out(s)%active_suction_sl(j) = check_layer_water( bc_in(s)%h2o_liqvol_sl(j),bc_in(s)%tempk_sl(j) ) end do else - bc_out(s)%active_suction_gl(:) = .false. + bc_out(s)%active_suction_sl(:) = .false. end if end do @@ -86,11 +85,11 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) ! --------------------------------------------------------------------------------- ! Calculate the transpiration wetness function (BTRAN) and the root uptake ! distribution (ROOTR). - ! Boundary conditions in: bc_in(s)%eff_porosity_gl(j) unfrozen porosity - ! bc_in(s)%watsat_gl(j) porosity - ! bc_in(s)%active_uptake_gl(j) frozen/not frozen - ! bc_in(s)%smp_gl(j) suction - ! Boundary conditions out: bc_out(s)%rootr_pagl root uptake distribution + ! Boundary conditions in: bc_in(s)%eff_porosity_sl(j) unfrozen porosity + ! bc_in(s)%watsat_sl(j) porosity + ! bc_in(s)%active_uptake_sl(j) frozen/not frozen + ! bc_in(s)%smp_sl(j) suction + ! Boundary conditions out: bc_out(s)%rootr_pasl root uptake distribution ! bc_out(s)%btran_pa wetness factor ! --------------------------------------------------------------------------------- @@ -124,7 +123,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) do s = 1,nsites - bc_out(s)%rootr_pagl(:,:) = 0._r8 + bc_out(s)%rootr_pasl(:,:) = 0._r8 ifp = 0 cpatch => sites(s)%oldest_patch @@ -135,16 +134,16 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) do ft = 1,numpft cpatch%btran_ft(ft) = 0.0_r8 - do j = 1,hlm_numlevgrnd + do j = 1,bc_in(s)%nlevsoil ! Calculations are only relevant where liquid water exists ! see clm_fates%wrap_btran for calculation with CLM/ALM - if ( check_layer_water(bc_in(s)%h2o_liqvol_gl(j),bc_in(s)%tempk_gl(j)) ) then + if ( check_layer_water(bc_in(s)%h2o_liqvol_sl(j),bc_in(s)%tempk_sl(j)) ) then - smp_node = max(smpsc(ft), bc_in(s)%smp_gl(j)) + smp_node = max(smpsc(ft), bc_in(s)%smp_sl(j)) - rresis = min( (bc_in(s)%eff_porosity_gl(j)/bc_in(s)%watsat_gl(j))* & + rresis = min( (bc_in(s)%eff_porosity_sl(j)/bc_in(s)%watsat_sl(j))* & (smp_node - smpsc(ft)) / (smpso(ft) - smpsc(ft)), 1._r8) cpatch%rootr_ft(ft,j) = cpatch%rootfr_ft(ft,j)*rresis @@ -152,7 +151,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) ! root water uptake is not linearly proportional to root density, ! to allow proper deep root funciton. Replace with equations from SPA/Newman. FIX(RF,032414) ! cpatch%rootr_ft(ft,j) = cpatch%rootfr_ft(ft,j)**0.3*rresis_ft(ft,j)/ & - ! sum(cpatch%rootfr_ft(ft,1:nlevgrnd)**0.3) + ! sum(cpatch%rootfr_ft(ft,1:nlevsoil)**0.3) cpatch%btran_ft(ft) = cpatch%btran_ft(ft) + cpatch%rootr_ft(ft,j) else @@ -162,7 +161,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) end do !j ! Normalize root resistances to get layer contribution to ET - do j = 1,hlm_numlevgrnd + do j = 1,bc_in(s)%nlevsoil if (cpatch%btran_ft(ft) > 0.0_r8) then cpatch%rootr_ft(ft,j) = cpatch%rootr_ft(ft,j)/cpatch%btran_ft(ft) else @@ -189,15 +188,15 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) ! distributed over the soil layers. sum_pftgs = sum(pftgs(1:numpft)) - do j = 1,hlm_numlevgrnd - bc_out(s)%rootr_pagl(ifp,j) = 0._r8 + do j = 1, bc_in(s)%nlevsoil + bc_out(s)%rootr_pasl(ifp,j) = 0._r8 do ft = 1,numpft if( sum_pftgs > 0._r8)then !prevent problem with the first timestep - might fail !bit-retart test as a result? FIX(RF,032414) - bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j) + & + bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + & cpatch%rootr_ft(ft,j) * pftgs(ft)/sum_pftgs else - bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j) + & + bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j) + & cpatch%rootr_ft(ft,j) * 1._r8/dble(numpft) end if enddo @@ -220,12 +219,12 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) enddo end if - temprootr = sum(bc_out(s)%rootr_pagl(ifp,1:hlm_numlevgrnd)) + temprootr = sum(bc_out(s)%rootr_pasl(ifp,1:bc_in(s)%nlevsoil)) if(abs(1.0_r8-temprootr) > 1.0e-10_r8 .and. temprootr > 1.0e-10_r8)then write(fates_log(),*) 'error with rootr in canopy fluxes',temprootr,sum_pftgs - do j = 1,hlm_numlevgrnd - bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j)/temprootr + do j = 1,bc_in(s)%nlevsoil + bc_out(s)%rootr_pasl(ifp,j) = bc_out(s)%rootr_pasl(ifp,j)/temprootr enddo end if diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 4b3fca7041..e0d68ec76b 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -70,7 +70,6 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) use EDTypesMod , only : ed_cohort_type use EDTypesMod , only : ed_site_type use EDTypesMod , only : maxpft - use FatesInterfaceMod , only : hlm_numlevsoil use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : bc_out_type use EDCanopyStructureMod, only : calc_areaindex @@ -549,8 +548,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! Fine Root MR (kgC/plant/s) ! ------------------------------------------------------------------ currentCohort%froot_mr = 0._r8 - do j = 1,hlm_numlevsoil - tcsoi = q10**((bc_in(s)%t_soisno_gl(j)-tfrz - 20.0_r8)/10.0_r8) + do j = 1,bc_in(s)%nlevsoil + tcsoi = q10**((bc_in(s)%t_soisno_sl(j)-tfrz - 20.0_r8)/10.0_r8) currentCohort%froot_mr = currentCohort%froot_mr + & froot_n * ED_val_base_mr_20 * tcsoi * currentPatch%rootfr_ft(ft,j) * maintresp_reduction_factor enddo @@ -559,9 +558,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! ------------------------------------------------------------------ if (woody(ft) == 1) then currentCohort%livecroot_mr = 0._r8 - do j = 1,hlm_numlevsoil + do j = 1,bc_in(s)%nlevsoil ! Soil temperature used to adjust base rate of MR - tcsoi = q10**((bc_in(s)%t_soisno_gl(j)-tfrz - 20.0_r8)/10.0_r8) + tcsoi = q10**((bc_in(s)%t_soisno_sl(j)-tfrz - 20.0_r8)/10.0_r8) currentCohort%livecroot_mr = currentCohort%livecroot_mr + & live_croot_n * ED_val_base_mr_20 * tcsoi * & currentPatch%rootfr_ft(ft,j) * maintresp_reduction_factor diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index bb34503463..894f954a2a 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -307,7 +307,7 @@ subroutine init_patches( nsites, sites, bc_in) ! make new patch... call create_patch(sites(s), newp, age, AREA, & cwd_ag_local, cwd_bg_local, leaf_litter_local, & - root_litter_local) + root_litter_local, bc_in(s)%nlevsoil ) call init_cohorts(sites(s), newp, bc_in(s)) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 502c52ada1..3fce4d3a24 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -62,7 +62,8 @@ module FatesInterfaceMod integer, protected :: hlm_numlevgrnd ! Number of ground layers - integer, protected :: hlm_numlevsoil ! Number of soil layers + ! NOTE! SOIL LAYERS ARE NOT A GLOBAL, THEY + ! ARE VARIABLE BY SITE integer, protected :: hlm_numlevdecomp_full ! Number of GROUND layers for the purposes @@ -253,10 +254,11 @@ module FatesInterfaceMod ! or alias other variables, ALLOCATABLES cannnot. According to S. Lionel ! (Intel-Forum Post), ALLOCATABLES are better perfomance wise as long as they point ! to contiguous memory spaces and do not alias other variables, the case here. - ! Naming conventions: _gl means ground layer dimensions - ! _si means site dimensions (scalar in that case) + ! Naming conventions: _si means site dimensions (scalar in that case) ! _pa means patch dimensions ! _rb means radiation band + ! _sl means soil layer + ! _sisl means site x soil layer ! ------------------------------------------------------------------------------------ type, public :: bc_in_type @@ -266,10 +268,15 @@ module FatesInterfaceMod ! Soil layer structure + + integer :: nlevsoil ! the number of soil layers in this column real(r8),allocatable :: zi_sisl(:) ! interface level below a "z" level (m) ! this contains a zero index for surface. real(r8),allocatable :: dz_sisl(:) ! layer thickness (m) - real(r8),allocatable :: z_sisl(:) ! layer depth (m) (1:hlm_nlevsoil) + real(r8),allocatable :: z_sisl(:) ! layer depth (m) + + + ! Decomposition Layer Structure real(r8), allocatable :: dz_decomp_sisl(:) @@ -348,7 +355,7 @@ module FatesInterfaceMod real(r8), allocatable :: tgcm_pa(:) ! soil temperature (Kelvin) - real(r8), allocatable :: t_soisno_gl(:) + real(r8), allocatable :: t_soisno_sl(:) ! Canopy Radiation Boundaries ! --------------------------------------------------------------------------------- @@ -389,19 +396,19 @@ module FatesInterfaceMod ! --------------------------------------------------------------------------------- ! Soil suction potential of layers in each site, negative, [mm] - real(r8), allocatable :: smp_gl(:) + real(r8), allocatable :: smp_sl(:) ! Effective porosity = porosity - vol_ic, of layers in each site [-] - real(r8), allocatable :: eff_porosity_gl(:) + real(r8), allocatable :: eff_porosity_sl(:) ! volumetric soil water at saturation (porosity) - real(r8), allocatable :: watsat_gl(:) + real(r8), allocatable :: watsat_sl(:) ! Temperature of ground layers [K] - real(r8), allocatable :: tempk_gl(:) + real(r8), allocatable :: tempk_sl(:) ! Liquid volume in ground layer (m3/m3) - real(r8), allocatable :: h2o_liqvol_gl(:) + real(r8), allocatable :: h2o_liqvol_sl(:) ! Site level filter for uptake response functions logical :: filter_btran @@ -439,10 +446,10 @@ module FatesInterfaceMod ! Logical stating whether a ground layer can have water uptake by plants ! The only condition right now is that liquid water exists ! The name (suction) is used to indicate that soil suction should be calculated - logical, allocatable :: active_suction_gl(:) + logical, allocatable :: active_suction_sl(:) ! Effective fraction of roots in each soil layer - real(r8), allocatable :: rootr_pagl(:,:) + real(r8), allocatable :: rootr_pasl(:,:) ! Integrated (vertically) transpiration wetness factor (0 to 1) ! (diagnostic, should not be used by HLM) @@ -612,7 +619,7 @@ end subroutine fates_clean ! ==================================================================================== - subroutine allocate_bcin(bc_in) + subroutine allocate_bcin(bc_in, nlevsoil_in) ! --------------------------------------------------------------------------------- ! Allocate and Initialze the FATES boundary condition vectors @@ -620,11 +627,16 @@ subroutine allocate_bcin(bc_in) implicit none type(bc_in_type), intent(inout) :: bc_in + integer,intent(in) :: nlevsoil_in ! Allocate input boundaries - allocate(bc_in%zi_sisl(0:hlm_numlevsoil)) - allocate(bc_in%dz_sisl(hlm_numlevsoil)) - allocate(bc_in%z_sisl(hlm_numlevsoil)) + + + bc_in%nlevsoil = nlevsoil_in + + allocate(bc_in%zi_sisl(0:nlevsoil_in)) + allocate(bc_in%dz_sisl(nlevsoil_in)) + allocate(bc_in%z_sisl(nlevsoil_in)) allocate(bc_in%dz_decomp_sisl(hlm_numlevdecomp_full)) @@ -640,11 +652,11 @@ subroutine allocate_bcin(bc_in) allocate(bc_in%solai_parb(maxPatchesPerSite,hlm_numSWb)) ! Hydrology - allocate(bc_in%smp_gl(hlm_numlevgrnd)) - allocate(bc_in%eff_porosity_gl(hlm_numlevgrnd)) - allocate(bc_in%watsat_gl(hlm_numlevgrnd)) - allocate(bc_in%tempk_gl(hlm_numlevgrnd)) - allocate(bc_in%h2o_liqvol_gl(hlm_numlevgrnd)) + allocate(bc_in%smp_sl(nlevsoil_in)) + allocate(bc_in%eff_porosity_sl(nlevsoil_in)) + allocate(bc_in%watsat_sl(nlevsoil_in)) + allocate(bc_in%tempk_sl(nlevsoil_in)) + allocate(bc_in%h2o_liqvol_sl(nlevsoil_in)) ! Photosynthesis allocate(bc_in%filter_photo_pa(maxPatchesPerSite)) @@ -656,7 +668,7 @@ subroutine allocate_bcin(bc_in) allocate(bc_in%rb_pa(maxPatchesPerSite)) allocate(bc_in%t_veg_pa(maxPatchesPerSite)) allocate(bc_in%tgcm_pa(maxPatchesPerSite)) - allocate(bc_in%t_soisno_gl(hlm_numlevgrnd)) + allocate(bc_in%t_soisno_sl(nlevsoil_in)) ! Canopy Radiation allocate(bc_in%filter_vegzen_pa(maxPatchesPerSite)) @@ -670,18 +682,18 @@ subroutine allocate_bcin(bc_in) allocate(bc_in%qflx_transp_pa(maxPatchesPerSite)) allocate(bc_in%swrad_net_pa(maxPatchesPerSite)) allocate(bc_in%lwrad_net_pa(maxPatchesPerSite)) - allocate(bc_in%watsat_sisl(hlm_numlevsoil)) - allocate(bc_in%watres_sisl(hlm_numlevsoil)) - allocate(bc_in%sucsat_sisl(hlm_numlevsoil)) - allocate(bc_in%bsw_sisl(hlm_numlevsoil)) - allocate(bc_in%hksat_sisl(hlm_numlevsoil)) - allocate(bc_in%h2o_liq_sisl(hlm_numlevsoil)); bc_in%h2o_liq_sisl = nan + allocate(bc_in%watsat_sisl(nlevsoil_in)) + allocate(bc_in%watres_sisl(nlevsoil_in)) + allocate(bc_in%sucsat_sisl(nlevsoil_in)) + allocate(bc_in%bsw_sisl(nlevsoil_in)) + allocate(bc_in%hksat_sisl(nlevsoil_in)) + allocate(bc_in%h2o_liq_sisl(nlevsoil_in)); bc_in%h2o_liq_sisl = nan end if return end subroutine allocate_bcin - subroutine allocate_bcout(bc_out) + subroutine allocate_bcout(bc_out, nlevsoil_in) ! --------------------------------------------------------------------------------- ! Allocate and Initialze the FATES boundary condition vectors @@ -689,7 +701,7 @@ subroutine allocate_bcout(bc_out) implicit none type(bc_out_type), intent(inout) :: bc_out - + integer,intent(in) :: nlevsoil_in ! Radiation allocate(bc_out%fsun_pa(maxPatchesPerSite)) @@ -697,8 +709,8 @@ subroutine allocate_bcout(bc_out) allocate(bc_out%laisha_pa(maxPatchesPerSite)) ! Hydrology - allocate(bc_out%active_suction_gl(hlm_numlevgrnd)) - allocate(bc_out%rootr_pagl(maxPatchesPerSite,hlm_numlevgrnd)) + allocate(bc_out%active_suction_sl(nlevsoil_in)) + allocate(bc_out%rootr_pasl(maxPatchesPerSite,nlevsoil_in)) allocate(bc_out%btran_pa(maxPatchesPerSite)) ! Photosynthesis @@ -737,7 +749,7 @@ subroutine allocate_bcout(bc_out) ! Plant-Hydro BC's if (hlm_use_planthydro.eq.itrue) then - allocate(bc_out%qflx_soil2root_sisl(hlm_numlevsoil)) + allocate(bc_out%qflx_soil2root_sisl(nlevsoil_in)) end if return @@ -765,11 +777,11 @@ subroutine zero_bcs(this,s) this%bc_in(s)%solad_parb(:,:) = 0.0_r8 this%bc_in(s)%solai_parb(:,:) = 0.0_r8 - this%bc_in(s)%smp_gl(:) = 0.0_r8 - this%bc_in(s)%eff_porosity_gl(:) = 0.0_r8 - this%bc_in(s)%watsat_gl(:) = 0.0_r8 - this%bc_in(s)%tempk_gl(:) = 0.0_r8 - this%bc_in(s)%h2o_liqvol_gl(:) = 0.0_r8 + this%bc_in(s)%smp_sl(:) = 0.0_r8 + this%bc_in(s)%eff_porosity_sl(:) = 0.0_r8 + this%bc_in(s)%watsat_sl(:) = 0.0_r8 + this%bc_in(s)%tempk_sl(:) = 0.0_r8 + this%bc_in(s)%h2o_liqvol_sl(:) = 0.0_r8 this%bc_in(s)%filter_vegzen_pa(:) = .false. this%bc_in(s)%coszen_pa(:) = 0.0_r8 this%bc_in(s)%albgr_dir_rb(:) = 0.0_r8 @@ -795,11 +807,11 @@ subroutine zero_bcs(this,s) ! Output boundaries - this%bc_out(s)%active_suction_gl(:) = .false. + this%bc_out(s)%active_suction_sl(:) = .false. this%bc_out(s)%fsun_pa(:) = 0.0_r8 this%bc_out(s)%laisun_pa(:) = 0.0_r8 this%bc_out(s)%laisha_pa(:) = 0.0_r8 - this%bc_out(s)%rootr_pagl(:,:) = 0.0_r8 + this%bc_out(s)%rootr_pasl(:,:) = 0.0_r8 this%bc_out(s)%btran_pa(:) = 0.0_r8 this%bc_out(s)%FATES_c_to_litr_lab_c_col(:) = 0.0_r8 @@ -1142,8 +1154,6 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_inir = unset_int hlm_ivis = unset_int hlm_is_restart = unset_int - hlm_numlevgrnd = unset_int - hlm_numlevsoil = unset_int hlm_numlevdecomp_full = unset_int hlm_numlevdecomp = unset_int hlm_name = 'unset' @@ -1269,13 +1279,6 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if(hlm_numlevsoil .eq. unset_int) then - if (fates_global_verbose()) then - write(fates_log(), *) 'FATES dimension/parameter unset: numlevground, exiting' - end if - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - if(hlm_numlevdecomp_full .eq. unset_int) then if (fates_global_verbose()) then write(fates_log(), *) 'FATES dimension/parameter unset: numlevdecomp_full, exiting' @@ -1386,12 +1389,6 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(),*) 'Transfering num_lev_ground = ',ival,' to FATES' end if - case('num_lev_soil') - hlm_numlevsoil = ival - if (fates_global_verbose()) then - write(fates_log(),*) 'Transfering num_lev_ground = ',ival,' to FATES' - end if - case('num_levdecomp_full') hlm_numlevdecomp_full = ival if (fates_global_verbose()) then diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 73a2e35338..e14dee607c 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -252,7 +252,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) call create_patch(sites(s), newpatch, age_init, area_init, & cwd_ag_init, cwd_bg_init, & - leaf_litter_init, root_litter_init ) + leaf_litter_init, root_litter_init, bc_in(s)%nlevsoil ) if( inv_format_list(invsite) == 1 ) then call set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 33c1990b00..cba58621a1 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -1482,7 +1482,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) ! make new patch call create_patch(sites(s), newp, patch_age, area, & cwd_ag_local, cwd_bg_local, & - leaf_litter_local, root_litter_local) + leaf_litter_local, root_litter_local,bc_in(s)%nlevsoil ) ! give this patch a unique patch number newp%patchno = idx_pa From a4600f70d4b2644c44793dd196619989f22ecd3f Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 9 May 2018 14:33:32 -0700 Subject: [PATCH 02/12] Keeping hlm_numlevgrnd for time being. --- main/FatesInterfaceMod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 3fce4d3a24..47399fb0b4 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1156,6 +1156,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_is_restart = unset_int hlm_numlevdecomp_full = unset_int hlm_numlevdecomp = unset_int + hlm_numlevgrnd = unset_int hlm_name = 'unset' hlm_hio_ignore_val = unset_double hlm_masterproc = unset_int From c6583510c7ede39acdd4ab62abcb10e4b6897db4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 10 May 2018 12:20:02 -0700 Subject: [PATCH 03/12] Added site specific nlevdecomp. Changing root profiles to be modular. --- biogeochem/EDCanopyStructureMod.F90 | 2 +- biogeochem/EDPatchDynamicsMod.F90 | 115 ++++++++++++++++++--- biogeochem/EDPhysiologyMod.F90 | 104 ++++++++++--------- biogeophys/FatesPlantHydraulicsMod.F90 | 1 - biogeophys/FatesPlantRespPhotosynthMod.F90 | 3 +- main/FatesInterfaceMod.F90 | 66 ++++-------- 6 files changed, 173 insertions(+), 118 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index af2d033c64..2ea1ce0179 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -920,7 +920,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) do while(associated(currentPatch)) - call set_root_fraction(currentPatch,bc_in(s)%zi_sisl) + call set_root_fraction(currentPatch%rootfr_ft,bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) !zero cohort-summed variables. currentPatch%total_canopy_area = 0.0_r8 diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 2364407bf7..6555369d54 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1143,8 +1143,8 @@ subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_l allocate(new_patch%fabi(hlm_numSWb)) allocate(new_patch%sabs_dir(hlm_numSWb)) allocate(new_patch%sabs_dif(hlm_numSWb)) - allocate(new_patch%rootfr_ft(numpft,nlevsoil) - allocate(new_patch%rootr_ft(numpft,nlevsoil) + allocate(new_patch%rootfr_ft(numpft,nlevsoil)) + allocate(new_patch%rootr_ft(numpft,nlevsoil)) call zero_patch(new_patch) !The nan value in here is not working?? @@ -1891,7 +1891,7 @@ end function countPatches ! ==================================================================================== - subroutine set_root_fraction( root_fraction , zi ) + subroutine set_root_fraction( root_fraction , zi, lowerb ) ! ! !DESCRIPTION: ! Calculates the fractions of the root biomass in each layer for each pft. @@ -1904,42 +1904,123 @@ subroutine set_root_fraction( root_fraction , zi ) ! ! !ARGUMENTS real(r8),intent(out) :: root_fraction(:,:) - real(r8),intent(in) :: zi(:) - ! - ! !LOCAL VARIABLES: - integer :: lev,p,c,ft - integer :: nlevsoil - real(r8) :: sum_rootfr + real(r8),intent(in) :: zi(lowerb:) + integer,intent(in) :: lowerb + + ! Parameters + integer, parameter :: exponential_1p_profile_type = 1 + integer, parameter :: jackson_beta_profile_type = 2 + integer, parameter :: exponential_2p_profile_type = 3 + + integer, parameter :: root_profile_type = exponential_2p_profile_type + !---------------------------------------------------------------------- if(lbound(zi,1).ne.0) then + write(fates_log(),*) 'lbound:',lbound(zi) + write(fates_log(),*) 'ubound:',ubound(zi) write(fates_log(),*) 'layer interface levels should have 0 index' call endrun(msg=errMsg(sourcefile, __LINE__)) end if + + select case(root_profile_type) + case ( exponential_1p_profile_type ) + call exponential_1p_root_profile(root_fraction , zi, lowerb) + case ( jackson_beta_profile_type ) + call jackson_beta_root_profile(root_fraction , zi, lowerb) + case ( exponential_2p_profile_type ) + call exponential_2p_root_profile(root_fraction , zi, lowerb) + case default + write(fates_log(),*) 'An undefined root profile type was specified' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + + end subroutine set_root_fraction + ! ===================================================================================== + + subroutine exponential_2p_root_profile(root_fraction , zi, lowerb ) + ! + ! !ARGUMENTS + real(r8),intent(out) :: root_fraction(:,:) + real(r8),intent(in) :: zi(lowerb:) + integer,intent(in) :: lowerb + + ! Locals + + integer :: nlevsoil ! Number of soil layers + integer :: ft ! pft index + integer :: lev ! soil layer index + real(r8) :: sum_rootfr ! sum of root fraction for normalization + nlevsoil = ubound(zi,1) + root_fraction = 0.0_r8 + do ft = 1,numpft sum_rootfr = 0.0_r8 - do lev = 1, nlevsoil root_fraction(ft,lev) = .5_r8*( & - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev-1)) & - + exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev-1)) & - - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev)) & - - exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev))) + exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev-1)) & + + exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev-1)) & + - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev)) & + - exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev))) - sum_rootfr = sum_rootfr + cpatch%rootfr_ft(ft,lev) + sum_rootfr = sum_rootfr + root_fraction(ft,lev) end do ! Normalize the root profile - root_fraction(ft,:) = cpatch%rootfr_ft(ft,:)/sum_rootfr + root_fraction(ft,:) = root_fraction(ft,:)/sum_rootfr end do + return + end subroutine exponential_2p_root_profile + + ! ===================================================================================== + + subroutine exponential_1p_root_profile(root_fraction , zi, lowerb) + + ! + ! !ARGUMENTS + real(r8),intent(out) :: root_fraction(:,:) + real(r8),intent(in) :: zi(lowerb:) + integer,intent(in) :: lowerb + + ! + ! LOCAL VARIABLES: + integer :: ft ! pft index + integer :: lev ! soil depth layer index + integer :: nlevsoil ! + + real(r8), parameter :: rootprof_exp = 3. ! how steep profile is + ! for root C inputs (1/ e-folding depth) (1/m) + + nlevsoil = ubound(zi,1) + + ! define rooting profile from exponential parameters + do ft = 1, numpft + depth = 0.0_r8 + do j = 1, nlevsoil + depth = depth + 0.5*(zi(j)+zi(j-1)) + root_fraction(ft,j) = exp(-rootprof_exp * depth) + end do + end do + + ! Normalize the root profile + root_fraction(ft,:) = root_fraction(ft,:)/sum_rootfr + + + return + end subroutine exponential_1p_root_profile - end subroutine set_root_fraction + + + + end subroutine jackson_beta_root_profile + end module EDPatchDynamicsMod diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 5bdc6c9001..f96bdead63 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -2005,11 +2005,11 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) use EDTypesMod, only : AREA - use FatesInterfaceMod, only : hlm_numlevdecomp use EDPftvarcon, only : EDPftvarcon_inst use FatesConstantsMod, only : sec_per_day use FatesInterfaceMod, only : bc_in_type, bc_out_type use FatesInterfaceMod, only : hlm_use_vertsoilc + use FatesInterfaceMod, only : hlm_numlevgrnd use FatesConstantsMod, only : itrue use FatesGlobals, only : endrun => fates_endrun use EDParamsMod , only : ED_val_cwd_flig, ED_val_cwd_fcel @@ -2045,12 +2045,12 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) real(r8) :: delta ! NOTE(bja, 201608) these were removed from clm in clm4_5_10_r187 - logical, parameter :: exponential_rooting_profile = .true. - logical, parameter :: pftspecific_rootingprofile = .true. +! logical, parameter :: exponential_rooting_profile = .true. +! logical, parameter :: pftspecific_rootingprofile = .true. ! NOTE(bja, 201608) as of clm4_5_10_r187 rootprof_exp is now a ! private function level parameter in RootBiophysMod.F90::exponential_rootfr() - real(r8), parameter :: rootprof_exp = 3. ! how steep profile is +! real(r8), parameter :: rootprof_exp = 3. ! how steep profile is ! for root C inputs (1/ e-folding depth) (1/m) ! NOTE(rgk, 201705) this parameter was brought over from SoilBiogeochemVerticalProfile @@ -2066,10 +2066,10 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! Doing so will be answer changing though so perhaps easiest to do this in steps. integer, parameter :: rooting_profile_varindex_water = 1 - real(r8) :: leaf_prof(1:nsites, 1:hlm_numlevdecomp) - real(r8) :: froot_prof(1:nsites, 1:maxpft, 1:hlm_numlevdecomp) - real(r8) :: croot_prof(1:nsites, 1:hlm_numlevdecomp) - real(r8) :: stem_prof(1:nsites, 1:hlm_numlevdecomp) + real(r8) :: leaf_prof(1:nsites, 1:hlm_numlevgrnd) + real(r8) :: froot_prof(1:nsites, 1:maxpft, 1:hlm_numlevgrnd) + real(r8) :: croot_prof(1:nsites, 1:hlm_numlevgrnd) + real(r8) :: stem_prof(1:nsites, 1:hlm_numlevgrnd) delta = 0.001_r8 !no of seconds in a year. @@ -2092,59 +2092,61 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) if (hlm_use_vertsoilc == itrue) then ! initialize profiles to zero - leaf_prof(1:nsites, :) = 0._r8 - froot_prof(1:nsites, 1:maxpft, :) = 0._r8 - stem_prof(1:nsites, :) = 0._r8 + leaf_prof(1:nsites, 1:hlm_numlevgrnd ) = 0._r8 + froot_prof(1:nsites, 1:maxpft, 1:hlm_numlevgrnd) = 0._r8 + stem_prof(1:nsites, 1:hlm_numlevgrnd) = 0._r8 do s = 1,nsites + ! define a single shallow surface profile for surface additions (leaves, stems, and N deposition) surface_prof(:) = 0._r8 - do j = 1, hlm_numlevdecomp + do j = 1, bc_in(s)%nlevdecomp surface_prof(j) = exp(-surfprof_exp * bc_in(s)%z_sisl(j)) / bc_in(s)%dz_decomp_sisl(j) end do + ! ----------------------------------------------------------------------------- + ! This is the rooting profile. cinput_rootfr + ! This array will calculate root + ! mass as far down as the soil column goes. It is possible + ! that the active layers are not as deep as the roots go. + ! That is ok, the roots in the active layers will be talied up and + ! normalized. + ! ----------------------------------------------------------------------------- + cinput_rootfr(:,:) = 0._r8 + + + ! calculate pft-specific rooting profiles in the absence of permafrost or bedrock limitations if ( exponential_rooting_profile ) then if ( .not. pftspecific_rootingprofile ) then ! define rooting profile from exponential parameters do ft = 1, numpft - do j = 1, hlm_numlevdecomp - cinput_rootfr(ft,j) = exp(-rootprof_exp * bc_in(s)%z_sisl(j)) / bc_in(s)%dz_decomp_sisl(j) + do j = 1, bc_in(s)%nlevdecomp + cinput_rootfr(ft,j) = exp(-rootprof_exp * bc_in(s)%z_sisl(j)) end do end do else ! use beta distribution parameter from Jackson et al., 1996 do ft = 1, numpft - do j = 1, hlm_numlevdecomp + do j = 1, bc_in(s)%nlevdecomp cinput_rootfr(ft,j) = & ( EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & (bc_in(s)%zi_sisl(j-1)*100._r8) - & EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & - (bc_in(s)%zi_sisl(j)*100._r8) ) & - / bc_in(s)%dz_decomp_sisl(j) + (bc_in(s)%zi_sisl(j)*100._r8) ) end do end do endif else + ! This generates a rooting profile over the whole soil column + call set_root_fraction(cinput_rootfr(:,1:bc_in(s)%nlevsoil), & + bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) - - - do ft = 1,numpft - do j = 1, hlm_numlevdecomp - ! use standard CLM root fraction profiles; - cinput_rootfr(ft,j) = ( .5_r8*( & - exp(-EDPftvarcon_inst%roota_par(ft) * bc_in(s)%zi_sisl(j-1)) & - + exp(-EDPftvarcon_inst%rootb_par(ft) * bc_in(s)%zi_sisl(j-1)) & - - exp(-EDPftvarcon_inst%roota_par(ft) * bc_in(s)%zi_sisl(j)) & - - exp(-EDPftvarcon_inst%rootb_par(ft) * bc_in(s)%zi_sisl(j)))) & - / bc_in(s)%dz_decomp_sisl(j) - end do - end do endif - + ! ! now add permafrost constraint: integrate rootfr over active layer of soil site, ! truncate below permafrost or bedrock table where present, and rescale so that integral = 1 @@ -2152,12 +2154,13 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) surface_prof_tot = 0._r8 ! - do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), hlm_numlevdecomp) + do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), bc_in(s)%nlevdecomp) surface_prof_tot = surface_prof_tot + surface_prof(j) * bc_in(s)%dz_decomp_sisl(j) end do + do ft = 1,numpft - do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), hlm_numlevdecomp) - rootfr_tot(ft) = rootfr_tot(ft) + cinput_rootfr(ft,j) * bc_in(s)%dz_decomp_sisl(j) + do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), bc_in(s)%nlevdecomp) + rootfr_tot(ft) = rootfr_tot(ft) + cinput_rootfr(ft,j) end do end do ! @@ -2166,7 +2169,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) if ( (bc_in(s)%max_rooting_depth_index_col > 0) .and. (rootfr_tot(ft) > 0._r8) ) then ! where there is not permafrost extending to the surface, integrate the profiles over the active layer ! this is equivalent to integrating over all soil layers outside of permafrost regions - do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), hlm_numlevdecomp) + do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), bc_in(s)%nlevdecomp) froot_prof(s,ft,j) = cinput_rootfr(ft,j) / rootfr_tot(ft) end do else @@ -2174,12 +2177,13 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) froot_prof(s,ft,1) = 1._r8/bc_in(s)%dz_decomp_sisl(1) endif end do + ! ! rescale the shallow profiles if ( (bc_in(s)%max_rooting_depth_index_col > 0) .and. (surface_prof_tot > 0._r8) ) then ! where there is not permafrost extending to the surface, integrate the profiles over the active layer ! this is equivalent to integrating over all soil layers outside of permafrost regions - do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), hlm_numlevdecomp) + do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), bc_in(s)%nlevdecomp) ! set all surface processes to shallower profile leaf_prof(s,j) = surface_prof(j)/ surface_prof_tot stem_prof(s,j) = surface_prof(j)/ surface_prof_tot @@ -2188,7 +2192,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! if fully frozen, or no roots, put everything in the top layer leaf_prof(s,1) = 1._r8/bc_in(s)%dz_decomp_sisl(1) stem_prof(s,1) = 1._r8/bc_in(s)%dz_decomp_sisl(1) - do j = 2, hlm_numlevdecomp + do j = 2, bc_in(s)%nlevdecomp leaf_prof(s,j) = 0._r8 stem_prof(s,j) = 0._r8 end do @@ -2209,7 +2213,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! check the leaf and stem profiles leaf_prof_sum = 0._r8 stem_prof_sum = 0._r8 - do j = 1, hlm_numlevdecomp + do j = 1, bc_in(s)%nlevdecomp leaf_prof_sum = leaf_prof_sum + leaf_prof(s,j) * bc_in(s)%dz_decomp_sisl(j) stem_prof_sum = stem_prof_sum + stem_prof(s,j) * bc_in(s)%dz_decomp_sisl(j) end do @@ -2226,7 +2230,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! now check each fine root profile do ft = 1,numpft froot_prof_sum = 0._r8 - do j = 1, hlm_numlevdecomp + do j = 1, bc_in(s)%nlevdecomp froot_prof_sum = froot_prof_sum + froot_prof(s,ft,j) * bc_in(s)%dz_decomp_sisl(j) end do if ( ( abs(froot_prof_sum - 1._r8) > delta ) ) then @@ -2238,7 +2242,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! zero the site-level C input variables do s = 1, nsites - do j = 1, hlm_numlevdecomp + do j = 1, bc_in(s)%nlevdecomp bc_out(s)%FATES_c_to_litr_lab_c_col(j) = 0._r8 bc_out(s)%FATES_c_to_litr_cel_c_col(j) = 0._r8 bc_out(s)%FATES_c_to_litr_lig_c_col(j) = 0._r8 @@ -2276,14 +2280,14 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) biomass_bg_tot = biomass_bg_tot + biomass_bg_ft(ft) end do ! - do j = 1, hlm_numlevdecomp - ! zero this for each patch - croot_prof_perpatch(j) = 0._r8 - end do + + ! zero this for each patch + croot_prof_perpatch(1:bc_in(s)%nlevdecomp) = 0._r8 + ! if ( biomass_bg_tot .gt. 0._r8) then do ft = 1,numpft - do j = 1, hlm_numlevdecomp + do j = 1, bc_in(s)%nlevdecomp croot_prof_perpatch(j) = croot_prof_perpatch(j) + froot_prof(s,ft,j) * biomass_bg_ft(ft) / biomass_bg_tot end do end do @@ -2293,7 +2297,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! ! add croot_prof as weighted average (weighted by patch area) of croot_prof_perpatch - do j = 1, hlm_numlevdecomp + do j = 1, bc_in(s)%nlevdecomp croot_prof(s, j) = croot_prof(s, j) + croot_prof_perpatch(j) * currentPatch%area / AREA end do ! @@ -2310,7 +2314,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! ! ! CWD pools fragmenting into decomposing litter pools. do ci = 1, ncwd - do j = 1, hlm_numlevdecomp + do j = 1, bc_in(s)%nlevdecomp bc_out(s)%FATES_c_to_litr_cel_c_col(j) = bc_out(s)%FATES_c_to_litr_cel_c_col(j) + & currentpatch%CWD_AG_out(ci) * ED_val_cwd_fcel * currentpatch%area/AREA * stem_prof(s,j) bc_out(s)%FATES_c_to_litr_lig_c_col(j) = bc_out(s)%FATES_c_to_litr_lig_c_col(j) + & @@ -2325,7 +2329,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! leaf and fine root pools. do ft = 1,numpft - do j = 1, hlm_numlevdecomp + do j = 1, bc_in(s)%nlevdecomp bc_out(s)%FATES_c_to_litr_lab_c_col(j) = bc_out(s)%FATES_c_to_litr_lab_c_col(j) + & currentpatch%leaf_litter_out(ft) * EDPftvarcon_inst%lf_flab(ft) * currentpatch%area/AREA * leaf_prof(s,j) bc_out(s)%FATES_c_to_litr_cel_c_col(j) = bc_out(s)%FATES_c_to_litr_cel_c_col(j) + & @@ -2357,7 +2361,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) end do ! do sites(s) do s = 1, nsites - do j = 1, hlm_numlevdecomp + do j = 1, bc_in(s)%nlevdecomp ! time unit conversion bc_out(s)%FATES_c_to_litr_lab_c_col(j)=bc_out(s)%FATES_c_to_litr_lab_c_col(j) * mass_convert / time_convert bc_out(s)%FATES_c_to_litr_cel_c_col(j)=bc_out(s)%FATES_c_to_litr_cel_c_col(j) * mass_convert / time_convert @@ -2368,7 +2372,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! write(fates_log(),*)'cdk FATES_c_to_litr_lab_c: ', FATES_c_to_litr_lab_c ! write_col(fates_log(),*)'cdk FATES_c_to_litr_cel_c: ', FATES_c_to_litr_cel_c ! write_col(fates_log(),*)'cdk FATES_c_to_litr_lig_c: ', FATES_c_to_litr_lig_c - ! write_col(fates_log(),*)'cdk hlm_numlevdecomp_full, bounds%begc, bounds%endc: ', hlm_numlevdecomp_full, bounds%begc, bounds%endc + ! write_col(fates_log(),*)'cdk bounds%begc, bounds%endc: ', bounds%begc, bounds%endc ! write(fates_log(),*)'cdk leaf_prof: ', leaf_prof ! write(fates_log(),*)'cdk stem_prof: ', stem_prof ! write(fates_log(),*)'cdk froot_prof: ', froot_prof diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index e68291057a..36336adc13 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -22,7 +22,6 @@ module FatesPlantHydraulicsMod use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : bc_out_type - use FatesInterfaceMod , only : hlm_numlevsoil use FatesHydraulicsMemMod, only: ed_site_hydr_type use FatesHydraulicsMemMod, only: ed_patch_hydr_type diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index e0d68ec76b..d96fbe8503 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -302,7 +302,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) end do !ft - call set_root_fraction(currentPatch,bc_in(s)%zi_sisl) + call set_root_fraction(currentPatch%rootfr_ft, & + bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) ! ------------------------------------------------------------------------ ! Part VI: Loop over all leaf layers. diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 47399fb0b4..4439127571 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -65,17 +65,6 @@ module FatesInterfaceMod ! NOTE! SOIL LAYERS ARE NOT A GLOBAL, THEY ! ARE VARIABLE BY SITE - - integer, protected :: hlm_numlevdecomp_full ! Number of GROUND layers for the purposes - ! of biogeochemistry; can be either 1 - ! or the total number of soil layers - ! (includes bedrock) - - - integer, protected :: hlm_numlevdecomp ! Number of SOIL layers for the purposes of - ! biogeochemistry; can be either 1 or the total - ! number of soil layers - integer, protected :: hlm_is_restart ! Is the HLM signalling that this is a restart ! type simulation? ! 1=TRUE, 0=FALSE @@ -270,6 +259,8 @@ module FatesInterfaceMod ! Soil layer structure integer :: nlevsoil ! the number of soil layers in this column + integer :: nlevdecomp ! the number of soil layers in the column + ! that are biogeochemically active real(r8),allocatable :: zi_sisl(:) ! interface level below a "z" level (m) ! this contains a zero index for surface. real(r8),allocatable :: dz_sisl(:) ! layer thickness (m) @@ -423,8 +414,8 @@ module FatesInterfaceMod real(r8),allocatable :: lwrad_net_pa(:) ! Net absorbed longwave radiation (W/m2) real(r8),allocatable :: watsat_sisl(:) ! volumetric soil water at saturation (porosity) real(r8),allocatable :: watres_sisl(:) ! volumetric residual soil water - real(r8),allocatable :: sucsat_sisl(:) ! minimum soil suction (mm) (hlm_nlevsoil) - real(r8),allocatable :: bsw_sisl(:) ! Clapp and Hornberger "b" (hlm_nlevsoil) + real(r8),allocatable :: sucsat_sisl(:) ! minimum soil suction (mm) + real(r8),allocatable :: bsw_sisl(:) ! Clapp and Hornberger "b" real(r8),allocatable :: hksat_sisl(:) ! hydraulic conductivity at saturation (mm H2O /s) real(r8),allocatable :: h2o_liq_sisl(:) ! Liquid water mass in each layer (kg/m2) real(r8) :: smpmin_si ! restriction for min of soil potential (mm) @@ -619,7 +610,7 @@ end subroutine fates_clean ! ==================================================================================== - subroutine allocate_bcin(bc_in, nlevsoil_in) + subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in) ! --------------------------------------------------------------------------------- ! Allocate and Initialze the FATES boundary condition vectors @@ -628,17 +619,19 @@ subroutine allocate_bcin(bc_in, nlevsoil_in) implicit none type(bc_in_type), intent(inout) :: bc_in integer,intent(in) :: nlevsoil_in + integer,intent(in) :: nlevdecomp_in ! Allocate input boundaries - bc_in%nlevsoil = nlevsoil_in + bc_in%nlevsoil = nlevsoil_in + bc_in%nlevdecomp = nlevdecomp_in allocate(bc_in%zi_sisl(0:nlevsoil_in)) allocate(bc_in%dz_sisl(nlevsoil_in)) allocate(bc_in%z_sisl(nlevsoil_in)) - allocate(bc_in%dz_decomp_sisl(hlm_numlevdecomp_full)) + allocate(bc_in%dz_decomp_sisl(nlevdecomp_in)) ! Vegetation Dynamics allocate(bc_in%t_veg24_pa(maxPatchesPerSite)) @@ -693,7 +686,7 @@ subroutine allocate_bcin(bc_in, nlevsoil_in) return end subroutine allocate_bcin - subroutine allocate_bcout(bc_out, nlevsoil_in) + subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in) ! --------------------------------------------------------------------------------- ! Allocate and Initialze the FATES boundary condition vectors @@ -702,6 +695,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in) implicit none type(bc_out_type), intent(inout) :: bc_out integer,intent(in) :: nlevsoil_in + integer,intent(in) :: nlevdecomp_in ! Radiation allocate(bc_out%fsun_pa(maxPatchesPerSite)) @@ -728,9 +722,9 @@ subroutine allocate_bcout(bc_out, nlevsoil_in) allocate(bc_out%ftii_parb(maxPatchesPerSite,hlm_numSWb)) ! biogeochemistry - allocate(bc_out%FATES_c_to_litr_lab_c_col(hlm_numlevdecomp_full)) - allocate(bc_out%FATES_c_to_litr_cel_c_col(hlm_numlevdecomp_full)) - allocate(bc_out%FATES_c_to_litr_lig_c_col(hlm_numlevdecomp_full)) + allocate(bc_out%FATES_c_to_litr_lab_c_col(nlevdecomp_in)) + allocate(bc_out%FATES_c_to_litr_cel_c_col(nlevdecomp_in)) + allocate(bc_out%FATES_c_to_litr_lig_c_col(nlevdecomp_in)) ! Canopy Structure allocate(bc_out%elai_pa(maxPatchesPerSite)) @@ -764,6 +758,10 @@ subroutine zero_bcs(this,s) integer, intent(in) :: s ! Input boundaries + ! Warning: these "z" type variables + ! are written only once at the beginning + ! so THIS ROUTINE SHOULD NOT BE CALLED AFTER + ! INITIALIZATION this%bc_in(s)%zi_sisl(:) = 0.0_r8 this%bc_in(s)%dz_sisl(:) = 0.0_r8 this%bc_in(s)%z_sisl(:) = 0.0_r8 @@ -1154,8 +1152,6 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_inir = unset_int hlm_ivis = unset_int hlm_is_restart = unset_int - hlm_numlevdecomp_full = unset_int - hlm_numlevdecomp = unset_int hlm_numlevgrnd = unset_int hlm_name = 'unset' hlm_hio_ignore_val = unset_double @@ -1280,20 +1276,6 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if(hlm_numlevdecomp_full .eq. unset_int) then - if (fates_global_verbose()) then - write(fates_log(), *) 'FATES dimension/parameter unset: numlevdecomp_full, exiting' - end if - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - if(hlm_numlevdecomp .eq. unset_int) then - if (fates_global_verbose()) then - write(fates_log(), *) 'FATES dimension/parameter unset: numlevdecomp, exiting' - end if - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - if(trim(hlm_name) .eq. 'unset') then if (fates_global_verbose()) then write(fates_log(),*) 'FATES dimension/parameter unset: hlm_name, exiting' @@ -1390,18 +1372,6 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(),*) 'Transfering num_lev_ground = ',ival,' to FATES' end if - case('num_levdecomp_full') - hlm_numlevdecomp_full = ival - if (fates_global_verbose()) then - write(fates_log(),*) 'Transfering num_levdecomp_full = ',ival,' to FATES' - end if - - case('num_levdecomp') - hlm_numlevdecomp = ival - if (fates_global_verbose()) then - write(fates_log(),*) 'Transfering num_levdecomp = ',ival,' to FATES' - end if - case('soilwater_ipedof') hlm_ipedof = ival if (fates_global_verbose()) then From 704f365e755482619af67b59956c4d7f5ec297ba Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 10 May 2018 13:03:02 -0700 Subject: [PATCH 04/12] Modularizing root depth profiles. --- biogeochem/EDCanopyStructureMod.F90 | 5 +- biogeochem/EDPatchDynamicsMod.F90 | 117 ++++++++++++--------- biogeochem/EDPhysiologyMod.F90 | 38 ++----- biogeophys/FatesPlantRespPhotosynthMod.F90 | 8 +- 4 files changed, 87 insertions(+), 81 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 2ea1ce0179..186d55b8b2 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -920,7 +920,10 @@ subroutine canopy_summarization( nsites, sites, bc_in ) do while(associated(currentPatch)) - call set_root_fraction(currentPatch%rootfr_ft,bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) + + + call set_root_fraction(currentPatch%rootfr_ft(ft,1:bc_in(s)%nlevsoil), ft, & + bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) !zero cohort-summed variables. currentPatch%total_canopy_area = 0.0_r8 diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 6555369d54..0429f7e77b 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1891,7 +1891,7 @@ end function countPatches ! ==================================================================================== - subroutine set_root_fraction( root_fraction , zi, lowerb ) + subroutine set_root_fraction(root_fraction, ft, zi, lowerb ) ! ! !DESCRIPTION: ! Calculates the fractions of the root biomass in each layer for each pft. @@ -1903,7 +1903,8 @@ subroutine set_root_fraction( root_fraction , zi, lowerb ) ! ! !ARGUMENTS - real(r8),intent(out) :: root_fraction(:,:) + real(r8),intent(inout) :: root_fraction(:) + integer, intent(in) :: ft real(r8),intent(in) :: zi(lowerb:) integer,intent(in) :: lowerb @@ -1925,76 +1926,69 @@ subroutine set_root_fraction( root_fraction , zi, lowerb ) select case(root_profile_type) case ( exponential_1p_profile_type ) - call exponential_1p_root_profile(root_fraction , zi, lowerb) + call exponential_1p_root_profile(root_fraction, ft, zi) case ( jackson_beta_profile_type ) - call jackson_beta_root_profile(root_fraction , zi, lowerb) + call jackson_beta_root_profile(root_fraction, ft, zi) case ( exponential_2p_profile_type ) - call exponential_2p_root_profile(root_fraction , zi, lowerb) + call exponential_2p_root_profile(root_fraction, ft, zi) case default write(fates_log(),*) 'An undefined root profile type was specified' write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end select - + + return end subroutine set_root_fraction ! ===================================================================================== - subroutine exponential_2p_root_profile(root_fraction , zi, lowerb ) + subroutine exponential_2p_root_profile(root_fraction, ft, zi ) ! ! !ARGUMENTS - real(r8),intent(out) :: root_fraction(:,:) - real(r8),intent(in) :: zi(lowerb:) - integer,intent(in) :: lowerb + real(r8),intent(out) :: root_fraction(:) + integer,intent(in) :: ft + real(r8),intent(in) :: zi(0:) ! Locals - integer :: nlevsoil ! Number of soil layers - integer :: ft ! pft index integer :: lev ! soil layer index real(r8) :: sum_rootfr ! sum of root fraction for normalization nlevsoil = ubound(zi,1) - root_fraction = 0.0_r8 - - do ft = 1,numpft - - sum_rootfr = 0.0_r8 - do lev = 1, nlevsoil - - root_fraction(ft,lev) = .5_r8*( & - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev-1)) & - + exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev-1)) & - - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev)) & - - exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev))) - - sum_rootfr = sum_rootfr + root_fraction(ft,lev) - - end do - - ! Normalize the root profile - root_fraction(ft,:) = root_fraction(ft,:)/sum_rootfr + sum_rootfr = 0.0_r8 + do lev = 1, nlevsoil + root_fraction(lev) = .5_r8*( & + exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev-1)) & + + exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev-1)) & + - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev)) & + - exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev))) + sum_rootfr = sum_rootfr + root_fraction(lev) end do + + ! Normalize the root profile + root_fraction(1:nlevsoil) = root_fraction(1:nlevsoil)/sum_rootfr + return end subroutine exponential_2p_root_profile ! ===================================================================================== - subroutine exponential_1p_root_profile(root_fraction , zi, lowerb) + subroutine exponential_1p_root_profile(root_fraction, ft, zi) ! ! !ARGUMENTS - real(r8),intent(out) :: root_fraction(:,:) - real(r8),intent(in) :: zi(lowerb:) - integer,intent(in) :: lowerb - + real(r8),intent(out) :: root_fraction(:) + integer,intent(in) :: ft + real(r8),intent(in) :: zi(0:) + ! ! LOCAL VARIABLES: - integer :: ft ! pft index - integer :: lev ! soil depth layer index - integer :: nlevsoil ! + integer :: lev ! soil depth layer index + integer :: nlevsoil ! number of soil layers + real(r8) :: depth ! Depth to middle of layer [m] + real(r8) :: sum_rootfr ! sum of rooting profile for normalization real(r8), parameter :: rootprof_exp = 3. ! how steep profile is ! for root C inputs (1/ e-folding depth) (1/m) @@ -2002,24 +1996,53 @@ subroutine exponential_1p_root_profile(root_fraction , zi, lowerb) nlevsoil = ubound(zi,1) ! define rooting profile from exponential parameters - do ft = 1, numpft - depth = 0.0_r8 - do j = 1, nlevsoil - depth = depth + 0.5*(zi(j)+zi(j-1)) - root_fraction(ft,j) = exp(-rootprof_exp * depth) - end do + sum_rootfr = 0.0_r8 + do lev = 1, nlevsoil + root_fraction(lev) = exp(-rootprof_exp * 0.5*(zi(lev)+zi(lev-1)) ) + sum_rootfr = sum_rootfr + root_fraction(lev) end do - + ! Normalize the root profile - root_fraction(ft,:) = root_fraction(ft,:)/sum_rootfr + root_fraction(1:nlevsoil) = root_fraction(1:nlevsoil)/sum_rootfr return end subroutine exponential_1p_root_profile + ! ===================================================================================== + subroutine jackson_beta_root_profile(root_fraction, ft, zi) + + ! !ARGUMENTS + real(r8),intent(out) :: root_fraction(:) ! fraction of root mass in each soil layer + integer,intent(in) :: ft ! functional type + real(r8),intent(in) :: zi(0:) ! depth of layer interfaces 0-nlevsoil + + ! + ! LOCAL VARIABLES: + integer :: lev ! soil depth layer index + integer :: nlevsoil ! number of soil layers + real(r8) :: sum_rootfr ! sum of rooting profile, for normalization + integer, parameter :: rooting_profile_varindex_water = 1 + + nlevsoil = ubound(zi,1) + ! use beta distribution parameter from Jackson et al., 1996 + sum_rootfr = 0.0_r8 + do lev = 1, nlevsoil + root_fraction(lev) = & + ( EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & + ( zi(lev-1)*100._r8) - & + EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & + ( zi(lev)*100._r8) ) + sum_rootfr = sum_rootfr + root_fraction(lev) + end do + + ! Normalize the root profile + root_fraction(1:nlevsoil) = root_fraction(1:nlevsoil)/sum_rootfr + + return end subroutine jackson_beta_root_profile diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index f96bdead63..02df97f47d 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -2013,6 +2013,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) use FatesConstantsMod, only : itrue use FatesGlobals, only : endrun => fates_endrun use EDParamsMod , only : ED_val_cwd_flig, ED_val_cwd_fcel + use EDPatchDynamicsMod, only : set_root_fraction implicit none @@ -2114,38 +2115,13 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! ----------------------------------------------------------------------------- cinput_rootfr(:,:) = 0._r8 - - - - - ! calculate pft-specific rooting profiles in the absence of permafrost or bedrock limitations - if ( exponential_rooting_profile ) then - if ( .not. pftspecific_rootingprofile ) then - ! define rooting profile from exponential parameters - do ft = 1, numpft - do j = 1, bc_in(s)%nlevdecomp - cinput_rootfr(ft,j) = exp(-rootprof_exp * bc_in(s)%z_sisl(j)) - end do - end do - else - ! use beta distribution parameter from Jackson et al., 1996 - do ft = 1, numpft - do j = 1, bc_in(s)%nlevdecomp - cinput_rootfr(ft,j) = & - ( EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & - (bc_in(s)%zi_sisl(j-1)*100._r8) - & - EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & - (bc_in(s)%zi_sisl(j)*100._r8) ) - end do - end do - endif - else - - ! This generates a rooting profile over the whole soil column - call set_root_fraction(cinput_rootfr(:,1:bc_in(s)%nlevsoil), & - bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) + do ft = 1, numpft - endif + ! This generates a rooting profile over the whole soil column + call set_root_fraction(cinput_rootfr(ft,1:bc_in(s)%nlevsoil), ft, & + bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) + + end do ! ! now add permafrost constraint: integrate rootfr over active layer of soil site, diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index d96fbe8503..9bcd36f846 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -299,11 +299,15 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) else kn(ft) = exp(0.00963_r8 * EDPftvarcon_inst%vcmax25top(ft) - 2.43_r8) end if + + ! This is probably unnecessary and already calculated (RGK 05-2018) + call set_root_fraction(currentPatch%rootfr_ft(ft,1:bc_in(s)%nlevsoil), ft, & + bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) + end do !ft - call set_root_fraction(currentPatch%rootfr_ft, & - bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) + ! ------------------------------------------------------------------------ ! Part VI: Loop over all leaf layers. From f3df7534d18006098c754488ac99416cedf930c1 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 11 May 2018 15:18:20 -0700 Subject: [PATCH 05/12] Moved root profiling to the allometry module. Added context switches --- biogeochem/EDCanopyStructureMod.F90 | 8 +- biogeochem/EDPatchDynamicsMod.F90 | 158 ----------------- biogeochem/EDPhysiologyMod.F90 | 8 +- biogeochem/FatesAllometryMod.F90 | 189 +++++++++++++++++++++ biogeophys/FatesPlantRespPhotosynthMod.F90 | 17 +- 5 files changed, 209 insertions(+), 171 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 186d55b8b2..d86a2b387a 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -882,7 +882,8 @@ subroutine canopy_summarization( nsites, sites, bc_in ) use FatesInterfaceMod , only : bc_in_type use EDPatchDynamicsMod , only : set_patchno - use EDPatchDynamicsMod , only : set_root_fraction + use FatesAllometryMod , only : set_root_fraction + use FatesAllometryMod , only : i_hydro_rootprof_context use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index use EDtypesMod , only : area use EDPftvarcon , only : EDPftvarcon_inst @@ -920,10 +921,9 @@ subroutine canopy_summarization( nsites, sites, bc_in ) do while(associated(currentPatch)) - - call set_root_fraction(currentPatch%rootfr_ft(ft,1:bc_in(s)%nlevsoil), ft, & - bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) + bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1), & + icontext=i_hydro_rootprof_context) !zero cohort-summed variables. currentPatch%total_canopy_area = 0.0_r8 diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 0429f7e77b..41f979e85d 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -54,7 +54,6 @@ module EDPatchDynamicsMod public :: disturbance_rates public :: check_patch_area public :: set_patchno - public :: set_root_fraction private:: fuse_2_patches character(len=*), parameter, private :: sourcefile = & @@ -1889,161 +1888,4 @@ function countPatches( nsites, sites ) result ( totNumPatches ) end function countPatches - ! ==================================================================================== - - subroutine set_root_fraction(root_fraction, ft, zi, lowerb ) - ! - ! !DESCRIPTION: - ! Calculates the fractions of the root biomass in each layer for each pft. - ! It assumes an exponential decay. If the soil depth is shallower than - ! then exponential attenuation function, then it will normalize - ! the profile and divide through. - ! - ! !USES: - - ! - ! !ARGUMENTS - real(r8),intent(inout) :: root_fraction(:) - integer, intent(in) :: ft - real(r8),intent(in) :: zi(lowerb:) - integer,intent(in) :: lowerb - - ! Parameters - integer, parameter :: exponential_1p_profile_type = 1 - integer, parameter :: jackson_beta_profile_type = 2 - integer, parameter :: exponential_2p_profile_type = 3 - - integer, parameter :: root_profile_type = exponential_2p_profile_type - - !---------------------------------------------------------------------- - - if(lbound(zi,1).ne.0) then - write(fates_log(),*) 'lbound:',lbound(zi) - write(fates_log(),*) 'ubound:',ubound(zi) - write(fates_log(),*) 'layer interface levels should have 0 index' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - select case(root_profile_type) - case ( exponential_1p_profile_type ) - call exponential_1p_root_profile(root_fraction, ft, zi) - case ( jackson_beta_profile_type ) - call jackson_beta_root_profile(root_fraction, ft, zi) - case ( exponential_2p_profile_type ) - call exponential_2p_root_profile(root_fraction, ft, zi) - case default - write(fates_log(),*) 'An undefined root profile type was specified' - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end select - - return - end subroutine set_root_fraction - - ! ===================================================================================== - - subroutine exponential_2p_root_profile(root_fraction, ft, zi ) - ! - ! !ARGUMENTS - real(r8),intent(out) :: root_fraction(:) - integer,intent(in) :: ft - real(r8),intent(in) :: zi(0:) - - ! Locals - integer :: nlevsoil ! Number of soil layers - integer :: lev ! soil layer index - real(r8) :: sum_rootfr ! sum of root fraction for normalization - - nlevsoil = ubound(zi,1) - - sum_rootfr = 0.0_r8 - do lev = 1, nlevsoil - root_fraction(lev) = .5_r8*( & - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev-1)) & - + exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev-1)) & - - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev)) & - - exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev))) - - sum_rootfr = sum_rootfr + root_fraction(lev) - end do - - ! Normalize the root profile - root_fraction(1:nlevsoil) = root_fraction(1:nlevsoil)/sum_rootfr - - return - end subroutine exponential_2p_root_profile - - ! ===================================================================================== - - subroutine exponential_1p_root_profile(root_fraction, ft, zi) - - ! - ! !ARGUMENTS - real(r8),intent(out) :: root_fraction(:) - integer,intent(in) :: ft - real(r8),intent(in) :: zi(0:) - - ! - ! LOCAL VARIABLES: - integer :: lev ! soil depth layer index - integer :: nlevsoil ! number of soil layers - real(r8) :: depth ! Depth to middle of layer [m] - real(r8) :: sum_rootfr ! sum of rooting profile for normalization - - real(r8), parameter :: rootprof_exp = 3. ! how steep profile is - ! for root C inputs (1/ e-folding depth) (1/m) - - nlevsoil = ubound(zi,1) - - ! define rooting profile from exponential parameters - sum_rootfr = 0.0_r8 - do lev = 1, nlevsoil - root_fraction(lev) = exp(-rootprof_exp * 0.5*(zi(lev)+zi(lev-1)) ) - sum_rootfr = sum_rootfr + root_fraction(lev) - end do - - ! Normalize the root profile - root_fraction(1:nlevsoil) = root_fraction(1:nlevsoil)/sum_rootfr - - - return - end subroutine exponential_1p_root_profile - - ! ===================================================================================== - - subroutine jackson_beta_root_profile(root_fraction, ft, zi) - - - ! !ARGUMENTS - real(r8),intent(out) :: root_fraction(:) ! fraction of root mass in each soil layer - integer,intent(in) :: ft ! functional type - real(r8),intent(in) :: zi(0:) ! depth of layer interfaces 0-nlevsoil - - ! - ! LOCAL VARIABLES: - integer :: lev ! soil depth layer index - integer :: nlevsoil ! number of soil layers - real(r8) :: sum_rootfr ! sum of rooting profile, for normalization - - integer, parameter :: rooting_profile_varindex_water = 1 - - nlevsoil = ubound(zi,1) - ! use beta distribution parameter from Jackson et al., 1996 - sum_rootfr = 0.0_r8 - do lev = 1, nlevsoil - root_fraction(lev) = & - ( EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & - ( zi(lev-1)*100._r8) - & - EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & - ( zi(lev)*100._r8) ) - sum_rootfr = sum_rootfr + root_fraction(lev) - end do - - ! Normalize the root profile - root_fraction(1:nlevsoil) = root_fraction(1:nlevsoil)/sum_rootfr - - return - end subroutine jackson_beta_root_profile - - end module EDPatchDynamicsMod diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 02df97f47d..475faf27ec 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -2013,8 +2013,9 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) use FatesConstantsMod, only : itrue use FatesGlobals, only : endrun => fates_endrun use EDParamsMod , only : ED_val_cwd_flig, ED_val_cwd_fcel - use EDPatchDynamicsMod, only : set_root_fraction - + use FatesAllometryMod, only : set_root_fraction + use FatesAllometryMod, only : i_biomass_rootprof_context + implicit none @@ -2119,7 +2120,8 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! This generates a rooting profile over the whole soil column call set_root_fraction(cinput_rootfr(ft,1:bc_in(s)%nlevsoil), ft, & - bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) + bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1), & + icontext=i_biomass_rootprof_context) end do diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index c4b6b9f0c6..6094ff5f85 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -113,10 +113,17 @@ module FatesAllometryMod public :: bstore_allom ! Generic maximum storage carbon wrapper public :: StructureResetOfDH ! Method to set DBH to sync with structure biomass public :: CheckIntegratedAllometries + public :: set_root_fraction ! Generic wrapper to calculate normalized + ! root profiles logical , parameter :: verbose_logging = .false. character(len=*), parameter :: sourcefile = __FILE__ + + integer, parameter, public :: i_hydro_rootprof__context = 1 + integer, parameter, public :: i_biomass_rootprof_context = 2 + + ! If testing b4b with older versions, do not remove sapwood ! Our old methods with saldarriaga did not remove sapwood from the ! bdead pool. But newer allometries are providing total agb @@ -1775,6 +1782,188 @@ end subroutine carea_2pwr ! ========================================================================= + subroutine set_root_fraction(root_fraction, ft, zi, lowerb, icontext ) + ! + ! !DESCRIPTION: + ! Calculates the fractions of the root biomass in each layer for each pft. + ! It assumes an exponential decay. If the soil depth is shallower than + ! then exponential attenuation function, then it will normalize + ! the profile and divide through. + ! + ! !USES: + + ! + ! !ARGUMENTS + real(r8),intent(inout) :: root_fraction(:) + integer, intent(in) :: ft + real(r8),intent(in) :: zi(lowerb:) + integer,intent(in) :: lowerb + integer,intent(in) :: icontext + + ! Parameters + ! + ! TO-DO: NEXT TIME WE ROLL OUT A NEW PARAMETER INTERFACE, ADD + ! PROFILE SWAPPING FLAGS. OR IF THERE IS NO DEMAND< LEAVE AS IS. + ! + ! + ! Two context exist 'hydraulic' and 'biomass' + ! These two contexts are allowed to have different profiles + + integer, parameter :: exponential_1p_profile_type = 1 + integer, parameter :: jackson_beta_profile_type = 2 + integer, parameter :: exponential_2p_profile_type = 3 + + integer :: root_profile_type + + !---------------------------------------------------------------------- + + if(lbound(zi,1).ne.0) then + write(fates_log(),*) 'lbound:',lbound(zi) + write(fates_log(),*) 'ubound:',ubound(zi) + write(fates_log(),*) 'layer interface levels should have 0 index' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if(icontext == i_hydro_rootprof__context) then + + root_profile_type = exponential_2p_profile_type + + else if(icontext == i_biomass_rootprof_context) + + root_profile_type = jackson_beta_profile_type + + else + write(fates_log(),*) 'An undefined context for calculating root profiles was provided' + write(fates_log(),*) 'There are only two contexts, hydraulic and biomass, pick one.' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + + select case(root_profile_type) + case ( exponential_1p_profile_type ) + call exponential_1p_root_profile(root_fraction, ft, zi) + case ( jackson_beta_profile_type ) + call jackson_beta_root_profile(root_fraction, ft, zi) + case ( exponential_2p_profile_type ) + call exponential_2p_root_profile(root_fraction, ft, zi) + case default + write(fates_log(),*) 'An undefined root profile type was specified' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + + return + end subroutine set_root_fraction + + ! ===================================================================================== + + subroutine exponential_2p_root_profile(root_fraction, ft, zi ) + ! + ! !ARGUMENTS + real(r8),intent(out) :: root_fraction(:) + integer,intent(in) :: ft + real(r8),intent(in) :: zi(0:) + + ! Locals + integer :: nlevsoil ! Number of soil layers + integer :: lev ! soil layer index + real(r8) :: sum_rootfr ! sum of root fraction for normalization + + nlevsoil = ubound(zi,1) + + sum_rootfr = 0.0_r8 + do lev = 1, nlevsoil + root_fraction(lev) = .5_r8*( & + exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev-1)) & + + exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev-1)) & + - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev)) & + - exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev))) + + sum_rootfr = sum_rootfr + root_fraction(lev) + end do + + ! Normalize the root profile + root_fraction(1:nlevsoil) = root_fraction(1:nlevsoil)/sum_rootfr + + return + end subroutine exponential_2p_root_profile + + ! ===================================================================================== + + subroutine exponential_1p_root_profile(root_fraction, ft, zi) + + ! + ! !ARGUMENTS + real(r8),intent(out) :: root_fraction(:) + integer,intent(in) :: ft + real(r8),intent(in) :: zi(0:) + + ! + ! LOCAL VARIABLES: + integer :: lev ! soil depth layer index + integer :: nlevsoil ! number of soil layers + real(r8) :: depth ! Depth to middle of layer [m] + real(r8) :: sum_rootfr ! sum of rooting profile for normalization + + real(r8), parameter :: rootprof_exp = 3. ! how steep profile is + ! for root C inputs (1/ e-folding depth) (1/m) + + nlevsoil = ubound(zi,1) + + ! define rooting profile from exponential parameters + sum_rootfr = 0.0_r8 + do lev = 1, nlevsoil + root_fraction(lev) = exp(-rootprof_exp * 0.5*(zi(lev)+zi(lev-1)) ) + sum_rootfr = sum_rootfr + root_fraction(lev) + end do + + ! Normalize the root profile + root_fraction(1:nlevsoil) = root_fraction(1:nlevsoil)/sum_rootfr + + + return + end subroutine exponential_1p_root_profile + + ! ===================================================================================== + + subroutine jackson_beta_root_profile(root_fraction, ft, zi) + + + ! !ARGUMENTS + real(r8),intent(out) :: root_fraction(:) ! fraction of root mass in each soil layer + integer,intent(in) :: ft ! functional type + real(r8),intent(in) :: zi(0:) ! depth of layer interfaces 0-nlevsoil + + ! + ! LOCAL VARIABLES: + integer :: lev ! soil depth layer index + integer :: nlevsoil ! number of soil layers + real(r8) :: sum_rootfr ! sum of rooting profile, for normalization + + integer, parameter :: rooting_profile_varindex_water = 1 + + nlevsoil = ubound(zi,1) + ! use beta distribution parameter from Jackson et al., 1996 + sum_rootfr = 0.0_r8 + do lev = 1, nlevsoil + root_fraction(lev) = & + ( EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & + ( zi(lev-1)*100._r8) - & + EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & + ( zi(lev)*100._r8) ) + sum_rootfr = sum_rootfr + root_fraction(lev) + end do + + ! Normalize the root profile + root_fraction(1:nlevsoil) = root_fraction(1:nlevsoil)/sum_rootfr + + return + end subroutine jackson_beta_root_profile + + + ! ===================================================================================== + subroutine StructureResetOfDH( bdead, ipft, canopy_trim, d, h ) ! ========================================================================= diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 9bcd36f846..a1bb88445c 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -79,10 +79,11 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) use FatesConstantsMod, only : rgas => rgas_J_K_kmol use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use FatesParameterDerivedMod, only : param_derived - use EDPatchDynamicsMod, only: set_root_fraction use EDParamsMod, only : ED_val_bbopt_c3, ED_val_bbopt_c4, ED_val_base_mr_20 - use FatesAllometryMod, only : bleaf, storage_fraction_of_target - + use FatesAllometryMod, only : bleaf + use FatesAllometryMod, only : storage_fraction_of_target + use FatesAllometryMod, only : set_root_fraction + use FatesAllometryMod, only : i_hydro_rootprof_context ! ARGUMENTS: ! ----------------------------------------------------------------------------------- @@ -300,10 +301,14 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) kn(ft) = exp(0.00963_r8 * EDPftvarcon_inst%vcmax25top(ft) - 2.43_r8) end if - ! This is probably unnecessary and already calculated (RGK 05-2018) + ! This is probably unnecessary and already calculated + ! ALSO, THIS ROOTING PROFILE IS USED TO CALCULATE RESPIRATION + ! YET IT USES THE PROFILE THAT IS CONSISTENT WITH WATER UPTAKE + ! AND NOT THE PROFILE WE USE FOR DECOMPOSITION + ! SEEMS LIKE THE LATTER WOULD BE MORE APPROPRIATE, RIGHT? (RGK 05-2018) call set_root_fraction(currentPatch%rootfr_ft(ft,1:bc_in(s)%nlevsoil), ft, & - bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1)) - + bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1), & + icontext = i_hydro_rootprof_context) end do !ft From 62d7c2b590015bfce687712e46c24f67ef666da9 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 11 May 2018 15:56:24 -0700 Subject: [PATCH 06/12] Clean-up and tweaks bringing rooting depth to fates allometry. --- biogeochem/EDCanopyStructureMod.F90 | 15 ++++-- biogeochem/EDPhysiologyMod.F90 | 56 +++++++++------------- biogeochem/FatesAllometryMod.F90 | 13 +++-- biogeophys/FatesPlantRespPhotosynthMod.F90 | 2 +- 4 files changed, 43 insertions(+), 43 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index d86a2b387a..4a6a20e5a4 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -886,7 +886,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) use FatesAllometryMod , only : i_hydro_rootprof_context use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index use EDtypesMod , only : area - use EDPftvarcon , only : EDPftvarcon_inst + use EDPftvarcon , only : EDPftvarcon_inst ! !ARGUMENTS integer , intent(in) :: nsites @@ -921,10 +921,14 @@ subroutine canopy_summarization( nsites, sites, bc_in ) do while(associated(currentPatch)) - call set_root_fraction(currentPatch%rootfr_ft(ft,1:bc_in(s)%nlevsoil), ft, & - bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1), & - icontext=i_hydro_rootprof_context) - + ! Calculate rooting depth fractions for the patch x pft + do ft = 1, numpft + call set_root_fraction(currentPatch%rootfr_ft(ft,1:bc_in(s)%nlevsoil), ft, & + bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1), & + icontext=i_hydro_rootprof_context) + end do + + !zero cohort-summed variables. currentPatch%total_canopy_area = 0.0_r8 currentPatch%total_tree_area = 0.0_r8 @@ -937,6 +941,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) ft = currentCohort%pft + ! Update the cohort's index within the size bin classes ! Update the cohort's index within the SCPF classification system call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 475faf27ec..7e2b2f10ac 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -2041,33 +2041,16 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) real(r8) :: croot_prof_perpatch(1:hlm_numlevgrnd) real(r8) :: surface_prof(1:hlm_numlevgrnd) integer :: ft + integer :: nlev_eff_decomp real(r8) :: rootfr_tot(1:maxpft) real(r8) :: biomass_bg_ft(1:maxpft) real(r8) :: surface_prof_tot, leaf_prof_sum, stem_prof_sum, froot_prof_sum, biomass_bg_tot real(r8) :: delta - ! NOTE(bja, 201608) these were removed from clm in clm4_5_10_r187 -! logical, parameter :: exponential_rooting_profile = .true. -! logical, parameter :: pftspecific_rootingprofile = .true. - - ! NOTE(bja, 201608) as of clm4_5_10_r187 rootprof_exp is now a - ! private function level parameter in RootBiophysMod.F90::exponential_rootfr() -! real(r8), parameter :: rootprof_exp = 3. ! how steep profile is - ! for root C inputs (1/ e-folding depth) (1/m) - ! NOTE(rgk, 201705) this parameter was brought over from SoilBiogeochemVerticalProfile ! how steep profile is for surface components (1/ e_folding depth) (1/m) real(r8), parameter :: surfprof_exp = 10. - ! NOTE(bja, 201608) as of clm4_5_10_r187 rootprof_beta is now a - ! two dimensional array with the second dimension being water,1, - ! or carbon,2,. These are currently hard coded, but may be - ! overwritten by the namelist. - - ! Note cdk 2016/08 we actually want to use the carbon index here rather than the water index. - ! Doing so will be answer changing though so perhaps easiest to do this in steps. - integer, parameter :: rooting_profile_varindex_water = 1 - real(r8) :: leaf_prof(1:nsites, 1:hlm_numlevgrnd) real(r8) :: froot_prof(1:nsites, 1:maxpft, 1:hlm_numlevgrnd) real(r8) :: croot_prof(1:nsites, 1:hlm_numlevgrnd) @@ -2100,6 +2083,11 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) do s = 1,nsites + ! Calculate the number of effective decomposition layers + ! This takes into account if vertical soil biogeochem is on, how deep the soil column + ! is, and also which layers may be frozen + nlev_eff_decomp = min(max(bc_in(s)%max_rooting_depth_index_col, 1), bc_in(s)%nlevdecomp) + ! define a single shallow surface profile for surface additions (leaves, stems, and N deposition) surface_prof(:) = 0._r8 do j = 1, bc_in(s)%nlevdecomp @@ -2118,10 +2106,14 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) cinput_rootfr(:,:) = 0._r8 do ft = 1, numpft - ! This generates a rooting profile over the whole soil column + ! This generates a rooting profile over the whole soil column for each pft call set_root_fraction(cinput_rootfr(ft,1:bc_in(s)%nlevsoil), ft, & - bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1), & - icontext=i_biomass_rootprof_context) + bc_in(s)%zi_sisl, lowerb=lbound(bc_in(s)%zi_sisl,1), & + icontext=i_biomass_rootprof_context) + + do j=1,nlev_eff_decomp + cinput_rootfr(ft,j) = cinput_rootfr(ft,j)/bc_in(s)%dz_decomp_sisl(j) + end do end do @@ -2132,13 +2124,13 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) surface_prof_tot = 0._r8 ! - do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), bc_in(s)%nlevdecomp) + do j = 1, nlev_eff_decomp surface_prof_tot = surface_prof_tot + surface_prof(j) * bc_in(s)%dz_decomp_sisl(j) end do do ft = 1,numpft - do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), bc_in(s)%nlevdecomp) - rootfr_tot(ft) = rootfr_tot(ft) + cinput_rootfr(ft,j) + do j = 1, nlev_eff_decomp + rootfr_tot(ft) = rootfr_tot(ft) + cinput_rootfr(ft,j)*bc_in(s)%dz_decomp_sisl(j) end do end do ! @@ -2147,7 +2139,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) if ( (bc_in(s)%max_rooting_depth_index_col > 0) .and. (rootfr_tot(ft) > 0._r8) ) then ! where there is not permafrost extending to the surface, integrate the profiles over the active layer ! this is equivalent to integrating over all soil layers outside of permafrost regions - do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), bc_in(s)%nlevdecomp) + do j = 1, nlev_eff_decomp froot_prof(s,ft,j) = cinput_rootfr(ft,j) / rootfr_tot(ft) end do else @@ -2161,7 +2153,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) if ( (bc_in(s)%max_rooting_depth_index_col > 0) .and. (surface_prof_tot > 0._r8) ) then ! where there is not permafrost extending to the surface, integrate the profiles over the active layer ! this is equivalent to integrating over all soil layers outside of permafrost regions - do j = 1, min(max(bc_in(s)%max_rooting_depth_index_col, 1), bc_in(s)%nlevdecomp) + do j = 1, nlev_eff_decomp ! set all surface processes to shallower profile leaf_prof(s,j) = surface_prof(j)/ surface_prof_tot stem_prof(s,j) = surface_prof(j)/ surface_prof_tot @@ -2228,16 +2220,14 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) end do end do - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! now disaggregate the inputs vertically, using the vertical profiles - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! now disaggregate the inputs vertically, using the vertical profiles + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - do s = 1,nsites + do s = 1,nsites - ! do g = bounds%begg,bounds%endg - ! if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + currentPatch => sites(s)%oldest_patch - do while(associated(currentPatch)) ! the CWD pools lose information about which PFT they came from; for the stems this doesn't matter as they all have the same profile, diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 6094ff5f85..ad6f91dbf7 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -120,7 +120,10 @@ module FatesAllometryMod character(len=*), parameter :: sourcefile = __FILE__ - integer, parameter, public :: i_hydro_rootprof__context = 1 + ! These are public contexts that other routines + ! should pass as arguments to the generic root profile + ! wrapper. + integer, parameter, public :: i_hydro_rootprof_context = 1 integer, parameter, public :: i_biomass_rootprof_context = 2 @@ -1824,11 +1827,11 @@ subroutine set_root_fraction(root_fraction, ft, zi, lowerb, icontext ) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if(icontext == i_hydro_rootprof__context) then + if(icontext == i_hydro_rootprof_context) then root_profile_type = exponential_2p_profile_type - else if(icontext == i_biomass_rootprof_context) + else if(icontext == i_biomass_rootprof_context) then root_profile_type = jackson_beta_profile_type @@ -1940,7 +1943,9 @@ subroutine jackson_beta_root_profile(root_fraction, ft, zi) integer :: lev ! soil depth layer index integer :: nlevsoil ! number of soil layers real(r8) :: sum_rootfr ! sum of rooting profile, for normalization - + + ! Note cdk 2016/08 we actually want to use the carbon index here rather than the water index. + ! Doing so will be answer changing though so perhaps easiest to do this in steps. integer, parameter :: rooting_profile_varindex_water = 1 nlevsoil = ubound(zi,1) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index a1bb88445c..c4bd555814 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -309,7 +309,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) call set_root_fraction(currentPatch%rootfr_ft(ft,1:bc_in(s)%nlevsoil), ft, & bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1), & icontext = i_hydro_rootprof_context) - + end do !ft From d5e433c139da513d4ebfde3cbed99196b12f3654 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 11 May 2018 17:22:14 -0700 Subject: [PATCH 07/12] Syntax updates to decomp profile stuff. --- biogeochem/EDPhysiologyMod.F90 | 135 +++++++++++++++++++++------------ 1 file changed, 88 insertions(+), 47 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 7e2b2f10ac..b687d4579d 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -1985,24 +1985,35 @@ subroutine cwd_out( currentSite, currentPatch, bc_in ) end subroutine cwd_out - + ! ===================================================================================== subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) - ! Created by Charlie Koven and Rosie Fisher, 2014-2015 - ! take the flux out of the fragmenting litter pools and port into the decomposing litter pools. - ! in this implementation, decomposing pools are assumed to be humus and non-flammable, whereas fragmenting pools - ! are assumed to be physically fragmenting but not respiring. This is a simplification, but allows us to - ! a) reconcile the need to track both chemical fractions (lignin, cellulose, labile) and size fractions (trunk, branch, etc.) - ! b) to impose a realistic delay on the surge of nutrients into the litter pools when large CWD is added to the system via mortality - - ! because of the different subgrid structure, this subroutine includes the functionality that in the big-leaf BGC model, is calculated in SoilBiogeochemVerticalProfileMod - - ! The ED code is resolved at a daily timestep, but all of the CN-BGC fluxes are passed in as derivatives per second, - ! and then accumulated in the CNStateUpdate routines. One way of doing this is to pass back the CN fluxes per second, - ! and keep them constant for the whole day (making sure they are not overwritten. - ! This means that the carbon gets passed back and forth between the photosynthesis code (fast timestepping) to the ED code (slow timestepping), back to the BGC code (fast timestepping). - ! This means that the state update for the litter pools and for the CWD pools occurs at different timescales. + ! ----------------------------------------------------------------------------------- + ! Created by Charlie Koven and Rosie Fisher, 2014-2015 + ! take the flux out of the fragmenting litter pools and port into the decomposing + ! litter pools. + ! in this implementation, decomposing pools are assumed to be humus and non-flammable, + ! whereas fragmenting pools are assumed to be physically fragmenting but not + ! respiring. This is a simplification, but allows us to + ! + ! a) reconcile the need to track both chemical fractions (lignin, cellulose, labile) + ! and size fractions (trunk, branch, etc.) + ! b) to impose a realistic delay on the surge of nutrients into the litter pools + ! when large CWD is added to the system via mortality + ! + ! Because of the different subgrid structure, this subroutine includes the functionality + ! that in the big-leaf BGC model, is calculated in SoilBiogeochemVerticalProfileMod + ! + ! The ED code is resolved at a daily timestep, but all of the CN-BGC fluxes are passed + ! in as derivatives per second, and then accumulated in the CNStateUpdate routines. + ! One way of doing this is to pass back the CN fluxes per second, and keep them + ! constant for the whole day (making sure they are not overwritten. This means that + ! the carbon gets passed back and forth between the photosynthesis code + ! (fast timestepping) to the ED code (slow timestepping), back to the BGC code + ! (fast timestepping). This means that the state update for the litter pools and + ! for the CWD pools occurs at different timescales. + ! ----------------------------------------------------------------------------------- use EDTypesMod, only : AREA use EDPftvarcon, only : EDPftvarcon_inst @@ -2034,10 +2045,12 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) real(r8) mass_convert ! ED uses kg, CLM uses g integer :: begp,endp integer :: begc,endc !bounds + !------------------------------------------------------------------------ ! The following scratch arrays are allocated for maximum possible ! pft and layer usage - real(r8) :: cinput_rootfr(1:maxpft, 1:hlm_numlevgrnd) ! column by pft root fraction used for calculating inputs + + real(r8) :: cinput_rootfr(1:maxpft, 1:hlm_numlevgrnd) real(r8) :: croot_prof_perpatch(1:hlm_numlevgrnd) real(r8) :: surface_prof(1:hlm_numlevgrnd) integer :: ft @@ -2137,8 +2150,9 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! rescale the fine root profile do ft = 1,numpft if ( (bc_in(s)%max_rooting_depth_index_col > 0) .and. (rootfr_tot(ft) > 0._r8) ) then - ! where there is not permafrost extending to the surface, integrate the profiles over the active layer - ! this is equivalent to integrating over all soil layers outside of permafrost regions + ! where there is not permafrost extending to the surface, integrate the profiles + ! over the active layer this is equivalent to integrating over all soil layers + ! outside of permafrost regions do j = 1, nlev_eff_decomp froot_prof(s,ft,j) = cinput_rootfr(ft,j) / rootfr_tot(ft) end do @@ -2151,8 +2165,9 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! ! rescale the shallow profiles if ( (bc_in(s)%max_rooting_depth_index_col > 0) .and. (surface_prof_tot > 0._r8) ) then - ! where there is not permafrost extending to the surface, integrate the profiles over the active layer - ! this is equivalent to integrating over all soil layers outside of permafrost regions + ! where there is not permafrost extending to the surface, integrate the profiles over + ! the active layer this is equivalent to integrating over all soil layers outside of + ! permafrost regions do j = 1, nlev_eff_decomp ! set all surface processes to shallower profile leaf_prof(s,j) = surface_prof(j)/ surface_prof_tot @@ -2230,14 +2245,18 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) currentPatch => sites(s)%oldest_patch do while(associated(currentPatch)) - ! the CWD pools lose information about which PFT they came from; for the stems this doesn't matter as they all have the same profile, - ! however for the coarse roots they may have different profiles. to approximately recover this information, loop over all cohorts in patch - ! to calculate the total root biomass in that patch of each pft, and then rescale the croot_prof as the weighted average of the froot_prof - biomass_bg_ft(:) = 0._r8 + ! the CWD pools lose information about which PFT they came from; + ! for the stems this doesn't matter as they all have the same profile, + ! however for the coarse roots they may have different profiles. + ! to approximately recover this information, loop over all cohorts in patch + ! to calculate the total root biomass in that patch of each pft, and then + ! rescale the croot_prof as the weighted average of the froot_prof + biomass_bg_ft(1:numpft) = 0._r8 currentCohort => currentPatch%tallest do while(associated(currentCohort)) biomass_bg_ft(currentCohort%pft) = biomass_bg_ft(currentCohort%pft) + & - ((currentCohort%bdead + currentCohort%bsw ) * (1.0_r8-EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)) + & + ((currentCohort%bdead + currentCohort%bsw ) * & + (1.0_r8-EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)) + & (currentCohort%br + currentCohort%bstore )) * & (currentCohort%n / currentPatch%area) currentCohort => currentCohort%shorter @@ -2256,7 +2275,8 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) if ( biomass_bg_tot .gt. 0._r8) then do ft = 1,numpft do j = 1, bc_in(s)%nlevdecomp - croot_prof_perpatch(j) = croot_prof_perpatch(j) + froot_prof(s,ft,j) * biomass_bg_ft(ft) / biomass_bg_tot + croot_prof_perpatch(j) = croot_prof_perpatch(j) + & + froot_prof(s,ft,j) * biomass_bg_ft(ft) / biomass_bg_tot end do end do else ! no biomass @@ -2269,29 +2289,38 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) croot_prof(s, j) = croot_prof(s, j) + croot_prof_perpatch(j) * currentPatch%area / AREA end do ! - ! now disaggregate, vertically and by decomposition substrate type, the actual fluxes from CWD and litter pools + ! now disaggregate, vertically and by decomposition substrate type, the + ! actual fluxes from CWD and litter pools ! ! do c = 1, ncwd - ! write(fates_log(),*)'cdk CWD_AG_out', c, currentpatch%CWD_AG_out(c), ED_val_cwd_fcel, currentpatch%area/AREA - ! write(fates_log(),*)'cdk CWD_BG_out', c, currentpatch%CWD_BG_out(c), ED_val_cwd_fcel, currentpatch%area/AREA + ! write(fates_log(),*)'cdk CWD_AG_out', c, currentpatch%CWD_AG_out(c), + ! ED_val_cwd_fcel, currentpatch%area/AREA + ! write(fates_log(),*)'cdk CWD_BG_out', c, currentpatch%CWD_BG_out(c), + ! ED_val_cwd_fcel, currentpatch%area/AREA ! end do ! do ft = 1,numpft - ! write(fates_log(),*)'cdk leaf_litter_out', ft, currentpatch%leaf_litter_out(ft), ED_val_cwd_fcel, currentpatch%area/AREA - ! write(fates_log(),*)'cdk root_litter_out', ft, currentpatch%root_litter_out(ft), ED_val_cwd_fcel, currentpatch%area/AREA + ! write(fates_log(),*)'cdk leaf_litter_out', ft, currentpatch%leaf_litter_out(ft), + ! ED_val_cwd_fcel, currentpatch%area/AREA + ! write(fates_log(),*)'cdk root_litter_out', ft, currentpatch%root_litter_out(ft), + ! ED_val_cwd_fcel, currentpatch%area/AREA ! end do ! ! ! CWD pools fragmenting into decomposing litter pools. do ci = 1, ncwd do j = 1, bc_in(s)%nlevdecomp bc_out(s)%FATES_c_to_litr_cel_c_col(j) = bc_out(s)%FATES_c_to_litr_cel_c_col(j) + & - currentpatch%CWD_AG_out(ci) * ED_val_cwd_fcel * currentpatch%area/AREA * stem_prof(s,j) + currentpatch%CWD_AG_out(ci) * ED_val_cwd_fcel * & + currentpatch%area/AREA * stem_prof(s,j) bc_out(s)%FATES_c_to_litr_lig_c_col(j) = bc_out(s)%FATES_c_to_litr_lig_c_col(j) + & - currentpatch%CWD_AG_out(ci) * ED_val_cwd_flig * currentpatch%area/AREA * stem_prof(s,j) + currentpatch%CWD_AG_out(ci) * ED_val_cwd_flig * & + currentpatch%area/AREA * stem_prof(s,j) ! bc_out(s)%FATES_c_to_litr_cel_c_col(j) = bc_out(s)%FATES_c_to_litr_cel_c_col(j) + & - currentpatch%CWD_BG_out(ci) * ED_val_cwd_fcel * currentpatch%area/AREA * croot_prof_perpatch(j) + currentpatch%CWD_BG_out(ci) * ED_val_cwd_fcel * & + currentpatch%area/AREA * croot_prof_perpatch(j) bc_out(s)%FATES_c_to_litr_lig_c_col(j) = bc_out(s)%FATES_c_to_litr_lig_c_col(j) + & - currentpatch%CWD_BG_out(ci) * ED_val_cwd_flig * currentpatch%area/AREA * croot_prof_perpatch(j) + currentpatch%CWD_BG_out(ci) * ED_val_cwd_flig * & + currentpatch%area/AREA * croot_prof_perpatch(j) end do end do @@ -2299,26 +2328,35 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) do ft = 1,numpft do j = 1, bc_in(s)%nlevdecomp bc_out(s)%FATES_c_to_litr_lab_c_col(j) = bc_out(s)%FATES_c_to_litr_lab_c_col(j) + & - currentpatch%leaf_litter_out(ft) * EDPftvarcon_inst%lf_flab(ft) * currentpatch%area/AREA * leaf_prof(s,j) + currentpatch%leaf_litter_out(ft) * EDPftvarcon_inst%lf_flab(ft) * & + currentpatch%area/AREA * leaf_prof(s,j) bc_out(s)%FATES_c_to_litr_cel_c_col(j) = bc_out(s)%FATES_c_to_litr_cel_c_col(j) + & - currentpatch%leaf_litter_out(ft) * EDPftvarcon_inst%lf_fcel(ft) * currentpatch%area/AREA * leaf_prof(s,j) + currentpatch%leaf_litter_out(ft) * EDPftvarcon_inst%lf_fcel(ft) * & + currentpatch%area/AREA * leaf_prof(s,j) bc_out(s)%FATES_c_to_litr_lig_c_col(j) = bc_out(s)%FATES_c_to_litr_lig_c_col(j) + & - currentpatch%leaf_litter_out(ft) * EDPftvarcon_inst%lf_flig(ft) * currentpatch%area/AREA * leaf_prof(s,j) + currentpatch%leaf_litter_out(ft) * EDPftvarcon_inst%lf_flig(ft) * & + currentpatch%area/AREA * leaf_prof(s,j) ! bc_out(s)%FATES_c_to_litr_lab_c_col(j) = bc_out(s)%FATES_c_to_litr_lab_c_col(j) + & - currentpatch%root_litter_out(ft) * EDPftvarcon_inst%fr_flab(ft) * currentpatch%area/AREA * froot_prof(s,ft,j) + currentpatch%root_litter_out(ft) * EDPftvarcon_inst%fr_flab(ft) * & + currentpatch%area/AREA * froot_prof(s,ft,j) bc_out(s)%FATES_c_to_litr_cel_c_col(j) = bc_out(s)%FATES_c_to_litr_cel_c_col(j) + & - currentpatch%root_litter_out(ft) * EDPftvarcon_inst%fr_fcel(ft) * currentpatch%area/AREA * froot_prof(s,ft,j) + currentpatch%root_litter_out(ft) * EDPftvarcon_inst%fr_fcel(ft) * & + currentpatch%area/AREA * froot_prof(s,ft,j) bc_out(s)%FATES_c_to_litr_lig_c_col(j) = bc_out(s)%FATES_c_to_litr_lig_c_col(j) + & - currentpatch%root_litter_out(ft) * EDPftvarcon_inst%fr_flig(ft) * currentpatch%area/AREA * froot_prof(s,ft,j) + currentpatch%root_litter_out(ft) * EDPftvarcon_inst%fr_flig(ft) * & + currentpatch%area/AREA * froot_prof(s,ft,j) ! !! and seed_decay too. for now, use the same lability fractions as for leaf litter bc_out(s)%FATES_c_to_litr_lab_c_col(j) = bc_out(s)%FATES_c_to_litr_lab_c_col(j) + & - currentpatch%seed_decay(ft) * EDPftvarcon_inst%lf_flab(ft) * currentpatch%area/AREA * leaf_prof(s,j) + currentpatch%seed_decay(ft) * EDPftvarcon_inst%lf_flab(ft) * & + currentpatch%area/AREA * leaf_prof(s,j) bc_out(s)%FATES_c_to_litr_cel_c_col(j) = bc_out(s)%FATES_c_to_litr_cel_c_col(j) + & - currentpatch%seed_decay(ft) * EDPftvarcon_inst%lf_fcel(ft) * currentpatch%area/AREA * leaf_prof(s,j) + currentpatch%seed_decay(ft) * EDPftvarcon_inst%lf_fcel(ft) * & + currentpatch%area/AREA * leaf_prof(s,j) bc_out(s)%FATES_c_to_litr_lig_c_col(j) = bc_out(s)%FATES_c_to_litr_lig_c_col(j) + & - currentpatch%seed_decay(ft) * EDPftvarcon_inst%lf_flig(ft) * currentpatch%area/AREA * leaf_prof(s,j) + currentpatch%seed_decay(ft) * EDPftvarcon_inst%lf_flig(ft) * & + currentpatch%area/AREA * leaf_prof(s,j) ! enddo end do @@ -2331,9 +2369,12 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) do s = 1, nsites do j = 1, bc_in(s)%nlevdecomp ! time unit conversion - bc_out(s)%FATES_c_to_litr_lab_c_col(j)=bc_out(s)%FATES_c_to_litr_lab_c_col(j) * mass_convert / time_convert - bc_out(s)%FATES_c_to_litr_cel_c_col(j)=bc_out(s)%FATES_c_to_litr_cel_c_col(j) * mass_convert / time_convert - bc_out(s)%FATES_c_to_litr_lig_c_col(j)=bc_out(s)%FATES_c_to_litr_lig_c_col(j) * mass_convert / time_convert + bc_out(s)%FATES_c_to_litr_lab_c_col(j)=bc_out(s)%FATES_c_to_litr_lab_c_col(j) * & + mass_convert / time_convert + bc_out(s)%FATES_c_to_litr_cel_c_col(j)=bc_out(s)%FATES_c_to_litr_cel_c_col(j) * & + mass_convert / time_convert + bc_out(s)%FATES_c_to_litr_lig_c_col(j)=bc_out(s)%FATES_c_to_litr_lig_c_col(j) * & + mass_convert / time_convert end do end do From 127e2fec4e6aabba98f78a3789b13228f6486acf Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 22 May 2018 16:05:37 -0700 Subject: [PATCH 08/12] Bringing variable soil depths to fates-hydro. --- biogeochem/EDCanopyStructureMod.F90 | 4 +- biogeochem/EDCohortDynamicsMod.F90 | 6 +- biogeochem/EDPatchDynamicsMod.F90 | 2 +- biogeophys/FatesPlantHydraulicsMod.F90 | 205 +++++++++++++------------ main/EDMainMod.F90 | 4 +- main/FatesHistoryInterfaceMod.F90 | 111 +++++++++++-- main/FatesHydraulicsMemMod.F90 | 128 ++++++++++----- 7 files changed, 301 insertions(+), 159 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index d8aef91aec..00b80830d5 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -393,7 +393,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) allocate(copyc) if( hlm_use_planthydro.eq.itrue ) then - call InitHydrCohort(copyc) + call InitHydrCohort(currentSite,copyc) endif call copy_cohort(currentCohort, copyc) ! @@ -721,7 +721,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) if(cc_gain < currentCohort%c_area)then allocate(copyc) if( hlm_use_planthydro.eq.itrue ) then - call InitHydrCohort(copyc) + call InitHydrCohort(CurrentSite,copyc) endif call copy_cohort(currentCohort, copyc) !makes an identical copy... diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 49f466c038..2c96085a7d 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -184,9 +184,9 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine new_cohort%isnew = .true. if( hlm_use_planthydro.eq.itrue ) then - call InitHydrCohort(new_cohort) - call updateSizeDepTreeHydProps(new_cohort, bc_in) - call initTreeHydStates(currentSite,new_cohort, bc_in) + call InitHydrCohort(CurrentSite,new_cohort) + call updateSizeDepTreeHydProps(CurrentSite,new_cohort, bc_in) + call initTreeHydStates(CurrentSite,new_cohort, bc_in) if(recruitstatus==1)then new_cohort%co_hydr%is_newly_recuited = .true. endif diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 404c85747a..5288718828 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -378,7 +378,7 @@ subroutine spawn_patches( currentSite, bc_in) do while(associated(currentCohort)) allocate(nc) - if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(nc) + if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(CurrentSite,nc) call zero_cohort(nc) ! nc is the new cohort that goes in the disturbed patch (new_patch)... currentCohort diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 47c5960317..a6a8f2a7fc 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -63,14 +63,14 @@ module FatesPlantHydraulicsMod use FatesHydraulicsMemMod, only: n_hypool_ag use FatesHydraulicsMemMod, only: n_hypool_troot use FatesHydraulicsMemMod, only: porous_media - use FatesHydraulicsMemMod, only: nlevsoi_hyd use FatesHydraulicsMemMod, only: cap_slp use FatesHydraulicsMemMod, only: cap_int use FatesHydraulicsMemMod, only: cap_corr use FatesHydraulicsMemMod, only: rwcft use FatesHydraulicsMemMod, only: C2B use FatesHydraulicsMemMod, only: InitHydraulicsDerived - + use FatesHydraulicsMemMod, only: nlevsoi_hyd_max + use clm_time_manager , only : get_step_size, get_nstep use FatesConstantsMod, only: cm2_per_m2 @@ -194,7 +194,7 @@ subroutine initTreeHydStates(site_p, cc_p, bc_in) !convert soil water contents to water potential in each soil layer and !assign it to the absorbing root (assume absorbing root water potential !in equlibrium w/ surrounding soil) - do j=1, nlevsoi_hyd + do j=1, site_p%si_hydr%nlevsoi_hyd !call swcVG_psi_from_th(waterstate_inst%h2osoi_liqvol_shell(c,j,1), & ! watsat(c,j), watres(c,j), alpha_VG(c,j), n_VG(c,j), m_VG(c,j), l_VG(c,j), & @@ -246,7 +246,8 @@ end subroutine initTreeHydStates ! ===================================================================================== - subroutine updateSizeDepTreeHydProps(cc_p,bc_in) + subroutine updateSizeDepTreeHydProps(currentSite,cc_p,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 @@ -257,6 +258,7 @@ subroutine updateSizeDepTreeHydProps(cc_p,bc_in) use shr_sys_mod , only : shr_sys_abort ! ! !ARGUMENTS: + type(ed_site_type) , intent(in) :: currentSite ! Site stuff type(ed_cohort_type) , intent(inout), target :: cc_p ! current cohort pointer type(bc_in_type) , intent(in) :: bc_in ! Boundary Conditions @@ -308,9 +310,13 @@ subroutine updateSizeDepTreeHydProps(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] + integer :: nlevsoi_hyd ! Number of soil hydraulic layers type(ed_cohort_hydr_type), pointer :: ccohort_hydr !----------------------------------------------------------------------- + + nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd + cCohort => cc_p ccohort_hydr => cc_p%co_hydr cPatch => cCohort%patchptr @@ -500,19 +506,19 @@ subroutine updateSizeDepTreeHydProps(cc_p,bc_in) end subroutine updateSizeDepTreeHydProps ! ===================================================================================== - subroutine updateSizeDepTreeHydStates(cc_p) + subroutine updateSizeDepTreeHydStates(currentSite,cc_p) ! ! !DESCRIPTION: ! ! !USES: ! !ARGUMENTS: + type(ed_site_type) , intent(in) :: currentSite ! Site stuff type(ed_cohort_type) , intent(inout), target :: cc_p ! current cohort pointer ! ! !LOCAL VARIABLES: type(ed_cohort_type), pointer :: cCohort type(ed_cohort_hydr_type), pointer :: ccohort_hydr integer :: j,k,FT ! indices - !----------------------------------------------------------------------- cCohort => cc_p @@ -532,7 +538,7 @@ subroutine updateSizeDepTreeHydStates(cc_p) ccohort_hydr%th_troot(k) = ccohort_hydr%th_troot(k) * & ccohort_hydr%v_troot_init(k) /ccohort_hydr%v_troot(k) enddo - do j=1,nlevsoi_hyd + do j=1,currentSite%si_hydr%nlevsoi_hyd ccohort_hydr%th_aroot(j) = ccohort_hydr%th_aroot(j) * & ccohort_hydr%v_aroot_layer_init(j)/ccohort_hydr%v_aroot_layer(j) enddo @@ -644,12 +650,12 @@ subroutine FuseCohortHydraulics(currentSite,currentCohort, nextCohort, bc_in, ne call psi_from_th(currentCohort%pft, 3, ccohort_hydr%th_troot(k-n_hypool_ag), & ccohort_hydr%psi_troot(k-n_hypool_ag), site_hydr, bc_in) end do - do j=1,nlevsoi_hyd + do j=1,site_hydr%nlevsoi_hyd call psi_from_th(currentCohort%pft, 4, ccohort_hydr%th_aroot(j), & ccohort_hydr%psi_aroot(j), site_hydr, bc_in) end do call flc_gs_from_psi(currentCohort, ccohort_hydr%psi_ag(1)) - call updateSizeDepTreeHydProps(currentCohort, bc_in) !hydraulics quantities that are functions of hite & biomasses + call updateSizeDepTreeHydProps(currentSite,currentCohort, bc_in) !hydraulics quantities that are functions of hite & biomasses ccohort_hydr%qtop_dt = (currentCohort%n*ccohort_hydr%qtop_dt + & nextCohort%n*ncohort_hydr%qtop_dt)/newn ccohort_hydr%dqtopdth_dthdt = (currentCohort%n*ccohort_hydr%dqtopdth_dthdt + & @@ -665,15 +671,17 @@ end subroutine FuseCohortHydraulics ! ===================================================================================== ! Initialization Routines ! ===================================================================================== - subroutine InitHydrCohort(currentCohort) + subroutine InitHydrCohort(currentSite,currentCohort) ! Arguments + type(ed_site_type), target :: currentSite type(ed_cohort_type), target :: currentCohort type(ed_cohort_hydr_type), pointer :: ccohort_hydr if ( hlm_use_planthydro.eq.ifalse ) return allocate(ccohort_hydr) currentCohort%co_hydr => ccohort_hydr + call ccohort_hydr%AllocateHydrCohortArrays(currentSite%si_hydr%nlevsoi_hyd) ccohort_hydr%is_newly_recuited = .false. end subroutine InitHydrCohort @@ -688,7 +696,7 @@ subroutine DeallocateHydrCohort(currentCohort) if ( hlm_use_planthydro.eq.ifalse ) return ccohort_hydr => currentCohort%co_hydr - + call ccohort_hydr%DeAllocateHydrCohortArrays() deallocate(ccohort_hydr) return @@ -700,6 +708,7 @@ subroutine InitHydrSites(sites,bc_in) ! Arguments type(ed_site_type),intent(inout),target :: sites(:) + type(bc_in_type),intent(in) :: bc_in(:) ! Locals integer :: nsites @@ -730,7 +739,6 @@ subroutine HydrSiteColdStart(sites, bc_in )! , bc_out) ! Arguments type(ed_site_type),intent(inout),target :: sites(:) type(bc_in_type),intent(in) :: bc_in(:) -! type(bc_out_type),intent(inout) :: bc_out(:) ! Local type(ed_site_hydr_type), pointer :: site_hydr @@ -739,24 +747,27 @@ subroutine HydrSiteColdStart(sites, bc_in )! , bc_out) integer :: s integer :: j integer :: nsites + integer :: nlevsoil ! Number of soil layers + integer :: nlevsoil_hyd ! Number of hydraulically relevant soil layers nsites = ubound(sites,1) do s = 1,nsites site_hydr => sites(s)%si_hydr - nlevsoil = bc_in(s)%nlevsoil + nlevsoil = bc_in(s)%nlevsoil + nlevsoil_hyd = site_hydr%nlevsoi_hyd - if (nlevsoi_hyd == 1) then + if ( nlevsoil_hyd == 1) then h2osoi_liqvol = min(bc_in(s)%eff_porosity_sl(nlevsoil), & bc_in(s)%h2o_liq_sisl(nlevsoil)/(bc_in(s)%dz_sisl(nlevsoil)*denh2o)) - site_hydr%h2osoi_liqvol_shell(:,1:nshell) = h2osoi_liqvol - site_hydr%h2osoi_liq_prev(:) = bc_in(s)%h2o_liq_sisl(nlevsoil) + site_hydr%h2osoi_liqvol_shell(nlevsoil_hyd,1:nshell) = h2osoi_liqvol + site_hydr%h2osoi_liq_prev(nlevsoil_hyd) = bc_in(s)%h2o_liq_sisl(nlevsoil) else - do j = 1, nlevsoi_hyd - h2osoi_liqvol = min(bc_in(s)%eff_porosity_gl(j), & + do j = 1,nlevsoil_hyd + h2osoi_liqvol = min(bc_in(s)%eff_porosity_sl(j), & bc_in(s)%h2o_liq_sisl(j)/(bc_in(s)%dz_sisl(j)*denh2o)) site_hydr%h2osoi_liqvol_shell(j,1:nshell) = h2osoi_liqvol @@ -764,7 +775,7 @@ subroutine HydrSiteColdStart(sites, bc_in )! , bc_out) end do end if - do j = 1, nlevsoi_hyd + do j = 1, nlevsoil_hyd ! Calculate the matric potential on the innner shell (this is used to initialize ! xylem and root pressures in new cohorts) call swcCampbell_psi_from_th(site_hydr%h2osoi_liqvol_shell(j,1), & @@ -774,7 +785,7 @@ subroutine HydrSiteColdStart(sites, bc_in )! , bc_out) site_hydr%psisoi_liq_innershell(j) = smp end do - site_hydr%l_aroot_layer(1:nlevsoi_hyd) = 0.0_r8 + site_hydr%l_aroot_layer(1:site_hydr%nlevsoi_hyd) = 0.0_r8 end do @@ -815,6 +826,7 @@ subroutine HydrSiteColdStart(sites, bc_in )! , bc_out) end subroutine HydrSiteColdStart ! ===================================================================================== + subroutine UpdateH2OVeg(nsites,sites,bc_out) ! ---------------------------------------------------------------------------------- @@ -838,7 +850,6 @@ subroutine UpdateH2OVeg(nsites,sites,bc_out) integer :: s real(r8) :: balive_patch - do s = 1,nsites bc_out(s)%plant_stored_h2o_si = 0.0_r8 end do @@ -943,17 +954,17 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) csite_hydr%l_aroot_1D = sum( csite_hydr%l_aroot_layer(:)) ! update outer radii of column-level rhizosphere shells (same across patches and cohorts) - do j = 1,nlevsoi_hyd + do j = 1,csite_hydr%nlevsoi_hyd ! proceed only if l_aroot_coh has changed if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then call shellGeom( csite_hydr%l_aroot_layer(j), csite_hydr%rs1(j), AREA, bc_in%dz_sisl(j), & csite_hydr%r_out_shell(j,:), csite_hydr%r_node_shell(j,:),csite_hydr%v_shell(j,:)) end if !has l_aroot_layer changed? enddo - call shellGeom( csite_hydr%l_aroot_1D, csite_hydr%rs1(1), AREA, sum(bc_in%dz_sisl(1:nlevsoi_hyd)), & + call shellGeom( csite_hydr%l_aroot_1D, csite_hydr%rs1(1), AREA, sum(bc_in%dz_sisl(1:csite_hydr%nlevsoi_hyd)), & csite_hydr%r_out_shell_1D(:), csite_hydr%r_node_shell_1D(:), csite_hydr%v_shell_1D(:)) - do j = 1,nlevsoi_hyd + do j = 1,csite_hydr%nlevsoi_hyd hksat_s = bc_in%hksat_sisl(j) * 1.e-3_r8 * 1/grav * 1.e6_r8 @@ -1016,7 +1027,8 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) end subroutine updateSizeDepRhizHydProps - ! ================================================================================= + ! ================================================================================= + subroutine updateSizeDepRhizHydStates(currentSite, bc_in) ! ! !DESCRIPTION: Updates size of 'representative' rhizosphere -- node radii, volumes. @@ -1031,23 +1043,23 @@ subroutine updateSizeDepRhizHydStates(currentSite, bc_in) type(bc_in_type), intent(in) :: bc_in ! ! !LOCAL VARIABLES: - real(r8) :: v_rhiz(nlevsoi_hyd) ! updated volume of all rhizosphere compartments [m3] - real(r8) :: r_delta ! change in radius of innermost rhizosphere compartment [m] - real(r8) :: dpsidr ! water potential gradient near root surface [MPa/m] - real(r8) :: w_shell_new(nlevsoi_hyd,nshell) ! updated water mass in rhizosphere compartment [kg] - real(r8) :: w_layer_init(nlevsoi_hyd) ! initial water mass by layer [kg] - real(r8) :: w_layer_interp(nlevsoi_hyd) ! water mass after interpolating to new rhizosphere [kg] - real(r8) :: w_layer_new(nlevsoi_hyd) ! water mass by layer after interpolation and fudging [kg] - real(r8) :: h2osoi_liq_col_new(nlevsoi_hyd) ! water mass per area after interpolating to new rhizosphere [kg/m2] - real(r8) :: s_shell_init(nlevsoi_hyd,nshell) ! initial saturation fraction in rhizosphere compartment [0-1] - real(r8) :: s_shell_interp(nlevsoi_hyd,nshell) ! interpolated saturation fraction in rhizosphere compartment [0-1] - real(r8) :: psi_shell_init(nlevsoi_hyd,nshell) ! initial water potential in rhizosphere compartment [MPa] - real(r8) :: psi_shell_interp(nlevsoi_hyd,nshell) ! interpolated psi_shell to new r_node_shell [MPa] - real(r8) :: delta_s(nlevsoi_hyd) ! change in saturation fraction needed to ensure water bal [0-1] - real(r8) :: errh2o(nlevsoi_hyd) ! water budget error after updating [kg/m2] - integer :: j,k ! gridcell, column, soil layer, rhizosphere shell indicies - integer :: indexc,indexj ! column and layer indices where there is a water balance error - logical :: found ! flag in search loop + real(r8) :: v_rhiz(nlevsoi_hyd_max) ! updated volume of all rhizosphere compartments [m3] + real(r8) :: r_delta ! change in radius of innermost rhizosphere compartment [m] + real(r8) :: dpsidr ! water potential gradient near root surface [MPa/m] + real(r8) :: w_shell_new ! updated water mass in rhizosphere compartment [kg] + real(r8) :: w_layer_init(nlevsoi_hyd_max) ! initial water mass by layer [kg] + real(r8) :: w_layer_interp(nlevsoi_hyd_max) ! water mass after interpolating to new rhizosphere [kg] + real(r8) :: w_layer_new(nlevsoi_hyd_max) ! water mass by layer after interpolation and fudging [kg] + real(r8) :: h2osoi_liq_col_new(nlevsoi_hyd_max) ! water mass per area after interpolating to new rhizosphere [kg/m2] + real(r8) :: s_shell_init(nlevsoi_hyd_max,nshell) ! initial saturation fraction in rhizosphere compartment [0-1] + real(r8) :: s_shell_interp(nlevsoi_hyd_max,nshell) ! interpolated saturation fraction in rhizosphere compartment [0-1] + real(r8) :: psi_shell_init(nlevsoi_hyd_max,nshell) ! initial water potential in rhizosphere compartment [MPa] + real(r8) :: psi_shell_interp(nlevsoi_hyd_max,nshell) ! interpolated psi_shell to new r_node_shell [MPa] + real(r8) :: delta_s(nlevsoi_hyd_max) ! change in saturation fraction needed to ensure water bal [0-1] + real(r8) :: errh2o(nlevsoi_hyd_max) ! water budget error after updating [kg/m2] + integer :: j,k ! gridcell, column, soil layer, rhizosphere shell indicies + integer :: indexc,indexj ! column and layer indices where there is a water balance error + logical :: found ! flag in search loop type(ed_site_hydr_type), pointer :: csite_hydr !----------------------------------------------------------------------- @@ -1060,7 +1072,7 @@ subroutine updateSizeDepRhizHydStates(currentSite, bc_in) if(.false.) then ! calculate initial s, psi by layer and shell - do j = 1,nlevsoi_hyd + do j = 1, csite_hydr%nlevsoi_hyd ! proceed only if l_aroot_coh has changed if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then select case (iswc) @@ -1089,7 +1101,7 @@ subroutine updateSizeDepRhizHydStates(currentSite, bc_in) ! interpolate initial psi values by layer and shell ! BOC...To-Do: need to constrain psi to be within realistic limits (i.e., < 0) - do j = 1,nlevsoi_hyd + do j = 1,csite_hydr%nlevsoi_hyd ! proceed only if l_aroot_coh has changed if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then @@ -1139,7 +1151,7 @@ subroutine updateSizeDepRhizHydStates(currentSite, bc_in) enddo ! 1st guess at new s based on interpolated psi - do j = 1,nlevsoi_hyd + do j = 1,csite_hydr%nlevsoi_hyd ! proceed only if l_aroot_coh has changed if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then select case (iswc) @@ -1166,7 +1178,7 @@ subroutine updateSizeDepRhizHydStates(currentSite, bc_in) enddo ! accumlate water across shells for each layer (initial and interpolated) - do j = 1,nlevsoi_hyd + do j = 1,csite_hydr%nlevsoi_hyd ! proceed only if l_aroot_coh has changed if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then w_layer_init(j) = 0._r8 @@ -1186,7 +1198,7 @@ subroutine updateSizeDepRhizHydStates(currentSite, bc_in) ! estimate delta_s across all shells needed to ensure total water in each layer doesn't change ! BOC...FIX: need to handle special cases where delta_s causes s_shell to go above or below 1 or 0, respectively. - do j = 1,nlevsoi_hyd + do j = 1,csite_hydr%nlevsoi_hyd ! proceed only if l_aroot_coh has changed if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then delta_s(j) = (( w_layer_init(j) - w_layer_interp(j) )/( v_rhiz(j) * & @@ -1196,7 +1208,7 @@ subroutine updateSizeDepRhizHydStates(currentSite, bc_in) enddo ! update h2osoi_liqvol_shell and h2osoi_liq_shell - do j = 1,nlevsoi_hyd + do j = 1,csite_hydr%nlevsoi_hyd ! proceed only if l_aroot_coh has changed if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then w_layer_new(j) = 0._r8 @@ -1204,16 +1216,16 @@ subroutine updateSizeDepRhizHydStates(currentSite, bc_in) s_shell_interp(j,k) = s_shell_interp(j,k) + delta_s(j) csite_hydr%h2osoi_liqvol_shell(j,k) = s_shell_interp(j,k) * & ( bc_in%watsat_sisl(j)-bc_in%watres_sisl(j) ) + bc_in%watres_sisl(j) - w_shell_new(j,k) = csite_hydr%h2osoi_liqvol_shell(j,k) * & + w_shell_new = csite_hydr%h2osoi_liqvol_shell(j,k) * & csite_hydr%v_shell(j,k) * denh2o - w_layer_new(j) = w_layer_new(j) + w_shell_new(j,k) + w_layer_new(j) = w_layer_new(j) + w_shell_new enddo h2osoi_liq_col_new(j) = w_layer_new(j)/( v_rhiz(j)/bc_in%dz_sisl(j) ) end if !has l_aroot_coh changed? enddo ! balance check - do j = 1,nlevsoi_hyd + do j = 1,csite_hydr%nlevsoi_hyd ! BOC: PLEASE CHECK UNITS ON h2o_liq_sisl(j) (RGK) errh2o(j) = h2osoi_liq_col_new(j) - bc_in%h2o_liq_sisl(j) if (abs(errh2o(j)) > 1.e-9_r8) then @@ -1318,8 +1330,8 @@ subroutine FillDrainRhizShells(nsites, sites, bc_in, bc_out) real(r8) :: dwat_kg ! water remaining to be distributed across shells [kg] real(r8) :: thdiff ! water content difference between ordered adjacent rhiz shells [m3 m-3] real(r8) :: wdiff ! mass of water represented by thdiff over previous k shells [kg] - real(r8) :: errh2o(nlevsoi_hyd) ! water budget error after updating [kg/m2] - real(r8) :: h2osoi_liq_shell(nlevsoi_hyd,nshell) ! + real(r8) :: errh2o(nlevsoi_hyd_max) ! water budget error after updating [kg/m2] + real(r8) :: h2osoi_liq_shell(nlevsoi_hyd_max,nshell) ! integer :: tmp ! temporary logical :: found ! flag in search loop !----------------------------------------------------------------------- @@ -1336,12 +1348,12 @@ subroutine FillDrainRhizShells(nsites, sites, bc_in, bc_out) ! BOC: This was previously in HydrologyDrainage: csite_hydr => sites(s)%si_hydr - - do j = 1,nlevsoi_hyd + + do j = 1,csite_hydr%nlevsoi_hyd - if(nlevsoi_hyd == 1) then - dwat_kgm2 = bc_in(s)%h2o_liq_sisl(hlm_numlevsoil) - csite_hydr%h2osoi_liq_prev(nlevsoi_hyd) - else if(nlevsoi_hyd == hlm_numlevsoil) then + if(csite_hydr%nlevsoi_hyd == 1) then + dwat_kgm2 = bc_in(s)%h2o_liq_sisl(bc_in(s)%nlevsoil) - csite_hydr%h2osoi_liq_prev(csite_hydr%nlevsoi_hyd) + else if(csite_hydr%nlevsoi_hyd == bc_in(s)%nlevsoil ) then dwat_kgm2 = bc_in(s)%h2o_liq_sisl(j) - csite_hydr%h2osoi_liq_prev(j) end if @@ -1409,8 +1421,8 @@ subroutine FillDrainRhizShells(nsites, sites, bc_in, bc_out) enddo ! balance check - if(nlevsoi_hyd == hlm_numlevsoil) then - do j = 1,nlevsoi_hyd + if(csite_hydr%nlevsoi_hyd == bc_in(s)%nlevsoil) then + do j = 1,csite_hydr%nlevsoi_hyd errh2o(j) = sum(h2osoi_liq_shell(j,:))/AREA - bc_in(s)%h2o_liq_sisl(j) if (abs(errh2o(j)) > 1.e-9_r8) then @@ -1423,8 +1435,8 @@ subroutine FillDrainRhizShells(nsites, sites, bc_in, bc_out) end if end if enddo - else if(nlevsoi_hyd == 1) then - errh2o(nlevsoi_hyd) = sum(h2osoi_liq_shell(nlevsoi_hyd,:))/AREA - sum( bc_in(s)%h2o_liq_sisl(:) ) + else if(csite_hydr%nlevsoi_hyd == 1) then + errh2o(csite_hydr%nlevsoi_hyd) = sum(h2osoi_liq_shell(csite_hydr%nlevsoi_hyd,:))/AREA - sum( bc_in(s)%h2o_liq_sisl(:) ) end if end do @@ -1450,7 +1462,6 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ! !DESCRIPTION: !s ! !USES: - use FatesInterfaceMod , only : hlm_numlevsoil use EDTypesMod , only : AREA ! ARGUMENTS: @@ -1527,7 +1538,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) real(r8) :: h2osoi_liqvol real(r8) :: dz_tot ! total soil depth (to bottom of bottom layer) [m] real(r8) :: l_aroot_tot_col ! total length of absorbing roots across all soil layers [m] - real(r8) :: dth_layershell_col(nlevsoi_hyd,nshell) ! accumulated water content change over all cohorts in a column [m3 m-3] + real(r8) :: dth_layershell_col(nlevsoi_hyd_max,nshell) ! accumulated water content change over all cohorts in a column [m3 m-3] real(r8) :: ths_shell_1D(nshell) ! saturated water content of rhizosphere compartment [m3 m-3] real(r8) :: thr_shell_1D(nshell) ! residual water content of rhizosphere compartment [m3 m-3] real(r8) :: kmax_bound_shell_1l(nshell) ! like kmax_bound_shell_1D(:) but for specific single soil layer [kg s-1 MPa-1] @@ -1538,17 +1549,17 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) real(r8) :: psi_node_aroot_1D ! water potential of absorbing root [MPa] ! hydraulics conductances - real(r8) :: kmax_bound_bylayershell(nlevsoi_hyd,nshell) ! maximum conductance at shell boundaries in each layer [kg s-1 MPa-1] + real(r8) :: kmax_bound_bylayershell(nlevsoi_hyd_max,nshell) ! maximum conductance at shell boundaries in each layer [kg s-1 MPa-1] real(r8) :: kmax_bound_aroot_soil1 ! maximum radial conductance of absorbing roots [kg s-1 MPa-1] real(r8) :: kmax_bound_aroot_soil2 ! maximum conductance to root surface from innermost rhiz shell [kg s-1 MPa-1] - real(r8) :: ksoil_bylayer(nlevsoi_hyd) ! total rhizosphere conductance (over all shells) by soil layer [MPa] + real(r8) :: ksoil_bylayer(nlevsoi_hyd_max) ! total rhizosphere conductance (over all shells) by soil layer [MPa] real(r8) :: ksoil_tot ! total rhizosphere conductance (over all shells and soil layers [MPa] - real(r8) :: kbg_layer(nlevsoi_hyd) ! total absorbing root & rhizosphere conductance (over all shells) by soil layer [MPa] + real(r8) :: kbg_layer(nlevsoi_hyd_max) ! total absorbing root & rhizosphere conductance (over all shells) by soil layer [MPa] real(r8) :: kbg_tot ! total absorbing root & rhizosphere conductance (over all shells and soil layers [MPa] real(r8) :: kmax_stem ! maximum whole-stem (above troot to leaf) conductance [kg s-1 MPa-1] ! hydraulics other - integer :: ordered(nlevsoi_hyd) = (/(j,j=1,nlevsoi_hyd,1)/) ! array of soil layer indices which have been ordered + integer :: ordered(nlevsoi_hyd_max) = (/(j,j=1,nlevsoi_hyd_max,1)/) ! array of soil layer indices which have been ordered real(r8) :: qflx_tran_veg_indiv ! individiual transpiration rate [kgh2o indiv-1 s-1] real(r8) :: qflx_tran_veg_patch_coh real(r8) :: gscan_patch ! sum of ccohort%gscan across all cohorts within a patch @@ -1688,11 +1699,11 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) sum(ccohort_hydr%th_troot(:)*ccohort_hydr%v_troot(:)) + & sum(ccohort_hydr%th_aroot(:)*ccohort_hydr%v_aroot_layer(:)))* & denh2o*ccohort%n/AREA/dtime - if(nlevsoi_hyd == 1) then + if( site_hydr%nlevsoi_hyd == 1) then site_hydr%recruit_w_uptake(1) = site_hydr%recruit_w_uptake(1)+ & recruitw else - do j=1,nlevsoi_hyd + do j=1,site_hydr%nlevsoi_hyd if(j == 1) then rootfr = zeng2001_crootfr(roota, rootb, bc_in(s)%zi_sisl(j)) else @@ -1719,7 +1730,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) - if(nlevsoi_hyd > 1) then + if(site_hydr%nlevsoi_hyd > 1) then ! BUCKET APPROXIMATION OF THE SOIL-ROOT HYDRAULIC GRADIENT (weighted average across layers) !call map2d_to_1d_shells(soilstate_inst, waterstate_inst, g, c, rs1(c,1), ccohort_hydr%l_aroot_layer*ccohort%n, & ! (2._r8 * ccohort_hydr%kmax_treebg_layer(:)), ths_shell_1D, & @@ -1730,7 +1741,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ! REPRESENTATIVE SINGLE FINE ROOT POOL (weighted average across layers) !call map2d_to_1d_aroot(ft, ccohort, kmax_bound_bylayershell, ths_aroot_1D, thr_aroot_1D, vtot_aroot_1D, psi_node_aroot_1D) !psi_node(n_hypool_ag+n_hypool_troot+1) = psi_node_aroot_1D - else if(nlevsoi_hyd == 1) then + else if(site_hydr%nlevsoi_hyd == 1) then write(fates_log(),*) 'Single layer hydraulics currently inoperative nlevsoi_hyd==1' call endrun(msg=errMsg(sourcefile, __LINE__)) !psi_node( (n_hypool_tot-nshell+1):n_hypool_tot) = psisoi_liq_shell(c,1,:) @@ -1744,14 +1755,14 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) z_node((n_hypool_ag+n_hypool_troot+1): n_hypool_tot) = ccohort_hydr%z_node_aroot(1) ! absorbing root and rhizosphere shells v_node( 1 : n_hypool_ag) = ccohort_hydr%v_ag(:) ! leaf and stem v_node( (n_hypool_ag+1):(n_hypool_ag+n_hypool_troot)) = ccohort_hydr%v_troot(:) ! transporting root - if(nlevsoi_hyd == 1) then + if(site_hydr%nlevsoi_hyd == 1) then v_node((n_hypool_ag+n_hypool_troot+1) ) = ccohort_hydr%v_aroot_tot ! absorbing root v_node( (n_hypool_tot-nshell+1): n_hypool_tot) = & site_hydr%v_shell_1D(:)*ccohort_hydr%l_aroot_tot/sum(bc_in(s)%dz_sisl(:)) ! rhizosphere shells end if ! SET SATURATED & RESIDUAL WATER CONTENTS - if(nlevsoi_hyd == 1) then + if(site_hydr%nlevsoi_hyd == 1) then ths_node( (n_hypool_tot-nshell+1):n_hypool_tot) = bc_in(s)%watsat_sisl(1) !! BOC... should the below code exist on HLM side? watres_col is a new SWC parameter 1 ! introduced for the van Genuchten, but does not exist for Campbell SWC. @@ -1787,7 +1798,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) kmax_upper( 1 : n_hypool_ag ) = ccohort_hydr%kmax_upper(:) kmax_lower( 1 : n_hypool_ag ) = ccohort_hydr%kmax_lower(:) kmax_upper(( n_hypool_ag+1) ) = ccohort_hydr%kmax_upper_troot - if(nlevsoi_hyd == 1) then + if(site_hydr%nlevsoi_hyd == 1) then !! estimate troot-aroot and aroot-radial components as a residual: !! 25% each of total (surface of aroots to leaves) resistance kmax_bound(( n_hypool_ag+1):(n_hypool_ag+2 )) = 2._r8 * ccohort_hydr%kmax_treebg_tot @@ -1807,7 +1818,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ccohort_hydr%l_aroot_tot / site_hydr%l_aroot_1D end if - if(nlevsoi_hyd == 1) then + if(site_hydr%nlevsoi_hyd == 1) then ! CONVERT WATER POTENTIALS TO WATER CONTENTS FOR THE NEW 'BUCKET' ! RHIZOSPHERE (fine roots and rhizosphere shells) do k = (n_hypool_tot - nshell), n_hypool_tot @@ -1829,7 +1840,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ! end if ! end do - if(nlevsoi_hyd == 1) then + if(site_hydr%nlevsoi_hyd == 1) then ! 1-D THETA-BASED SOLUTION TO RICHARDS' EQUATION call Hydraulics_1DSolve(ccohort, ft, z_node, v_node, ths_node, & thr_node, kmax_bound, kmax_upper, kmax_lower, & @@ -1867,11 +1878,11 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) site_hydr%errh2o_hyd = site_hydr%errh2o_hyd + ccohort_hydr%errh2o*ccohort%c_area /AREA !*patch_wgt ! ACCUMULATE CHANGE IN SOIL WATER CONTENT OF EACH COHORT TO COLUMN-LEVEL - dth_layershell_col(nlevsoi_hyd,:) = dth_layershell_col(nlevsoi_hyd,:) + & + dth_layershell_col(site_hydr%nlevsoi_hyd,:) = dth_layershell_col(site_hydr%nlevsoi_hyd,:) + & dth_node((n_hypool_tot-nshell+1):n_hypool_tot) * & ccohort_hydr%l_aroot_tot * ccohort%n / site_hydr%l_aroot_1D! * & !patch_wgt - else if(nlevsoi_hyd > 1) then + else if(site_hydr%nlevsoi_hyd > 1) then ! VERTICAL LAYER CONTRIBUTION TO TOTAL ROOT WATER UPTAKE OR LOSS ! _____ ! | | @@ -1912,7 +1923,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ! layers have transporting-to-absorbing root water potential gradients of opposite sign kbg_layer(:) = 0._r8 kbg_tot = 0._r8 - do j=1,nlevsoi_hyd + do j=1,site_hydr%nlevsoi_hyd z_node_1l((n_hypool_ag+n_hypool_troot+1):(n_hypool_tot)) = bc_in(s)%z_sisl(j) v_node_1l((n_hypool_ag+n_hypool_troot+1) ) = ccohort_hydr%v_aroot_layer(j) v_node_1l((n_hypool_tot-nshell+1):(n_hypool_tot)) = site_hydr%v_shell(j,:) * & @@ -1988,7 +1999,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ! order soil layers in terms of decreasing volumetric water content ! algorithm same as that used in histFileMod.F90 to alphabetize history tape contents - do j = nlevsoi_hyd-1,1,-1 + do j = site_hydr%nlevsoi_hyd-1,1,-1 do jj = 1,j if (kbg_layer(ordered(jj)) <= kbg_layer(ordered(jj+1))) then tmp = ordered(jj) @@ -2010,7 +2021,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ccohort_hydr%errh2o = 0._r8 ! do j=1,nlevsoi_hyd ! replace j with ordered(jj) in order ! to go through soil layers in order of decreasing total root-soil conductance - do jj=1,size(ordered) + do jj=1,site_hydr%nlevsoi_hyd !initialize state variables in absorbing roots and rhizosphere shells in each soil layer !z_node_1l( ( n_hypool_ag+1):(n_hypool_troot )) = -bc_in(s)%z_sisl(ordered(jj)) @@ -2155,7 +2166,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ! UPDATE WATER CONTENT & POTENTIAL IN LEAVES, STEM, AND TROOT (COHORT-LEVEL) do k=1,n_hypool_ag - if(nlevsoi_hyd == 1) then + if(site_hydr%nlevsoi_hyd == 1) then ccohort_hydr%th_ag(k) = th_node(k) else ccohort_hydr%th_ag(k) = th_node_1l(k) @@ -2166,7 +2177,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ccohort_hydr%flc_ag(k), site_hydr, bc_in(s) ) enddo do k=(n_hypool_ag+1),(n_hypool_ag+n_hypool_troot) - if(nlevsoi_hyd == 1) then + if(site_hydr%nlevsoi_hyd == 1) then ccohort_hydr%th_troot(k-n_hypool_ag) = th_node(k) else ccohort_hydr%th_troot(k-n_hypool_ag) = th_node_1l(k) @@ -2202,7 +2213,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) (ccohort_hydr%flc_troot(k) - ccohort_hydr%flc_min_troot(k))*exp(-refill_rate*dtime) end if end do - do j=1,nlevsoi_hyd + do j=1,site_hydr%nlevsoi_hyd ccohort_hydr%flc_min_aroot(j) = min(ccohort_hydr%flc_min_aroot(j), ccohort_hydr%flc_aroot(j)) if(ccohort_hydr%psi_aroot(j) >= ccohort_hydr%refill_thresh .and. & ccohort_hydr%flc_aroot(j) > ccohort_hydr%flc_min_aroot(j)) then ! then refilling @@ -2219,7 +2230,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) ! UPDATE THE COLUMN-LEVEL SOIL WATER CONTENT (LAYER x SHELL) site_hydr%supsub_flag(:) = 999 - do j=1,nlevsoi_hyd + do j=1,site_hydr%nlevsoi_hyd !! BOC... should the below code exist on HLM side? watres_col is a new SWC parameter ! introduced for the van Genuchten, but does not exist for Campbell SWC. select case (iswc) @@ -2260,23 +2271,23 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) site_hydr%psisoi_liq_innershell(j) = smp - if(nlevsoi_hyd == 1) then - bc_out(s)%qflx_soil2root_sisl(1:hlm_numlevsoil-1) = 0._r8 + if(site_hydr%nlevsoi_hyd == 1) then + bc_out(s)%qflx_soil2root_sisl(1:bc_in(s)%nlevsoil-1) = 0._r8 - ! qflx_rootsoi(c,hlm_numlevsoil) = + ! qflx_rootsoi(c,bc_in(s)%nlevsoil) = ! -(sum(dth_layershell_col(j,:))*bc_in(s)%dz_sisl(j)*denh2o/dtime) - bc_out(s)%qflx_soil2root_sisl(hlm_numlevsoil) = & + bc_out(s)%qflx_soil2root_sisl(bc_in(s)%nlevsoil) = & -(sum(dth_layershell_col(j,:)*site_hydr%v_shell_1D(:)) * & site_hydr%l_aroot_1D/bc_in(s)%dz_sisl(j)/AREA*denh2o/dtime)- & - site_hydr%recruit_w_uptake(nlevsoi_hyd) + site_hydr%recruit_w_uptake(site_hydr%nlevsoi_hyd) -! h2osoi_liqvol = min(bc_in(s)%eff_porosity_gl(hlm_numlevsoil), & -! bc_in(s)%h2o_liq_sisl(hlm_numlevsoil)/(bc_in(s)%dz_sisl(hlm_numlevsoil)*denh2o)) +! h2osoi_liqvol = min(bc_in(s)%eff_porosity_sl(bc_in(s)%nlevsoil), & +! bc_in(s)%h2o_liq_sisl(bc_in(s)%nlevsoil)/(bc_in(s)%dz_sisl(bc_in(s)%nlevsoil)*denh2o)) ! Save the amount of liquid soil water known to the model after root uptake - site_hydr%h2osoi_liq_prev(nlevsoi_hyd) = bc_in(s)%h2o_liq_sisl(hlm_numlevsoil) - & - dtime*bc_out(s)%qflx_soil2root_sisl(hlm_numlevsoil) + site_hydr%h2osoi_liq_prev(site_hydr%nlevsoi_hyd) = bc_in(s)%h2o_liq_sisl(bc_in(s)%nlevsoil) - & + dtime*bc_out(s)%qflx_soil2root_sisl(bc_in(s)%nlevsoil) else !qflx_rootsoi(c,j) = -(sum(dth_layershell_col(j,:))*bc_in(s)%dz_sisl(j)*denh2o/dtime) @@ -2285,7 +2296,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) site_hydr%l_aroot_layer(j)/bc_in(s)%dz_sisl(j)/AREA*denh2o/dtime)- & site_hydr%recruit_w_uptake(j) - ! h2osoi_liqvol = min(bc_in(s)%eff_porosity_gl(j), & + ! h2osoi_liqvol = min(bc_in(s)%eff_porosity_sl(j), & ! bc_in(s)%h2o_liq_sisl(j)/(bc_in(s)%dz_sisl(j)*denh2o)) ! Save the amount of liquid soil water known to the model after root uptake @@ -2295,7 +2306,7 @@ subroutine hydraulics_bc ( nsites, sites,bc_in,bc_out,dtime ) end if - enddo !nlevsoi_hyd + enddo !site_hydr%nlevsoi_hyd !----------------------------------------------------------------------- ! mass balance check and pass the total stored vegetation water to HLM ! in order for it to fill its balance checks diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 3c23a41d8a..817adf0850 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -283,8 +283,8 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) ! (size --> heights of elements --> hydraulic path lengths --> ! maximum node-to-node conductances) if( (hlm_use_planthydro.eq.itrue) .and. do_growthrecruiteffects) then - call updateSizeDepTreeHydProps(currentCohort, bc_in) - call updateSizeDepTreeHydStates(currentCohort) + call updateSizeDepTreeHydProps(currentSite,currentCohort, bc_in) + call updateSizeDepTreeHydStates(currentSite,currentCohort) end if currentCohort => currentCohort%taller diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 5139c4a266..e05255ef6a 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -2408,7 +2408,6 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) AREA use FatesHydraulicsMemMod, only : ed_cohort_hydr_type - use FatesHydraulicsMemMod, only : nlevsoi_hyd use EDTypesMod , only : maxpft @@ -2443,6 +2442,17 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) real(r8), parameter :: daysecs = 86400.0_r8 ! What modeler doesn't recognize 86400? real(r8), parameter :: yeardays = 365.0_r8 ! Should this be 365.25? + + logical :: layer1_present + logical :: layer2_present + logical :: layer3_present + logical :: layer4_present + logical :: layer5_present + logical :: layer6_present + logical :: layer7_present + logical :: layer8_present + logical :: layer9_present + logical :: layer10_present if(hlm_use_planthydro.eq.ifalse) return @@ -2490,6 +2500,69 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) ! This is actually used as a check ! on hio_nplant_si_scpf + + ! Determine which hydraulic soil layers are present + if( sites(s)%si_hydr%nlevsoi_hyd >=1 ) then + layer1_present = .true. + else + layer1_present = .false. + end if + ! Determine which hydraulic soil layers are present + if( sites(s)%si_hydr%nlevsoi_hyd >=2 ) then + layer2_present = .true. + else + layer2_present = .false. + end if + ! Determine which hydraulic soil layers are present + if( sites(s)%si_hydr%nlevsoi_hyd >=3 ) then + layer3_present = .true. + else + layer3_present = .false. + end if + ! Determine which hydraulic soil layers are present + if( sites(s)%si_hydr%nlevsoi_hyd >=4 ) then + layer4_present = .true. + else + layer4_present = .false. + end if + ! Determine which hydraulic soil layers are present + if( sites(s)%si_hydr%nlevsoi_hyd >=5 ) then + layer5_present = .true. + else + layer5_present = .false. + end if + ! Determine which hydraulic soil layers are present + if( sites(s)%si_hydr%nlevsoi_hyd >=6 ) then + layer6_present = .true. + else + layer6_present = .false. + end if + ! Determine which hydraulic soil layers are present + if( sites(s)%si_hydr%nlevsoi_hyd >=7 ) then + layer7_present = .true. + else + layer7_present = .false. + end if + ! Determine which hydraulic soil layers are present + if( sites(s)%si_hydr%nlevsoi_hyd >=8 ) then + layer8_present = .true. + else + layer8_present = .false. + end if + ! Determine which hydraulic soil layers are present + if( sites(s)%si_hydr%nlevsoi_hyd >=9 ) then + layer9_present = .true. + else + layer9_present = .false. + end if + ! Determine which hydraulic soil layers are present + if( sites(s)%si_hydr%nlevsoi_hyd >=10 ) then + layer10_present = .true. + else + layer10_present = .false. + end if + + cpatch => sites(s)%oldest_patch do while(associated(cpatch)) ccohort => cpatch%shortest @@ -2545,37 +2618,49 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) hio_rootuptake_scpf(io_si,iscpf) = hio_rootuptake_scpf(io_si,iscpf) + & ccohort_hydr%rootuptake * number_fraction_rate ! [kg/indiv/s] - if(nlevsoi_hyd == 10) then + + ! Not sure how to simplify this + ! All of these if's inside a cohort loop is not good.... + + if (layer1_present)then hio_rootuptake01_scpf(io_si,iscpf) = hio_rootuptake01_scpf(io_si,iscpf) + & ccohort_hydr%rootuptake01 * number_fraction_rate ! [kg/indiv/s] - + end if + if (layer2_present) then hio_rootuptake02_scpf(io_si,iscpf) = hio_rootuptake02_scpf(io_si,iscpf) + & ccohort_hydr%rootuptake02 * number_fraction_rate ! [kg/indiv/s] - + end if + if (layer3_present) then hio_rootuptake03_scpf(io_si,iscpf) = hio_rootuptake03_scpf(io_si,iscpf) + & ccohort_hydr%rootuptake03 * number_fraction_rate ! [kg/indiv/s] - + end if + if (layer4_present) then hio_rootuptake04_scpf(io_si,iscpf) = hio_rootuptake04_scpf(io_si,iscpf) + & ccohort_hydr%rootuptake04 * number_fraction_rate ! [kg/indiv/s] - + end if + if (layer5_present) then hio_rootuptake05_scpf(io_si,iscpf) = hio_rootuptake05_scpf(io_si,iscpf) + & ccohort_hydr%rootuptake05 * number_fraction_rate ! [kg/indiv/s] - + end if + if (layer6_present) then hio_rootuptake06_scpf(io_si,iscpf) = hio_rootuptake06_scpf(io_si,iscpf) + & ccohort_hydr%rootuptake06 * number_fraction_rate ! [kg/indiv/s] - + end if + if (layer7_present) then hio_rootuptake07_scpf(io_si,iscpf) = hio_rootuptake07_scpf(io_si,iscpf) + & ccohort_hydr%rootuptake07 * number_fraction_rate ! [kg/indiv/s] - + end if + if (layer8_present) then hio_rootuptake08_scpf(io_si,iscpf) = hio_rootuptake08_scpf(io_si,iscpf) + & ccohort_hydr%rootuptake08 * number_fraction_rate ! [kg/indiv/s] - + end if + if (layer9_present) then hio_rootuptake09_scpf(io_si,iscpf) = hio_rootuptake09_scpf(io_si,iscpf) + & - ccohort_hydr%rootuptake09 * number_fraction_rate ! [kg/indiv/s] - + ccohort_hydr%rootuptake09 * number_fraction_rate ! [kg/indiv/s] + end if + if (layer10_present) then hio_rootuptake10_scpf(io_si,iscpf) = hio_rootuptake10_scpf(io_si,iscpf) + & ccohort_hydr%rootuptake10 * number_fraction_rate ! [kg/indiv/s] - end if hio_sapflow_scpf(io_si,iscpf) = hio_sapflow_scpf(io_si,iscpf) + & diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90 index 469ff386bd..fa386fde10 100644 --- a/main/FatesHydraulicsMemMod.F90 +++ b/main/FatesHydraulicsMemMod.F90 @@ -164,7 +164,6 @@ module FatesHydraulicsMemMod ! Heights are referenced to soil surface (+ = above; - = below) real(r8) :: z_node_ag(n_hypool_ag) ! nodal height of aboveground water storage compartments [m] real(r8) :: z_node_troot(n_hypool_troot) ! nodal height of belowground water storage compartments [m] - real(r8) :: z_node_aroot(nlevsoi_hyd) ! nodal height of absorbing root water storage compartments [m] real(r8) :: z_upper_ag(n_hypool_ag) ! upper boundary height of aboveground water storage compartments [m] real(r8) :: z_upper_troot(n_hypool_troot) ! upper boundary height of belowground water storage compartments [m] real(r8) :: z_lower_ag(n_hypool_ag) ! lower boundary height of aboveground water storage compartments [m] @@ -174,39 +173,34 @@ module FatesHydraulicsMemMod real(r8) :: kmax_upper_troot ! maximum hydraulic conductance from troot node to upper boundary [kg s-1 MPa-1] real(r8) :: kmax_bound(n_hypool_ag) ! maximum hydraulic conductance at lower boundary (canopy to troot) [kg s-1 MPa-1] real(r8) :: kmax_treebg_tot ! total belowground tree kmax (troot to surface of absorbing roots) [kg s-1 MPa-1] - real(r8) :: kmax_treebg_layer(nlevsoi_hyd) ! total belowground tree kmax partitioned by soil layer [kg s-1 MPa-1] real(r8) :: v_ag_init(n_hypool_ag) ! previous day's volume of aboveground water storage compartments [m3] real(r8) :: v_ag(n_hypool_ag) ! volume of aboveground water storage compartments [m3] real(r8) :: v_troot_init(n_hypool_troot) ! previous day's volume of belowground water storage compartments [m3] real(r8) :: v_troot(n_hypool_troot) ! volume of belowground water storage compartments [m3] real(r8) :: v_aroot_tot ! total volume of absorbing roots [m3] - real(r8) :: v_aroot_layer_init(nlevsoi_hyd) ! previous day's volume of absorbing roots by soil layer [m3] - real(r8) :: v_aroot_layer(nlevsoi_hyd) ! volume of absorbing roots by soil layer [m3] real(r8) :: l_aroot_tot ! total length of absorbing roots [m] - real(r8) :: l_aroot_layer(nlevsoi_hyd) ! length of absorbing roots by soil layer [m] - - real(r8),allocatable :: z_node_aroot(nlevsoi_hyd) ! nodal height of absorbing root water storage compartments [m] + + real(r8),allocatable :: z_node_aroot(:) ! nodal height of absorbing root water storage compartments [m] + real(r8),allocatable :: kmax_treebg_layer(:) ! total belowground tree kmax partitioned by soil layer [kg s-1 MPa-1] + real(r8),allocatable :: v_aroot_layer_init(:) ! previous day's volume of absorbing roots by soil layer [m3] + real(r8),allocatable :: v_aroot_layer(:) ! volume of absorbing roots by soil layer [m3] + real(r8),allocatable :: l_aroot_layer(:) ! length of absorbing roots by soil layer [m] ! BC PLANT HYDRAULICS - state variables real(r8) :: th_ag(n_hypool_ag) ! water in aboveground compartments [kgh2o/indiv] real(r8) :: th_troot(n_hypool_troot) ! water in belowground compartments [kgh2o/indiv] - real(r8) :: th_aroot(nlevsoi_hyd) ! water in absorbing roots [kgh2o/indiv] real(r8) :: lwp_mem(numLWPmem) ! leaf water potential over the previous numLWPmem timesteps [MPa] real(r8) :: lwp_stable ! leaf water potential just before it became unstable [MPa] logical :: lwp_is_unstable ! flag for instability of leaf water potential over previous timesteps real(r8) :: psi_ag(n_hypool_ag) ! water potential in aboveground compartments [MPa] real(r8) :: psi_troot(n_hypool_troot) ! water potential in belowground compartments [MPa] - real(r8) :: psi_aroot(nlevsoi_hyd) ! water potential in absorbing roots [MPa] real(r8) :: flc_ag(n_hypool_ag) ! fractional loss of conductivity in aboveground compartments [-] real(r8) :: flc_troot(n_hypool_troot) ! fractional loss of conductivity in belowground compartments [-] - real(r8) :: flc_aroot(nlevsoi_hyd) ! fractional loss of conductivity in absorbing roots [-] real(r8) :: flc_min_ag(n_hypool_ag) ! min attained fractional loss of conductivity in ! aboveground compartments (for tracking xylem refilling dynamics) [-] real(r8) :: flc_min_troot(n_hypool_troot) ! min attained fractional loss of conductivity in ! belowground compartments (for tracking xylem refilling dynamics) [-] - real(r8) :: flc_min_aroot(nlevsoi_hyd) ! min attained fractional loss of conductivity in absorbing roots - ! (for tracking xylem refilling dynamics) [-] real(r8) :: refill_thresh ! water potential threshold for xylem refilling to occur [MPa] real(r8) :: refill_days ! number of days required for 50% of xylem refilling to occur [days] real(r8) :: btran(nlevcan_hyd) ! leaf water potential limitation on gs [0-1] @@ -215,7 +209,13 @@ module FatesHydraulicsMemMod real(r8) :: iterh1 ! number of iterations required to achieve tolerable water balance error real(r8) :: iterh2 ! number of inner iterations real(r8) :: errh2o ! total water balance error per unit crown area [kgh2o/m2] - + + real(r8),allocatable :: th_aroot(:) ! water in absorbing roots [kgh2o/indiv] + real(r8),allocatable :: psi_aroot(:) ! water potential in absorbing roots [MPa] + real(r8),allocatable :: flc_aroot(:) ! fractional loss of conductivity in absorbing roots [-] + real(r8),allocatable :: flc_min_aroot(:) ! min attained fractional loss of conductivity in absorbing roots + ! (for tracking xylem refilling dynamics) [-] + ! BC PLANT HYDRAULICS - fluxes real(r8) :: qtop_dt ! transpiration boundary condition (+ to atm) [kg/indiv/timestep] real(r8) :: dqtopdth_dthdt ! transpiration tendency term (+ to atm) [kg/indiv/timestep] @@ -234,9 +234,52 @@ module FatesHydraulicsMemMod real(r8) :: rootuptake10 ! net flow into roots (+ into roots), soil layer 10 [kg/indiv/timestep] ! BC PLANT HYDRAULICS - flags logical :: is_newly_recuited !whether the new cohort is newly recuited + + contains + + procedure :: AllocateHydrCohortArrays + procedure :: DeallocateHydrCohortArrays + end type ed_cohort_hydr_type contains + + subroutine AllocateHydrCohortArrays(this,nlevsoil_hydr) + + ! Arguments + class(ed_cohort_hydr_type),intent(inout) :: this + integer, intent(in) :: nlevsoil_hydr + + allocate(this%z_node_aroot(1:nlevsoil_hydr)) + allocate(this%kmax_treebg_layer(1:nlevsoil_hydr)) + allocate(this%v_aroot_layer_init(1:nlevsoil_hydr)) + allocate(this%v_aroot_layer(1:nlevsoil_hydr)) + allocate(this%l_aroot_layer(1:nlevsoil_hydr)) + allocate(this%th_aroot(1:nlevsoil_hydr)) + allocate(this%psi_aroot(1:nlevsoil_hydr)) + allocate(this%flc_aroot(1:nlevsoil_hydr)) + allocate(this%flc_min_aroot(1:nlevsoil_hydr)) + + return + end subroutine AllocateHydrCohortArrays + + ! =================================================================================== + + subroutine DeallocateHydrCohortArrays(this) + + class(ed_cohort_hydr_type),intent(inout) :: this + deallocate(this%z_node_aroot) + deallocate(this%kmax_treebg_layer) + deallocate(this%v_aroot_layer_init) + deallocate(this%v_aroot_layer) + deallocate(this%l_aroot_layer) + deallocate(this%th_aroot) + deallocate(this%psi_aroot) + deallocate(this%flc_aroot) + deallocate(this%flc_min_aroot) + + return + end subroutine DeallocateHydrCohortArrays ! =================================================================================== @@ -246,35 +289,38 @@ subroutine InitHydrSite(this) class(ed_site_hydr_type),intent(inout) :: this associate( nlevsoil_hyd => this%nlevsoi_hyd ) + + allocate(this%v_shell(1:nlevsoil_hyd,1:nshell)) ; this%v_shell = nan + allocate(this%v_shell_init(1:nlevsoil_hyd,1:nshell)) ; this%v_shell_init = nan + allocate(this%v_shell_1D(1:nshell)) ; this%v_shell_1D = nan + allocate(this%r_node_shell(1:nlevsoil_hyd,1:nshell)) ; this%r_node_shell = nan + allocate(this%r_node_shell_init(1:nlevsoil_hyd,1:nshell)); this%r_node_shell_init = nan + allocate(this%r_out_shell(1:nlevsoil_hyd,1:nshell)) ; this%r_out_shell = nan + allocate(this%l_aroot_layer(1:nlevsoil_hyd)) ; this%l_aroot_layer = nan + allocate(this%l_aroot_layer_init(1:nlevsoil_hyd)) ; this%l_aroot_layer_init = nan + allocate(this%kmax_upper_shell(1:nlevsoil_hyd,1:nshell)); this%kmax_upper_shell = nan + allocate(this%kmax_bound_shell(1:nlevsoil_hyd,1:nshell)); this%kmax_bound_shell = nan + allocate(this%kmax_lower_shell(1:nlevsoil_hyd,1:nshell)); this%kmax_lower_shell = nan + allocate(this%r_out_shell_1D(1:nshell)) ; this%r_out_shell_1D = nan + allocate(this%r_node_shell_1D(1:nshell)) ; this%r_node_shell_1D = nan + allocate(this%kmax_upper_shell_1D(1:nshell)) ; this%kmax_upper_shell_1D = nan + allocate(this%kmax_bound_shell_1D(1:nshell)) ; this%kmax_bound_shell_1D = nan + allocate(this%kmax_lower_shell_1D(1:nshell)) ; this%kmax_lower_shell_1D = nan + allocate(this%supsub_flag(1:nlevsoil_hyd)) ; this%supsub_flag = -999 + allocate(this%h2osoi_liqvol_shell(1:nlevsoil_hyd,1:nshell)) ; this%h2osoi_liqvol_shell = nan + allocate(this%h2osoi_liq_prev(1:nlevsoil_hyd)) ; this%h2osoi_liq_prev = nan + allocate(this%psisoi_liq_innershell(1:nlevsoil_hyd)); this%psisoi_liq_innershell = nan + allocate(this%rs1(1:nlevsoil_hyd)); this%rs1(:) = fine_root_radius_const + allocate(this%recruit_w_uptake(1:nlevsoil_hyd)); this%recruit_w_uptake = nan + + this%l_aroot_1D = nan + this%errh2o_hyd = nan + this%dwat_veg = nan + this%h2oveg = nan + this%h2oveg_dead = 0.0_r8 + + end associate - allocate(this%v_shell(1:nlevsoil_hyd,1:nshell)) ; this%v_shell = nan - allocate(this%v_shell_init(1:nlevsoil_hyd,1:nshell)) ; this%v_shell_init = nan - allocate(this%v_shell_1D(1:nshell)) ; this%v_shell_1D = nan - allocate(this%r_node_shell(1:nlevsoil_hyd,1:nshell)) ; this%r_node_shell = nan - allocate(this%r_node_shell_init(1:nlevsoil_hyd,1:nshell)); this%r_node_shell_init = nan - allocate(this%r_out_shell(1:nlevsoil_hyd,1:nshell)) ; this%r_out_shell = nan - allocate(this%l_aroot_layer(1:nlevsoil_hyd)) ; this%l_aroot_layer = nan - allocate(this%l_aroot_layer_init(1:nlevsoil_hyd)) ; this%l_aroot_layer_init = nan - allocate(this%kmax_upper_shell(1:nlevsoil_hyd,1:nshell)); this%kmax_upper_shell = nan - allocate(this%kmax_bound_shell(1:nlevsoil_hyd,1:nshell)); this%kmax_bound_shell = nan - allocate(this%kmax_lower_shell(1:nlevsoil_hyd,1:nshell)); this%kmax_lower_shell = nan - allocate(this%r_out_shell_1D(1:nshell)) ; this%r_out_shell_1D = nan - allocate(this%r_node_shell_1D(1:nshell)) ; this%r_node_shell_1D = nan - allocate(this%kmax_upper_shell_1D(1:nshell)) ; this%kmax_upper_shell_1D = nan - allocate(this%kmax_bound_shell_1D(1:nshell)) ; this%kmax_bound_shell_1D = nan - allocate(this%kmax_lower_shell_1D(1:nshell)) ; this%kmax_lower_shell_1D = nan - allocate(this%supsub_flag(1:nlevsoil_hyd)) ; this%supsub_flag = -999 - allocate(this%h2osoi_liqvol_shell(1:nlevsoil_hyd,1:nshell)) ; this%h2osoi_liqvol_shell = nan - allocate(this%h2osoi_liq_prev(1:nlevsoil_hyd)) ; this%h2osoi_liq_prev = nan - allocate(this%psisoi_liq_innershell(1:nlevsoil_hyd)); this%psisoi_liq_innershell = nan - allocate(this%rs1(1:nlevsoil_hyd)); this%rs1(:) = fine_root_radius_const - allocate(this%recruit_w_uptake(1:nlevsoil_hyd)); this%recruit_w_uptake = nan - - this%l_aroot_1D = nan - this%errh2o_hyd = nan - this%dwat_veg = nan - this%h2oveg = nan - this%h2oveg_dead = 0.0_r8 return end subroutine InitHydrSite From bc2ffeafb8fecaee7c24ee0bf0c67149b03b4470 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 23 May 2018 10:43:24 -0700 Subject: [PATCH 09/12] Updated descriptions about the calling context to root fractions. --- biogeochem/EDCanopyStructureMod.F90 | 4 ++++ biogeochem/EDPhysiologyMod.F90 | 5 ++++- biogeochem/FatesAllometryMod.F90 | 15 ++++++++++++--- 3 files changed, 20 insertions(+), 4 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 00b80830d5..be6d04eb70 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -930,6 +930,10 @@ subroutine canopy_summarization( nsites, sites, bc_in ) do while(associated(currentPatch)) ! Calculate rooting depth fractions for the patch x pft + ! Note that we are calling for the root fractions in the hydrologic context. + ! See explanation in FatesAllometryMod. In other locations, this + ! function is called to return the profile of biomass as used for litter + do ft = 1, numpft call set_root_fraction(currentPatch%rootfr_ft(ft,1:bc_in(s)%nlevsoil), ft, & bc_in(s)%zi_sisl,lowerb=lbound(bc_in(s)%zi_sisl,1), & diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 5353728ed0..50e2c9548e 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -2034,7 +2034,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) use FatesGlobals, only : endrun => fates_endrun use EDParamsMod , only : ED_val_cwd_flig, ED_val_cwd_fcel use FatesAllometryMod, only : set_root_fraction - use FatesAllometryMod, only : i_biomass_rootprof_context + use FatesAllometryMod, only : i_biomass_rootprof_context implicit none @@ -2129,6 +2129,9 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) do ft = 1, numpft ! This generates a rooting profile over the whole soil column for each pft + ! Note that we are calling for the root fractions in the biomass + ! for litter context, and not the hydrologic uptake context. + call set_root_fraction(cinput_rootfr(ft,1:bc_in(s)%nlevsoil), ft, & bc_in(s)%zi_sisl, lowerb=lbound(bc_in(s)%zi_sisl,1), & icontext=i_biomass_rootprof_context) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index ad6f91dbf7..229cb4966b 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -120,9 +120,18 @@ module FatesAllometryMod character(len=*), parameter :: sourcefile = __FILE__ - ! These are public contexts that other routines - ! should pass as arguments to the generic root profile - ! wrapper. + ! The code will call the wrapper routine "set_root_fraction" + ! in at least two different context. In one context it will query + ! set_root_fraction to describe the depth profile of hydraulicly + ! active roots. In the other context, it will ask the wrapper + ! to define the profile of roots as litter. We allow these + ! two contexts to differ. While not fully implemented, the use + ! will have control parameters to choose from different relationships + ! in these two contexts. The calling function, therefore + ! has to tell the wrapper function which context (water or biomass) + ! is being querried. So that we don't have to do messy string + ! parsing, we have two pre-defined flags. + integer, parameter, public :: i_hydro_rootprof_context = 1 integer, parameter, public :: i_biomass_rootprof_context = 2 From 3b83faa6550a606c93403b8497c4197eff50fd27 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 23 May 2018 11:20:05 -0700 Subject: [PATCH 10/12] Added descriptions of root profile contexts in the wrapper function. --- biogeochem/FatesAllometryMod.F90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 229cb4966b..c09274b985 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -1818,8 +1818,14 @@ subroutine set_root_fraction(root_fraction, ft, zi, lowerb, icontext ) ! PROFILE SWAPPING FLAGS. OR IF THERE IS NO DEMAND< LEAVE AS IS. ! ! - ! Two context exist 'hydraulic' and 'biomass' - ! These two contexts are allowed to have different profiles + ! Two context exist 'hydraulic' and 'biomass'. This allows us to + ! allow different profiles for how water is drawn from the soil + ! and different profiles to define the biomass for litter flux. + ! These two context can currently choose 1 of the following three + ! methods of defining the profile: 1) A 1 parameter exponential, 2) + ! a beta profile defined by Jackson et al. and 3) a 2 parameter + ! exponential. + ! All methods return a normalized profile. integer, parameter :: exponential_1p_profile_type = 1 integer, parameter :: jackson_beta_profile_type = 2 From f74a4167cbd0c9671dbd6f1fd25b230ffa7ff87d Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 23 May 2018 17:33:34 -0700 Subject: [PATCH 11/12] Potentially temporary change that tries to rectify setting ccohort_hydr%z_node_aroot in plant hydraulics with new soil depth system --- biogeophys/FatesPlantHydraulicsMod.F90 | 35 ++++++++++++++++---------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index a6a8f2a7fc..4f25c08daa 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -61,7 +61,6 @@ module FatesPlantHydraulicsMod use FatesHydraulicsMemMod, only: n_porous_media use FatesHydraulicsMemMod, only: nshell use FatesHydraulicsMemMod, only: n_hypool_ag - use FatesHydraulicsMemMod, only: n_hypool_troot use FatesHydraulicsMemMod, only: porous_media use FatesHydraulicsMemMod, only: cap_slp use FatesHydraulicsMemMod, only: cap_int @@ -311,12 +310,13 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) 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] 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 - + nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd + nlevsoil = bc_in%nlevsoil cCohort => cc_p ccohort_hydr => cc_p%co_hydr cPatch => cCohort%patchptr @@ -382,7 +382,8 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) ! 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/n_hypool_troot + dcumul_rf = 1._r8/dble(n_hypool_troot) + do k=1,n_hypool_troot cumul_rf = dcumul_rf*k call bisect_rootfr(roota, rootb, 0._r8, 1.E10_r8, & @@ -410,7 +411,12 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) ! ABSORBING ROOT DEPTH, LENGTH & VOLUME - ccohort_hydr%z_node_aroot(:) = -bc_in%z_sisl(:) + if ( nlevsoi_hyd == 1) then + ccohort_hydr%z_node_aroot(nlevsoi_hyd) = -bc_in%z_sisl(nlevsoi_hyd) + else + ccohort_hydr%z_node_aroot(1:nlevsoi_hyd) = -bc_in%z_sisl(1:nlevsoi_hyd) + end if + ccohort_hydr%l_aroot_tot = cCohort%br*C2B*EDPftvarcon_inst%hydr_srl(FT) !ccohort_hydr%v_aroot_tot = cCohort%br/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 @@ -433,7 +439,7 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) end do end if if(nlevsoi_hyd == 1) then - ccohort_hydr%z_node_troot(:) = ccohort_hydr%z_node_aroot(1) + ccohort_hydr%z_node_troot(:) = ccohort_hydr%z_node_aroot(nlevsoi_hyd) end if ! MAXIMUM (SIZE-DEPENDENT) HYDRAULIC CONDUCTANCES @@ -786,7 +792,7 @@ subroutine HydrSiteColdStart(sites, bc_in )! , bc_out) end do site_hydr%l_aroot_layer(1:site_hydr%nlevsoi_hyd) = 0.0_r8 - + end do ! @@ -929,10 +935,13 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) ! (kg water/m2 root area/Mpa/s) ! 1.e-5_r8 from Rudinger et al 1994 real(r8) :: kmax_root_surf_total !maximum conducitivity for total root surface(kg water/Mpa/s) - real(r8) :: kmax_soil_total !maximum conducitivity for total root surface(kg water/Mpa/s) + real(r8) :: kmax_soil_total !maximum conducitivity for total root surface(kg water/Mpa/s) + 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(:,:) @@ -954,14 +963,14 @@ subroutine updateSizeDepRhizHydProps(currentSite, bc_in ) csite_hydr%l_aroot_1D = sum( csite_hydr%l_aroot_layer(:)) ! update outer radii of column-level rhizosphere shells (same across patches and cohorts) - do j = 1,csite_hydr%nlevsoi_hyd + do j = 1,nlevsoi_hyd ! proceed only if l_aroot_coh has changed if( csite_hydr%l_aroot_layer(j) /= csite_hydr%l_aroot_layer_init(j) ) then call shellGeom( csite_hydr%l_aroot_layer(j), csite_hydr%rs1(j), AREA, bc_in%dz_sisl(j), & csite_hydr%r_out_shell(j,:), csite_hydr%r_node_shell(j,:),csite_hydr%v_shell(j,:)) end if !has l_aroot_layer changed? enddo - call shellGeom( csite_hydr%l_aroot_1D, csite_hydr%rs1(1), AREA, sum(bc_in%dz_sisl(1:csite_hydr%nlevsoi_hyd)), & + call shellGeom( csite_hydr%l_aroot_1D, csite_hydr%rs1(1), AREA, sum(bc_in%dz_sisl(1:nlevsoi_hyd)), & csite_hydr%r_out_shell_1D(:), csite_hydr%r_node_shell_1D(:), csite_hydr%v_shell_1D(:)) do j = 1,csite_hydr%nlevsoi_hyd @@ -1353,7 +1362,7 @@ subroutine FillDrainRhizShells(nsites, sites, bc_in, bc_out) if(csite_hydr%nlevsoi_hyd == 1) then dwat_kgm2 = bc_in(s)%h2o_liq_sisl(bc_in(s)%nlevsoil) - csite_hydr%h2osoi_liq_prev(csite_hydr%nlevsoi_hyd) - else if(csite_hydr%nlevsoi_hyd == bc_in(s)%nlevsoil ) then + else ! if(csite_hydr%nlevsoi_hyd == bc_in(s)%nlevsoil ) then dwat_kgm2 = bc_in(s)%h2o_liq_sisl(j) - csite_hydr%h2osoi_liq_prev(j) end if @@ -1421,7 +1430,7 @@ subroutine FillDrainRhizShells(nsites, sites, bc_in, bc_out) enddo ! balance check - if(csite_hydr%nlevsoi_hyd == bc_in(s)%nlevsoil) then + if(csite_hydr%nlevsoi_hyd .ne. 1) then do j = 1,csite_hydr%nlevsoi_hyd errh2o(j) = sum(h2osoi_liq_shell(j,:))/AREA - bc_in(s)%h2o_liq_sisl(j) @@ -1435,7 +1444,7 @@ subroutine FillDrainRhizShells(nsites, sites, bc_in, bc_out) end if end if enddo - else if(csite_hydr%nlevsoi_hyd == 1) then + else errh2o(csite_hydr%nlevsoi_hyd) = sum(h2osoi_liq_shell(csite_hydr%nlevsoi_hyd,:))/AREA - sum( bc_in(s)%h2o_liq_sisl(:) ) end if From 02792d0b2b9c102a704ecbf77b78e50a49aadda4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 24 May 2018 16:19:20 -0700 Subject: [PATCH 12/12] hydraulics: changed maximum hydraulics layers to a relatively large value, added code that checks to make sure it is larger than number of soil layers. --- biogeophys/FatesPlantHydraulicsMod.F90 | 19 +++++++++++-------- main/FatesHydraulicsMemMod.F90 | 9 ++++++++- 2 files changed, 19 insertions(+), 9 deletions(-) diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 4f25c08daa..98ff027ad6 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -410,13 +410,9 @@ subroutine updateSizeDepTreeHydProps(currentSite,cc_p,bc_in) 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) - if ( nlevsoi_hyd == 1) then - ccohort_hydr%z_node_aroot(nlevsoi_hyd) = -bc_in%z_sisl(nlevsoi_hyd) - else - ccohort_hydr%z_node_aroot(1:nlevsoi_hyd) = -bc_in%z_sisl(1:nlevsoi_hyd) - end if - + ccohort_hydr%l_aroot_tot = cCohort%br*C2B*EDPftvarcon_inst%hydr_srl(FT) !ccohort_hydr%v_aroot_tot = cCohort%br/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 @@ -731,8 +727,15 @@ subroutine InitHydrSites(sites,bc_in) do s=1,nsites allocate(csite_hydr) sites(s)%si_hydr => csite_hydr - sites(s)%si_hydr%nlevsoi_hyd = min(bc_in(s)%nlevsoil,nlevsoi_hyd_max) - + if ( bc_in(s)%nlevsoil > nlevsoi_hyd_max ) then + write(fates_log(),*) 'The host land model has defined soil with' + write(fates_log(),*) bc_in(s)%nlevsoil,' layers, for one of its columns.' + write(fates_log(),*) 'Fates-hydro temporary array spaces with size' + write(fates_log(),*) 'nlevsoi_hyd_max = ',nlevsoi_hyd_max,' must be larger' + write(fates_log(),*) 'see main/FatesHydraulicsMemMod.F90' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + sites(s)%si_hydr%nlevsoi_hyd = bc_in(s)%nlevsoil call sites(s)%si_hydr%InitHydrSite() end do diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90 index fa386fde10..5934aaa1af 100644 --- a/main/FatesHydraulicsMemMod.F90 +++ b/main/FatesHydraulicsMemMod.F90 @@ -9,7 +9,14 @@ module FatesHydraulicsMemMod implicit none ! Number of soil layers for indexing cohort fine root quanitities - integer, parameter :: nlevsoi_hyd_max = 10 + ! NOTE: The hydraulics code does have some capacity to run a single soil + ! layer that was developed for comparisons with TFS. However, this has + ! not maintained support through the integration with FATES and its + ! communications with the LSM. Please do not set nlevsoi_hyd_max + ! to 1 unless you are developing and testing. + + + integer, parameter :: nlevsoi_hyd_max = 40 ! number of distinct types of plant porous media (leaf, stem, troot, aroot) integer, parameter :: n_porous_media = 4