Merge pull request #1005 from YanlanLiu/cross_grid_seed

Cross grid seed dispersal

biogeochem/EDPhysiologyMod.F90

@@ -156,7 +156,7 @@ module EDPhysiologyMod
  public :: PreDisturbanceLitterFluxes
  public :: PreDisturbanceIntegrateLitter
  public :: GenerateDamageAndLitterFluxes
- public :: SeedIn
+ public :: SeedUpdate
  public :: UpdateRecruitL2FR
  public :: UpdateRecruitStoicH
  public :: SetRecruitL2FR
@@ -555,7 +555,6 @@ subroutine PreDisturbanceIntegrateLitter(currentPatch)
  litt%seed_germ_in(pft) - &
  litt%seed_germ_decay(pft)
 
-
  enddo
 
  ! Update the Coarse Woody Debris pools (above and below)
@@ -2022,8 +2021,8 @@ subroutine assign_cohort_SP_properties(currentCohort, htop, tlai, tsai, parea, i
  end subroutine assign_cohort_SP_properties
 
  ! =====================================================================================
-
- subroutine SeedIn( currentSite, bc_in )
+ 
+ subroutine SeedUpdate( currentSite )
 
  ! -----------------------------------------------------------------------------------
  ! Flux from plants into the seed pool.
@@ -2040,10 +2039,11 @@ subroutine SeedIn( currentSite, bc_in )
  ! !USES:
  use EDTypesMod, only : area
  use EDTypesMod, only : homogenize_seed_pfts
+ use FatesInterfaceTypesMod, only : hlm_seeddisp_cadence
+ use FatesInterfaceTypesMod, only : fates_dispersal_cadence_none
  !
  ! !ARGUMENTS
  type(ed_site_type), intent(inout), target :: currentSite
- type(bc_in_type), intent(in) :: bc_in
 
  type(fates_patch_type), pointer :: currentPatch
  type(litter_type), pointer :: litt
@@ -2052,26 +2052,32 @@ subroutine SeedIn( currentSite, bc_in )
 
  integer :: pft
  real(r8) :: store_m_to_repro ! mass sent from storage to reproduction upon death [kg/plant]
- real(r8) :: site_seed_rain(maxpft) ! This is the sum of seed-rain for the site [kg/site/day]
+ real(r8) :: site_seed_rain(numpft) ! This is the sum of seed-rain for the site [kg/site/day]
+ real(r8) :: site_disp_frac(numpft) ! Fraction of seeds from prodeced in current grid cell to
+ ! disperse out to other gridcells
  real(r8) :: seed_in_external ! Mass of externally generated seeds [kg/m2/day]
  real(r8) :: seed_stoich ! Mass ratio of nutrient per C12 in seeds [kg/kg]
  real(r8) :: seed_prod ! Seed produced in this dynamics step [kg/day]
  integer :: n_litt_types ! number of litter element types (c,n,p, etc)
  integer :: el ! loop counter for litter element types
  integer :: element_id ! element id consistent with parteh/PRTGenericMod.F90
- !------------------------------------------------------------------------------------
 
- do el = 1, num_elements
+ ! If the dispersal kernel is not turned on, keep the dispersal fraction at zero
+ site_disp_frac(:) = 0._r8
+ if (hlm_seeddisp_cadence .ne. fates_dispersal_cadence_none) then
+ site_disp_frac(:) = EDPftvarcon_inst%seed_dispersal_fraction(:)
+ end if
 
-  site_seed_rain(:) = 0._r8
+ el_loop: do el = 1, num_elements
 
+ site_seed_rain(:) = 0._r8
  element_id = element_list(el)
 
  site_mass => currentSite%mass_balance(el)
 
  ! Loop over all patches and sum up the seed input for each PFT
  currentPatch => currentSite%oldest_patch
- do while (associated(currentPatch))
+ seed_rain_loop: do while (associated(currentPatch))
 
  currentCohort => currentPatch%tallest
  do while (associated(currentCohort))
@@ -2100,14 +2106,15 @@ subroutine SeedIn( currentSite, bc_in )
  currentcohort%seed_prod = seed_prod
  end if
 
+
  site_seed_rain(pft) = site_seed_rain(pft) + &
- (seed_prod * currentCohort%n + store_m_to_repro)
+ (seed_prod * currentCohort%n + store_m_to_repro) ![kg/site/day, kg/ha/day]
 
  currentCohort => currentCohort%shorter
  enddo !cohort loop
 
  currentPatch => currentPatch%younger
- enddo
+ enddo seed_rain_loop
 
  ! We can choose to homogenize seeds. This is simple, we just
  ! add up all the seed from each pft at the site level, and then
@@ -2116,22 +2123,21 @@ subroutine SeedIn( currentSite, bc_in )
  site_seed_rain(1:numpft) = sum(site_seed_rain(:))/real(numpft,r8)
  end if
 
-
  ! Loop over all patches again and disperse the mixed seeds into the input flux
  ! arrays
-
  ! Loop over all patches and sum up the seed input for each PFT
  currentPatch => currentSite%oldest_patch
- do while (associated(currentPatch))
+ seed_in_loop: do while (associated(currentPatch))
 
  litt => currentPatch%litter(el)
  do pft = 1,numpft
 
  if(currentSite%use_this_pft(pft).eq.itrue)then
-
- ! Seed input from local sources (within site)
- litt%seed_in_local(pft) = litt%seed_in_local(pft) + site_seed_rain(pft)/area
-
+
+ ! Seed input from local sources (within site). Note that a fraction of the
+ ! internal seed rain is sent out to neighboring gridcells.
+ litt%seed_in_local(pft) = litt%seed_in_local(pft) + site_seed_rain(pft)*(1-site_disp_frac(pft))/area ![kg/m2/day]
+
  ! If we are using the Tree Recruitment Scheme (TRS) with or w/o seedling dynamics
  if ( any(regeneration_model == [TRS_regeneration, TRS_no_seedling_dyn]) .and. &
  prt_params%allom_dbh_maxheight(pft) > min_max_dbh_for_trees) then
@@ -2161,22 +2167,29 @@ subroutine SeedIn( currentSite, bc_in )
  end select
 
  ! Seed input from external sources (user param seed rain, or dispersal model)
- seed_in_external = seed_stoich*EDPftvarcon_inst%seed_suppl(pft)*years_per_day
+ ! Include both prescribed seed_suppl and seed_in dispersed from neighbouring gridcells
+ seed_in_external = seed_stoich*(currentSite%seed_in(pft)/area + EDPftvarcon_inst%seed_suppl(pft)*years_per_day) ![kg/m2/day]
  litt%seed_in_extern(pft) = litt%seed_in_extern(pft) + seed_in_external
 
  ! Seeds entering externally [kg/site/day]
  site_mass%seed_in = site_mass%seed_in + seed_in_external*currentPatch%area
  end if !use this pft 
  enddo
-
 
  currentPatch => currentPatch%younger
- enddo
+ enddo seed_in_loop
 
- end do
+ ! Determine the total site-level seed output for the current element and update the seed_out mass
+ ! for each element loop since the site_seed_rain is resent and updated for each element loop iteration
+ do pft = 1,numpft
+ site_mass%seed_out = site_mass%seed_out + site_seed_rain(pft)*site_disp_frac(pft) ![kg/site/day]
+ currentSite%seed_out(pft) = currentSite%seed_out(pft) + site_seed_rain(pft)*site_disp_frac(pft) ![kg/site/day]
+ end do
+
+ end do el_loop
 
  return
- end subroutine SeedIn
+ end subroutine SeedUpdate
 
  ! ============================================================================
 
@@ -2404,6 +2417,7 @@ subroutine SeedGermination( litt, cold_stat, drought_stat, bc_in, currentPatch )
 
  if ((prt_params%season_decid(pft) == itrue ) .and. &
  (any(cold_stat == [phen_cstat_nevercold,phen_cstat_iscold]))) then
+ ! no germination for all PFTs when cold
  litt%seed_germ_in(pft) = 0.0_r8
  endif
 
@@ -2422,7 +2436,6 @@ end subroutine SeedGermination
 
  ! =====================================================================================
 
-
  subroutine recruitment(currentSite, currentPatch, bc_in)
  !
  ! DESCRIPTION:

main/ChecksBalancesMod.F90

@@ -105,6 +105,7 @@ subroutine PatchMassStock(currentPatch,el,live_stock,seed_stock,litter_stock)
  seed_stock = currentPatch%area * &
  (sum(litt%seed) + sum(litt%seed_germ))
 
+
  ! Total mass on living plants
  live_stock = 0._r8
  currentCohort => currentPatch%tallest

main/EDInitMod.F90

@@ -208,12 +208,14 @@ subroutine init_site_vars( site_in, bc_in, bc_out )
 
  ! Initialize the static soil
  ! arrays from the boundary (initial) condition
-
  site_in%zi_soil(:) = bc_in%zi_sisl(:)
  site_in%dz_soil(:) = bc_in%dz_sisl(:)
  site_in%z_soil(:) = bc_in%z_sisl(:)
 
- !
+ ! Seed dispersal
+ allocate(site_in%seed_in(1:numpft))
+ allocate(site_in%seed_out(1:numpft))
+
  end subroutine init_site_vars
 
  ! ============================================================================
@@ -339,6 +341,10 @@ subroutine zero_site( site_in )
  ! canopy spread
  site_in%spread = 0._r8
 
+ ! Seed dispersal
+ site_in%seed_in(:) = 0.0_r8
+ site_in%seed_out(:) = 0.0_r8
+
  site_in%area_pft(:) = 0._r8
  site_in%use_this_pft(:) = fates_unset_int
  site_in%area_by_age(:) = 0._r8

main/EDMainMod.F90

@@ -45,7 +45,7 @@ module EDMainMod
  use EDPhysiologyMod , only : satellite_phenology
  use EDPhysiologyMod , only : recruitment
  use EDPhysiologyMod , only : trim_canopy
- use EDPhysiologyMod , only : SeedIn
+ use EDPhysiologyMod , only : SeedUpdate
  use EDPhysiologyMod , only : ZeroAllocationRates
  use EDPhysiologyMod , only : ZeroLitterFluxes
  use EDPhysiologyMod , only : PreDisturbanceLitterFluxes
@@ -246,6 +246,8 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out)
 
  ! adds small cohort of each PFT
  call recruitment(currentSite, currentPatch, bc_in)
+ !YL --------------
+ ! call recruitment(currentSite, currentPatch, bc_in, bc_out)
 
  currentPatch => currentPatch%younger
  enddo
@@ -717,8 +719,8 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out )
  ! With growth and mortality rates now calculated we can determine the seed rain
  ! fluxes. However, because this is potentially a cross-patch mixing model
  ! we will calculate this as a group
-
- call SeedIn(currentSite,bc_in)
+ 
+ call SeedUpdate(currentSite)
 
  ! Calculate all other litter fluxes
  ! -----------------------------------------------------------------------------------

main/EDPftvarcon.F90

@@ -110,7 +110,11 @@ module EDPftvarcon
  real(r8), allocatable :: germination_rate(:) ! Fraction of seed mass germinating per year (yr-1)
  real(r8), allocatable :: seed_decay_rate(:) ! Fraction of seed mass (both germinated and
  ! ungerminated), decaying per year (yr-1)
-
+ real(r8), allocatable :: seed_dispersal_pdf_scale(:) ! Seed dispersal scale parameter, Bullock et al. (2017)
+ real(r8), allocatable :: seed_dispersal_pdf_shape(:) ! Seed dispersal shape parameter, Bullock et al. (2017)
+ real(r8), allocatable :: seed_dispersal_max_dist(:) ! Maximum seed dispersal distance parameter (m)
+ real(r8), allocatable :: seed_dispersal_fraction(:) ! Fraction of seed rain to disperse, per pft
+
  real(r8), allocatable :: repro_frac_seed(:) ! fraciton of reproductive carbon that is seed
  real(r8), allocatable :: a_emerg(:) ! mean fraction of seed bank emerging [day-1]
  real(r8), allocatable :: b_emerg(:) ! seedling emergence sensitivity to soil moisture
@@ -676,7 +680,23 @@ subroutine Register_PFT(this, fates_params)
  name = 'fates_frag_seed_decay_rate'
  call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
  dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
+ name = 'fates_seed_dispersal_pdf_scale'
+ call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+ dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
+ name = 'fates_seed_dispersal_pdf_shape'
+ call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+ dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
+ name = 'fates_seed_dispersal_max_dist'
+ call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+ dimension_names=dim_names, lower_bounds=dim_lower_bound) 
 
+ name = 'fates_seed_dispersal_fraction'
+ call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+ dimension_names=dim_names, lower_bounds=dim_lower_bound) 
+
  name = 'fates_trim_limit'
  call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
  dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -1116,6 +1136,22 @@ subroutine Receive_PFT(this, fates_params)
  call fates_params%RetrieveParameterAllocate(name=name, &
  data=this%seed_decay_rate)
 
+ name = 'fates_seed_dispersal_pdf_scale'
+ call fates_params%RetrieveParameterAllocate(name=name, &
+ data=this%seed_dispersal_pdf_scale)
+
+ name = 'fates_seed_dispersal_pdf_shape'
+ call fates_params%RetrieveParameterAllocate(name=name, &
+ data=this%seed_dispersal_pdf_shape)
+
+ name = 'fates_seed_dispersal_max_dist'
+ call fates_params%RetrieveParameterAllocate(name=name, &
+ data=this%seed_dispersal_max_dist)
+
+ name = 'fates_seed_dispersal_fraction'
+ call fates_params%RetrieveParameterAllocate(name=name, &
+ data=this%seed_dispersal_fraction) 
+
  name = 'fates_trim_limit'
  call fates_params%RetrieveParameterAllocate(name=name, &
  data=this%trim_limit)
@@ -1651,6 +1687,10 @@ subroutine FatesReportPFTParams(is_master)
  write(fates_log(),fmt0) 'jmaxse = ',EDPftvarcon_inst%jmaxse
  write(fates_log(),fmt0) 'germination_timescale = ',EDPftvarcon_inst%germination_rate
  write(fates_log(),fmt0) 'seed_decay_turnover = ',EDPftvarcon_inst%seed_decay_rate
+ write(fates_log(),fmt0) 'seed_dispersal_pdf_scale = ',EDPftvarcon_inst%seed_dispersal_pdf_scale
+ write(fates_log(),fmt0) 'seed_dispersal_pdf_shape = ',EDPftvarcon_inst%seed_dispersal_pdf_shape
+ write(fates_log(),fmt0) 'seed_dispersal_max_dist = ',EDPftvarcon_inst%seed_dispersal_max_dist
+ write(fates_log(),fmt0) 'seed_dispersal_fraction = ',EDPftvarcon_inst%seed_dispersal_fraction 
  write(fates_log(),fmt0) 'repro_frac_seed = ',EDPftvarcon_inst%repro_frac_seed
  write(fates_log(),fmt0) 'a_emerg = ',EDPftvarcon_inst%a_emerg
  write(fates_log(),fmt0) 'b_emerg = ',EDPftvarcon_inst%b_emerg
@@ -1883,6 +1923,68 @@ subroutine FatesCheckParams(is_master)
  call endrun(msg=errMsg(sourcefile, __LINE__))
  end if 
 
+ ! Check that the seed dispersal parameters are all set if one of them is set
+ !-----------------------------------------------------------------------------------
+ if (( EDPftvarcon_inst%seed_dispersal_pdf_scale(ipft) < fates_check_param_set ) .and. &
+ (( EDPftvarcon_inst%seed_dispersal_max_dist(ipft) > fates_check_param_set ) .or. &
+ ( EDPftvarcon_inst%seed_dispersal_pdf_shape(ipft) > fates_check_param_set ) .or. &
+ ( EDPftvarcon_inst%seed_dispersal_fraction(ipft) > fates_check_param_set ))) then
+
+ write(fates_log(),*) 'Seed dispersal is on per fates_seed_dispersal_pdf_scale being set'
+ write(fates_log(),*) 'Please provide values for all other seed_dispersal parameters'
+ write(fates_log(),*) 'See Bullock et al. (2017) for reasonable values'
+ write(fates_log(),*) 'Aborting'
+ call endrun(msg=errMsg(sourcefile, __LINE__))
+ end if
+
+ if (( EDPftvarcon_inst%seed_dispersal_pdf_shape(ipft) < fates_check_param_set ) .and. &
+ (( EDPftvarcon_inst%seed_dispersal_max_dist(ipft) > fates_check_param_set ) .or. &
+ ( EDPftvarcon_inst%seed_dispersal_pdf_scale(ipft) > fates_check_param_set ) .or. &
+ ( EDPftvarcon_inst%seed_dispersal_fraction(ipft) > fates_check_param_set ))) then
+
+ write(fates_log(),*) 'Seed dispersal is on per fates_seed_dispersal_pdf_shape being set'
+ write(fates_log(),*) 'Please provide values for all other seed_dispersal parameters'
+ write(fates_log(),*) 'See Bullock et al. (2017) for reasonable values'
+ write(fates_log(),*) 'Aborting'
+ call endrun(msg=errMsg(sourcefile, __LINE__))
+ end if
+
+ if (( EDPftvarcon_inst%seed_dispersal_max_dist(ipft) < fates_check_param_set ) .and. &
+ (( EDPftvarcon_inst%seed_dispersal_pdf_shape(ipft) > fates_check_param_set ) .or. &
+ ( EDPftvarcon_inst%seed_dispersal_pdf_scale(ipft) > fates_check_param_set ) .or. &
+ ( EDPftvarcon_inst%seed_dispersal_fraction(ipft) > fates_check_param_set ))) then
+
+ write(fates_log(),*) 'Seed dispersal is on per seed_dispersal_max_dist being set'
+ write(fates_log(),*) 'Please provide values for all other seed_dispersal parameters'
+ write(fates_log(),*) 'See Bullock et al. (2017) for reasonable values'
+ write(fates_log(),*) 'Aborting'
+ call endrun(msg=errMsg(sourcefile, __LINE__))
+ end if
+
+ if (( EDPftvarcon_inst%seed_dispersal_fraction(ipft) < fates_check_param_set ) .and. &
+ (( EDPftvarcon_inst%seed_dispersal_pdf_shape(ipft) > fates_check_param_set ) .or. &
+ ( EDPftvarcon_inst%seed_dispersal_pdf_scale(ipft) > fates_check_param_set ) .or. &
+ ( EDPftvarcon_inst%seed_dispersal_max_dist(ipft) > fates_check_param_set ))) then
+
+ write(fates_log(),*) 'Seed dispersal is on per seed_dispersal_fraction being set'
+ write(fates_log(),*) 'Please provide values for all other seed_dispersal parameters'
+ write(fates_log(),*) 'See Bullock et al. (2017) for reasonable values'
+ write(fates_log(),*) 'Aborting'
+ call endrun(msg=errMsg(sourcefile, __LINE__))
+ end if
+
+ ! Check that parameter ranges for the seed dispersal fraction make sense
+ !-----------------------------------------------------------------------------------
+ if (( EDPftvarcon_inst%seed_dispersal_fraction(ipft) < fates_check_param_set ) .and. &
+ (( EDPftvarcon_inst%seed_dispersal_fraction(ipft) > 1.0_r8 ) .or. &
+ ( EDPftvarcon_inst%seed_dispersal_fraction(ipft) < 0.0_r8 ))) then
+
+ write(fates_log(),*) 'Seed dispersal is on per seed_dispersal_fraction being set'
+ write(fates_log(),*) 'Please make sure the fraction value is between 0 and 1'
+ write(fates_log(),*) 'Aborting'
+ call endrun(msg=errMsg(sourcefile, __LINE__))
+ end if
+
  ! Check that parameter ranges for age-dependent mortality make sense
  !-----------------------------------------------------------------------------------
  if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < fates_check_param_set ) .and. &

main/EDTypesMod.F90

@@ -420,6 +420,10 @@ module EDTypesMod
  ! Canopy Spread
  real(r8) :: spread ! dynamic canopy allometric term [unitless]
 
+ ! Seed dispersal
+ real(r8), allocatable :: seed_out(:) ! amount of seed leaving the site [kg/site/day]
+ real(r8), allocatable :: seed_in(:) ! amount of seed dispersed into the site from neighbouring cells [kg/site/day]
+
  ! site-level variables to keep track of the disturbance rates, both actual and "potential"
  real(r8) :: disturbance_rates_primary_to_primary(N_DIST_TYPES) ! actual disturbance rates from primary patches to primary patches [m2/m2/day]
  real(r8) :: disturbance_rates_primary_to_secondary(N_DIST_TYPES) ! actual disturbance rates from primary patches to secondary patches [m2/m2/day]