Merge pull request #3 from rgknox/SF_fuels_refactor

fix CWD issue, indexing and parameter file
rgknox · Apr 29, 2019 · 12a6c7e · 12a6c7e
2 parents 7e6d3dd + 1015c4b
commit 12a6c7e
Show file tree

Hide file tree

Showing 6 changed files with 57 additions and 29 deletions.
diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90
@@ -59,6 +59,9 @@ module EDLoggingMortalityMod
 
    character(len=*), parameter, private :: sourcefile = &
          __FILE__
+
+
+   real(r8), public, parameter :: logging_export_frac = 0.8_r8
 
    public :: LoggingMortality_frac
    public :: logging_litter_fluxes
@@ -170,6 +173,8 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, &
       real(r8), parameter   :: adjustment = 1.0 ! adjustment for mortality rates
 
       if (logging_time) then 
+
+
          if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees
 
             ! Pass logging rates to cohort level 
@@ -194,6 +199,9 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, &
             ! will be applied via "understory_death" via the disturbance algorithm
             if (canopy_layer .eq. 1) then
                lmort_collateral = logging_collateral_frac * adjustment
+            else
+               lmort_collateral = 0._r8
+               l_degrad         = l_degrad + logging_collateral_frac * adjustment
             endif
 
          else
@@ -202,6 +210,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, &
             lmort_infra      = 0.0_r8
             l_degrad         = 0.0_r8
          end if
+
       else 
          lmort_direct    = 0.0_r8
          lmort_collateral = 0.0_r8

diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90
@@ -200,15 +200,18 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in)
                                currentCohort%lmort_infra,                        &
                                currentCohort%l_degrad)
 
+
+
+
     if (currentCohort%canopy_layer > 1)then 
-
        ! Include understory logging mortality rates not associated with disturbance
        dndt_logging = (currentCohort%lmort_direct     + &
                        currentCohort%lmort_collateral + &
                        currentCohort%lmort_infra)/hlm_freq_day
-
        currentCohort%dndt = -1.0_r8 * (cmort+hmort+bmort+frmort+dndt_logging) * currentCohort%n
     else
+       ! Mortality from logging in the canopy is ONLY disturbance generating, don't
+       ! update number densities via non-disturbance inducing death
        currentCohort%dndt = -(1.0_r8 - fates_mortality_disturbance_fraction) &
             * (cmort+hmort+bmort+frmort) * currentCohort%n
     endif

diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90
@@ -20,6 +20,7 @@ module EDPatchDynamicsMod
   use EDTypesMod           , only : dtype_ilog
   use EDTypesMod           , only : dtype_ifire
   use EDTypesMod           , only : ican_upper
+  use EDTypesMod           , only : lg_sf
   use FatesInterfaceMod    , only : hlm_use_planthydro
   use FatesInterfaceMod    , only : hlm_numSWb
   use FatesInterfaceMod    , only : bc_in_type
@@ -477,6 +478,8 @@ subroutine spawn_patches( currentSite, bc_in)
 
                 call logging_litter_fluxes(currentSite, currentPatch, new_patch, patch_site_areadis)
 
+                if(debug) write(fates_log(),*) "Logging disturbance generated:",patch_site_areadis
+
              elseif ((currentPatch%disturbance_rates(dtype_ifire) > &
                       currentPatch%disturbance_rates(dtype_ifall)) .and. &
                      (currentPatch%disturbance_rates(dtype_ifire) > &
@@ -1085,7 +1088,7 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si
        !************************************/ 
        do c = 1,ncwd
           burned_litter = new_patch%cwd_ag(c) * patch_site_areadis/new_patch%area * &
-                currentPatch%burnt_frac_litter(c+1) !kG/m2/day
+                currentPatch%burnt_frac_litter(c) !kG/m2/day
           new_patch%cwd_ag(c) = new_patch%cwd_ag(c) - burned_litter
           currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/day
           currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + &
@@ -1240,7 +1243,7 @@ subroutine fire_litter_fluxes(currentSite, cp_target, new_patch_target, patch_si
           if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then
              burned_leaves = leaf_c * currentCohort%fraction_crown_burned
           else
-             burned_leaves = leaf_c * currentPatch%burnt_frac_litter(6)
+             burned_leaves = leaf_c * currentPatch%burnt_frac_litter(lg_sf)
           endif
 
           if (burned_leaves > 0.0_r8) then

diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90
@@ -42,7 +42,8 @@ module EDPhysiologyMod
   use FatesGlobals          , only : fates_log
   use FatesGlobals          , only : endrun => fates_endrun
   use EDParamsMod           , only : fates_mortality_disturbance_fraction
-  use EDParamsMod           , only : logging_export_frac
+  use EDLoggingMortalityMod , only : logging_export_frac
+  !use EDParamsMod           , only : logging_export_frac
 
   use FatesPlantHydraulicsMod  , only : AccumulateMortalityWaterStorage
 
@@ -1167,26 +1168,38 @@ subroutine CWD_Input( currentSite, currentPatch)
               (1.0_r8-EDPftvarcon_inst%allom_agb_frac(currentCohort%pft))
       enddo
 
-      !if (currentCohort%canopy_layer > 1)then   
-
           ! ================================================        
           ! Litter fluxes for understorey  mortality. KgC/m2/year
           ! ================================================
 
           ! Total number of dead understory (n/m2)
           dead_n = -1.0_r8 * currentCohort%dndt / currentPatch%area
 
-          ! Total number of dead understory from direct logging
-          ! (it is possible that large harvestable trees are in the understory)
-          dead_n_dlogging = ( currentCohort%lmort_direct) * &
-               currentCohort%n/hlm_freq_day/currentPatch%area
-
-          ! Total number of dead understory from indirect logging
-          dead_n_ilogging = ( currentCohort%lmort_collateral + currentCohort%lmort_infra) * &
-                currentCohort%n/hlm_freq_day/currentPatch%area
-
-          dead_n_natural = dead_n - dead_n_dlogging - dead_n_ilogging
+          if(currentCohort%canopy_layer > 1)then   
+
+             ! Total number of dead understory from direct logging
+             ! (it is possible that large harvestable trees are in the understory)
+             dead_n_dlogging = ( currentCohort%lmort_direct) * &
+                  currentCohort%n/hlm_freq_day/currentPatch%area
+
+             ! Total number of dead understory from indirect logging
+             dead_n_ilogging = ( currentCohort%lmort_collateral + currentCohort%lmort_infra) * &
+                  currentCohort%n/hlm_freq_day/currentPatch%area
+
+          else
+
+             ! All mortality from logging in the canopy is
+             ! is disturbance generating
+
+             dead_n_dlogging = 0._r8
+             dead_n_ilogging = 0._r8
+
+          end if
+
+
 
+          dead_n_natural = dead_n - dead_n_dlogging - dead_n_ilogging
+
 
           currentPatch%leaf_litter_in(pft) = currentPatch%leaf_litter_in(pft) + &
                (leaf_c)* dead_n
@@ -1274,7 +1287,8 @@ subroutine CWD_Input( currentSite, currentPatch)
 
              if (currentPatch%cwd_AG_in(c) < 0.0_r8)then
                 write(fates_log(),*) 'negative CWD in flux',currentPatch%cwd_AG_in(c), &
-                      (struct_c + sapw_c), dead_n
+                      (struct_c + sapw_c), dead_n, dead_n_natural, dead_n_ilogging, dead_n_dlogging
+                call endrun(msg=errMsg(sourcefile, __LINE__))
              endif
 
           end do
@@ -1283,7 +1297,6 @@ subroutine CWD_Input( currentSite, currentPatch)
                 currentSite%resources_management%delta_individual + &
                 (dead_n_dlogging+dead_n_ilogging) * hlm_freq_day * currentPatch%area
 
-       !endif !canopy layer
 
        currentCohort => currentCohort%taller
     enddo  ! end loop over cohorts 

diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90
@@ -750,8 +750,8 @@ subroutine bypass_dynamics(currentSite)
           currentCohort%frmort = 0.0_r8
 
           currentCohort%dndt      = 0.0_r8
-	  currentCohort%dhdt      = 0.0_r8
-	  currentCohort%ddbhdt    = 0.0_r8
+          currentCohort%dhdt      = 0.0_r8
+          currentCohort%ddbhdt    = 0.0_r8
 
           currentCohort => currentCohort%taller
        enddo

diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90
@@ -118,8 +118,8 @@ module EDParamsMod
    real(r8),protected :: logging_dbhmax_infra        ! "Tree diameter, above which infrastructure from logging does not impact damage or mortality.
    character(len=param_string_length),parameter :: logging_name_dbhmax_infra = "fates_logging_dbhmax_infra"
 
-   real(r8),protected :: logging_export_frac        ! "fraction of trunk product being shipped offsite, the leftovers will be left onsite as large CWD
-   character(len=param_string_length),parameter :: logging_name_export_frac ="fates_logging_export_frac"
+!   real(r8),protected :: logging_export_frac        ! "fraction of trunk product being shipped offsite, the leftovers will be left onsite as large CWD
+!   character(len=param_string_length),parameter :: logging_name_export_frac ="fates_logging_export_frac"
 
    public :: FatesParamsInit
    public :: FatesRegisterParams
@@ -172,7 +172,7 @@ subroutine FatesParamsInit()
     logging_mechanical_frac               = nan
     logging_event_code                    = nan
     logging_dbhmax_infra                  = nan
-    logging_export_frac                   = nan
+!    logging_export_frac                   = nan
 
   end subroutine FatesParamsInit
 
@@ -296,8 +296,8 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=logging_name_dbhmax_infra, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
-    call fates_params%RegisterParameter(name=logging_export_frac, dimension_shape=dimension_shape_1d, &
-         dimension_names=dim_names)
+!    call fates_params%RegisterParameter(name=logging_export_frac, dimension_shape=dimension_shape_1d, &
+!         dimension_names=dim_names)
 
     ! non-scalar parameters
     call fates_params%RegisterParameter(name=ED_name_history_sizeclass_bin_edges, dimension_shape=dimension_shape_1d, &
@@ -420,8 +420,8 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetreiveParameter(name=logging_name_dbhmax_infra, &
           data=logging_dbhmax_infra)
 
-    call fates_params%RetreiveParameter(name=logging_export_frac, &
-          data=logging_export_frac)
+!    call fates_params%RetreiveParameter(name=logging_name_export_frac, &
+!          data=logging_export_frac)
 
     ! parameters that are arrays of size defined within the params file and thus need allocating as well
     call fates_params%RetreiveParameterAllocate(name=ED_name_history_sizeclass_bin_edges, &
@@ -481,7 +481,7 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),fmt0) 'logging_mechanical_frac = ',logging_mechanical_frac
         write(fates_log(),fmt0) 'logging_event_code = ',logging_event_code
         write(fates_log(),fmt0) 'logging_dbhmax_infra = ',logging_dbhmax_infra
-        write(fates_log(),fmt0) 'logging_export_frac = ',logging_export_frac
+!        write(fates_log(),fmt0) 'logging_export_frac = ',logging_export_frac
         write(fates_log(),*) '------------------------------------------------------'
 
      end if