Merge pull request ESCOMP#7 from ekluzek/master

Bring in changes from FATES sci 1.40 to 1.4.1

biogeochem/EDCohortDynamicsMod.F90

@@ -863,7 +863,6 @@ subroutine fuse_cohorts(patchptr, bc_in)
                                    currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn
                                    currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn
                                    currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn
-                                   currentCohort%imort = (currentCohort%n*currentCohort%imort + nextc%n*nextc%imort)/newn
                                    currentCohort%fmort = (currentCohort%n*currentCohort%fmort + nextc%n*nextc%fmort)/newn
 
                                    ! logging mortality, Yi Xu
@@ -1254,7 +1253,6 @@ subroutine copy_cohort( currentCohort,copyc )
     ! Mortality diagnostics
     n%cmort = o%cmort
     n%bmort = o%bmort
-    n%imort = o%imort
     n%fmort = o%fmort
     n%hmort = o%hmort
 

biogeochem/EDPatchDynamicsMod.F90

@@ -31,6 +31,10 @@ module EDPatchDynamicsMod
   use EDLoggingMortalityMod, only : logging_litter_fluxes 
   use EDLoggingMortalityMod, only : logging_time
   use EDParamsMod          , only : fates_mortality_disturbance_fraction
+  use FatesConstantsMod    , only : g_per_kg
+  use FatesConstantsMod    , only : ha_per_m2
+  use FatesConstantsMod    , only : days_per_sec
+  use FatesConstantsMod    , only : years_per_day
 
 
   ! CIME globals
@@ -117,7 +121,6 @@ subroutine disturbance_rates( site_in)
           currentCohort%cmort = cmort
           currentCohort%bmort = bmort
           currentCohort%hmort = hmort
-          currentCohort%imort = 0.0_r8 ! Impact mortality is always zero except in new patches
           currentCohort%fmort = 0.0_r8 ! Fire mortality is initialized as zero, but may be changed
 
           call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, &
@@ -199,7 +202,6 @@ subroutine disturbance_rates( site_in)
                 currentCohort%hmort = currentCohort%hmort*(1.0_r8 - fates_mortality_disturbance_fraction)
                 currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction)
                 currentCohort%dmort = currentCohort%dmort*(1.0_r8 - fates_mortality_disturbance_fraction)
-                ! currentCohort%imort will likely exist with logging
              end if
              currentCohort => currentCohort%taller
           enddo !currentCohort
@@ -399,7 +401,6 @@ subroutine spawn_patches( currentSite, bc_in)
                    nc%hmort = nan
                    nc%bmort = nan
                    nc%fmort = nan
-                   nc%imort = nan
                    nc%lmort_logging    = nan
                    nc%lmort_collateral = nan
                    nc%lmort_infra      = nan
@@ -414,7 +415,18 @@ subroutine spawn_patches( currentSite, bc_in)
                       !          The number density per square are doesn't change, but since the patch is smaller
                       !          and cohort counts are absolute, reduce this number.
                       nc%n = currentCohort%n * patch_site_areadis/currentPatch%area
-
+
+                      ! because the mortality rate due to impact for the cohorts which had been in the understory and are now in the newly-
+                      ! disturbed patch is very high, passing the imort directly to history results in large numerical errors, on account
+                      ! of the sharply reduced number densities.  so instead pass this info via a site-level diagnostic variable before reducing 
+                      ! the number density.
+                      currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = &
+                           currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + &
+                           nc%n * ED_val_understorey_death / hlm_freq_day
+                      currentSite%imort_carbonflux = currentSite%imort_carbonflux + &
+                           (nc%n * ED_val_understorey_death / hlm_freq_day ) * &
+                           currentCohort%b * g_per_kg * days_per_sec * years_per_day * ha_per_m2
+
                       ! Step 2:  Apply survivor ship function based on the understory death fraction
                       ! remaining of understory plants of those that are knocked over by the overstorey trees dying...  
                       nc%n = nc%n * (1.0_r8 - ED_val_understorey_death)
@@ -428,7 +440,6 @@ subroutine spawn_patches( currentSite, bc_in)
                       ! so with the number density must come the effective mortality rates.
 
                       nc%fmort            = 0.0_r8               ! Should had also been zero in the donor
-                      nc%imort            = ED_val_understorey_death/hlm_freq_day  ! This was zero in the donor
                       nc%cmort            = currentCohort%cmort
                       nc%hmort            = currentCohort%hmort
                       nc%bmort            = currentCohort%bmort
@@ -443,6 +454,7 @@ subroutine spawn_patches( currentSite, bc_in)
                       ! Besides, the current and newly created patch sum to unity                      
 
                       currentCohort%n = currentCohort%n * (1._r8 -  patch_site_areadis/currentPatch%area)
+
                    else 
                       ! grass is not killed by mortality disturbance events. Just move it into the new patch area. 
                       ! Just split the grass into the existing and new patch structures
@@ -452,7 +464,6 @@ subroutine spawn_patches( currentSite, bc_in)
                       currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area)
 
                       nc%fmort            = 0.0_r8
-                      nc%imort            = 0.0_r8
                       nc%cmort            = currentCohort%cmort
                       nc%hmort            = currentCohort%hmort
                       nc%bmort            = currentCohort%bmort
@@ -478,7 +489,6 @@ subroutine spawn_patches( currentSite, bc_in)
                 nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) 
 
                 nc%fmort            = currentCohort%fire_mort/hlm_freq_day
-                nc%imort            = 0.0_r8
 
                 nc%cmort            = currentCohort%cmort
                 nc%hmort            = currentCohort%hmort
@@ -508,7 +518,6 @@ subroutine spawn_patches( currentSite, bc_in)
                    nc%hmort            = nan
                    nc%bmort            = nan
                    nc%fmort            = nan
-                   nc%imort            = nan
                    nc%lmort_logging    = nan
                    nc%lmort_collateral = nan
                    nc%lmort_infra      = nan
@@ -526,6 +535,17 @@ subroutine spawn_patches( currentSite, bc_in)
                       !          The number density per square are doesn't change, but since the patch is smaller
                       !          and cohort counts are absolute, reduce this number.
                       nc%n = currentCohort%n * patch_site_areadis/currentPatch%area
+
+                      ! because the mortality rate due to impact for the cohorts which had been in the understory and are now in the newly-
+                      ! disturbed patch is very high, passing the imort directly to history results in large numerical errors, on account
+                      ! of the sharply reduced number densities.  so instead pass this info via a site-level diagnostic variable before reducing 
+                      ! the number density.
+                      currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = &
+                           currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + &
+                           nc%n * ED_val_understorey_death / hlm_freq_day
+                      currentSite%imort_carbonflux = currentSite%imort_carbonflux + &
+                           (nc%n * ED_val_understorey_death / hlm_freq_day ) * &
+                           currentCohort%b * g_per_kg * days_per_sec * years_per_day * ha_per_m2
 
                       ! Step 2:  Apply survivor ship function based on the understory death fraction
 
@@ -540,7 +560,6 @@ subroutine spawn_patches( currentSite, bc_in)
 
 
                       nc%fmort = 0.0_r8
-                      nc%imort = ED_val_understorey_death/hlm_freq_day
                       nc%cmort            = currentCohort%cmort
                       nc%hmort            = currentCohort%hmort
                       nc%bmort            = currentCohort%bmort
@@ -560,7 +579,6 @@ subroutine spawn_patches( currentSite, bc_in)
 
                       ! No grass impact mortality imposed on the newly created patch
                       nc%fmort            = 0.0_r8
-                      nc%imort            = 0.0_r8
                       nc%cmort            = currentCohort%cmort
                       nc%hmort            = currentCohort%hmort
                       nc%bmort            = currentCohort%bmort
@@ -992,7 +1010,7 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat
              !currentCohort%dmort = mortality_rates(currentCohort) 
              !the disturbance calculations are done with the previous n, c_area and d_mort. So it's probably &
              !not right to recalcualte dmort here.
-             canopy_dead = currentCohort%n * min(1.0_r8,currentCohort%dmort * hlm_freq_day)
+             canopy_dead = currentCohort%n * min(1.0_r8,currentCohort%dmort * hlm_freq_day * fates_mortality_disturbance_fraction)
 
              canopy_mortality_woody_litter   = canopy_mortality_woody_litter  + &
                   canopy_dead*(currentCohort%bdead+currentCohort%bsw)

biogeochem/EDPhysiologyMod.F90

@@ -845,9 +845,11 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in)
        if (currentCohort%canopy_layer .eq. 1) then
           currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_canopy(currentCohort%pft) &
                                        * currentCohort%c_area / currentCohort%n
+          currentCohort%npp_acc = currentCohort%npp_acc_hold / hlm_days_per_year ! add these for balance checking purposes
        else
           currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_understory(currentCohort%pft) &
                                        * currentCohort%c_area / currentCohort%n
+          currentCohort%npp_acc = currentCohort%npp_acc_hold / hlm_days_per_year ! add these for balance checking purposes
        endif
     endif
 
@@ -1081,22 +1083,27 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
        temp_cohort%bstore      = EDPftvarcon_inst%cushion(ft)*(temp_cohort%balive/ (1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) &
             + EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))
 
-       if (hlm_use_ed_prescribed_phys .eq. ifalse) then
+       if (hlm_use_ed_prescribed_phys .eq. ifalse .or. EDPftvarcon_inst%prescribed_recruitment(ft) .lt. 0. ) then
           temp_cohort%n           = currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day &
                / (temp_cohort%bdead+temp_cohort%balive+temp_cohort%bstore)
        else
           ! prescribed recruitment rates. number per sq. meter per year
           temp_cohort%n        = currentPatch%area * EDPftvarcon_inst%prescribed_recruitment(ft) * hlm_freq_day
+          ! modify the carbon balance accumulators to take into account the different way of defining recruitment
+          ! add prescribed rates as an input C flux, and the recruitment that would have otherwise occured as an output flux
+          ! (since the carbon associated with them effectively vanishes)
+          currentSite%flux_in = currentSite%flux_in + temp_cohort%n * (temp_cohort%bstore + temp_cohort%balive + temp_cohort%bdead)
+          currentSite%flux_out = currentSite%flux_out + currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day
        endif
 
        temp_cohort%laimemory = 0.0_r8     
        if (EDPftvarcon_inst%season_decid(temp_cohort%pft) == 1.and.currentSite%status == 1)then
-         temp_cohort%laimemory = (1.0_r8/(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) + &
-              EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))*temp_cohort%balive
+          temp_cohort%laimemory = (1.0_r8/(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) + &
+               EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))*temp_cohort%balive
        endif
        if (EDPftvarcon_inst%stress_decid(temp_cohort%pft) == 1.and.currentSite%dstatus == 1)then
-         temp_cohort%laimemory = (1.0_r8/(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) + &
-            EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))*temp_cohort%balive
+          temp_cohort%laimemory = (1.0_r8/(1.0_r8 + EDPftvarcon_inst%allom_l2fr(ft) + &
+               EDpftvarcon_inst%allom_latosa_int(ft)*temp_cohort%hite))*temp_cohort%balive
        endif
 
        cohortstatus = currentSite%status
@@ -1105,17 +1112,16 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
        endif
 
        if (temp_cohort%n > 0.0_r8 )then
-           if ( DEBUG ) write(fates_log(),*) 'EDPhysiologyMod.F90 call create_cohort '
-           call create_cohort(currentPatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, &
-                temp_cohort%balive, temp_cohort%bdead, temp_cohort%bstore,  &
-                temp_cohort%laimemory, cohortstatus, temp_cohort%canopy_trim, currentPatch%NCL_p, &
-                bc_in)
+          if ( DEBUG ) write(fates_log(),*) 'EDPhysiologyMod.F90 call create_cohort '
+          call create_cohort(currentPatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, &
+               temp_cohort%balive, temp_cohort%bdead, temp_cohort%bstore,  &
+               temp_cohort%laimemory, cohortstatus, temp_cohort%canopy_trim, currentPatch%NCL_p, &
+               bc_in)
 
-           ! keep track of how many individuals were recruited for passing to history
-           currentSite%recruitment_rate(ft) = currentSite%recruitment_rate(ft) + temp_cohort%n
+          ! keep track of how many individuals were recruited for passing to history
+          currentSite%recruitment_rate(ft) = currentSite%recruitment_rate(ft) + temp_cohort%n
 
        endif
-
     enddo  !pft loop
 
     deallocate(temp_cohort) ! delete temporary cohort
@@ -1176,7 +1182,7 @@ subroutine CWD_Input( currentSite, currentPatch)
               SF_val_CWD_frac(c) * currentCohort%n/currentPatch%area *(1.0_r8-EDPftvarcon_inst%allom_agb_frac(currentCohort%pft))
       enddo
 
-      if (currentCohort%canopy_layer > 1)then   
+      !if (currentCohort%canopy_layer > 1)then   
 
           ! ================================================        
           ! Litter fluxes for understorey  mortality. KgC/m2/year
@@ -1275,7 +1281,7 @@ subroutine CWD_Input( currentSite, currentPatch)
                 currentSite%resources_management%delta_individual + &
                 (dead_n_dlogging+dead_n_ilogging) * hlm_freq_day * currentPatch%area
 
-       endif !canopy layer
+       !endif !canopy layer
 
        currentCohort => currentCohort%taller
     enddo  ! end loop over cohorts 

main/EDInitMod.F90

@@ -64,6 +64,7 @@ subroutine init_site_vars( site_in )
     allocate(site_in%terminated_nindivs(1:nlevsclass,1:numpft,2))
     allocate(site_in%demotion_rate(1:nlevsclass))
     allocate(site_in%promotion_rate(1:nlevsclass))
+    allocate(site_in%imort_rate(1:nlevsclass,1:numpft))
     !
     end subroutine init_site_vars
 
@@ -115,6 +116,8 @@ subroutine zero_site( site_in )
     site_in%terminated_nindivs(:,:,:) = 0._r8
     site_in%termination_carbonflux(:) = 0._r8
     site_in%recruitment_rate(:) = 0._r8
+    site_in%imort_rate(:,:) = 0._r8
+    site_in%imort_carbonflux = 0._r8
 
     ! demotion/promotion info
     site_in%demotion_rate(:) = 0._r8

main/EDMainMod.F90

@@ -609,7 +609,6 @@ subroutine bypass_dynamics(currentSite)
           currentCohort%bmort = 0.0_r8
           currentCohort%hmort = 0.0_r8
           currentCohort%cmort = 0.0_r8
-          currentCohort%imort = 0.0_r8
           currentCohort%fmort = 0.0_r8
 
           currentCohort => currentCohort%taller

main/EDTypesMod.F90

@@ -226,7 +226,6 @@ module EDTypesMod
      real(r8) ::  bmort                                  ! background mortality rate        n/year
      real(r8) ::  cmort                                  ! carbon starvation mortality rate n/year
      real(r8) ::  hmort                                  ! hydraulic failure mortality rate n/year
-     real(r8) ::  imort                                  ! mortality from impacts by others n/year
      real(r8) ::  fmort                                  ! fire mortality                   n/year
 
       ! Logging Mortality Rate 
@@ -540,6 +539,8 @@ module EDTypesMod
      real(r8) :: demotion_carbonflux                           ! biomass of demoted individuals from canopy to understory [kgC/ha/day]
      real(r8), allocatable :: promotion_rate(:)               ! rate of individuals promoted from understory to canopy per FATES timestep
      real(r8) :: promotion_carbonflux                          ! biomass of promoted individuals from understory to canopy [kgC/ha/day]
+     real(r8), allocatable :: imort_rate(:,:)                    ! rate of individuals killed due to impact mortality per year.  on size x pft array
+     real(r8) :: imort_carbonflux                                ! biomass of individuals killed due to impact mortality per year. [kgC/ha/day]
 
      ! some diagnostic-only (i.e. not resolved by ODE solver) flux of carbon to CWD and litter pools from termination and canopy mortality
      real(r8) :: CWD_AG_diagnostic_input_carbonflux(1:ncwd)       ! diagnostic flux to AG CWD [kg C / m2 / yr]
@@ -743,7 +744,6 @@ subroutine dump_cohort(ccohort)
      write(fates_log(),*) 'co%woody_turnover         = ', ccohort%woody_turnover
      write(fates_log(),*) 'co%cmort                  = ', ccohort%cmort
      write(fates_log(),*) 'co%bmort                  = ', ccohort%bmort
-     write(fates_log(),*) 'co%imort                  = ', ccohort%imort
      write(fates_log(),*) 'co%fmort                  = ', ccohort%fmort
      write(fates_log(),*) 'co%hmort                  = ', ccohort%hmort
      write(fates_log(),*) 'co%isnew                  = ', ccohort%isnew