Skip to content

Commit

Permalink
Merge pull request #561 from jkshuman/ignition-fixes
Browse files Browse the repository at this point in the history 
Ignition fixes
  • Loading branch information
@rgknox
rgknox authored Sep 13, 2019
2 parents 3f6749d + bd688a3 commit 393b905
Show file tree
Hide file tree
Showing 5 changed files with 61 additions and 81 deletions.
18 changes: 7 additions & 11 deletions biogeochem/EDPatchDynamicsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,14 +252,14 @@ subroutine disturbance_rates( site_in, bc_in)
endif

! Fire Disturbance Rate
! Fudge - fires can't burn the whole patch, as this causes /0 errors.
! This is accumulating the daily fires over the whole 30 day patch generation phase.
currentPatch%disturbance_rates(dtype_ifire) = &
min(0.99_r8,currentPatch%disturbance_rates(dtype_ifire) + currentPatch%frac_burnt)
! Fires can't burn the whole patch, as this causes /0 errors.
currentPatch%disturbance_rates(dtype_ifire) = currentPatch%frac_burnt

if (currentPatch%disturbance_rates(dtype_ifire) > 0.98_r8)then
if (debug) then
if (currentPatch%disturbance_rates(dtype_ifire) > 0.98_r8)then
write(fates_log(),*) 'very high fire areas', &
currentPatch%disturbance_rates(dtype_ifire),currentPatch%frac_burnt
currentPatch%disturbance_rates(dtype_ifire),currentPatch%frac_burnt
endif
endif


Expand Down Expand Up @@ -1980,10 +1980,8 @@ subroutine zero_patch(cp_p)
currentPatch%fire = 999 ! sr decide_fire.1=fire hot enough to proceed. 0=stop everything- no fires today
currentPatch%fd = 0.0_r8 ! fire duration (mins)
currentPatch%ros_back = 0.0_r8 ! backward ros (m/min)
currentPatch%ab = 0.0_r8 ! area burnt daily m2
currentPatch%nf = 0.0_r8 ! number of fires initiated daily
currentPatch%sh = 0.0_r8 ! average scorch height for the patch(m)
currentPatch%frac_burnt = 0.0_r8 ! fraction burnt in each timestep.
currentPatch%frac_burnt = 0.0_r8 ! fraction burnt daily
currentPatch%burnt_frac_litter(:) = 0.0_r8
currentPatch%btran_ft(:) = 0.0_r8

Expand Down Expand Up @@ -2298,8 +2296,6 @@ subroutine fuse_2_patches(csite, dp, rp)
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
rp%ab = (dp%ab*dp%area + rp%ab*rp%area) * inv_sum_area
rp%nf = (dp%nf*dp%area + rp%nf*rp%area) * inv_sum_area
rp%sh = (dp%sh*dp%area + rp%sh*rp%area) * inv_sum_area
rp%frac_burnt = (dp%frac_burnt*dp%area + rp%frac_burnt*rp%area) * inv_sum_area
rp%burnt_frac_litter(:) = (dp%burnt_frac_litter(:)*dp%area + rp%burnt_frac_litter(:)*rp%area) * inv_sum_area
Expand Down
106 changes: 46 additions & 60 deletions fire/SFMainMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,6 @@ subroutine fire_model( currentSite, bc_in)
currentPatch => currentSite%youngest_patch
do while(associated(currentPatch))
currentPatch%frac_burnt = 0.0_r8
currentPatch%AB = 0.0_r8
currentPatch%fire = 0
currentPatch => currentPatch%older
enddo
Expand Down Expand Up @@ -121,14 +120,14 @@ subroutine fire_danger_index ( currentSite, bc_in)
type(ed_site_type) , intent(inout), target :: currentSite
type(bc_in_type) , intent(in) :: bc_in

real(r8) :: temp_in_C ! daily averaged temperature in celcius
real(r8) :: rainfall ! daily precip in mm/day
real(r8) :: rh ! daily rh
real(r8) :: temp_in_C ! daily averaged temperature in celcius
real(r8) :: rainfall ! daily precip in mm/day
real(r8) :: rh ! daily rh

real yipsolon; !intermediate varable for dewpoint calculation
real dewpoint; !dewpoint in K
real d_NI; !daily change in Nesterov Index. C^2
integer :: iofp ! index of oldest the fates patch
real(r8) :: yipsolon !intermediate varable for dewpoint calculation
real(r8) :: dewpoint !dewpoint in K
real(r8) :: d_NI !daily change in Nesterov Index. C^2
integer :: iofp ! index of oldest the fates patch

! NOTE that the boundary conditions of temperature, precipitation and relative humidity
! are available at the patch level. We are currently using a simplification where the whole site
Expand Down Expand Up @@ -169,7 +168,6 @@ subroutine charecteristics_of_fuel ( currentSite )
type(ed_cohort_type), pointer :: currentCohort
type(litter_type), pointer :: litt_c

real(r8) timeav_swc
real(r8) alpha_FMC(nfsc) ! Relative fuel moisture adjusted per drying ratio
real(r8) fuel_moisture(nfsc) ! Scaled moisture content of small litter fuels.
real(r8) MEF(nfsc) ! Moisture extinction factor of fuels integer n
Expand Down Expand Up @@ -204,7 +202,7 @@ subroutine charecteristics_of_fuel ( currentSite )

! zero fire arrays.
currentPatch%fuel_eff_moist = 0.0_r8
currentPatch%fuel_bulkd = 0.0_r8 !this is kgBiomass/m2 for use in rate of spread equations
currentPatch%fuel_bulkd = 0.0_r8 !this is kgBiomass/m3 for use in rate of spread equations
currentPatch%fuel_sav = 0.0_r8
currentPatch%fuel_frac(:) = 0.0_r8
currentPatch%fuel_mef = 0.0_r8
Expand Down Expand Up @@ -241,13 +239,13 @@ subroutine charecteristics_of_fuel ( currentSite )
endif

currentPatch%fuel_frac(lg_sf) = currentPatch%livegrass / currentPatch%sum_fuel
MEF(1:nfsc) = 0.524_r8 - 0.066_r8 * log10(SF_val_SAV(1:nfsc))
MEF(1:nfsc) = 0.524_r8 - 0.066_r8 * log10(SF_val_SAV(1:nfsc))

!--- weighted average of relative moisture content---
! Equation 6 in Thonicke et al. 2010. across twig, small branch, large branch, and dead leaves
! dead leaves and twigs included in 1hr pool per Thonicke (2010)
! Calculate fuel moisture for trunks to hold value for fuel consumption
alpha_FMC(tw_sf:dl_sf) = SF_val_SAV(tw_sf:dl_sf)/SF_val_drying_ratio
alpha_FMC(tw_sf:dl_sf) = SF_val_SAV(tw_sf:dl_sf)/SF_val_drying_ratio

fuel_moisture(tw_sf:dl_sf) = exp(-1.0_r8 * alpha_FMC(tw_sf:dl_sf) * currentSite%acc_NI)

Expand All @@ -257,12 +255,11 @@ subroutine charecteristics_of_fuel ( currentSite )
if ( hlm_masterproc == itrue ) write(fates_log(),*) 'csa ',currentSite%acc_NI
if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sfv ',alpha_FMC
endif
! FIX(RF,032414): needs refactoring.
! average water content !is this the correct metric?
timeav_swc = sum(currentSite%water_memory(1:numWaterMem)) / real(numWaterMem,r8)
! Equation B2 in Thonicke et al. 2010
! live grass moisture content depends on upper soil layer
fuel_moisture(lg_sf) = max(0.0_r8, 10.0_r8/9._r8 * timeav_swc - 1.0_r8/9.0_r8)

! live grass moisture is a function of SAV and changes via Nesterov Index
! along the same relationship as the 1 hour fuels (live grass has same SAV as dead grass,
! but retains more moisture with this calculation.)
fuel_moisture(lg_sf) = exp(-1.0_r8 * ((SF_val_SAV(tw_sf)/SF_val_drying_ratio) * currentSite%acc_NI))

! Average properties over the first three litter pools (twigs, s branches, l branches)
currentPatch%fuel_bulkd = sum(currentPatch%fuel_frac(tw_sf:lb_sf) * SF_val_FBD(tw_sf:lb_sf))
Expand Down Expand Up @@ -481,7 +478,8 @@ subroutine rate_of_spread ( currentSite )

! ---heat of pre-ignition---
! Equation A4 in Thonicke et al. 2010
! conversion of Rohtermal (1972) equation 12 in BTU/lb to current kJ/kg
! Rothermal EQ12= 250 Btu/lb + 1116 Btu/lb * fuel_eff_moist
! conversion of Rothermal (1972) EQ12 in BTU/lb to current kJ/kg
! q_ig in kJ/kg
q_ig = 581.0_r8 +2594.0_r8 * currentPatch%fuel_eff_moist

Expand Down Expand Up @@ -678,7 +676,10 @@ subroutine fire_intensity ( currentSite )
do while(associated(currentPatch))
ROS = currentPatch%ROS_front / 60.0_r8 !m/min to m/sec
W = currentPatch%TFC_ROS / 0.45_r8 !kgC/m2 to kgbiomass/m2

!units of fire intensity = (kJ/kg)*(kgBiomass/m2)*(m/min)
currentPatch%FI = SF_val_fuel_energy * W * ROS !kj/m/s, or kW/m

if(write_sf == itrue)then
if( hlm_masterproc == itrue ) write(fates_log(),*) 'fire_intensity',currentPatch%fi,W,currentPatch%ROS_front
endif
Expand All @@ -688,7 +689,8 @@ subroutine fire_intensity ( currentSite )

! Equation 7 from Venevsky et al GCB 2002 (modification of equation 8 in Thonicke et al. 2010)
! FDI 0.1 = low, 0.3 moderate, 0.75 high, and 1 = extreme ignition potential for alpha 0.000337
currentSite%FDI = 1.0_r8 - exp(-SF_val_fdi_alpha*currentSite%acc_NI)
currentSite%FDI = 1.0_r8 - exp(-SF_val_fdi_alpha*currentSite%acc_NI)

! Equation 14 in Thonicke et al. 2010
! fire duration in minutes

Expand Down Expand Up @@ -719,30 +721,35 @@ end subroutine fire_intensity
!*****************************************************************
subroutine area_burnt ( currentSite )
!*****************************************************************
!currentPatch%AB !daily area burnt (m2)
!currentPatch%NF !Daily number of ignitions (lightning and human-caused), adjusted for size of patch.

use EDParamsMod, only : ED_val_nignitions
use EDParamsMod, only : ED_val_nignitions
use FatesConstantsMod, only : years_per_day

type(ed_site_type), intent(inout), target :: currentSite
type(ed_patch_type), pointer :: currentPatch

real lb !length to breadth ratio of fire ellipse
real df !distance fire has travelled forward
real db !distance fire has travelled backward
real patch_area_in_m2 !'actual' patch area as applied to whole grid cell
real(r8) gridarea
real(r8) lb !length to breadth ratio of fire ellipse (unitless)
real(r8) df !distance fire has travelled forward in m
real(r8) db !distance fire has travelled backward in m
real(r8) NF !number of lightning strikes per day per km2
real(r8) AB !daily area burnt in m2 per km2
real(r8) size_of_fire !in m2
real(r8),parameter :: km2_to_m2 = 1000000.0_r8 !area conversion for square km to square m
integer g, p

currentSite%frac_burnt = 0.0_r8
! ---initialize site parameters to zero---
currentSite%frac_burnt = 0.0_r8

!NF = number of lighting strikes per day per km2
NF = ED_val_nignitions * years_per_day

! If there are 15 lightning strikes per year, per km2. (approx from NASA product for S.A.)
! then there are 15 * 1/365 strikes/km2 each day.

currentPatch => currentSite%oldest_patch;
do while(associated(currentPatch))
currentPatch%AB = 0.0_r8
! ---initialize patch parameters to zero---
currentPatch%frac_burnt = 0.0_r8
lb = 0.0_r8; db = 0.0_r8; df = 0.0_r8

if (currentPatch%fire == 1) then
! The feedback between vegetation structure and ellipse size if turned off for now,
Expand All @@ -767,50 +774,29 @@ subroutine area_burnt ( currentSite )

! --- calculate area burnt---
if(lb > 0.0_r8) then

! INTERF-TODO:
! THIS SHOULD HAVE THE COLUMN AND LU AREA WEIGHT ALSO, NO?

gridarea = km2_to_m2 ! 1M m2 in a km2

! NF = number of lighting strikes per day per km2
currentPatch%NF = ED_val_nignitions * currentPatch%area/area /365

! If there are 15 lightening strickes per year, per km2. (approx from NASA product)
! then there are 15/365 s/km2 each day.

! Equation 1 in Thonicke et al. 2010
! To Do: Connect here with the Li & Levis GDP fire suppression algorithm.
! Equation 16 in arora and boer model JGR 2005
!currentPatch%AB = currentPatch%AB *3.0_r8
! AB = AB *3.0_r8

!size of fire = equation 14 Arora and Boer JGR 2005
size_of_fire = ((3.1416_r8/(4.0_r8*lb))*((df+db)**2.0_r8))

!AB = daily area burnt = size fires in m2 * num ignitions * prob ignition starts fire
! m2 per km2 per day
currentPatch%AB = size_of_fire * currentPatch%NF * currentSite%FDI

patch_area_in_m2 = gridarea *currentPatch%area/area
!AB = daily area burnt = size fires in m2 * num ignitions per day per km2 * prob ignition starts fire
!AB = m2 per km2 per day
AB = size_of_fire * NF * currentSite%FDI

!frac_burnt in units of m2 here.
currentPatch%frac_burnt = min(0.99_r8, AB / km2_to_m2)

currentPatch%frac_burnt = currentPatch%AB / patch_area_in_m2
if(write_SF == itrue)then
if ( hlm_masterproc == itrue ) write(fates_log(),*) 'frac_burnt',currentPatch%frac_burnt
endif

if (currentPatch%frac_burnt > 1.0_r8 ) then !all of patch burnt.

currentPatch%frac_burnt = 1.0_r8 ! capping at 1 same as %AB/patch_area_in_m2

if ( hlm_masterproc == itrue ) write(fates_log(),*) 'burnt all of patch',currentPatch%patchno
if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ros',currentPatch%ROS_front,currentPatch%FD, &
currentPatch%NF,currentPatch%FI,size_of_fire

endif


endif
endif! fire
! convert frac_burnt to % prior to accumulating at site level
currentSite%frac_burnt = currentSite%frac_burnt + currentPatch%frac_burnt * currentPatch%area/area

currentPatch => currentPatch%younger
Expand Down
2 changes: 0 additions & 2 deletions main/EDTypesMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -519,11 +519,9 @@ module EDTypesMod
real(r8) :: fi ! average fire intensity of flaming front: kj/m/s or kw/m
integer :: fire ! Is there a fire? 1=yes 0=no
real(r8) :: fd ! fire duration: mins
real(r8) :: nf ! number of fires initiated daily: n/gridcell/day
real(r8) :: sh ! average scorch height: m

! FIRE EFFECTS
real(r8) :: ab ! area burnt: m2/day
real(r8) :: frac_burnt ! fraction burnt: frac gridcell/day
real(r8) :: tfc_ros ! total fuel consumed - no trunks. KgC/m2/day
real(r8) :: burnt_frac_litter(nfsc) ! fraction of each litter pool burned:-
Expand Down
10 changes: 5 additions & 5 deletions main/FatesHistoryInterfaceMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3891,7 +3891,7 @@ subroutine define_history_vars(this, initialize_variables)
avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, &
ivar=ivar, initialize=initialize_variables, index = ih_effect_wspeed_pa )

call this%set_history_var(vname='FIRE_TFC_ROS', units='none', &
call this%set_history_var(vname='FIRE_TFC_ROS', units='kgC/m2', &
long ='total fuel consumed', use_default='active', &
avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, &
ivar=ivar, initialize=initialize_variables, index = ih_TFC_ROS_pa )
Expand All @@ -3902,12 +3902,12 @@ subroutine define_history_vars(this, initialize_variables)
ivar=ivar, initialize=initialize_variables, index = ih_fire_intensity_pa )

call this%set_history_var(vname='FIRE_AREA', units='fraction', &
long='spitfire fire area:m2', use_default='active', &
long='spitfire fire area burn fraction', use_default='active', &
avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, &
ivar=ivar, initialize=initialize_variables, index = ih_fire_area_pa )

call this%set_history_var(vname='SCORCH_HEIGHT', units='m', &
long='spitfire fire area:m2', use_default='active', &
long='spitfire flame height:m', use_default='active', &
avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, &
ivar=ivar, initialize=initialize_variables, index = ih_scorch_height_pa )

Expand All @@ -3916,7 +3916,7 @@ subroutine define_history_vars(this, initialize_variables)
avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, &
ivar=ivar, initialize=initialize_variables, index = ih_fire_fuel_mef_pa )

call this%set_history_var(vname='fire_fuel_bulkd', units='m', &
call this%set_history_var(vname='fire_fuel_bulkd', units='kg biomass/m3', &
long='spitfire fuel bulk density', use_default='active', &
avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, &
ivar=ivar, initialize=initialize_variables, index = ih_fire_fuel_bulkd_pa )
Expand All @@ -3926,7 +3926,7 @@ subroutine define_history_vars(this, initialize_variables)
avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, &
ivar=ivar, initialize=initialize_variables, index = ih_fire_fuel_eff_moist_pa )

call this%set_history_var(vname='fire_fuel_sav', units='m', &
call this%set_history_var(vname='fire_fuel_sav', units='per m', &
long='spitfire fuel surface/volume ', use_default='active', &
avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, &
ivar=ivar, initialize=initialize_variables, index = ih_fire_fuel_sav_pa )
Expand Down
6 changes: 3 additions & 3 deletions parameter_files/fates_params_default.cdl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -502,8 +502,8 @@ variables:
fates_fire_miner_total:units = "fraction" ;
fates_fire_miner_total:long_name = "spitfire parameter, total mineral content, Table A1 Thonicke et al 2010" ;
double fates_fire_nignitions ;
fates_fire_nignitions:units = "ignitions per day" ;
fates_fire_nignitions:long_name = "number of daily ignitions" ;
fates_fire_nignitions:units = "ignitions per year per km2" ;
fates_fire_nignitions:long_name = "number of annual ignitions per square km" ;
double fates_fire_part_dens ;
fates_fire_part_dens:units = "kg/m2" ;
fates_fire_part_dens:long_name = "spitfire parameter, oven dry particle density, Table A1 Thonicke et al 2010" ;
Expand Down Expand Up @@ -1104,7 +1104,7 @@ data:

fates_cwd_flig = 0.24 ;

fates_fire_drying_ratio = 13000 ;
fates_fire_drying_ratio = 66000 ;

fates_fire_durat_slope = -11.06 ;

Expand Down

0 comments on commit 393b905

Please sign in to comment.