Fire emissions history output.

biogeochem/EDCanopyStructureMod.F90

@@ -1877,6 +1877,8 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
     use EDTypesMod        , only : ed_site_type, AREA
     use FatesPatchMod,      only : fates_patch_type
     use FatesInterfaceTypesMod , only : bc_out_type
+    use EDParamsMod, only : num_emission_compounds
+    use FatesConstantsMod   , only : sec_per_day
 
     !
     ! !ARGUMENTS
@@ -1903,6 +1905,8 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
        bc_out(s)%dleaf_pa(:) = 0._r8
        bc_out(s)%z0m_pa(:) = 0._r8
        bc_out(s)%displa_pa(:) = 0._r8
+       bc_out(s)%fire_emissions_pa(:,:) = 0._r8
+       bc_out(s)%fire_emission_height_pa(:) = 0._r8
 
        currentPatch => sites(s)%oldest_patch
        c = fcolumn(s)
@@ -1956,6 +1960,13 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
                    currentCohort => currentCohort%taller
                 end do
 
+                do c = 1, num_emission_compounds
+                   !is this in the right place in the code?
+                   !what else do we need to do with restarts to not mess up the first day of emissions?
+                   bc_out(s)%fire_emissions_pa(ifp,c) = currentPatch%fire_emissions(c) / sec_per_day
+                end do
+                bc_out(s)%fire_emission_height_pa(ifp) = currentPatch%fire_emission_height
+
                 ! make sure there is some leaf and stem area
                 if (total_patch_leaf_stem_area > nearzero) then
                    currentCohort => currentPatch%shortest

biogeochem/EDPatchDynamicsMod.F90

@@ -3136,6 +3136,7 @@ subroutine fuse_2_patches(csite, dp, rp)
     !
     ! !USES:
     use FatesSizeAgeTypeIndicesMod, only: get_age_class_index
+    use EDParamsMod, only: num_emission_compounds
     !
     ! !ARGUMENTS:
     type (ed_site_type), intent(inout),target :: csite  ! Current site 
@@ -3203,6 +3204,10 @@ subroutine fuse_2_patches(csite, dp, rp)
     rp%ros_front               = (dp%ros_front*dp%area + rp%ros_front*rp%area) * inv_sum_area
     rp%tau_l                   = (dp%tau_l*dp%area + rp%tau_l*rp%area) * inv_sum_area
     rp%tfc_ros              = (dp%tfc_ros*dp%area + rp%tfc_ros*rp%area) * inv_sum_area
+    rp%fire_emission_height = (dp%fire_emission_height*dp%area + rp%fire_emission_height*rp%area) * inv_sum_area
+    do c = 1, num_emission_compounds
+       rp%fire_emissions(c) = (dp%fire_emissions(c)*dp%area + rp%fire_emissions(c)*rp%area) * inv_sum_area
+    enddo 
     rp%fi                   = (dp%fi*dp%area + rp%fi*rp%area) * inv_sum_area
     rp%fd                   = (dp%fd*dp%area + rp%fd*rp%area) * inv_sum_area
     rp%ros_back             = (dp%ros_back*dp%area + rp%ros_back*rp%area) * inv_sum_area

biogeochem/FatesPatchMod.F90

@@ -22,6 +22,7 @@ module FatesPatchMod
   use PRTParametersMod,       only : prt_params
   use FatesConstantsMod,      only : nocomp_bareground
   use EDParamsMod,            only : nlevleaf, nclmax, maxpft
+  use EDParamsMod,            only : num_emission_compounds
   use FatesConstantsMod,      only : n_dbh_bins, n_dist_types
   use FatesConstantsMod,      only : t_water_freeze_k_1atm
   use FatesRunningMeanMod,    only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm
@@ -216,6 +217,9 @@ module FatesPatchMod
     ! fire effects      
     real(r8)              :: scorch_ht(maxpft)       ! scorch height [m] 
     real(r8)              :: tfc_ros                 ! total intensity-relevant fuel consumed - no trunks [kgC/m2 of burned ground/day]
+    real(r8)              :: fire_emission_height    ! Height of fire emissions into the atmsophere. m.
+    real(r8)              :: fire_emissions(num_emission_compounds) ! Emissions from fires. g emissions/m2/day
+
     !---------------------------------------------------------------------------
 
     ! PLANT HYDRAULICS (not currently used in hydraulics RGK 03-2018)  
@@ -507,6 +511,9 @@ subroutine NanValues(this)
       this%scorch_ht(:)                 = nan 
       this%tfc_ros                      = nan
       this%frac_burnt                   = nan
+      this%fire_emissions(:)            = nan
+      this%fire_emission_height         = nan
+
 
     end subroutine NanValues
 
@@ -592,6 +599,8 @@ subroutine ZeroValues(this)
       this%fd                                = 0.0_r8
       this%scorch_ht(:)                      = 0.0_r8  
       this%tfc_ros                           = 0.0_r8
+      this%fire_emissions(:)                 = 0.0_r8
+      this%fire_emission_height              = 0.0_r8     
       this%frac_burnt                        = 0.0_r8
 
     end subroutine ZeroValues

fire/SFMainMod.F90

@@ -37,6 +37,24 @@ module SFMainMod
   implicit none
   private
 
+
+  public :: fire_model
+  public :: fire_danger_index 
+  public :: charecteristics_of_fuel
+  public :: rate_of_spread
+  public :: ground_fuel_consumption
+  public :: wind_effect
+  public :: area_burnt_intensity
+  public :: crown_scorching
+  public :: crown_damage
+  public :: cambial_damage_kill
+  public :: post_fire_mortality
+  public :: fire_emissions
+
+  ! The following parameter represents one of the values of hlm_spitfire_mode
+  ! and more of these appear in subroutine area_burnt_intensity below
+  ! NB. The same parameters are set in /src/biogeochem/CNFireFactoryMod
+
   public :: DailyFireModel
   public :: UpdateFuelCharacteristics
 
@@ -70,6 +88,8 @@ subroutine DailyFireModel(currentSite, bc_in)
       call crown_damage(currentSite)
       call cambial_damage_kill(currentSite)
       call post_fire_mortality(currentSite)
+      call fire_emissions(currentSite)
+
     end if
 
   end subroutine DailyFireModel
@@ -693,5 +713,46 @@ subroutine  post_fire_mortality ( currentSite )
 
   end subroutine post_fire_mortality
 
+  !*****************************************************************                                   
+  subroutine  fire_emissions ( currentSite )
+  !*****************************************************************                                   
+
+    use EDParamsMod                , only : num_emission_compounds
+    use FatesInterfaceTypesMod        , only : hlm_use_nocomp
+
+    type(ed_site_type), intent(in), target :: currentSite
+    type(fates_patch_type),  pointer :: currentPatch
+    type(fates_cohort_type), pointer :: currentCohort
+    real(r8) biomass_burned  ! Local biomass burned variable
+    real(r8) emission_factor ! Local emission factor variable
+    integer c
+
+    currentPatch => currentSite%oldest_patch
+    if(hlm_use_nocomp == itrue)then
+       ! Do not do fire emissions if we are not in nocomp mode
+       !this capability has not been added yet. 
+
+       do while(associated(currentPatch))
+          if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then
+             biomass_burned = currentPatch%TFC_ROS / 0.45_r8 ! kg biomass/m2/day
+             !n.b. does this also need to include the amount of tree canopy consumed?
+
+             do c = 1, num_emission_compounds
+                emission_factor = EDPftvarcon_inst%fire_emission_factors(currentPatch%nocomp_pft_label,c)
+                currentPatch%fire_emissions(c) = biomass_burned * emission_factor
+             enddo 
+
+             currentPatch%fire_emission_height =  EDPftvarcon_inst%fire_emission_heights(currentPatch%nocomp_pft_label)
+             if(currentPatch%fire_emissions(c).gt.0.0_r8)then
+                write(*,*) 'postivie emissionssf',currentPatch%fire_emissions(1),biomass_burned
+             endif
+          endif ! bare ground
+          currentPatch => currentPatch%younger
+
+       enddo !end patch loop
+    end if ! is nocomp 
+
+  end subroutine fire_emissions   
+
   ! ============================================================================
 end module SFMainMod

main/EDInitMod.F90

@@ -1008,6 +1008,9 @@ subroutine init_patches( nsites, sites, bc_in)
           currentPatch%ros_back                   = 0._r8
           currentPatch%scorch_ht(:)               = 0._r8
           currentPatch%frac_burnt                 = 0._r8
+          currentPatch%fire_emissions(:)          = 0._r8
+          currentPatch%fire_emission_height       = 0._r8
+
           currentPatch => currentPatch%older
        enddo
     enddo

main/EDParamsMod.F90

@@ -93,6 +93,8 @@ module EDParamsMod
 
    integer, parameter, public :: maxpft = 16      ! maximum number of PFTs allowed
 
+   integer, parameter, public :: num_emission_compounds = 13 ! number of fire emission fields calculated
+
    real(r8),protected,public  :: q10_mr     ! Q10 for respiration rate (for soil fragmenation and plant respiration)    (unitless)
    real(r8),protected,public  :: q10_froz   ! Q10 for frozen-soil respiration rates (for soil fragmentation)            (unitless)
 
@@ -342,7 +344,7 @@ subroutine FatesRegisterParams(fates_params)
     use FatesParametersInterface, only : dimension_name_history_size_bins, dimension_name_history_age_bins
     use FatesParametersInterface, only : dimension_name_history_height_bins, dimension_name_hydr_organs
     use FatesParametersInterface, only : dimension_name_history_coage_bins, dimension_name_history_damage_bins
-    use FatesParametersInterface, only : dimension_shape_scalar, dimension_name_landuse
+    use FatesParametersInterface, only : dimension_shape_scalar, dimension_name_landuse, dimension_name_nemission_compounds
 
 
     implicit none
@@ -357,7 +359,8 @@ subroutine FatesRegisterParams(fates_params)
     character(len=param_string_length), parameter :: dim_names_hydro_organs(1) = (/dimension_name_hydr_organs/)
     character(len=param_string_length), parameter :: dim_names_damageclass(1)= (/dimension_name_history_damage_bins/)
     character(len=param_string_length), parameter :: dim_names_landuse(1)= (/dimension_name_landuse/)
-
+    character(len=param_string_length), parameter :: dim_names_nemission_compounds(1)= (/dimension_name_nemission_compounds/)
+
     call FatesParamsInit()
 
     call fates_params%RegisterParameter(name=ED_name_photo_temp_acclim_timescale, dimension_shape=dimension_shape_scalar, &
@@ -550,6 +553,7 @@ subroutine FatesRegisterParams(fates_params)
 
     call fates_params%RegisterParameter(name=name_landuse_grazing_rate, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names_landuse)
+
   end subroutine FatesRegisterParams
 
 

main/EDPftvarcon.F90

@@ -129,7 +129,9 @@ module EDPftvarcon
      ! -------------------------------------------------------------------------------------------
      real(r8), allocatable :: fire_alpha_SH(:)      ! spitfire parameter, alpha scorch height
                                                     ! Equation 16 Thonicke et al 2010
-
+     real(r8), allocatable :: fire_emission_factors(:,:) ! emission factorss indexed by PFT and emission type/species. g emissions/m2/s per kg biomass burned. 
+     real(r8), allocatable :: fire_emission_heights(:)   ! heights that fire emissions are injected into the atmosphere (m)
+
      ! Non-PARTEH Allometry Parameters
      ! --------------------------------------------------------------------------------------------
 
@@ -342,7 +344,8 @@ subroutine Register_PFT(this, fates_params)
 
     use FatesParametersInterface, only : fates_parameters_type, param_string_length
     use FatesParametersInterface, only : dimension_name_pft, dimension_shape_1d
-    use FatesParametersInterface, only : dimension_name_hlm_pftno, dimension_shape_2d
+    use FatesParametersInterface, only : dimension_name_hlm_pftno, dimension_name_nemission_compounds
+    use FatesParametersInterface, only : dimension_shape_2d
 
     implicit none
 
@@ -351,7 +354,8 @@ subroutine Register_PFT(this, fates_params)
 
     character(len=param_string_length), parameter :: dim_names(1) = (/dimension_name_pft/)
     character(len=param_string_length) :: pftmap_dim_names(2)
-
+    character(len=param_string_length) :: emission_factors_dim_names(2)    
+
     integer, parameter :: dim_lower_bound(1) = (/ lower_bound_pft /)
 
 
@@ -473,6 +477,18 @@ subroutine Register_PFT(this, fates_params)
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
           dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
+    ! adding the hlm_pft_map variable with two dimensions - FATES PFTno and N emission factors
+    emission_factors_dim_names(2) = dimension_name_pft
+    emission_factors_dim_names(1) = dimension_name_nemission_compounds
+
+    name = 'fates_fire_emission_factors'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_2d, &
+         dimension_names=emission_factors_dim_names, lower_bounds=dim_lower_bound)
+
+    name = 'fates_fire_emission_heights'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+          dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
     name = 'fates_allom_frbstor_repro'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -882,6 +898,14 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%fire_alpha_SH)
 
+    name = 'fates_fire_emission_factors'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=this%fire_emission_factors)
+
+    name = 'fates_fire_emission_heights'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=this%fire_emission_heights)
+
     name = 'fates_allom_frbstor_repro'
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%allom_frbstor_repro)
@@ -1177,7 +1201,7 @@ subroutine Register_PFT_numrad(this, fates_params)
     ! arrays. We have to register the parameters as 1-d arrays as they
     ! are on the parameter file. We store them as 2-d in the receive step.
     use FatesParametersInterface, only : fates_parameters_type, param_string_length
-    use FatesParametersInterface, only : dimension_name_pft, dimension_shape_1d
+    use FatesParametersInterface, only : dimension_name_pft, dimension_name_nemission_compounds, dimension_shape_1d
 
     implicit none
 
@@ -1634,6 +1658,8 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'phen_flush_fraction',EDpftvarcon_inst%phenflush_fraction
         write(fates_log(),fmt0) 'phen_cold_size_threshold = ',EDPftvarcon_inst%phen_cold_size_threshold
         write(fates_log(),fmt0) 'fire_alpha_SH = ',EDPftvarcon_inst%fire_alpha_SH
+        write(fates_log(),fmt0) 'fire_emission_factors = ',EDPftvarcon_inst%fire_emission_factors
+        write(fates_log(),fmt0) 'fire_emission_factors = ',EDPftvarcon_inst%fire_emission_heights
         write(fates_log(),fmt0) 'allom_frbstor_repro = ',EDPftvarcon_inst%allom_frbstor_repro
         write(fates_log(),fmt0) 'hydro_p_taper = ',EDPftvarcon_inst%hydr_p_taper
         write(fates_log(),fmt0) 'hydro_rs2 = ',EDPftvarcon_inst%hydr_rs2

main/EDTypesMod.F90

@@ -88,8 +88,6 @@ module EDTypesMod
                                                                    ! The actual number of soil layers should not exceed this
 
 
-
-
   ! BIOLOGY/BIOGEOCHEMISTRY        
   integer , parameter, public :: num_vegtemp_mem      = 10         ! Window of time over which we track temp for cold sensecence (days)
 
@@ -457,6 +455,7 @@ module EDTypesMod
      real(r8) ::  NF_successful                                ! daily ignitions in km2 that actually lead to fire
      class(fire_weather), pointer :: fireWeather               ! fire weather object
 
+
      ! PLANT HYDRAULICS
      type(ed_site_hydr_type), pointer :: si_hydr