Fix bugs in the Tree Recruitment Scheme with seedling dynamics

biogeochem/EDPhysiologyMod.F90

@@ -541,6 +541,7 @@ subroutine PreDisturbanceIntegrateLitter(currentPatch)
     integer :: nlevsoil    ! number of soil layers
     integer :: ilyr        ! soil layer loop counter
     integer :: dcmpy       ! decomposability index
+    real(r8) :: net_seed_available  ! Available seed after inputs and decay
 
     do el = 1, num_elements
 
@@ -550,15 +551,20 @@ subroutine PreDisturbanceIntegrateLitter(currentPatch)
        ! -----------------------------------------------------------------------------------
 
        do pft = 1,numpft
-          litt%seed(pft) = litt%seed(pft) + &
-               litt%seed_in_local(pft) +   &
-               litt%seed_in_extern(pft) -  &
-               litt%seed_decay(pft) -      &
-               litt%seed_germ_in(pft)
-
-          ! Note that the recruitment scheme will use seed_germ
-          ! for its construction costs.
-          litt%seed_germ(pft) = litt%seed_germ(pft) + &
+
+         ! Calculate net available seed after inputs and decay
+         net_seed_available = litt%seed(pft) + litt%seed_in_local(pft) + &
+            litt%seed_in_extern(pft) - litt%seed_decay(pft)
+
+         ! Cap germination at available seed to prevent negative seed pool
+         litt%seed_germ_in(pft) = min(litt%seed_germ_in(pft), net_seed_available)
+
+         ! Update pools 
+         litt%seed(pft) = net_seed_available - litt%seed_germ_in(pft)
+
+         ! Note that the recruitment scheme will use seed_germ
+         ! for its construction costs.
+         litt%seed_germ(pft) = litt%seed_germ(pft) + &
                litt%seed_germ_in(pft) - &
                litt%seed_germ_decay(pft)
 
@@ -2253,6 +2259,7 @@ subroutine SeedDecay( litt , currentPatch )
     real(r8) ::  seedling_light_mort_rate    ! daily seedling mortality rate from light stress
     real(r8) ::  seedling_h2o_mort_rate      ! daily seedling mortality rate from moisture stress
     real(r8) ::  seedling_mdds               ! moisture deficit days accumulated in the seedling layer
+    real(r8) ::  total_seedling_mort_rate    ! summed seedling mortality rates
 
     !----------------------------------------------------------------------
 
@@ -2317,25 +2324,42 @@ subroutine SeedDecay( litt , currentPatch )
 
           ! Get the current seedling moisture deficit days (tracked as a pft-specific exponential
           ! average)
-          seedling_mdds = currentPatch%sdlng_mdd(pft)%p%GetMean()     
+          seedling_mdds = currentPatch%sdlng_mdd(pft)%p%GetMean()
 
           ! Calculate seedling mortality as a function of moisture deficit days (mdd)
           ! If the seedling mmd value is below a critical threshold then moisture-based mortality is zero
+
           if (seedling_mdds < EDPftvarcon_inst%seedling_mdd_crit(pft)) then
              seedling_h2o_mort_rate = 0.0_r8
           else
              seedling_h2o_mort_rate = EDPftvarcon_inst%seedling_h2o_mort_a(pft) * seedling_mdds**2 + &
                   EDPftvarcon_inst%seedling_h2o_mort_b(pft) * seedling_mdds + &
                   EDPftvarcon_inst%seedling_h2o_mort_c(pft)
+             ! Cap h2o mortality rate at 1, quadratic can produce values >1 for large moisture
+             ! deficit days, if rate exceeds 1 write to fates log 
+
+             if (seedling_h2o_mort_rate > 1.0_r8) then
+               write(fates_log(), *) 'TRS seedling_h2o_mort_rate > 1', &
+                     'pft=', pft, 'uncapped=', seedling_h2o_mort_rate, 'mdds=', seedling_mdds
+             end if 
+             seedling_h2o_mort_rate = min(1.0_r8, seedling_h2o_mort_rate)
           end if ! mdd threshold check
 
           ! Step 3. Sum modes of mortality (including background mortality) and send dead seedlings
           ! to litter        
-          litt%seed_germ_decay(pft) = (litt%seed_germ(pft) * seedling_light_mort_rate) + &
-               (litt%seed_germ(pft) * seedling_h2o_mort_rate) + &
-               (litt%seed_germ(pft) * EDPftvarcon_inst%background_seedling_mort(pft) &
-               * years_per_day)
-
+          ! Cap total mortality rate at 1 to prevent negative seedling pool
+          total_seedling_mort_rate = seedling_light_mort_rate + &
+                        seedling_h2o_mort_rate + &
+                        (EDPftvarcon_inst%background_seedling_mort(pft) * years_per_day)
+          if (total_seedling_mort_rate > 1.0_r8) then 
+            write(fates_log(), *) 'TRS total seedling mortality rate > 1', &
+                  'pft=', pft, 'uncapped=', total_seedling_mort_rate, &
+                  'seedling_light_mort_rate=', seedling_light_mort_rate, & 
+                  'seedling_h2o_mort_rate=', seedling_h2o_mort_rate
+          end if
+          litt%seed_germ_decay(pft) = litt%seed_germ(pft) * &
+            min(1.0_r8, total_seedling_mort_rate)
+
        else
 
           litt%seed_germ_decay(pft) = litt%seed_germ(pft) * &
@@ -2440,6 +2464,15 @@ subroutine SeedGermination( litt, cold_stat, drought_stat, currentPatch )
           if ( seedling_layer_smp .GE. EDPftvarcon_inst%seedling_psi_emerg(pft) ) then
              seedling_emerg_rate = photoblastic_germ_modifier * EDPftvarcon_inst%a_emerg(pft) * &
                   wetness_index**EDPftvarcon_inst%b_emerg(pft)
+             ! Cap emergence rate at 1, rate can exceed 1 for some parameter combinations,
+             ! leading to negative seed bank, write to log if this is the case 
+             if (seedling_emerg_rate > 1.0_r8) then
+               write(fates_log(),*) 'TRS seedling_emerg_rate > 1', &
+                     'pft=', pft, 'uncapped=', seedling_emerg_rate, &
+                     'photoblastic_germ_modifier=', photoblastic_germ_modifier, &
+                     'wetness_index=', wetness_index
+             end if 
+             seedling_emerg_rate = min(1.0_r8, seedling_emerg_rate )
           else 
 
              seedling_emerg_rate = 0.0_r8
@@ -2522,6 +2555,7 @@ subroutine recruitment(currentSite, currentPatch, bc_in)
       real(r8)                          :: stem_drop_fraction ! 
       real(r8)                          :: fnrt_drop_fraction ! 
       real(r8)                          :: sdlng2sap_par      ! running mean of PAR at the seedling layer [MJ/m2/day]
+      real(r8)                          :: sdlng2sap_rate     ! seedling-to-sapling transition rate factor (unitless)
       real(r8)                          :: seedling_layer_smp ! soil matric potential at seedling rooting depth [mm H2O suction]
       integer, parameter                :: recruitstatus = 1  ! whether the newly created cohorts are recruited or initialized
       integer                           :: ilayer_seedling_root ! the soil layer at seedling rooting depth
@@ -2667,11 +2701,22 @@ subroutine recruitment(currentSite, currentPatch, bc_in)
                      sdlng2sap_par = currentPatch%sdlng2sap_par%GetMean()*     &
                         sec_per_day*megajoules_per_joule
 
+                     sdlng2sap_rate = EDPftvarcon_inst%seedling_light_rec_a(ft)*             &
+                        sdlng2sap_par**EDPftvarcon_inst%seedling_light_rec_b(ft)
+
+                     ! If the seedling to sapling transition rate exceeds 1, 
+                     ! cap at 1 (prevent seed_germ from going negative) and write to fates log 
+
+                     if (sdlng2sap_rate > 1.0_r8) then
+                        write(fates_log(), *) 'TRS seedling-to-sapling transition rate > 1', &
+                              'pft=', ft, 'uncapped=', sdlng2sap_rate,         &
+                              'sdlng2sap_par=', sdlng2sap_par
+                     end if 
+
                      mass_avail = currentPatch%area*                           &
                         currentPatch%litter(el)%seed_germ(ft)*                 & 
-                        EDPftvarcon_inst%seedling_light_rec_a(ft)*             &
-                        sdlng2sap_par**EDPftvarcon_inst%seedling_light_rec_b(ft) 
-
+                        min(1.0_r8, sdlng2sap_rate)
+
                      ! If soil moisture is below pft-specific seedling  moisture stress threshold the 
                      ! recruitment does not occur.
                      ilayer_seedling_root = minloc(abs(bc_in%z_sisl(:) -       &
@@ -2780,6 +2825,7 @@ subroutine recruitment(currentSite, currentPatch, bc_in)
                      currentPatch%litter(el)%seed_germ(ft) =                   &
                      currentPatch%litter(el)%seed_germ(ft) - cohort_n / currentPatch%area *   &
                      (m_struct + m_leaf + m_fnrt + m_sapw + m_store + m_repro)
+
                   end if
 
                end do

biogeochem/FatesPatchMod.F90

@@ -666,7 +666,7 @@ subroutine InitRunningMeans(this, current_tod, regeneration_model, numpft)
         call this%seedling_layer_par24%InitRMean(fixed_24hr,                   &
           init_value=init_seedling_par, init_offset=real(current_tod, r8))
         call this%sdlng_mort_par%InitRMean(ema_sdlng_mort_par,                 &
-          init_value=temp_init_veg)
+          init_value=init_seedling_par)
         call this%sdlng2sap_par%InitRMean(ema_sdlng2sap_par,                   &
           init_value=init_seedling_par)
 

main/FatesInterfaceMod.F90

@@ -2335,8 +2335,12 @@ subroutine UpdateFatesRMeansTStep(sites,bc_in, bc_out)
                  endif
 
                  ! Update the seedling layer smp and mdd running means
-                 call cpatch%sdlng_emerg_smp(pft)%p%UpdateRMean(new_seedling_layer_smp)
-                 call cpatch%sdlng_mdd(pft)%p%UpdateRMean(new_seedling_mdd)
+                 ! Only update after first model day to avoid adding
+                 ! initialization spike in smp to running means
+                 if (hlm_model_day > 2.0_r8) then 
+                   call cpatch%sdlng_mdd(pft)%p%UpdateRMean(new_seedling_mdd)
+                   call cpatch%sdlng_emerg_smp(pft)%p%UpdateRMean(new_seedling_layer_smp)
+                 endif
 
               enddo !end pft loop
 
@@ -2410,8 +2414,8 @@ subroutine SeedlingParPatch(cpatch, &
 
   ! Start with the assumption that there is a single canopy layer
   seedling_par_high = atm_par_dir+atm_par_dif
-  par_high_frac     = 1._r8-cpatch%total_canopy_area
-  par_low_frac      = cpatch%total_canopy_area
+  par_high_frac     = 1._r8 - (cpatch%total_canopy_area / cpatch%area)
+  par_low_frac      = cpatch%total_canopy_area / cpatch%area
 
   ! Work up through the canopy layers from the bottom layer
   do cl = cpatch%NCL_p,max(1,cpatch%NCL_p-1),-1