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]