diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 696a88514c..9b02c9e14a 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -38,8 +38,7 @@ module EDPatchDynamicsMod use FatesConstantsMod , only : ican_upper use PRTGenericMod , only : num_elements use PRTGenericMod , only : element_list - use FatesLitterMod , only : lg_sf - use FatesLitterMod , only : dl_sf + use FatesFuelClassesMod , only : fuel_classes use FatesConstantsMod , only : N_DIST_TYPES use EDTypesMod , only : AREA_INV use EDTypesMod , only : dump_site @@ -765,11 +764,11 @@ subroutine spawn_patches( currentSite, bc_in) ! Transfer the litter existing already in the donor patch to the new patch ! This call will only transfer non-burned litter to new patch - ! and burned litter to atmosphere. Thus it is important to zero burnt_frac_litter when + ! and burned litter to atmosphere. Thus it is important to zero fuel%frac_burnt when ! fire is not the current disturbance regime. if(i_disturbance_type .ne. dtype_ifire) then - currentPatch%burnt_frac_litter(:) = 0._r8 + currentPatch%fuel%frac_burnt(:) = 0._r8 end if call CopyPatchMeansTimers(currentPatch, newPatch) @@ -1053,7 +1052,7 @@ subroutine spawn_patches( currentSite, bc_in) ! Grasses determine their fraction of leaves burned here - leaf_burn_frac = currentPatch%burnt_frac_litter(lg_sf) + leaf_burn_frac = currentPatch%fuel%frac_burnt(fuel_classes%live_grass()) endif ! Perform a check to make sure that spitfire gave @@ -1727,7 +1726,7 @@ subroutine split_patch(currentSite, currentPatch, new_patch, fraction_to_keep, a call TransLitterNewPatch( currentSite, currentPatch, new_patch, temp_area) - currentPatch%burnt_frac_litter(:) = 0._r8 + currentPatch%fuel%frac_burnt(:) = 0._r8 ! Next, we loop through the cohorts in the donor patch, copy them with ! area modified number density into the new-patch, and apply survivorship. @@ -2082,10 +2081,10 @@ subroutine TransLitterNewPatch(currentSite, & ! Transfer above ground CWD donatable_mass = curr_litt%ag_cwd(c) * patch_site_areadis * & - (1._r8 - currentPatch%burnt_frac_litter(c)) + (1._r8 - currentPatch%fuel%frac_burnt(c)) burned_mass = curr_litt%ag_cwd(c) * patch_site_areadis * & - currentPatch%burnt_frac_litter(c) + currentPatch%fuel%frac_burnt(c) new_litt%ag_cwd(c) = new_litt%ag_cwd(c) + donatable_mass*donate_m2 curr_litt%ag_cwd(c) = curr_litt%ag_cwd(c) + donatable_mass*retain_m2 @@ -2106,10 +2105,10 @@ subroutine TransLitterNewPatch(currentSite, & ! Transfer leaf fines donatable_mass = curr_litt%leaf_fines(dcmpy) * patch_site_areadis * & - (1._r8 - currentPatch%burnt_frac_litter(dl_sf)) + (1._r8 - currentPatch%fuel%frac_burnt(fuel_classes%dead_leaves())) burned_mass = curr_litt%leaf_fines(dcmpy) * patch_site_areadis * & - currentPatch%burnt_frac_litter(dl_sf) + currentPatch%fuel%frac_burnt(fuel_classes%dead_leaves()) new_litt%leaf_fines(dcmpy) = new_litt%leaf_fines(dcmpy) + donatable_mass*donate_m2 curr_litt%leaf_fines(dcmpy) = curr_litt%leaf_fines(dcmpy) + donatable_mass*retain_m2 @@ -3266,6 +3265,8 @@ subroutine fuse_2_patches(csite, dp, rp) do el = 1,num_elements call rp%litter(el)%FuseLitter(rp%area,dp%area,dp%litter(el)) end do + + call rp%fuel%Fuse(rp%area, dp%area, dp%fuel) if ( rp%land_use_label .ne. dp%land_use_label) then write(fates_log(),*) 'trying to fuse patches with different land_use_label values' @@ -3296,22 +3297,15 @@ subroutine fuse_2_patches(csite, dp, rp) call rp%tveg_longterm%FuseRMean(dp%tveg_longterm,rp%area*inv_sum_area) - rp%fuel_eff_moist = (dp%fuel_eff_moist*dp%area + rp%fuel_eff_moist*rp%area) * inv_sum_area - rp%livegrass = (dp%livegrass*dp%area + rp%livegrass*rp%area) * inv_sum_area - rp%sum_fuel = (dp%sum_fuel*dp%area + rp%sum_fuel*rp%area) * inv_sum_area - rp%fuel_bulkd = (dp%fuel_bulkd*dp%area + rp%fuel_bulkd*rp%area) * inv_sum_area - rp%fuel_sav = (dp%fuel_sav*dp%area + rp%fuel_sav*rp%area) * inv_sum_area - rp%fuel_mef = (dp%fuel_mef*dp%area + rp%fuel_mef*rp%area) * inv_sum_area - 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%fuel_frac(:) = (dp%fuel_frac(:)*dp%area + rp%fuel_frac(:)*rp%area) * inv_sum_area + rp%livegrass = (dp%livegrass*dp%area + rp%livegrass*rp%area) * inv_sum_area + 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%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%scorch_ht(:) = (dp%scorch_ht(:)*dp%area + rp%scorch_ht(:)*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 rp%btran_ft(:) = (dp%btran_ft(:)*dp%area + rp%btran_ft(:)*rp%area) * inv_sum_area rp%zstar = (dp%zstar*dp%area + rp%zstar*rp%area) * inv_sum_area rp%c_stomata = (dp%c_stomata*dp%area + rp%c_stomata*rp%area) * inv_sum_area @@ -3406,7 +3400,6 @@ subroutine fuse_2_patches(csite, dp, rp) olderp%younger => null() end if - if(associated(olderp))then ! Update the older patch's new younger patch (becuase it isn't dp anymore) olderp%younger => youngerp diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index dddfd83113..8d689a759a 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -53,8 +53,8 @@ module EDPhysiologyMod use EDTypesMod , only : site_massbal_type use EDTypesMod , only : numlevsoil_max use EDTypesMod , only : numWaterMem + use FatesFuelClassesMod , only : fuel_classes use EDTypesMod , only : elem_diag_type - use FatesLitterMod , only : dl_sf use EDParamsMod , only : dinc_vai, dlower_vai use EDTypesMod , only : area_inv use EDTypesMod , only : AREA @@ -3269,11 +3269,11 @@ subroutine CWDOut( litt, fragmentation_scaler, nlev_eff_decomp ) do dcmpy = 1,ndcmpy litt%leaf_fines_frag(dcmpy) = litt%leaf_fines(dcmpy) * & - years_per_day * SF_val_max_decomp(dl_sf) * fragmentation_scaler(soil_layer_index) + years_per_day * SF_val_max_decomp(fuel_classes%dead_leaves()) * fragmentation_scaler(soil_layer_index) do ilyr = 1,nlev_eff_decomp litt%root_fines_frag(dcmpy,ilyr) = litt%root_fines(dcmpy,ilyr) * & - years_per_day * SF_val_max_decomp(dl_sf) * fragmentation_scaler(ilyr) + years_per_day * SF_val_max_decomp(fuel_classes%dead_leaves()) * fragmentation_scaler(ilyr) end do enddo diff --git a/biogeochem/FatesLitterMod.F90 b/biogeochem/FatesLitterMod.F90 index 7d55dc9aab..9a240d45e5 100644 --- a/biogeochem/FatesLitterMod.F90 +++ b/biogeochem/FatesLitterMod.F90 @@ -54,18 +54,6 @@ module FatesLitterMod integer, public, parameter :: ilabile = 1 ! Array index for labile portion integer, public, parameter :: icellulose = 2 ! Array index for cellulose portion integer, public, parameter :: ilignin = 3 ! Array index for the lignin portion - - ! SPITFIRE - - integer, parameter, public :: NFSC = NCWD+2 ! number fuel size classes (4 cwd size classes, leaf litter, and grass) - integer, parameter, public :: tw_sf = 1 ! array index of twig pool for spitfire - integer, parameter, public :: lb_sf = 3 ! array index of large branch pool for spitfire - integer, parameter, public :: tr_sf = 4 ! array index of dead trunk pool for spitfire - integer, parameter, public :: dl_sf = 5 ! array index of dead leaf pool for spitfire (dead grass and dead leaves) - integer, parameter, public :: lg_sf = 6 ! array index of live grass pool for spitfire - - - type, public :: litter_type ! This object is allocated for each element (C, N, P, etc) that we wish to track. diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 999bd73228..f8afc711db 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -6,17 +6,19 @@ module FatesPatchMod use FatesConstantsMod, only : primaryland, secondaryland use FatesConstantsMod, only : n_landuse_cats use FatesConstantsMod, only : TRS_regeneration - use FatesConstantsMod, only : itrue + use FatesConstantsMod, only : itrue, ifalse use FatesGlobals, only : fates_log use FatesGlobals, only : endrun => fates_endrun use FatesUtilsMod, only : check_hlm_list use FatesUtilsMod, only : check_var_real use FatesCohortMod, only : fates_cohort_type use FatesRunningMeanMod, only : rmean_type, rmean_arr_type - use FatesLitterMod, only : nfsc use FatesLitterMod, only : litter_type + use FatesFuelMod, only : fuel_type use PRTGenericMod, only : num_elements use PRTGenericMod, only : element_list + use PRTGenericMod, only : carbon12_element + use PRTGenericMod, only : struct_organ, leaf_organ, sapw_organ use PRTParametersMod, only : prt_params use FatesConstantsMod, only : nocomp_bareground use EDParamsMod, only : nlevleaf, nclmax, maxpft @@ -195,24 +197,14 @@ module FatesPatchMod ! LITTER AND COARSE WOODY DEBRIS type(litter_type), pointer :: litter(:) ! litter (leaf,fnrt,CWD and seeds) for different elements + type(fuel_type), pointer :: fuel ! fuel class real(r8), allocatable :: fragmentation_scaler(:) ! scale rate of litter fragmentation based on soil layer [0-1] !--------------------------------------------------------------------------- ! FUELS AND FIRE ! fuel characteristics - real(r8) :: sum_fuel ! total ground fuel related to ROS (omits 1000 hr fuels) [kgC/m2] - real(r8) :: fuel_frac(nfsc) ! fraction of each litter class in the ros_fuel [0-1] real(r8) :: livegrass ! total aboveground grass biomass in patch [kgC/m2] - real(r8) :: fuel_bulkd ! average fuel bulk density of the ground fuel. [kg/m3] - ! (incl. live grasses, omits 1000hr fuels) - real(r8) :: fuel_sav ! average surface area to volume ratio of the ground fuel [cm-1] - ! (incl. live grasses, omits 1000hr fuels) - real(r8) :: fuel_mef ! average moisture of extinction factor - ! of the ground fuel (incl. live grasses, omits 1000hr fuels) - real(r8) :: fuel_eff_moist ! effective avearage fuel moisture content of the ground fuel - ! (incl. live grasses. omits 1000hr fuels) - real(r8) :: litter_moisture(nfsc) ! moisture of litter [m3/m3] ! fire spread real(r8) :: ros_front ! rate of forward spread of fire [m/min] @@ -221,13 +213,11 @@ module FatesPatchMod 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 [min] + real(r8) :: frac_burnt ! fraction of patch burnt by fire ! fire effects real(r8) :: scorch_ht(maxpft) ! scorch height [m] - real(r8) :: frac_burnt ! fraction burnt [0-1/day] real(r8) :: tfc_ros ! total intensity-relevant fuel consumed - no trunks [kgC/m2 of burned ground/day] - real(r8) :: burnt_frac_litter(nfsc) ! fraction of each litter pool burned, conditional on it being burned [0-1] - !--------------------------------------------------------------------------- ! PLANT HYDRAULICS (not currently used in hydraulics RGK 03-2018) @@ -245,6 +235,7 @@ module FatesPatchMod procedure :: InitLitter procedure :: Create procedure :: UpdateTreeGrassArea + procedure :: UpdateLiveGrass procedure :: FreeMemory procedure :: Dump procedure :: CheckVars @@ -506,14 +497,7 @@ subroutine NanValues(this) this%fragmentation_scaler(:) = nan ! FUELS AND FIRE - this%sum_fuel = nan - this%fuel_frac(:) = nan this%livegrass = nan - this%fuel_bulkd = nan - this%fuel_sav = nan - this%fuel_mef = nan - this%fuel_eff_moist = nan - this%litter_moisture(:) = nan this%ros_front = nan this%ros_back = nan this%tau_l = nan @@ -521,10 +505,9 @@ subroutine NanValues(this) this%fire = fates_unset_int this%fd = nan this%scorch_ht(:) = nan + this%tfc_ros = nan this%frac_burnt = nan - this%tfc_ros = nan - this%burnt_frac_litter(:) = nan - + end subroutine NanValues !=========================================================================== @@ -600,23 +583,15 @@ subroutine ZeroValues(this) this%fragmentation_scaler(:) = 0.0_r8 ! FIRE - this%sum_fuel = 0.0_r8 - this%fuel_frac(:) = 0.0_r8 this%livegrass = 0.0_r8 - this%fuel_bulkd = 0.0_r8 - this%fuel_sav = 0.0_r8 - this%fuel_mef = 0.0_r8 - this%fuel_eff_moist = 0.0_r8 - this%litter_moisture(:) = 0.0_r8 this%ros_front = 0.0_r8 this%ros_back = 0.0_r8 this%tau_l = 0.0_r8 this%fi = 0.0_r8 this%fd = 0.0_r8 this%scorch_ht(:) = 0.0_r8 - this%frac_burnt = 0.0_r8 this%tfc_ros = 0.0_r8 - this%burnt_frac_litter(:) = 0.0_r8 + this%frac_burnt = 0.0_r8 end subroutine ZeroValues @@ -738,6 +713,10 @@ subroutine Create(this, age, area, land_use_label, nocomp_pft, num_swb, num_pft, ! initialize litter call this%InitLitter(num_pft, num_levsoil) + ! initialize fuel + allocate(this%fuel) + call this%fuel%Init() + this%twostr%scelg => null() ! The radiation module will check if this ! is associated, since it is not, it will then ! initialize and allocate @@ -776,7 +755,7 @@ subroutine UpdateTreeGrassArea(this) class(fates_patch_type), intent(inout) :: this ! patch object ! LOCALS: - type(fates_cohort_Type), pointer :: currentCohort ! cohort object + type(fates_cohort_type), pointer :: currentCohort ! cohort object real(r8) :: tree_area ! treed area of patch [m2] real(r8) :: grass_area ! grass area of patch [m2] @@ -802,6 +781,38 @@ end subroutine UpdateTreeGrassArea !=========================================================================== + subroutine UpdateLiveGrass(this) + ! + ! DESCRIPTION: + ! Calculates the sum of live grass biomass [kgC/m2] on a patch + + ! ARGUMENTS: + class(fates_patch_type), intent(inout) :: this ! patch + + ! LOCALS: + real(r8) :: live_grass ! live grass [kgC/m2] + type(fates_cohort_type), pointer :: currentCohort ! cohort type + + live_grass = 0.0_r8 + currentCohort => this%tallest + do while(associated(currentCohort)) + ! for grasses sum all aboveground tissues + if (prt_params%woody(currentCohort%pft) == ifalse) then + live_grass = live_grass + & + (currentCohort%prt%GetState(leaf_organ, carbon12_element) + & + currentCohort%prt%GetState(sapw_organ, carbon12_element) + & + currentCohort%prt%GetState(struct_organ, carbon12_element))* & + currentCohort%n/this%area + endif + currentCohort => currentCohort%shorter + enddo + + this%livegrass = live_grass + + end subroutine UpdateLiveGrass + + !=========================================================================== + subroutine FreeMemory(this, regeneration_model, numpft) ! ! DESCRIPTION: @@ -849,7 +860,13 @@ subroutine FreeMemory(this, regeneration_model, numpft) write(fates_log(),*) 'dealloc008: fail on deallocate(this%litter):'//trim(smsg) call endrun(msg=errMsg(sourcefile, __LINE__)) endif - + + deallocate(this%fuel, stat=istat, errmsg=smsg) + if (istat/=0) then + write(fates_log(),*) 'dealloc009: fail on deallocate patch fuel:'//trim(smsg) + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif + ! deallocate the allocatable arrays deallocate(this%tr_soil_dir, & this%tr_soil_dif, & @@ -889,24 +906,24 @@ subroutine FreeMemory(this, regeneration_model, numpft) end if if (istat/=0) then - write(fates_log(),*) 'dealloc009: fail on deallocate patch vectors:'//trim(smsg) + write(fates_log(),*) 'dealloc010: fail on deallocate patch vectors:'//trim(smsg) call endrun(msg=errMsg(sourcefile, __LINE__)) endif ! deallocate running means deallocate(this%tveg24, stat=istat, errmsg=smsg) if (istat/=0) then - write(fates_log(),*) 'dealloc010: fail on deallocate(this%tveg24):'//trim(smsg) + write(fates_log(),*) 'dealloc011: fail on deallocate(this%tveg24):'//trim(smsg) call endrun(msg=errMsg(sourcefile, __LINE__)) endif deallocate(this%tveg_lpa, stat=istat, errmsg=smsg) if (istat/=0) then - write(fates_log(),*) 'dealloc011: fail on deallocate(this%tveg_lpa):'//trim(smsg) + write(fates_log(),*) 'dealloc012: fail on deallocate(this%tveg_lpa):'//trim(smsg) call endrun(msg=errMsg(sourcefile, __LINE__)) endif deallocate(this%tveg_longterm, stat=istat, errmsg=smsg) if (istat/=0) then - write(fates_log(),*) 'dealloc012: fail on deallocate(this%tveg_longterm):'//trim(smsg) + write(fates_log(),*) 'dealloc013: fail on deallocate(this%tveg_longterm):'//trim(smsg) call endrun(msg=errMsg(sourcefile, __LINE__)) endif diff --git a/fire/CMakeLists.txt b/fire/CMakeLists.txt index 341065c3c3..f3b0f84ca5 100644 --- a/fire/CMakeLists.txt +++ b/fire/CMakeLists.txt @@ -3,6 +3,8 @@ list(APPEND fates_sources SFParamsMod.F90 SFFireWeatherMod.F90 SFNesterovMod.F90 + FatesFuelMod.F90 + FatesFuelClassesMod.F90 ) sourcelist_to_parent(fates_sources) diff --git a/fire/FatesFuelClassesMod.F90 b/fire/FatesFuelClassesMod.F90 new file mode 100644 index 0000000000..257c822a9f --- /dev/null +++ b/fire/FatesFuelClassesMod.F90 @@ -0,0 +1,63 @@ +module FatesFuelClassesMod + + use FatesLitterMod, only : ncwd + + implicit none + private + + integer, parameter, public :: num_fuel_classes = 6 ! number of total fuel classes + + type :: fuel_classes_type + ! There are six fuel classes: + ! 1) twigs, 2) small branches, 3) large branches 4) trunks + ! 5) dead leaves, 6) live grass + integer, private :: twigs_i = 1 ! array index for twigs pool + integer, private :: small_branches_i = 2 ! array index for small branches pool + integer, private :: large_branches_i = 3 ! array index for large branches pool + integer, private :: dead_leaves_i = 5 ! array index for dead leaves pool + integer, private :: live_grass_i = 6 ! array index for live grass pool + integer, private :: trunks_i = 4 ! array index for trunks pool + + contains + + procedure :: twigs, small_branches, large_branches, trunks + procedure :: dead_leaves, live_grass + + end type fuel_classes_type + + ! actual type we can pass around + type(fuel_classes_type), public :: fuel_classes + + contains + + integer function twigs(this) + class(fuel_classes_type), intent(in) :: this + twigs = this%twigs_i + end function twigs + + integer function small_branches(this) + class(fuel_classes_type), intent(in) :: this + small_branches = this%small_branches_i + end function small_branches + + integer function large_branches(this) + class(fuel_classes_type), intent(in) :: this + large_branches = this%large_branches_i + end function large_branches + + integer function trunks(this) + class(fuel_classes_type), intent(in) :: this + trunks = this%trunks_i + end function trunks + + integer function dead_leaves(this) + class(fuel_classes_type), intent(in) :: this + dead_leaves = this%dead_leaves_i + end function dead_leaves + + integer function live_grass(this) + class(fuel_classes_type), intent(in) :: this + live_grass = this%live_grass_i + end function live_grass + +end module FatesFuelClassesMod diff --git a/fire/FatesFuelMod.F90 b/fire/FatesFuelMod.F90 new file mode 100644 index 0000000000..8aa8d00d74 --- /dev/null +++ b/fire/FatesFuelMod.F90 @@ -0,0 +1,379 @@ +module FatesFuelMod + + use FatesFuelClassesMod, only : num_fuel_classes, fuel_classes + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : nearzero + use SFNesterovMod, only : nesterov_index + use SFFireWeatherMod, only : fire_weather + use FatesGlobals, only : fates_log + use FatesGlobals, only : endrun => fates_endrun + use shr_log_mod, only : errMsg => shr_log_errMsg + + implicit none + private + + type, public :: fuel_type + + real(r8) :: loading(num_fuel_classes) ! fuel loading of each fuel class [kgC/m2] + real(r8) :: effective_moisture(num_fuel_classes) ! fuel effective moisture all fuel class (moisture/MEF) [m3/m3] + real(r8) :: frac_loading(num_fuel_classes) ! fractional loading of all fuel classes [0-1] + real(r8) :: frac_burnt(num_fuel_classes) ! fraction of litter burnt by fire [0-1] + real(r8) :: non_trunk_loading ! total fuel loading excluding trunks [kgC/m2] + real(r8) :: average_moisture_notrunks ! weighted average of fuel moisture across non-trunk fuel classes [m3/m3] + real(r8) :: bulk_density_notrunks ! weighted average of bulk density across non-trunk fuel classes [kg/m3] + real(r8) :: SAV_notrunks ! weighted average of surface area to volume ratio across non-trunk fuel classes [/cm] + real(r8) :: MEF_notrunks ! weighted average of moisture of extinction across non-trunk fuel classes [m3/m3] + + contains + + procedure :: Init + procedure :: Fuse + procedure :: UpdateLoading + procedure :: SumLoading + procedure :: CalculateFractionalLoading + procedure :: UpdateFuelMoisture + procedure :: AverageBulkDensity_NoTrunks + procedure :: AverageSAV_NoTrunks + + end type fuel_type + + contains + + subroutine Init(this) + ! DESCRIPTION: + ! Initialize fuel class + + ! ARGUMENTS: + class(fuel_type), intent(inout) :: this ! fuel class + + ! just zero everything + this%loading(1:num_fuel_classes) = 0.0_r8 + this%frac_loading(1:num_fuel_classes) = 0.0_r8 + this%frac_burnt(1:num_fuel_classes) = 0.0_r8 + this%effective_moisture(1:num_fuel_classes) = 0.0_r8 + this%non_trunk_loading = 0.0_r8 + this%average_moisture_notrunks = 0.0_r8 + this%bulk_density_notrunks = 0.0_r8 + this%SAV_notrunks = 0.0_r8 + this%MEF_notrunks = 0.0_r8 + + end subroutine Init + + !------------------------------------------------------------------------------------- + + subroutine Fuse(this, self_area, donor_area, donor_fuel) + ! DESCRIPTION: + ! Fuse attributes of this object with another + + ! ARGUMENTS: + class(fuel_type), intent(inout) :: this ! fuel class + real(r8), intent(in) :: self_area ! area of this fuel class's patch [m2] + real(r8), intent(in) :: donor_area ! area of donor fuel class's patch [m2] + type(fuel_type), intent(in) :: donor_fuel ! donor fuel class + + ! LOCALS: + integer :: i ! looping index + real(r8) :: self_weight ! weighting of the receiving fuel class + real(r8) :: donor_weight ! weighting of the donor fuel class + + self_weight = self_area/(donor_area + self_area) + donor_weight = 1.0_r8 - self_weight + + do i = 1, num_fuel_classes + this%loading(i) = this%loading(i)*self_weight + & + donor_fuel%loading(i)*donor_weight + this%frac_loading(i) = this%frac_loading(i)*self_weight + & + donor_fuel%frac_loading(i)*donor_weight + this%frac_burnt(i) = this%frac_burnt(i)*self_weight + & + donor_fuel%frac_burnt(i)*donor_weight + this%effective_moisture(i) = this%effective_moisture(i)*self_weight + & + donor_fuel%effective_moisture(i)*donor_weight + end do + + this%non_trunk_loading = this%non_trunk_loading*self_weight + & + donor_fuel%non_trunk_loading*donor_weight + this%average_moisture_notrunks = this%average_moisture_notrunks*self_weight + & + donor_fuel%average_moisture_notrunks*donor_weight + this%bulk_density_notrunks = this%bulk_density_notrunks*self_weight + & + donor_fuel%bulk_density_notrunks*donor_weight + this%SAV_notrunks = this%SAV_notrunks*self_weight + donor_fuel%SAV_notrunks*donor_weight + this%MEF_notrunks = this%MEF_notrunks*self_weight + donor_fuel%MEF_notrunks*donor_weight + + end subroutine Fuse + + !------------------------------------------------------------------------------------- + + subroutine UpdateLoading(this, leaf_litter, twig_litter, small_branch_litter, & + large_branch_litter, trunk_litter, live_grass) + ! DESCRIPTION: + ! Updates loading for each fuel type + + ! ARGUMENTS: + class(fuel_type), intent(inout) :: this ! fuel class + real(r8), intent(in) :: leaf_litter ! input leaf litter [kgC/m2] + real(r8), intent(in) :: twig_litter ! input twig litter [kgC/m2] + real(r8), intent(in) :: small_branch_litter ! input small branch litter [kgC/m2] + real(r8), intent(in) :: large_branch_litter ! input leaf litter [kgC/m2] + real(r8), intent(in) :: trunk_litter ! input leaf litter [kgC/m2] + real(r8), intent(in) :: live_grass ! input live grass [kgC/m2] + + this%loading(fuel_classes%dead_leaves()) = leaf_litter + this%loading(fuel_classes%twigs()) = twig_litter + this%loading(fuel_classes%small_branches()) = small_branch_litter + this%loading(fuel_classes%large_branches()) = large_branch_litter + this%loading(fuel_classes%live_grass()) = live_grass + this%loading(fuel_classes%trunks()) = trunk_litter + + end subroutine UpdateLoading + + !------------------------------------------------------------------------------------- + + subroutine SumLoading(this) + ! DESCRIPTION: + ! Sums up the loading - excludes trunks + ! + ! Only the 1-h, 10-h and 100-h fuel classes influence fire spread + ! Rothermel, 1972 (USDA FS GTR INT-115) + ! Wilson, 1982 (UTINT-289) + ! Pyne et al., 1996 (Introduction to wildland fire) + + ! ARGUMENTS: + class(fuel_type), intent(inout) :: this ! fuel class + + ! LOCALS: + integer :: i ! looping index + + this%non_trunk_loading = 0.0_r8 + do i = 1, num_fuel_classes + if (i /= fuel_classes%trunks()) then + this%non_trunk_loading = this%non_trunk_loading + this%loading(i) + end if + end do + + end subroutine SumLoading + + !------------------------------------------------------------------------------------- + + subroutine CalculateFractionalLoading(this) + ! DESCRIPTION: + ! Calculates fractional loading for fuel + + ! ARGUMENTS: + class(fuel_type), intent(inout) :: this ! fuel class + + ! LOCALS: + integer :: i ! looping index + + ! sum up loading just in case + call this%SumLoading() + + if (this%non_trunk_loading > nearzero) then + do i = 1, num_fuel_classes + if (i /= fuel_classes%trunks()) then + this%frac_loading(i) = this%loading(i)/this%non_trunk_loading + else + this%frac_loading(i) = 0.0_r8 + end if + end do + else + this%frac_loading(1:num_fuel_classes) = 0.0_r8 + this%non_trunk_loading = 0.0_r8 + end if + + end subroutine CalculateFractionalLoading + + !------------------------------------------------------------------------------------- + + subroutine UpdateFuelMoisture(this, sav_fuel, drying_ratio, fireWeatherClass) + ! DESCRIPTION: + ! Updates fuel moisture depending on what fire weather class is in use + + ! ARGUMENTS: + class(fuel_type), intent(inout) :: this ! fuel class + real(r8), intent(in) :: sav_fuel(num_fuel_classes) ! surface area to volume ratio of all fuel types [/cm] + real(r8), intent(in) :: drying_ratio ! drying ratio + class(fire_weather), intent(in) :: fireWeatherClass ! fireWeatherClass + + real(r8) :: moisture(num_fuel_classes) ! fuel moisture [m3/m3] + real(r8) :: moisture_of_extinction(num_fuel_classes) ! fuel moisture of extinction [m3/m3] + integer :: i ! looping index + + if (this%non_trunk_loading + this%loading(fuel_classes%trunks()) > nearzero) then + ! calculate fuel moisture [m3/m3] for each fuel class depending on what + ! fire weather class is in use + select type (fireWeatherClass) + class is (nesterov_index) + call CalculateFuelMoistureNesterov(sav_fuel, drying_ratio, & + fireWeatherClass%fire_weather_index, moisture) + class default + write(fates_log(), *) 'Unknown fire weather class selected.' + write(fates_log(), *) 'Choose a different fire weather class or upate this subroutine.' + call endrun(msg=errMsg( __FILE__, __LINE__)) + end select + + this%average_moisture_notrunks = 0.0_r8 + this%MEF_notrunks = 0.0_r8 + do i = 1, num_fuel_classes + ! calculate moisture of extinction and fuel effective moisture + moisture_of_extinction(i) = MoistureOfExtinction(sav_fuel(i)) + this%effective_moisture(i) = moisture(i)/moisture_of_extinction(i) + + ! average fuel moisture and MEF across all fuel types except trunks [m3/m3] + if (i /= fuel_classes%trunks()) then + this%average_moisture_notrunks = this%average_moisture_notrunks + this%frac_loading(i)*moisture(i) + this%MEF_notrunks = this%MEF_notrunks + this%frac_loading(i)*moisture_of_extinction(i) + end if + end do + + else + this%effective_moisture(1:num_fuel_classes) = 0.0_r8 + this%average_moisture_notrunks = 0.0_r8 + this%MEF_notrunks = 0.0_r8 + end if + + end subroutine UpdateFuelMoisture + + !------------------------------------------------------------------------------------- + + subroutine CalculateFuelMoistureNesterov(sav_fuel, drying_ratio, NI, moisture) + ! + ! DESCRIPTION: + ! Updates fuel moisture + + ! ARGUMENTS: + real(r8), intent(in) :: sav_fuel(num_fuel_classes) ! surface area to volume ratio of all fuel types [/cm] + real(r8), intent(in) :: drying_ratio ! drying ratio + real(r8), intent(in) :: NI ! Nesterov Index + real(r8), intent(out) :: moisture(num_fuel_classes) ! moisture of litter [m3/m3] + + ! LOCALS + integer :: i ! looping index + real(r8) :: alpha_FMC ! intermediate variable for calculating fuel moisture + + do i = 1, num_fuel_classes + if (i == fuel_classes%live_grass()) then + ! 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 + alpha_FMC = sav_fuel(fuel_classes%twigs())/drying_ratio + else + alpha_FMC = sav_fuel(i)/drying_ratio + end if + moisture(i) = exp(-1.0_r8*alpha_FMC*NI) + end do + + end subroutine CalculateFuelMoistureNesterov + + !------------------------------------------------------------------------------------- + + real(r8) function MoistureOfExtinction(sav) + ! + ! DESCRIPTION: + ! Calculates moisture of extinction based on input surface area to volume ratio + + ! MEF (moisure of extinction) depends on compactness of fuel, depth, particle size, + ! wind, and slope + ! Equation here is Eq. 27 from Peterson and Ryan (1986) "Modeling Postfire Conifer + ! Mortality for Long-Range Planning" + ! + ! Example MEFs: + ! pine needles = 0.30 (Rothermel 1972) + ! short grass = 0.12 (Rothermel 1983; Gen. Tech. Rep. INT-143; Table II-1) + ! tall grass = 0.24 (Rothermel 1983) + ! chaparral = 0.20 (Rothermel 1983) + ! closed timber litter = 0.30 (Rothermel 1983) + ! hardwood litter = 0.25 (Rothermel 1983) + ! grass = 0.2 (Lasslop 2014; Table 1) + ! shrubs = 0.3 (Lasslop 2014; Table 1) + ! tropical evergreen trees = 0.2 (Lasslop 2014; Table 1) + ! tropical deciduous trees = 0.3 (Lasslop 2014; Table 1) + ! extratropical trees = 0.3 (Lasslop 2014; Table 1) + ! + ! SAV values from Thonicke 2010 give: + ! twigs = 0.355, small branches = 0.44, large branches = 0.525, trunks = 0.63 + ! dead leaves = 0.248, live grass = 0.248 + ! + + ! ARGUMENTS: + real(r8), intent(in) :: sav ! fuel surface area to volume ratio [/cm] + + ! CONSTANTS: + real(r8), parameter :: MEF_a = 0.524_r8 + real(r8), parameter :: MEF_b = 0.066_r8 + + if (sav <= 0.0_r8) then + write(fates_log(), *) 'SAV cannot be negative - SAV' + call endrun(msg=errMsg(__FILE__, __LINE__)) + else + MoistureOfExtinction = MEF_a - MEF_b*log(sav) + end if + + end function MoistureOfExtinction + + !------------------------------------------------------------------------------------- + + subroutine AverageBulkDensity_NoTrunks(this, bulk_density) + ! DESCRIPTION: + ! Calculates average bulk density (not including trunks) + ! + ! Only the 1-h, 10-h and 100-h fuel classes influence fire spread + ! Rothermel, 1972 (USDA FS GTR INT-115) + ! Wilson, 1982 (UTINT-289) + ! Pyne et al., 1996 (Introduction to wildland fire) + + ! ARGUMENTS: + class(fuel_type), intent(inout) :: this ! fuel class + real(r8), intent(in) :: bulk_density(num_fuel_classes) ! bulk density of all fuel types [kg/m2] + + ! LOCALS: + integer :: i ! looping index + + if (this%non_trunk_loading > nearzero) then + this%bulk_density_notrunks = 0.0_r8 + do i = 1, num_fuel_classes + ! average bulk density across all fuel types except trunks + if (i /= fuel_classes%trunks()) then + this%bulk_density_notrunks = this%bulk_density_notrunks + this%frac_loading(i)*bulk_density(i) + end if + end do + else + this%bulk_density_notrunks = sum(bulk_density(1:num_fuel_classes))/num_fuel_classes + end if + + end subroutine AverageBulkDensity_NoTrunks + + !------------------------------------------------------------------------------------- + + subroutine AverageSAV_NoTrunks(this, sav_fuel) + ! DESCRIPTION: + ! Calculates average surface area to volume ratio (not including trunks) + ! + ! Only the 1-h, 10-h and 100-h fuel classes influence fire spread + ! Rothermel, 1972 (USDA FS GTR INT-115) + ! Wilson, 1982 (UTINT-289) + ! Pyne et al., 1996 (Introduction to wildland fire) + + ! ARGUMENTS: + class(fuel_type), intent(inout) :: this ! fuel class + real(r8), intent(in) :: sav_fuel(num_fuel_classes) ! surface area to volume ratio of all fuel types [/cm] + + ! LOCALS: + integer :: i ! looping index + + if (this%non_trunk_loading > nearzero) then + this%SAV_notrunks = 0.0_r8 + do i = 1, num_fuel_classes + ! average bulk density across all fuel types except trunks + if (i /= fuel_classes%trunks()) then + this%SAV_notrunks = this%SAV_notrunks + this%frac_loading(i)*sav_fuel(i) + end if + end do + else + this%SAV_notrunks = sum(sav_fuel(1:num_fuel_classes))/num_fuel_classes + end if + + end subroutine AverageSAV_NoTrunks + + !--------------------------------------------------------------------------------------- + +end module FatesFuelMod diff --git a/fire/SFFireWeatherMod.F90 b/fire/SFFireWeatherMod.F90 index 11b8602fa4..3191b460b1 100644 --- a/fire/SFFireWeatherMod.F90 +++ b/fire/SFFireWeatherMod.F90 @@ -40,7 +40,8 @@ subroutine update_fire_weather(this, temp_C, precip, rh, wind) real(r8), intent(in) :: wind end subroutine update_fire_weather - end interface + + end interface contains diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index eeb37e79f7..97bbcb6d54 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -1,54 +1,44 @@ - module SFMainMod +module SFMainMod ! ============================================================================ ! All subroutines related to the SPITFIRE fire routine. ! Code originally developed by Allan Spessa & Rosie Fisher as part of the NERC-QUEST project. ! ============================================================================ - use FatesConstantsMod , only : r8 => fates_r8 - use FatesConstantsMod , only : itrue, ifalse - use FatesConstantsMod , only : pi_const - use FatesConstantsMod , only : nocomp_bareground - use FatesInterfaceTypesMod, only : hlm_masterproc ! 1= master process, 0=not master process - use FatesGlobals , only : fates_log + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : itrue, ifalse + use FatesConstantsMod, only : pi_const + use FatesConstantsMod, only : nocomp_bareground, nearzero + use FatesGlobals, only : fates_log + use FatesInterfaceTypesMod, only : hlm_masterproc use FatesInterfaceTypesMod, only : hlm_spitfire_mode use FatesInterfaceTypesMod, only : hlm_sf_nofire_def use FatesInterfaceTypesMod, only : hlm_sf_scalar_lightning_def use FatesInterfaceTypesMod, only : hlm_sf_successful_ignitions_def use FatesInterfaceTypesMod, only : hlm_sf_anthro_ignitions_def use FatesInterfaceTypesMod, only : bc_in_type - - use EDPftvarcon , only : EDPftvarcon_inst - use PRTParametersMod , only : prt_params - - use PRTGenericMod , only : element_pos - use EDtypesMod , only : ed_site_type - use FatesPatchMod , only : fates_patch_type - use FatesCohortMod , only : fates_cohort_type - use EDtypesMod , only : AREA - use FatesLitterMod , only : DL_SF - use FatesLitterMod , only : TW_SF - use FatesLitterMod , only : LB_SF - use FatesLitterMod , only : LG_SF - use FatesLitterMod , only : ncwd - use FatesLitterMod , only : NFSC - use FatesLitterMod , only : TR_SF - use FatesLitterMod , only : litter_type - + use EDPftvarcon, only : EDPftvarcon_inst + use PRTParametersMod, only : prt_params + use PRTGenericMod, only : element_pos + use EDtypesMod, only : ed_site_type + use FatesPatchMod, only : fates_patch_type + use FatesCohortMod, only : fates_cohort_type + use EDtypesMod, only : AREA + use FatesLitterMod, only : litter_type + use FatesFuelClassesMod, only : num_fuel_classes use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : carbon12_element - use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : sapw_organ use PRTGenericMod, only : struct_organ use FatesInterfaceTypesMod, only : numpft use FatesAllometryMod, only : CrownDepth - use FatesConstantsMod, only : nearzero - + use FatesFuelClassesMod, only : fuel_classes + implicit none private public :: fire_model - public :: charecteristics_of_fuel + public :: UpdateFuelCharacteristics public :: rate_of_spread public :: ground_fuel_consumption public :: area_burnt_intensity @@ -57,9 +47,6 @@ module SFMainMod public :: cambial_damage_kill public :: post_fire_mortality - ! 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 integer :: write_SF = ifalse ! for debugging logical :: debug = .false. ! for debugging @@ -75,7 +62,7 @@ subroutine fire_model(currentSite, bc_in) ! ARGUMENTS: type(ed_site_type), intent(inout), target :: currentSite ! site object type(bc_in_type), intent(in) :: bc_in ! BC in object - + ! LOCALS: type (fates_patch_type), pointer :: currentPatch ! patch object @@ -86,10 +73,10 @@ subroutine fire_model(currentSite, bc_in) currentPatch%fire = 0 currentPatch => currentPatch%older end do - + if (hlm_spitfire_mode > hlm_sf_nofire_def) then call UpdateFireWeather(currentSite, bc_in) - call charecteristics_of_fuel(currentSite) + call UpdateFuelCharacteristics(currentSite) call rate_of_spread(currentSite) call ground_fuel_consumption(currentSite) call area_burnt_intensity(currentSite, bc_in) @@ -108,6 +95,7 @@ subroutine UpdateFireWeather(currentSite, bc_in) ! DESCRIPTION: ! Updates the site's fire weather index and calculates effective windspeed based on ! vegetation characteristics + ! ! Currently we use tree and grass fraction averaged over whole grid (site) to ! prevent extreme divergence @@ -115,7 +103,6 @@ subroutine UpdateFireWeather(currentSite, bc_in) use FatesConstantsMod, only : sec_per_day, sec_per_min use EDTypesMod, only : CalculateTreeGrassAreaSite - ! ARGUMENTS: type(ed_site_type), intent(inout), target :: currentSite type(bc_in_type), intent(in) :: bc_in @@ -130,12 +117,12 @@ subroutine UpdateFireWeather(currentSite, bc_in) real(r8) :: grass_fraction ! site-level grass fraction [0-1] real(r8) :: bare_fraction ! site-level bare ground fraction [0-1] 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 ! is simply using the values associated with the first patch. ! which probably won't have much impact, unless we decide to ever calculated fire weather for each patch. - + currentPatch => currentSite%oldest_patch ! If the oldest patch is a bareground patch (i.e. nocomp mode is on) use the first vegetated patch @@ -151,7 +138,7 @@ subroutine UpdateFireWeather(currentSite, bc_in) wind = bc_in%wind24_pa(iofp) ! convert to m/min - currentSite%wind = wind*sec_per_min + currentSite%wind = wind*sec_per_min ! update fire weather index call currentSite%fireWeather%UpdateIndex(temp_C, precip, rh, wind) @@ -167,191 +154,57 @@ end subroutine UpdateFireWeather !--------------------------------------------------------------------------------------- - subroutine charecteristics_of_fuel ( currentSite ) - - use SFParamsMod, only : SF_val_drying_ratio, SF_val_SAV, SF_val_FBD - - type(ed_site_type), intent(in), target :: currentSite + subroutine UpdateFuelCharacteristics(currentSite) + ! + ! DESCRIPTION: + ! Updates fuel characteristics on each patch of the site + ! - type(fates_patch_type), pointer :: currentPatch - type(fates_cohort_type), pointer :: currentCohort - type(litter_type), pointer :: litt_c + use SFParamsMod, only : SF_val_drying_ratio, SF_val_SAV, SF_val_FBD - 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 + ! ARGUMENTS: + type(ed_site_type), intent(in), target :: currentSite ! site object - fuel_moisture(:) = 0.0_r8 + ! LOCALS: + type(fates_patch_type), pointer :: currentPatch ! FATES patch + type(litter_type), pointer :: litter ! pointer to patch litter class + real(r8) :: MEF_trunks, fuel_moisture_trunks - currentPatch => currentSite%oldest_patch; + currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - - litt_c => currentPatch%litter(element_pos(carbon12_element)) - - ! How much live grass is there? - currentPatch%livegrass = 0.0_r8 - currentCohort => currentPatch%tallest - do while(associated(currentCohort)) - ! for grasses sum all aboveground tissues - if( prt_params%woody(currentCohort%pft) == ifalse)then - - currentPatch%livegrass = currentPatch%livegrass + & - ( currentCohort%prt%GetState(leaf_organ, carbon12_element) + & - currentCohort%prt%GetState(sapw_organ, carbon12_element) + & - currentCohort%prt%GetState(struct_organ, carbon12_element) ) * & - currentCohort%n/currentPatch%area + if (currentPatch%nocomp_pft_label /= nocomp_bareground) then - endif - currentCohort => currentCohort%shorter - enddo - - ! There are SIX fuel classes - ! 1:4) four CWD_AG pools (twig, s branch, l branch, trunk), 5) dead leaves and 6) live grass - ! NCWD =4 NFSC = 6 - ! tw_sf = 1, lb_sf = 3, tr_sf = 4, dl_sf = 5, lg_sf = 6, - + ! calculate live grass [kgC/m2] + call currentPatch%UpdateLiveGrass() - if(write_sf == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) ' leaf_litter1 ',sum(litt_c%leaf_fines(:)) - if ( hlm_masterproc == itrue ) write(fates_log(),*) ' leaf_litter2 ',sum(litt_c%ag_cwd(:)) - if ( hlm_masterproc == itrue ) write(fates_log(),*) ' leaf_litter3 ',currentPatch%livegrass - endif - - currentPatch%sum_fuel = sum(litt_c%leaf_fines(:)) + & - sum(litt_c%ag_cwd(:)) + & - currentPatch%livegrass - if(write_SF == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sum fuel', currentPatch%sum_fuel,currentPatch%area - endif - ! =============================================== - ! Average moisture, bulk density, surface area-volume and moisture extinction of fuel - ! ================================================ - - if (currentPatch%sum_fuel > 0.0) then - ! Fraction of fuel in litter classes - currentPatch%fuel_frac(dl_sf) = sum(litt_c%leaf_fines(:))/ currentPatch%sum_fuel - currentPatch%fuel_frac(tw_sf:tr_sf) = litt_c%ag_cwd(:) / currentPatch%sum_fuel - - if(write_sf == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff2a ', & - lg_sf,currentPatch%livegrass,currentPatch%sum_fuel - endif - - currentPatch%fuel_frac(lg_sf) = currentPatch%livegrass / currentPatch%sum_fuel - - ! MEF (moisure of extinction) depends on compactness of fuel, depth, particle size, wind, slope - ! Eqn here is eqn 27 from Peterson and Ryan (1986) "Modeling Postfire Conifer Mortality for Long-Range Planning" - ! but lots of other approaches in use out there... - ! MEF: pine needles=0.30 (text near EQ 28 Rothermal 1972) - ! Table II-1 NFFL mixed fuels models from Rothermal 1983 Gen. Tech. Rep. INT-143 - ! MEF: short grass=0.12,tall grass=0.25,chaparral=0.20,closed timber litter=0.30,hardwood litter=0.25 - ! Thonicke 2010 SAV values propagated thru P&R86 eqn below gives MEF:tw=0.355, sb=0.44, lb=0.525, tr=0.63, dg=0.248, lg=0.248 - ! Lasslop 2014 Table 1 MEF PFT level:grass=0.2,shrubs=0.3,TropEverGrnTree=0.2,TropDecid Tree=0.3, Extra-trop Tree=0.3 - MEF(1:nfsc) = 0.524_r8 - 0.066_r8 * log(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 + ! update fuel loading [kgC/m2] + litter => currentPatch%litter(element_pos(carbon12_element)) + call currentPatch%fuel%UpdateLoading(sum(litter%leaf_fines(:)), & + litter%ag_cwd(1), litter%ag_cwd(2), litter%ag_cwd(3), litter%ag_cwd(4), & + currentPatch%livegrass) + + ! sum up fuel classes and calculate fractional loading for each + call currentPatch%fuel%SumLoading() + call currentPatch%fuel%CalculateFractionalLoading() - fuel_moisture(tw_sf:dl_sf) = exp(-1.0_r8 * alpha_FMC(tw_sf:dl_sf) * & - currentSite%fireWeather%fire_weather_index) - - if(write_SF == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff3 ',currentPatch%fuel_frac - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'fm ',fuel_moisture - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'csa ',currentSite%fireWeather%fire_weather_index - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sfv ',alpha_FMC - endif - - ! 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%fireWeather%fire_weather_index)) - - ! 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)) - currentPatch%fuel_sav = sum(currentPatch%fuel_frac(tw_sf:lb_sf) * SF_val_SAV(tw_sf:lb_sf)) - currentPatch%fuel_mef = sum(currentPatch%fuel_frac(tw_sf:lb_sf) * MEF(tw_sf:lb_sf)) - currentPatch%fuel_eff_moist = sum(currentPatch%fuel_frac(tw_sf:lb_sf) * fuel_moisture(tw_sf:lb_sf)) - if(write_sf == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff4 ',currentPatch%fuel_eff_moist - endif - ! Add on properties of dead leaves and live grass pools (5 & 6) - currentPatch%fuel_bulkd = currentPatch%fuel_bulkd + sum(currentPatch%fuel_frac(dl_sf:lg_sf) * SF_val_FBD(dl_sf:lg_sf)) - currentPatch%fuel_sav = currentPatch%fuel_sav + sum(currentPatch%fuel_frac(dl_sf:lg_sf) * SF_val_SAV(dl_sf:lg_sf)) - currentPatch%fuel_mef = currentPatch%fuel_mef + sum(currentPatch%fuel_frac(dl_sf:lg_sf) * MEF(dl_sf:lg_sf)) - currentPatch%fuel_eff_moist = currentPatch%fuel_eff_moist+ sum(currentPatch%fuel_frac(dl_sf:lg_sf) * fuel_moisture(dl_sf:lg_sf)) - - ! Correct averaging for the fact that we are not using the trunks pool for fire ROS and intensity (5) - ! Consumption of fuel in trunk pool does not influence fire ROS or intensity (Pyne 1996) - if ( (1.0_r8-currentPatch%fuel_frac(tr_sf)) .gt. nearzero ) then - currentPatch%fuel_bulkd = currentPatch%fuel_bulkd * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf))) - currentPatch%fuel_sav = currentPatch%fuel_sav * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf))) - currentPatch%fuel_mef = currentPatch%fuel_mef * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf))) - currentPatch%fuel_eff_moist = currentPatch%fuel_eff_moist * (1.0_r8/(1.0_r8-currentPatch%fuel_frac(tr_sf))) - else - ! somehow the fuel is all trunk. put dummy values from large branches so as not to break things later in code. - currentPatch%fuel_bulkd = SF_val_FBD(lb_sf) - currentPatch%fuel_sav = SF_val_SAV(lb_sf) - currentPatch%fuel_mef = MEF(lb_sf) - currentPatch%fuel_eff_moist = fuel_moisture(lb_sf) - endif - - ! Pass litter moisture into the fuel burning routine (all fuels: twigs,s branch,l branch,trunk,dead leaves,live grass) - ! (wo/me term in Thonicke et al. 2010) - currentPatch%litter_moisture(tw_sf:lb_sf) = fuel_moisture(tw_sf:lb_sf)/MEF(tw_sf:lb_sf) - currentPatch%litter_moisture(tr_sf) = fuel_moisture(tr_sf)/MEF(tr_sf) - currentPatch%litter_moisture(dl_sf) = fuel_moisture(dl_sf)/MEF(dl_sf) - currentPatch%litter_moisture(lg_sf) = fuel_moisture(lg_sf)/MEF(lg_sf) - - ! remove trunks from patch%sum_fuel because they should not be included in fire equations - ! NOTE: ACF will update this soon to be more clean/bug-proof - currentPatch%sum_fuel = currentPatch%sum_fuel - litt_c%ag_cwd(tr_sf) - - else - - if(write_SF == itrue)then - - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'no litter fuel at all',currentPatch%patchno, & - currentPatch%sum_fuel,sum(litt_c%ag_cwd(:)),sum(litt_c%leaf_fines(:)) - - endif - currentPatch%fuel_sav = sum(SF_val_SAV(1:nfsc))/(nfsc) ! make average sav to avoid crashing code. - - if ( hlm_masterproc == itrue .and. write_SF == itrue)then - write(fates_log(),*) 'problem with spitfire fuel averaging' - end if - - ! FIX(SPM,032414) refactor...should not have 0 fuel unless everything is burnt - ! off. - currentPatch%fuel_eff_moist = 0.0000000001_r8 - currentPatch%fuel_bulkd = 0.0000000001_r8 - currentPatch%fuel_frac(:) = 0.0000000001_r8 - currentPatch%fuel_mef = 0.0000000001_r8 - currentPatch%sum_fuel = 0.0000000001_r8 - - endif - ! check values. - ! FIX(SPM,032414) refactor... - if(write_SF == itrue.and.currentPatch%fuel_sav <= 0.0_r8.or.currentPatch%fuel_bulkd <= & - 0.0_r8.or.currentPatch%fuel_mef <= 0.0_r8.or.currentPatch%fuel_eff_moist <= 0.0_r8)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'problem with spitfire fuel averaging' - endif - endif !nocomp_pft_label check - currentPatch => currentPatch%younger - - enddo !end patch loop - - end subroutine charecteristics_of_fuel + ! calculate fuel moisture [m3/m3] + call currentPatch%fuel%UpdateFuelMoisture(SF_val_SAV, SF_val_drying_ratio, & + currentSite%fireWeather) + + ! calculate geometric properties + call currentPatch%fuel%AverageBulkDensity_NoTrunks(SF_val_FBD) + call currentPatch%fuel%AverageSAV_NoTrunks(SF_val_SAV) + + end if + currentPatch => currentPatch%younger + end do + end subroutine UpdateFuelCharacteristics + !--------------------------------------------------------------------------------------- - subroutine rate_of_spread ( currentSite ) + subroutine rate_of_spread (currentSite) !*****************************************************************. !Routine called daily from within ED within a site loop. !Returns the updated currentPatch%ROS_front value for each patch. @@ -360,12 +213,12 @@ subroutine rate_of_spread ( currentSite ) SF_val_part_dens, & SF_val_miner_damp, & SF_val_fuel_energy - + use FatesConstantsMod, only : nearzero type(ed_site_type), intent(in), target :: currentSite type(fates_patch_type), pointer :: currentPatch - ! Rothermal fire spread model parameters. + ! Rothermel fire spread model parameters. real(r8) beta,beta_op ! weighted average of packing ratio (unitless) real(r8) ir ! reaction intensity (kJ/m2/min) real(r8) xi,eps,phi_wind ! all are unitless @@ -383,50 +236,58 @@ subroutine rate_of_spread ( currentSite ) do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - + if(currentPatch%nocomp_pft_label .ne. nocomp_bareground .and. currentPatch%fuel%non_trunk_loading > nearzero) then + ! remove mineral content from net fuel load per Thonicke 2010 for ir calculation - currentPatch%sum_fuel = currentPatch%sum_fuel * (1.0_r8 - SF_val_miner_total) !net of minerals + currentPatch%fuel%non_trunk_loading = currentPatch%fuel%non_trunk_loading * (1.0_r8 - SF_val_miner_total) !net of minerals ! ----start spreading--- if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) & - 'SF - currentPatch%fuel_bulkd ',currentPatch%fuel_bulkd + 'SF - currentPatch%fuel%bulk_density_notrunks',currentPatch%fuel%bulk_density_notrunks if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) & 'SF - SF_val_part_dens ',SF_val_part_dens ! beta = packing ratio (unitless) - ! fraction of fuel array volume occupied by fuel or compactness of fuel bed - beta = currentPatch%fuel_bulkd / SF_val_part_dens + ! fraction of fuel array volume occupied by fuel or compactness of fuel bed + beta = currentPatch%fuel%bulk_density_notrunks/SF_val_part_dens ! Equation A6 in Thonicke et al. 2010 - ! packing ratio (unitless) - beta_op = 0.200395_r8 *(currentPatch%fuel_sav**(-0.8189_r8)) + ! packing ratio (unitless) + if (currentPatch%fuel%SAV_notrunks < nearzero) then + beta_op = 0.0_r8 + else + beta_op = 0.200395_r8 *(currentPatch%fuel%SAV_notrunks**(-0.8189_r8)) + end if if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) 'SF - beta ',beta if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) 'SF - beta_op ',beta_op - beta_ratio = beta/beta_op !unitless + if (beta_op < nearzero) then + beta_ratio = 0.0_r8 + else + beta_ratio = beta/beta_op !unitless + end if if(write_sf == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'esf ',currentPatch%fuel_eff_moist + if ( hlm_masterproc == itrue ) write(fates_log(),*) 'average moisture',currentPatch%fuel%average_moisture_notrunks endif ! ---heat of pre-ignition--- ! Equation A4 in Thonicke et al. 2010 - ! Rothermal EQ12= 250 Btu/lb + 1116 Btu/lb * fuel_eff_moist - ! conversion of Rothermal (1972) EQ12 in BTU/lb to current kJ/kg + ! Rothermel EQ12= 250 Btu/lb + 1116 Btu/lb * average_moisture + ! conversion of Rothermel (1972) EQ12 in BTU/lb to current kJ/kg ! q_ig in kJ/kg - q_ig = q_dry +2594.0_r8 * currentPatch%fuel_eff_moist + q_ig = q_dry +2594.0_r8 * currentPatch%fuel%average_moisture_notrunks ! ---effective heating number--- ! Equation A3 in Thonicke et al. 2010. - eps = exp(-4.528_r8 / currentPatch%fuel_sav) + eps = exp(-4.528_r8 / currentPatch%fuel%SAV_notrunks) ! Equation A7 in Thonicke et al. 2010 per eqn 49 from Rothermel 1972 - b = 0.15988_r8 * (currentPatch%fuel_sav**0.54_r8) + b = 0.15988_r8 * (currentPatch%fuel%SAV_notrunks**0.54_r8) ! Equation A8 in Thonicke et al. 2010 per eqn 48 from Rothermel 1972 - c = 7.47_r8 * (exp(-0.8711_r8 * (currentPatch%fuel_sav**0.55_r8))) + c = 7.47_r8 * (exp(-0.8711_r8 * (currentPatch%fuel%SAV_notrunks**0.55_r8))) ! Equation A9 in Thonicke et al. 2010. (appears to have typo, using coefficient eqn.50 Rothermel 1972) - e = 0.715_r8 * (exp(-0.01094_r8 * currentPatch%fuel_sav)) + e = 0.715_r8 * (exp(-0.01094_r8 * currentPatch%fuel%SAV_notrunks)) if (debug) then if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) 'SF - c ',c @@ -442,26 +303,26 @@ subroutine rate_of_spread ( currentSite ) ! ---propagating flux---- - ! Equation A2 in Thonicke et al.2010 and Eq. 42 Rothermal 1972 + ! Equation A2 in Thonicke et al.2010 and Eq. 42 Rothermel 1972 ! xi (unitless) - xi = (exp((0.792_r8 + 3.7597_r8 * (currentPatch%fuel_sav**0.5_r8)) * (beta+0.1_r8))) / & - (192_r8+7.9095_r8 * currentPatch%fuel_sav) + xi = (exp((0.792_r8 + 3.7597_r8 * (currentPatch%fuel%SAV_notrunks**0.5_r8)) * (beta+0.1_r8))) / & + (192_r8+7.9095_r8 * currentPatch%fuel%SAV_notrunks) ! ---reaction intensity---- ! Equation in table A1 Thonicke et al. 2010. - a = 8.9033_r8 * (currentPatch%fuel_sav**(-0.7913_r8)) + a = 8.9033_r8 * (currentPatch%fuel%SAV_notrunks**(-0.7913_r8)) a_beta = exp(a*(1.0_r8-beta_ratio)) !dummy variable for reaction_v_opt equation ! Equation in table A1 Thonicke et al. 2010. ! reaction_v_max and reaction_v_opt = reaction velocity in units of per min - ! reaction_v_max = Equation 36 in Rothermal 1972 and Fig 12 - reaction_v_max = 1.0_r8 / (0.0591_r8 + 2.926_r8* (currentPatch%fuel_sav**(-1.5_r8))) - ! reaction_v_opt = Equation 38 in Rothermal 1972 and Fig 11 + ! reaction_v_max = Equation 36 in Rothermel 1972 and Fig 12 + reaction_v_max = 1.0_r8 / (0.0591_r8 + 2.926_r8* (currentPatch%fuel%SAV_notrunks**(-1.5_r8))) + ! reaction_v_opt = Equation 38 in Rothermel 1972 and Fig 11 reaction_v_opt = reaction_v_max*(beta_ratio**a)*a_beta ! mw_weight = relative fuel moisture/fuel moisture of extinction ! average values for litter pools (dead leaves, twigs, small and large branches) plus grass - mw_weight = currentPatch%fuel_eff_moist/currentPatch%fuel_mef + mw_weight = currentPatch%fuel%average_moisture_notrunks/currentPatch%fuel%MEF_notrunks ! Equation in table A1 Thonicke et al. 2010. ! moist_damp is unitless @@ -469,23 +330,22 @@ subroutine rate_of_spread ( currentSite ) (3.52_r8*(mw_weight**3.0_r8)))) ! ir = reaction intenisty in kJ/m2/min - ! currentPatch%sum_fuel converted from kgC/m2 to kgBiomass/m2 for ir calculation - ir = reaction_v_opt*(currentPatch%sum_fuel/0.45_r8)*SF_val_fuel_energy*moist_damp*SF_val_miner_damp + ! currentPatch%fuel%non_trunk_loading converted from kgC/m2 to kgBiomass/m2 for ir calculation + ir = reaction_v_opt*(currentPatch%fuel%non_trunk_loading/0.45_r8)*SF_val_fuel_energy*moist_damp*SF_val_miner_damp ! write(fates_log(),*) 'ir',gamma_aptr,moist_damp,SF_val_fuel_energy,SF_val_miner_damp - if (((currentPatch%fuel_bulkd) <= 0.0_r8).or.(eps <= 0.0_r8).or.(q_ig <= 0.0_r8)) then + if (((currentPatch%fuel%bulk_density_notrunks) <= 0.0_r8).or.(eps <= 0.0_r8).or.(q_ig <= 0.0_r8)) then currentPatch%ROS_front = 0.0_r8 else ! Equation 9. Thonicke et al. 2010. ! forward ROS in m/min - currentPatch%ROS_front = (ir*xi*(1.0_r8+phi_wind)) / (currentPatch%fuel_bulkd*eps*q_ig) - ! write(fates_log(),*) 'ros calcs',currentPatch%fuel_bulkd,ir,xi,eps,q_ig + currentPatch%ROS_front = (ir*xi*(1.0_r8+phi_wind)) / (currentPatch%fuel%bulk_density_notrunks*eps*q_ig) endif ! Equation 10 in Thonicke et al. 2010 ! backward ROS from Can FBP System (1992) in m/min ! backward ROS wind not changed by vegetation currentPatch%ROS_back = currentPatch%ROS_front*exp(-0.012_r8*currentSite%wind) - + end if ! nocomp_pft_label check currentPatch => currentPatch%younger @@ -497,66 +357,72 @@ end subroutine rate_of_spread subroutine ground_fuel_consumption ( currentSite ) !***************************************************************** !returns the the hypothetic fuel consumed by the fire - - use SFParamsMod, only : SF_val_miner_total, SF_val_min_moisture, & - SF_val_mid_moisture, SF_val_low_moisture_Coeff, SF_val_low_moisture_Slope, & - SF_val_mid_moisture_Coeff, SF_val_mid_moisture_Slope + use SFParamsMod, only: SF_val_mid_moisture, SF_val_mid_moisture_Coeff, SF_val_mid_moisture_Slope + use SFParamsMod, only : SF_val_min_moisture, SF_val_low_moisture_Coeff, SF_val_low_moisture_Slope + use SFParamsMod, only : SF_val_miner_total type(ed_site_type) , intent(in), target :: currentSite type(fates_patch_type), pointer :: currentPatch type(litter_type), pointer :: litt_c ! carbon 12 litter pool real(r8) :: moist !effective fuel moisture - real(r8) :: tau_b(nfsc) !lethal heating rates for each fuel class (min) - real(r8) :: fc_ground(nfsc) !total amount of fuel consumed per area of burned ground (kg C / m2 of burned area) - + real(r8) :: tau_b(num_fuel_classes) !lethal heating rates for each fuel class (min) + real(r8) :: fc_ground(num_fuel_classes) !total amount of fuel consumed per area of burned ground (kg C / m2 of burned area) + integer :: tr_sf, tw_sf, dl_sf, lg_sf integer :: c + + tr_sf = fuel_classes%trunks() + tw_sf = fuel_classes%twigs() + dl_sf = fuel_classes%dead_leaves() + lg_sf = fuel_classes%live_grass() currentPatch => currentSite%oldest_patch; do while(associated(currentPatch)) if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - - currentPatch%burnt_frac_litter(:) = 1.0_r8 - ! Calculate fraction of litter is burnt for all classes. - ! Equation B1 in Thonicke et al. 2010--- - do c = 1, nfsc !work out the burnt fraction for all pools, even if those pools dont exist. - moist = currentPatch%litter_moisture(c) - ! 1. Very dry litter - if (moist <= SF_val_min_moisture(c)) then - currentPatch%burnt_frac_litter(c) = 1.0_r8 - endif - ! 2. Low to medium moistures - if (moist > SF_val_min_moisture(c).and.moist <= SF_val_mid_moisture(c)) then - currentPatch%burnt_frac_litter(c) = max(0.0_r8,min(1.0_r8,SF_val_low_moisture_Coeff(c)- & - SF_val_low_moisture_Slope(c)*moist)) - else - ! For medium to high moistures. - if (moist > SF_val_mid_moisture(c).and.moist <= 1.0_r8) then - currentPatch%burnt_frac_litter(c) = max(0.0_r8,min(1.0_r8,SF_val_mid_moisture_Coeff(c)- & - SF_val_mid_moisture_Slope(c)*moist)) - endif - - endif - ! Very wet litter - if (moist >= 1.0_r8) then !this shouldn't happen? - currentPatch%burnt_frac_litter(c) = 0.0_r8 - endif - enddo !c - - ! we can't ever kill -all- of the grass. - currentPatch%burnt_frac_litter(lg_sf) = min(0.8_r8,currentPatch%burnt_frac_litter(lg_sf )) - - ! reduce burnt amount for mineral content. - currentPatch%burnt_frac_litter(:) = currentPatch%burnt_frac_litter(:) * (1.0_r8-SF_val_miner_total) - - !---Calculate amount of fuel burnt.--- - - litt_c => currentPatch%litter(element_pos(carbon12_element)) - FC_ground(tw_sf:tr_sf) = currentPatch%burnt_frac_litter(tw_sf:tr_sf) * litt_c%ag_cwd(tw_sf:tr_sf) - FC_ground(dl_sf) = currentPatch%burnt_frac_litter(dl_sf) * sum(litt_c%leaf_fines(:)) - FC_ground(lg_sf) = currentPatch%burnt_frac_litter(lg_sf) * currentPatch%livegrass + + currentPatch%fuel%frac_burnt(:) = 1.0_r8 + ! Calculate fraction of litter is burnt for all classes. + ! Equation B1 in Thonicke et al. 2010--- + do c = 1, num_fuel_classes !work out the burnt fraction for all pools, even if those pools dont exist. + moist = currentPatch%fuel%effective_moisture(c) + ! 1. Very dry litter + if (moist <= SF_val_min_moisture(c)) then + currentPatch%fuel%frac_burnt(c) = 1.0_r8 + endif + ! 2. Low to medium moistures + if (moist > SF_val_min_moisture(c).and.moist <= SF_val_mid_moisture(c)) then + currentPatch%fuel%frac_burnt(c) = max(0.0_r8,min(1.0_r8,SF_val_low_moisture_Coeff(c)- & + SF_val_low_moisture_Slope(c)*moist)) + else + ! For medium to high moistures. + if (moist > SF_val_mid_moisture(c).and.moist <= 1.0_r8) then + currentPatch%fuel%frac_burnt(c) = max(0.0_r8,min(1.0_r8,SF_val_mid_moisture_Coeff(c)- & + SF_val_mid_moisture_Slope(c)*moist)) + endif + + endif + ! Very wet litter + if (moist >= 1.0_r8) then !this shouldn't happen? + currentPatch%fuel%frac_burnt(c) = 0.0_r8 + endif + enddo !c + + ! we can't ever kill -all- of the grass. + currentPatch%fuel%frac_burnt(lg_sf) = min(0.8_r8,currentPatch%fuel%frac_burnt(lg_sf )) + + ! reduce burnt amount for mineral content. + currentPatch%fuel%frac_burnt(:) = currentPatch%fuel%frac_burnt(:) * (1.0_r8-SF_val_miner_total) + + !---Calculate amount of fuel burnt.--- + + litt_c => currentPatch%litter(element_pos(carbon12_element)) + FC_ground(tw_sf:tr_sf) = currentPatch%fuel%frac_burnt(tw_sf:tr_sf) * litt_c%ag_cwd(tw_sf:tr_sf) + FC_ground(dl_sf) = currentPatch%fuel%frac_burnt(dl_sf) * sum(litt_c%leaf_fines(:)) + FC_ground(lg_sf) = currentPatch%fuel%frac_burnt(lg_sf) * currentPatch%livegrass + + !call currentPatch%fuel%BurnFuel(fc_ground) ! Following used for determination of cambial kill follows from Peterson & Ryan (1986) scheme ! less empirical cf current scheme used in SPITFIRE which attempts to mesh Rothermel @@ -565,9 +431,9 @@ subroutine ground_fuel_consumption ( currentSite ) ! taul is the duration of the lethal heating. ! The /10 is to convert from kgC/m2 into gC/cm2, as in the Peterson and Ryan paper #Rosie,Jun 2013 - do c = 1,nfsc - tau_b(c) = 39.4_r8 *(currentPatch%fuel_frac(c)*currentPatch%sum_fuel/0.45_r8/10._r8)* & - (1.0_r8-((1.0_r8-currentPatch%burnt_frac_litter(c))**0.5_r8)) + do c = 1,num_fuel_classes + tau_b(c) = 39.4_r8 *(currentPatch%fuel%frac_loading(c)*currentPatch%fuel%non_trunk_loading/0.45_r8/10._r8)* & + (1.0_r8-((1.0_r8-currentPatch%fuel%frac_burnt(c))**0.5_r8)) enddo tau_b(tr_sf) = 0.0_r8 ! Cap the residence time to 8mins, as suggested by literature survey by P&R (1986). diff --git a/fire/SFNesterovMod.F90 b/fire/SFNesterovMod.F90 index f8001c2552..23128c880f 100644 --- a/fire/SFNesterovMod.F90 +++ b/fire/SFNesterovMod.F90 @@ -15,7 +15,7 @@ module SFNesterovMod procedure, public :: Init => init_nesterov_fire_weather procedure, public :: UpdateIndex => update_nesterov_index - + end type nesterov_index real(r8), parameter :: min_precip_thresh = 3.0_r8 ! threshold for precipitation above which to zero NI [mm/day] @@ -104,5 +104,6 @@ real(r8) function dewpoint(temp_C, rh) dewpoint = (dewpoint_b*yipsolon)/(dewpoint_a - yipsolon) end function dewpoint - + + !------------------------------------------------------------------------------------- end module SFNesterovMod \ No newline at end of file diff --git a/fire/SFParamsMod.F90 b/fire/SFParamsMod.F90 index c2dbc3fcd6..65d87e5c6d 100644 --- a/fire/SFParamsMod.F90 +++ b/fire/SFParamsMod.F90 @@ -2,17 +2,17 @@ module SFParamsMod ! ! module that deals with reading the SF parameter file ! - use FatesConstantsMod , only: r8 => fates_r8 - use FatesConstantsMod , only: fates_check_param_set - use FatesLitterMod , only: NFSC - use FatesLitterMod , only: ncwd + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : fates_check_param_set + use FatesFuelClassesMod, only : num_fuel_classes + use FatesLitterMod, only : ncwd use FatesParametersInterface, only : param_string_length - use FatesGlobals, only : fates_log - use FatesGlobals, only : endrun => fates_endrun - use shr_log_mod , only : errMsg => shr_log_errMsg + use FatesGlobals, only : fates_log + use FatesGlobals, only : endrun => fates_endrun + use shr_log_mod, only : errMsg => shr_log_errMsg implicit none - private ! Modules are private by default + private save ! @@ -28,15 +28,15 @@ module SFParamsMod real(r8),protected, public :: SF_val_drying_ratio real(r8),protected, public :: SF_val_fire_threshold ! threshold for fires that spread or go out. kW/m (Pyne 1996) real(r8),protected, public :: SF_val_CWD_frac(ncwd) - real(r8),protected, public :: SF_val_max_decomp(NFSC) - real(r8),protected, public :: SF_val_SAV(NFSC) - real(r8),protected, public :: SF_val_FBD(NFSC) - real(r8),protected, public :: SF_val_min_moisture(NFSC) - real(r8),protected, public :: SF_val_mid_moisture(NFSC) - real(r8),protected, public :: SF_val_low_moisture_Coeff(NFSC) - real(r8),protected, public :: SF_val_low_moisture_Slope(NFSC) - real(r8),protected, public :: SF_val_mid_moisture_Coeff(NFSC) - real(r8),protected, public :: SF_val_mid_moisture_Slope(NFSC) + real(r8),protected, public :: SF_val_max_decomp(num_fuel_classes) + real(r8),protected, public :: SF_val_SAV(num_fuel_classes) + real(r8),protected, public :: SF_val_FBD(num_fuel_classes) + real(r8),protected, public :: SF_val_min_moisture(num_fuel_classes) + real(r8),protected, public :: SF_val_mid_moisture(num_fuel_classes) + real(r8),protected, public :: SF_val_low_moisture_Coeff(num_fuel_classes) + real(r8),protected, public :: SF_val_low_moisture_Slope(num_fuel_classes) + real(r8),protected, public :: SF_val_mid_moisture_Coeff(num_fuel_classes) + real(r8),protected, public :: SF_val_mid_moisture_Slope(num_fuel_classes) character(len=param_string_length),parameter :: SF_name_fdi_alpha = "fates_fire_fdi_alpha" character(len=param_string_length),parameter :: SF_name_miner_total = "fates_fire_miner_total" @@ -58,18 +58,13 @@ module SFParamsMod character(len=param_string_length),parameter :: SF_name_mid_moisture_Coeff = "fates_fire_mid_moisture_Coeff" character(len=param_string_length),parameter :: SF_name_mid_moisture_Slope = "fates_fire_mid_moisture_Slope" - character(len=*), parameter, private :: sourcefile = & - __FILE__ - - - real(r8), parameter,private :: min_fire_threshold = 0.0001_r8 ! The minimum reasonable fire intensity threshold [kW/m] - + character(len=*), parameter, private :: sourcefile = __FILE__ + real(r8), parameter, private :: min_fire_threshold = 0.0001_r8 ! The minimum reasonable fire intensity threshold [kW/m] public :: SpitFireRegisterParams public :: SpitFireReceiveParams public :: SpitFireCheckParams - contains ! ===================================================================================== @@ -97,7 +92,7 @@ subroutine SpitFireCheckParams(is_master) if(.not.is_master) return ! Move these checks to initialization - do c = 1,nfsc + do c = 1,num_fuel_classes if ( SF_val_max_decomp(c) < 0._r8) then write(fates_log(),*) 'Decomposition rates should be >0' write(fates_log(),*) 'c = ',c,' SF_val_max_decomp(c) = ',SF_val_max_decomp(c) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index f56e3e0dd2..9fc11491c1 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -997,16 +997,9 @@ subroutine init_patches( nsites, sites, bc_in) currentPatch => sites(s)%youngest_patch do while(associated(currentPatch)) - currentPatch%litter_moisture(:) = 0._r8 - currentPatch%fuel_eff_moist = 0._r8 currentPatch%livegrass = 0._r8 - currentPatch%sum_fuel = 0._r8 - currentPatch%fuel_bulkd = 0._r8 - currentPatch%fuel_sav = 0._r8 - currentPatch%fuel_mef = 0._r8 currentPatch%ros_front = 0._r8 currentPatch%tau_l = 0._r8 - currentPatch%fuel_frac(:) = 0._r8 currentPatch%tfc_ros = 0._r8 currentPatch%fi = 0._r8 currentPatch%fire = 0 @@ -1014,8 +1007,6 @@ subroutine init_patches( nsites, sites, bc_in) currentPatch%ros_back = 0._r8 currentPatch%scorch_ht(:) = 0._r8 currentPatch%frac_burnt = 0._r8 - currentPatch%burnt_frac_litter(:) = 0._r8 - currentPatch => currentPatch%older enddo enddo diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 61878db1ef..20881c6e34 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -22,7 +22,7 @@ module EDTypesMod use PRTGenericMod, only : num_element_types use PRTGenericMod, only : carbon12_element use FatesLitterMod, only : litter_type - use FatesLitterMod, only : ncwd, NFSC + use FatesLitterMod, only : ncwd use FatesConstantsMod, only : days_per_year use FatesRunningMeanMod, only : rmean_type,rmean_arr_type use FatesConstantsMod, only : fates_unset_r8 diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 364ad816cf..1215f1745a 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -123,7 +123,7 @@ module FatesHistoryInterfaceMod use FatesSizeAgeTypeIndicesMod, only : get_layersizetype_class_index use FatesSizeAgeTypeIndicesMod, only : get_age_class_index - use FatesLitterMod , only : nfsc + use FatesFuelClassesMod , only : num_fuel_classes use FatesLitterMod , only : ncwd use FatesConstantsMod , only : ican_upper use FatesConstantsMod , only : ican_ustory @@ -2684,11 +2684,11 @@ subroutine update_history_dyn1(this,nc,nsites,sites,bc_in) hio_tfc_ros_si(io_si) = hio_tfc_ros_si(io_si) + cpatch%TFC_ROS * cpatch%area * AREA_INV hio_fire_intensity_si(io_si) = hio_fire_intensity_si(io_si) + cpatch%FI * cpatch%area * AREA_INV * J_per_kJ hio_fire_area_si(io_si) = hio_fire_area_si(io_si) + cpatch%frac_burnt * cpatch%area * AREA_INV / sec_per_day - hio_fire_fuel_bulkd_si(io_si) = hio_fire_fuel_bulkd_si(io_si) + cpatch%fuel_bulkd * cpatch%area * AREA_INV - hio_fire_fuel_eff_moist_si(io_si) = hio_fire_fuel_eff_moist_si(io_si) + cpatch%fuel_eff_moist * cpatch%area * AREA_INV - hio_fire_fuel_sav_si(io_si) = hio_fire_fuel_sav_si(io_si) + cpatch%fuel_sav * cpatch%area * AREA_INV / m_per_cm - hio_fire_fuel_mef_si(io_si) = hio_fire_fuel_mef_si(io_si) + cpatch%fuel_mef * cpatch%area * AREA_INV - hio_sum_fuel_si(io_si) = hio_sum_fuel_si(io_si) + cpatch%sum_fuel * cpatch%area * AREA_INV + hio_fire_fuel_bulkd_si(io_si) = hio_fire_fuel_bulkd_si(io_si) + cpatch%fuel%bulk_density_notrunks * cpatch%area * AREA_INV + hio_fire_fuel_eff_moist_si(io_si) = hio_fire_fuel_eff_moist_si(io_si) + cpatch%fuel%average_moisture_notrunks * cpatch%area * AREA_INV + hio_fire_fuel_sav_si(io_si) = hio_fire_fuel_sav_si(io_si) + cpatch%fuel%SAV_notrunks * cpatch%area * AREA_INV / m_per_cm + hio_fire_fuel_mef_si(io_si) = hio_fire_fuel_mef_si(io_si) + cpatch%fuel%MEF_notrunks * cpatch%area * AREA_INV + hio_sum_fuel_si(io_si) = hio_sum_fuel_si(io_si) + cpatch%fuel%non_trunk_loading * cpatch%area * AREA_INV hio_fire_intensity_area_product_si(io_si) = hio_fire_intensity_area_product_si(io_si) + & cpatch%FI * cpatch%frac_burnt * cpatch%area * AREA_INV * J_per_kJ @@ -3490,7 +3490,7 @@ subroutine update_history_dyn2(this,nc,nsites,sites,bc_in) ! Fuel sum [kg/m2] hio_fire_sum_fuel_si_age(io_si, cpatch%age_class) = hio_fire_sum_fuel_si_age(io_si, cpatch%age_class) + & - cpatch%sum_fuel * cpatch%area * AREA_INV + cpatch%fuel%non_trunk_loading * cpatch%area * AREA_INV @@ -4270,20 +4270,20 @@ subroutine update_history_dyn2(this,nc,nsites,sites,bc_in) hio_fragmentation_scaler_sl(io_si,ilyr) = hio_fragmentation_scaler_sl(io_si,ilyr) + cpatch%fragmentation_scaler(ilyr) * cpatch%area * AREA_INV end do - do i_fuel = 1,nfsc + do i_fuel = 1, num_fuel_classes i_agefuel = get_agefuel_class_index(cpatch%age,i_fuel) hio_fuel_amount_age_fuel(io_si,i_agefuel) = hio_fuel_amount_age_fuel(io_si,i_agefuel) + & - cpatch%fuel_frac(i_fuel) * cpatch%sum_fuel * cpatch%area * AREA_INV + cpatch%fuel%frac_loading(i_fuel) * cpatch%fuel%non_trunk_loading * cpatch%area * AREA_INV hio_litter_moisture_si_fuel(io_si, i_fuel) = hio_litter_moisture_si_fuel(io_si, i_fuel) + & - cpatch%litter_moisture(i_fuel) * cpatch%area * AREA_INV + cpatch%fuel%effective_moisture(i_fuel) * cpatch%area * AREA_INV hio_fuel_amount_si_fuel(io_si, i_fuel) = hio_fuel_amount_si_fuel(io_si, i_fuel) + & - cpatch%fuel_frac(i_fuel) * cpatch%sum_fuel * cpatch%area * AREA_INV + cpatch%fuel%frac_loading(i_fuel) * cpatch%fuel%non_trunk_loading * cpatch%area * AREA_INV hio_burnt_frac_litter_si_fuel(io_si, i_fuel) = hio_burnt_frac_litter_si_fuel(io_si, i_fuel) + & - cpatch%burnt_frac_litter(i_fuel) * cpatch%frac_burnt * cpatch%area * AREA_INV + cpatch%fuel%frac_burnt(i_fuel) * cpatch%frac_burnt * cpatch%area * AREA_INV end do diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 12eb83bcc0..e77b3e34bc 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1122,7 +1122,7 @@ end subroutine InitPARTEHGlobals subroutine fates_history_maps - use FatesLitterMod, only : NFSC + use FatesFuelClassesMod, only : num_fuel_classes use EDParamsMod, only : nclmax use EDParamsMod, only : nlevleaf use EDParamsMod, only : ED_val_history_sizeclass_bin_edges @@ -1157,7 +1157,7 @@ subroutine fates_history_maps allocate( fates_hdim_scmap_levscpf(1:nlevsclass*numpft)) allocate( fates_hdim_levpft(1:numpft )) allocate( fates_hdim_levlanduse(1:n_landuse_cats)) - allocate( fates_hdim_levfuel(1:NFSC )) + allocate( fates_hdim_levfuel(1:num_fuel_classes )) allocate( fates_hdim_levcwdsc(1:NCWD )) allocate( fates_hdim_levage(1:nlevage )) allocate( fates_hdim_levheight(1:nlevheight )) @@ -1180,8 +1180,8 @@ subroutine fates_history_maps allocate( fates_hdim_pftmap_levscagpft(nlevsclass * nlevage * numpft)) allocate( fates_hdim_agmap_levagepft(nlevage * numpft)) allocate( fates_hdim_pftmap_levagepft(nlevage * numpft)) - allocate( fates_hdim_agmap_levagefuel(nlevage * nfsc)) - allocate( fates_hdim_fscmap_levagefuel(nlevage * nfsc)) + allocate( fates_hdim_agmap_levagefuel(nlevage * num_fuel_classes)) + allocate( fates_hdim_fscmap_levagefuel(nlevage * num_fuel_classes)) allocate( fates_hdim_elmap_levelpft(num_elements*numpft)) allocate( fates_hdim_elmap_levelcwd(num_elements*ncwd)) @@ -1211,7 +1211,7 @@ subroutine fates_history_maps end do ! make fuel array - do ifuel=1,NFSC + do ifuel=1,num_fuel_classes fates_hdim_levfuel(ifuel) = ifuel end do @@ -1356,7 +1356,7 @@ subroutine fates_history_maps i=0 do iage=1,nlevage - do ifuel=1,NFSC + do ifuel=1,num_fuel_classes i=i+1 fates_hdim_agmap_levagefuel(i) = iage fates_hdim_fscmap_levagefuel(i) = ifuel diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 6236f4ef29..7f890e5ecc 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -40,7 +40,8 @@ module FatesRestartInterfaceMod use FatesInterfaceTypesMod, only : nlevsclass use FatesInterfaceTypesMod, only : nlevdamage use FatesLitterMod, only : litter_type - use FatesLitterMod, only : ncwd, nfsc + use FatesLitterMod, only : ncwd + use FatesFuelClassesMod, only : num_fuel_classes use FatesLitterMod, only : ndcmpy use EDTypesMod, only : area use EDParamsMod, only : nlevleaf @@ -2572,8 +2573,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) end do io_idx_pa_cwd = io_idx_co_1st - do i = 1,nfsc - this%rvars(ir_litter_moisture_pa_nfsc)%r81d(io_idx_pa_cwd) = cpatch%litter_moisture(i) + do i = 1,num_fuel_classes + this%rvars(ir_litter_moisture_pa_nfsc)%r81d(io_idx_pa_cwd) = cpatch%fuel%effective_moisture(i) io_idx_pa_cwd = io_idx_pa_cwd + 1 end do @@ -3548,8 +3549,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) end do io_idx_pa_cwd = io_idx_co_1st - do i = 1,nfsc - cpatch%litter_moisture(i) = this%rvars(ir_litter_moisture_pa_nfsc)%r81d(io_idx_pa_cwd) + do i = 1,num_fuel_classes + cpatch%fuel%effective_moisture(i) = this%rvars(ir_litter_moisture_pa_nfsc)%r81d(io_idx_pa_cwd) io_idx_pa_cwd = io_idx_pa_cwd + 1 end do diff --git a/testing/CMakeLists.txt b/testing/CMakeLists.txt index d59509c73d..7c19f7cc99 100644 --- a/testing/CMakeLists.txt +++ b/testing/CMakeLists.txt @@ -3,6 +3,8 @@ ## Functional tests add_subdirectory(functional_testing/allometry fates_allom_ftest) add_subdirectory(functional_testing/math_utils fates_math_ftest) +add_subdirectory(functional_testing/fire fates_fuel_ftest) ## Unit tests -add_subdirectory(unit_testing/fire_weather_test fates_fire_weather_utest) \ No newline at end of file +add_subdirectory(unit_testing/fire_weather_test fates_fire_weather_utest) +add_subdirectory(unit_testing/fire_fuel_test fates_fire_fuel_utest) diff --git a/testing/functional_testing/allometry/FatesTestAllometry.F90 b/testing/functional_testing/allometry/FatesTestAllometry.F90 index 91ee4df2cf..d4e9d28e7a 100644 --- a/testing/functional_testing/allometry/FatesTestAllometry.F90 +++ b/testing/functional_testing/allometry/FatesTestAllometry.F90 @@ -6,6 +6,7 @@ program FatesTestAllometry use FatesAllometryMod, only : bfineroot, bstore_allom, bdead_allom use PRTParametersMod, only : prt_params use FatesUnitTestParamReaderMod, only : fates_unit_test_param_reader + use FatesArgumentUtils, only : command_line_arg implicit none @@ -13,10 +14,8 @@ program FatesTestAllometry type(fates_unit_test_param_reader) :: param_reader ! param reader instance character(len=:), allocatable :: param_file ! input parameter file integer :: numpft ! number of pfts (from parameter file) - integer :: arglen ! length of command line argument integer :: i, j ! looping indices integer :: numdbh ! size of dbh array - integer :: nargs ! number of command line arguments real(r8), allocatable :: dbh(:) ! diameter at breast height [cm] real(r8), allocatable :: height(:, :) ! height [m] real(r8), allocatable :: bagw(:, :) ! aboveground woody biomass [kgC] @@ -76,17 +75,9 @@ end subroutine WriteAllometryData end interface - ! get parameter file from command-line argument - nargs = command_argument_count() - if (nargs /= 1) then - write(*, '(a, i2, a)') "Incorrect number of arguments: ", nargs, ". Should be 1." - stop - else - call get_command_argument(1, length=arglen) - allocate(character(arglen) :: param_file) - call get_command_argument(1, value=param_file) - endif - + ! read in parameter file name from command line + param_file = command_line_arg(1) + ! read in parameter file call param_reader%Init(param_file) call param_reader%RetrieveParameters() diff --git a/testing/functional_testing/fire/CMakeLists.txt b/testing/functional_testing/fire/CMakeLists.txt new file mode 100644 index 0000000000..cf4e55508d --- /dev/null +++ b/testing/functional_testing/fire/CMakeLists.txt @@ -0,0 +1,28 @@ +set(fire_test_sources + FatesTestFuel.F90 + FatesTestFireMod.F90 + SyntheticFuelModels.F90) + +set(NETCDF_C_DIR ${NETCDF_C_PATH}) +set(NETCDF_FORTRAN_DIR ${NETCDF_F_PATH}) + +FIND_PATH(NETCDFC_FOUND libnetcdf.a ${NETCDF_C_DIR}/lib) +FIND_PATH(NETCDFF_FOUND libnetcdff.a ${NETCDF_FORTRAN_DIR}/lib) + + +include_directories(${NETCDF_C_DIR}/include + ${NETCDF_FORTRAN_DIR}/include) + +link_directories(${NETCDF_C_DIR}/lib + ${NETCDF_FORTRAN_DIR}/lib + ${PFUNIT_TOP_DIR}/lib) + +add_executable(FATES_fuel_exe ${fire_test_sources}) + +target_link_libraries(FATES_fuel_exe + netcdf + netcdff + fates + csm_share + funit) + \ No newline at end of file diff --git a/testing/functional_testing/fire/FatesTestFireMod.F90 b/testing/functional_testing/fire/FatesTestFireMod.F90 new file mode 100644 index 0000000000..c37982039b --- /dev/null +++ b/testing/functional_testing/fire/FatesTestFireMod.F90 @@ -0,0 +1,274 @@ +module FatesTestFireMod + ! + ! DESCRIPTION: + ! Module to support testing the FATES SPIFTIRE model + ! + + use FatesConstantsMod, only : r8 => fates_r8 + use EDTypesMod, only : ed_site_type + use FatesPatchMod, only : fates_patch_type + use SFNesterovMod, only : nesterov_index + use FatesUnitTestIOMod, only : OpenNCFile, GetVar, CloseNCFile, RegisterNCDims + use FatesUnitTestIOMod, only : RegisterVar, EndNCDef, WriteVar + use FatesUnitTestIOMod, only : type_double, type_int, type_char + use FatesFuelClassesMod, only : num_fuel_classes + use SyntheticFuelModels, only : fuel_models_array_class + use SFParamsMod, only : SF_val_CWD_frac + use FatesFuelMod, only : fuel_type + + implicit none + private + + public :: SetUpFuel, ReadDatmData, WriteFireData + + contains + + !===================================================================================== + + subroutine SetUpFuel(fuel, fuel_model_array, fuel_model_index, fuel_name, fuel_carrier) + ! + ! DESCRIPTION: + ! Sets up fuel loading + ! + + ! ARGUMENTS: + type(fuel_type), intent(inout) :: fuel ! fuel object + type(fuel_models_array_class), intent(in) :: fuel_model_array ! array of fuel models + integer, intent(in) :: fuel_model_index ! fuel model index + character(len=100), intent(out) :: fuel_name ! name of fuel model + character(len=2), intent(out) :: fuel_carrier ! fuel carrier for fuel model + + ! LOCALS: + integer :: i ! position of fuel model in array + real(r8) :: leaf_litter ! leaf litter [kg/m2] + real(r8) :: twig_litter ! twig litter [kg/m2] + real(r8) :: small_branch_litter ! small branch litter [kg/m2] + real(r8) :: large_branch_litter ! large branch litter [kg/m2] + real(r8) :: grass_litter ! grass litter [kg/m2] + + + ! get fuel model position in array + i = fuel_model_array%FuelModelPosition(fuel_model_index) + + ! fuel model data + leaf_litter = fuel_model_array%fuel_models(i)%hr1_loading + twig_litter = fuel_model_array%fuel_models(i)%hr10_loading + + ! small vs. large branches based on input parameter file + small_branch_litter = fuel_model_array%fuel_models(i)%hr100_loading*SF_val_CWD_frac(2)/ & + (SF_val_CWD_frac(2) + SF_val_CWD_frac(3)) + large_branch_litter = fuel_model_array%fuel_models(i)%hr100_loading*SF_val_CWD_frac(3)/ & + (SF_val_CWD_frac(2) + SF_val_CWD_frac(3)) + + grass_litter = fuel_model_array%fuel_models(i)%live_herb_loading + + fuel_name = fuel_model_array%fuel_models(i)%fuel_model_name + fuel_carrier = fuel_model_array%fuel_models(i)%carrier + + call fuel%UpdateLoading(leaf_litter, twig_litter, small_branch_litter, & + large_branch_litter, 0.0_r8, grass_litter) + + end subroutine SetUpFuel + + !===================================================================================== + + subroutine ReadDatmData(nc_file, temp_degC, precip, rh, wind) + ! + ! DESCRIPTION: + ! Reads and returns DATM data + ! + + ! ARGUMENTS: + character(len=*), intent(in) :: nc_file ! netcdf file with DATM data + real(r8), allocatable, intent(out) :: temp_degC(:) ! daily air temperature [degC] + real(r8), allocatable, intent(out) :: precip(:) ! daily precipitation [mm] + real(r8), allocatable, intent(out) :: rh(:) ! daily relative humidity [%] + real(r8), allocatable, intent(out) :: wind(:) ! daily wind speed [m/s] + + ! LOCALS: + integer :: ncid ! netcdf file unit number + + ! open file + call OpenNCFile(trim(nc_file), ncid, 'read') + + ! read in data + call GetVar(ncid, 'temp_degC', temp_degC) + call GetVar(ncid, 'precip', precip) + call GetVar(ncid, 'RH', rh) + call GetVar(ncid, 'wind', wind) + + ! close file + call CloseNCFile(ncid) + + end subroutine ReadDatmData + + !===================================================================================== + + subroutine WriteFireData(out_file, nsteps, nfuelmods, temp_degC, precip, rh, NI, & + loading, frac_loading, fuel_BD, fuel_SAV, non_trunk_loading, fuel_moisture, & + fuel_models, carriers) + ! + ! DESCRIPTION: + ! writes out data from the unit test + ! + + ! ARGUMENTS: + character(len=*), intent(in) :: out_file + integer, intent(in) :: nsteps + integer, intent(in) :: nfuelmods + real(r8), intent(in) :: temp_degC(:) + real(r8), intent(in) :: precip(:) + real(r8), intent(in) :: rh(:) + real(r8), intent(in) :: NI(:) + real(r8), intent(in) :: loading(:,:) + real(r8), intent(in) :: frac_loading(:,:) + real(r8), intent(in) :: non_trunk_loading(:) + real(r8), intent(in) :: fuel_moisture(:,:) + real(r8), intent(in) :: fuel_BD(:) + real(r8), intent(in) :: fuel_SAV(:) + integer, intent(in) :: fuel_models(:) + character(len=2), intent(in) :: carriers(:) + + ! LOCALS: + integer, allocatable :: time_index(:) ! array of time index + integer :: ncid ! netcdf id + integer :: i ! looping index + character(len=20) :: dim_names(3) ! dimension names + integer :: dimIDs(3) ! dimension IDs + ! variable IDS + integer :: timeID, litterID + integer :: modID + integer :: tempID, precipID + integer :: rhID, NIID, loadingID + integer :: frac_loadingID + integer :: tot_loadingID + integer :: BDID, SAVID + integer :: moistID + integer :: cID + + ! create pft indices + allocate(time_index(nsteps)) + do i = 1, nsteps + time_index(i) = i + end do + + ! dimension names + dim_names = [character(len=20) :: 'time', 'litter_class', 'fuel_model'] + + ! open file + call OpenNCFile(trim(out_file), ncid, 'readwrite') + + ! register dimensions + call RegisterNCDims(ncid, dim_names, (/nsteps, num_fuel_classes, nfuelmods/), 3, dimIDs) + + ! first register dimension variables + + ! register time + call RegisterVar(ncid, 'time', dimIDs(1:1), type_int, & + [character(len=20) :: 'time_origin', 'units', 'calendar', 'long_name'], & + [character(len=150) :: '2018-01-01 00:00:00', 'days since 2018-01-01 00:00:00', & + 'gregorian', 'time'], & + 4, timeID) + + ! register litter class + call RegisterVar(ncid, 'litter_class', dimIDs(2:2), type_int, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: '', 'fuel class'], 2, litterID) + + ! register fuel models + call RegisterVar(ncid, 'fuel_model', dimIDs(3:3), type_int, & + [character(len=20) :: 'units', 'long_name'], & + [character(len=150) :: '', 'fuel model index'], 2, modID) + + ! then register actual variables + + ! register fuel carriers + call RegisterVar(ncid, 'carrier', dimIDs(3:3), type_char, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'fuel_model_index', '', 'carrier of fuel'], & + 3, cID) + + ! register temperature + call RegisterVar(ncid, 'temp_degC', dimIDs(1:1), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'time', 'degrees C', 'air temperature'], & + 3, tempID) + + ! register precipitation + call RegisterVar(ncid, 'precip', dimIDs(1:1), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'time', 'mm', 'precipitation'], & + 3, precipID) + + ! register relative humidity + call RegisterVar(ncid, 'RH', dimIDs(1:1), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'time', '%', 'relative humidity'], & + 3, rhID) + + ! register Nesterov Index + call RegisterVar(ncid, 'NI', dimIDs(1:1), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'time', '', 'Nesterov Index'], & + 3, NIID) + + ! register fuel moisture + call RegisterVar(ncid, 'fuel_moisture', (/dimIDs(1), dimIDs(3)/), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'time fuel_model', 'm3 m-3', 'average fuel moisture'], & + 3, moistID) + + ! register fuel loading + call RegisterVar(ncid, 'fuel_loading', dimIDs(2:3), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'litter_class fuel_model', 'kgC m-2', 'fuel loading'], & + 3, loadingID) + + ! register fractional fuel loading + call RegisterVar(ncid, 'frac_loading', dimIDs(2:3), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'litter_class fuel_model', '', 'fractional loading'], & + 3, frac_loadingID) + + ! register non-trunk fuel loading + call RegisterVar(ncid, 'non_trunk_loading', dimIDs(3:3), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'fuel_model', 'kgC m-2', 'total loading'], & + 3, tot_loadingID) + + ! register fuel bulk density + call RegisterVar(ncid, 'bulk_density', dimIDs(3:3), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'fuel_model', 'kg m-3', 'fuel bulk density'], & + 3, BDID) + + ! register fuel SAV + call RegisterVar(ncid, 'SAV', dimIDs(3:3), type_double, & + [character(len=20) :: 'coordinates', 'units', 'long_name'], & + [character(len=150) :: 'fuel_model', 'cm-1', 'fuel surface area to volume ratio'], & + 3, SAVID) + + ! finish defining variables + call EndNCDef(ncid) + + ! write out data + call WriteVar(ncid, timeID, time_index) + call WriteVar(ncid, litterID, (/1, 2, 3, 4, 5, 6/)) + call WriteVar(ncid, modID, fuel_models(:)) + call WriteVar(ncid, cID, carriers(:)) + call WriteVar(ncid, tempID, temp_degC(:)) + call WriteVar(ncid, precipID, precip(:)) + call WriteVar(ncid, rhID, rh(:)) + call WriteVar(ncid, NIID, NI(:)) + call WriteVar(ncid, loadingID, loading(:,:)) + call WriteVar(ncid, frac_loadingID, frac_loading(:,:)) + call WriteVar(ncid, tot_loadingID, non_trunk_loading(:)) + call WriteVar(ncid, moistiD, fuel_moisture(:,:)) + call WriteVar(ncid, BDID, fuel_BD(:)) + call WriteVar(ncid, SAVID, fuel_SAV(:)) + + call CloseNCFile(ncid) + + end subroutine WriteFireData + +end module FatesTestFireMod diff --git a/testing/functional_testing/fire/FatesTestFuel.F90 b/testing/functional_testing/fire/FatesTestFuel.F90 new file mode 100644 index 0000000000..42237e095e --- /dev/null +++ b/testing/functional_testing/fire/FatesTestFuel.F90 @@ -0,0 +1,123 @@ +program FatesTestFuel + + use FatesConstantsMod, only : r8 => fates_r8 + use EDTypesMod, only : ed_site_type + use FatesTestFireMod, only : SetUpFuel, ReadDatmData, WriteFireData + use FatesArgumentUtils, only : command_line_arg + use FatesUnitTestParamReaderMod, only : fates_unit_test_param_reader + use SyntheticFuelModels, only : fuel_models_array_class + use SFFireWeatherMod, only : fire_weather + use SFNesterovMod, only : nesterov_index + use FatesFuelMod, only : fuel_type + use FatesFuelClassesMod, only : num_fuel_classes + use SFParamsMod, only : SF_val_SAV, SF_val_drying_ratio + use SFParamsMod, only : SF_val_FBD + + implicit none + + ! LOCALS: + type(fates_unit_test_param_reader) :: param_reader ! param reader instance + type(fuel_models_array_class) :: fuel_models_array ! array of fuel models + class(fire_weather), pointer :: fireWeather ! fire weather object + type(fuel_type), allocatable :: fuel(:) ! fuel objects + character(len=:), allocatable :: param_file ! input parameter file + character(len=:), allocatable :: datm_file ! input DATM driver file + real(r8), allocatable :: temp_degC(:) ! daily air temperature [degC] + real(r8), allocatable :: precip(:) ! daily precipitation [mm] + real(r8), allocatable :: rh(:) ! daily relative humidity [%] + real(r8), allocatable :: wind(:) ! daily wind speed [m/s] + real(r8), allocatable :: NI(:) ! Nesterov index + real(r8), allocatable :: fuel_loading(:,:) ! fuel loading [kgC/m2] + real(r8), allocatable :: non_trunk_loading(:) ! non-trunk fuel loading [kgC/m2] + real(r8), allocatable :: frac_loading(:,:) ! fractional fuel loading [0-1] + real(r8), allocatable :: fuel_BD(:) ! bulk density of fuel [kg/m3] + real(r8), allocatable :: fuel_SAV(:) ! fuel surface area to volume ratio [/cm] + real(r8), allocatable :: fuel_moisture(:,:) ! fuel moisture [m3/m3] + character(len=100), allocatable :: fuel_names(:) ! names of fuel models + character(len=2), allocatable :: carriers(:) ! carriers of fuel models + integer :: i, f ! looping indices + integer :: num_fuel_models ! number of fuel models to test + + ! CONSTANTS: + integer, parameter :: n_days = 365 ! number of days to run simulation + character(len=*), parameter :: out_file = 'fuel_out.nc' ! output file + + ! fuel models to test + integer, parameter, dimension(3) :: fuel_models = (/102, 183, 164/) + + ! number of fuel models to test + num_fuel_models = size(fuel_models) + + ! allocate arrays + allocate(temp_degC(n_days)) + allocate(precip(n_days)) + allocate(rh(n_days)) + allocate(wind(n_days)) + allocate(NI(n_days)) + allocate(fuel_moisture(n_days, num_fuel_models)) + allocate(fuel_loading(num_fuel_classes, num_fuel_models)) + allocate(frac_loading(num_fuel_classes, num_fuel_models)) + allocate(fuel_BD(num_fuel_models)) + allocate(fuel_SAV(num_fuel_models)) + allocate(non_trunk_loading(num_fuel_models)) + allocate(fuel_names(num_fuel_models)) + allocate(carriers(num_fuel_models)) + + ! read in parameter file name and DATM file from command line + param_file = command_line_arg(1) + datm_file = command_line_arg(2) + + ! read in parameter file + call param_reader%Init(param_file) + call param_reader%RetrieveParameters() + + ! read in DATM data + call ReadDatmData(datm_file, temp_degC, precip, rh, wind) + + ! set up fire weather class + allocate(nesterov_index :: fireWeather) + call fireWeather%Init() + + ! set up fuel objects and calculate loading + allocate(fuel(num_fuel_models)) + call fuel_models_array%GetFuelModels() + do f = 1, num_fuel_models + + ! uses data from fuel_models to initialize fuel + call SetUpFuel(fuel(f), fuel_models_array, fuel_models(f), fuel_names(f), carriers(f)) + + ! sum up fuel and calculate loading + call fuel(f)%SumLoading() + call fuel(f)%CalculateFractionalLoading() + + ! calculate geometric properties + call fuel(f)%AverageBulkDensity_NoTrunks(SF_val_FBD) + call fuel(f)%AverageSAV_NoTrunks(SF_val_SAV) + + ! save values + fuel_loading(:,f) = fuel(f)%loading(:) + non_trunk_loading(f) = fuel(f)%non_trunk_loading + frac_loading(:,f) = fuel(f)%frac_loading(:) + fuel_BD(f) = fuel(f)%bulk_density_notrunks + fuel_SAV(f) = fuel(f)%SAV_notrunks + + end do + + ! run on time steps + do i = 1, n_days + call fireWeather%UpdateIndex(temp_degC(i), precip(i), rh(i), wind(i)) + NI(i) = fireWeather%fire_weather_index + + ! calculate fuel moisture [m3/m3] + do f = 1, num_fuel_models + call fuel(f)%UpdateFuelMoisture(SF_val_SAV, SF_val_drying_ratio, fireWeather) + fuel_moisture(i, f) = fuel(f)%average_moisture_notrunks + end do + end do + + ! write out data + call WriteFireData(out_file, n_days, num_fuel_models, temp_degC, precip, rh, NI, & + fuel_loading, frac_loading, fuel_BD, fuel_SAV, non_trunk_loading, fuel_moisture, & + fuel_models, carriers) + +end program FatesTestFuel diff --git a/testing/functional_testing/fire/SyntheticFuelModels.F90 b/testing/functional_testing/fire/SyntheticFuelModels.F90 new file mode 100644 index 0000000000..ce1c8e85e2 --- /dev/null +++ b/testing/functional_testing/fire/SyntheticFuelModels.F90 @@ -0,0 +1,413 @@ +module SyntheticFuelModels + + use FatesConstantsMod, only : r8 => fates_r8 + + implicit none + private + + ! Fuel model numbers come from Scott and Burgen (2005) RMRS-GTR-153 + + integer, parameter, public, dimension(52) :: all_fuel_models = (/1, 2, 101, 102, 104, & + 107, 121, 122, 3, 103, 105, 106, & + 108, 109, 123, 124, 4, 5, 6, 141, & + 142, 145, 147, 161, 164, 10, 7, & + 143, 144, 146, 148, 149, 162, & + 163, 8, 9, 181, 182, 183, 184, & + 185, 186, 187, 188, 189, 11, 12, & + 13, 201, 202, 203, 204/) + + integer, parameter :: chunk_size = 10 + real(r8), parameter :: ustons_to_kg = 907.185_r8 + real(r8), parameter :: acres_to_m2 = 4046.86_r8 + real(r8), parameter :: ustons_acre_to_kgC_m2 = ustons_to_kg/acres_to_m2*0.45_r8 + real(r8), parameter :: ft_to_m = 0.3048_r8 + + ! holds data for fake fuel models that can be used for functional + ! testing of the FATES fire model + ! these are taken from the fire behavior fuel models in Scott & Burgan 2005 + type, public :: synthetic_fuel_model + + integer :: fuel_model_index ! fuel model index + character(len=2) :: carrier ! carrier ('GR', 'GS', etc.) + character(len=5) :: fuel_model_code ! carrier plus fuel model + character(len=100) :: fuel_model_name ! long name of fuel model + real(r8) :: wind_adj_factor ! wind adjustment factor + real(r8) :: hr1_loading ! fuel loading for 1 hour fuels [kg/m2] + real(r8) :: hr10_loading ! fuel loading for 10 hour fuels [kg/m2] + real(r8) :: hr100_loading ! fuel loading for 100 hour fuels [kg/m2] + real(r8) :: live_herb_loading ! fuel loading for live herbacious fuels [kg/m2] + real(r8) :: live_woody_loading ! fuel loading for live woody fuels [kg/m2] + real(r8) :: fuel_depth ! fuel bed depth [m] + contains + + procedure :: InitFuelModel + + end type synthetic_fuel_model + + ! -------------------------------------------------------------------------------------- + + ! a class to just hold an array of these fuel models + type, public :: fuel_models_array_class + + type(synthetic_fuel_model), allocatable :: fuel_models(:) ! array of fuel models + integer :: num_fuel_models ! number of total fuel models + + contains + + procedure :: AddFuelModel + procedure :: GetFuelModels + procedure :: FuelModelPosition + + end type fuel_models_array_class + + ! -------------------------------------------------------------------------------------- + + contains + + subroutine InitFuelModel(this, fuel_model_index, carrier, fuel_model_name, & + wind_adj_factor, hr1_loading, hr10_loading, hr100_loading, live_herb_loading, & + live_woody_loading, fuel_depth) + ! + ! DESCRIPTION: + ! Initializes the fuel model with input characteristics + ! Also converts units as needed + ! + ! NOTE THE UNITS ON INPUTS + ! + + ! ARGUMENTS: + class(synthetic_fuel_model), intent(inout) :: this + integer, intent(in) :: fuel_model_index ! fuel model index + character(len=2), intent(in) :: carrier ! main carrier + character(len=*), intent(in) :: fuel_model_name ! fuel model long name + real(r8), intent(in) :: wind_adj_factor ! wind adjustment factor + real(r8), intent(in) :: hr1_loading ! loading for 1-hr fuels [US tons/acre] + real(r8), intent(in) :: hr10_loading ! loading for 10-hr fuels [US tons/acre] + real(r8), intent(in) :: hr100_loading ! loading for 100-hr fuels [US tons/acre] + real(r8), intent(in) :: live_herb_loading ! loading for live herbacious fuels [US tons/acre] + real(r8), intent(in) :: live_woody_loading ! loading for live woody fuels [US tons/acre] + real(r8), intent(in) :: fuel_depth ! fuel bed depth [ft] + + this%fuel_model_index = fuel_model_index + this%carrier = carrier + this%fuel_model_name = fuel_model_name + this%wind_adj_factor = wind_adj_factor + this%hr1_loading = hr1_loading*ustons_acre_to_kgC_m2 ! convert to kgC/m2 + this%hr10_loading = hr10_loading*ustons_acre_to_kgC_m2 ! convert to kgC/m2 + this%hr100_loading = hr100_loading*ustons_acre_to_kgC_m2 ! convert to kgC/m2 + this%live_herb_loading = live_herb_loading*ustons_acre_to_kgC_m2 ! convert to kgC/m2 + this%live_woody_loading = live_woody_loading*ustons_acre_to_kgC_m2 ! convert to kgC/m2 + this%fuel_depth = fuel_depth*ft_to_m ! convert to m + + end subroutine InitFuelModel + + ! -------------------------------------------------------------------------------------- + + subroutine AddFuelModel(this, fuel_model_index, carrier, fuel_model_name, & + wind_adj_factor, hr1_loading, hr10_loading, hr100_loading, live_herb_loading, & + live_woody_loading, fuel_depth) + ! + ! DESCRIPTION: + ! Adds a fuel model to the dynamic array + ! + ! NOTE THE UNITS ON INPUTS + ! + + ! ARGUMENTS: + class(fuel_models_array_class), intent(inout) :: this ! array of fuel models + integer, intent(in) :: fuel_model_index ! fuel model index + character(len=2), intent(in) :: carrier ! main carrier + character(len=*), intent(in) :: fuel_model_name ! fuel model long name + real(r8), intent(in) :: wind_adj_factor ! wind adjustment factor + real(r8), intent(in) :: hr1_loading ! loading for 1-hr fuels [US tons/acre] + real(r8), intent(in) :: hr10_loading ! loading for 10-hr fuels [US tons/acre] + real(r8), intent(in) :: hr100_loading ! loading for 100-hr fuels [US tons/acre] + real(r8), intent(in) :: live_herb_loading ! loading for live herbacious fuels [US tons/acre] + real(r8), intent(in) :: live_woody_loading ! loading for live woody fuels [US tons/acre] + real(r8), intent(in) :: fuel_depth ! fuel bed depth [ft] + + ! LOCALS: + type(synthetic_fuel_model) :: fuel_model ! fuel model + type(synthetic_fuel_model), allocatable :: temporary_array(:) ! temporary array to hold data while re-allocating + + ! first make sure we have enough space in the array + if (allocated(this%fuel_models)) then + ! already allocated to some size + if (this%num_fuel_models == size(this%fuel_models)) then + ! need to add more space + allocate(temporary_array(size(this%fuel_models) + chunk_size)) + temporary_array(1:size(this%fuel_models)) = this%fuel_models + call move_alloc(temporary_array, this%fuel_models) + end if + + this%num_fuel_models = this%num_fuel_models + 1 + + else + ! first element in array + allocate(this%fuel_models(chunk_size)) + this%num_fuel_models = 1 + end if + + call fuel_model%InitFuelModel(fuel_model_index, carrier, fuel_model_name, & + wind_adj_factor, hr1_loading, hr10_loading, hr100_loading, live_herb_loading, & + live_woody_loading, fuel_depth) + + this%fuel_models(this%num_fuel_models) = fuel_model + + end subroutine AddFuelModel + + ! -------------------------------------------------------------------------------------- + + integer function FuelModelPosition(this, fuel_model_index) + ! + ! DESCRIPTION: + ! Returns the index of a desired fuel model + ! + + ! ARGUMENTS: + class(fuel_models_array_class), intent(in) :: this ! array of fuel models + integer, intent(in) :: fuel_model_index ! desired fuel model index + + ! LOCALS: + integer :: i ! looping index + + do i = 1, this%num_fuel_models + if (this%fuel_models(i)%fuel_model_index == fuel_model_index) then + FuelModelPosition = i + return + end if + end do + write(*, '(a, i2, a)') "Cannot find the fuel model index ", fuel_model_index, "." + stop + + end function FuelModelPosition + + ! -------------------------------------------------------------------------------------- + + subroutine GetFuelModels(this) + ! + ! DESCRIPTION: + ! Returns an array of hard-coded fuel models + ! these are taken from the fire behavior fuel models in Scott & Burgan 2005 + ! + + ! ARGUMENTS: + class(fuel_models_array_class), intent(inout) :: this ! array of fuel models + + call this%AddFuelModel(fuel_model_index=1, carrier='GR', fuel_model_name='short grass', & + wind_adj_factor=0.36_r8, hr1_loading=0.7_r8, hr10_loading=0.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=2, carrier='GR', fuel_model_name='timber and grass understory', & + wind_adj_factor=0.36_r8, hr1_loading=2.0_r8, hr10_loading=1.0_r8, hr100_loading=0.5_r8, & + live_herb_loading=0.5_r8, live_woody_loading=0.0_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=3, carrier='GR', fuel_model_name='tall grass', & + wind_adj_factor=0.44_r8, hr1_loading=3.0_r8, hr10_loading=0.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=2.5_r8) + + call this%AddFuelModel(fuel_model_index=4, carrier='SH', fuel_model_name='chapparal', & + wind_adj_factor=0.55_r8, hr1_loading=5.0_r8, hr10_loading=4.0_r8, hr100_loading=2.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=5.0_r8, fuel_depth=6.0_r8) + + call this%AddFuelModel(fuel_model_index=5, carrier='SH', fuel_model_name='brush', & + wind_adj_factor=0.42_r8, hr1_loading=1.0_r8, hr10_loading=0.5_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=2.0_r8, fuel_depth=2.0_r8) + + call this%AddFuelModel(fuel_model_index=6, carrier='SH', fuel_model_name='dormant brush', & + wind_adj_factor=0.44_r8, hr1_loading=1.5_r8, hr10_loading=2.5_r8, hr100_loading=2.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=2.5_r8) + + call this%AddFuelModel(fuel_model_index=7, carrier='SH', fuel_model_name='southern rough', & + wind_adj_factor=0.44_r8, hr1_loading=1.1_r8, hr10_loading=1.9_r8, hr100_loading=1.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.4_r8, fuel_depth=2.5_r8) + + call this%AddFuelModel(fuel_model_index=8, carrier='TL', fuel_model_name='compact timber litter', & + wind_adj_factor=0.28_r8, hr1_loading=1.5_r8, hr10_loading=1.0_r8, hr100_loading=2.5_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.2_r8) + + call this%AddFuelModel(fuel_model_index=9, carrier='TL', fuel_model_name='hardwood litter', & + wind_adj_factor=0.28_r8, hr1_loading=2.9_r8, hr10_loading=0.4_r8, hr100_loading=0.2_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.2_r8) + + call this%AddFuelModel(fuel_model_index=10, carrier='TU', fuel_model_name='timber and litter understorey', & + wind_adj_factor=0.46_r8, hr1_loading=3.0_r8, hr10_loading=2.0_r8, hr100_loading=5.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=2.0_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=11, carrier='SB', fuel_model_name='light slash', & + wind_adj_factor=0.36_r8, hr1_loading=1.5_r8, hr10_loading=4.5_r8, hr100_loading=5.5_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=12, carrier='SB', fuel_model_name='medium slash', & + wind_adj_factor=0.43_r8, hr1_loading=4.0_r8, hr10_loading=14.0_r8, hr100_loading=16.5_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=2.3_r8) + + call this%AddFuelModel(fuel_model_index=13, carrier='SB', fuel_model_name='heavy slash', & + wind_adj_factor=0.46_r8, hr1_loading=7.0_r8, hr10_loading=23.0_r8, hr100_loading=28.1_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=3.0_r8) + + call this%AddFuelModel(fuel_model_index=101, carrier='GR', fuel_model_name='short, sparse dry climate grass', & + wind_adj_factor=0.31_r8, hr1_loading=0.1_r8, hr10_loading=0.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.3_r8, live_woody_loading=0.0_r8, fuel_depth=0.4_r8) + + call this%AddFuelModel(fuel_model_index=102, carrier='GR', fuel_model_name='low load dry climate grass', & + wind_adj_factor=0.36_r8, hr1_loading=0.1_r8, hr10_loading=0.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=1.0_r8, live_woody_loading=0.0_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=103, carrier='GR', fuel_model_name='low load very coarse humid climate grass', & + wind_adj_factor=0.42_r8, hr1_loading=0.1_r8, hr10_loading=0.4_r8, hr100_loading=0.0_r8, & + live_herb_loading=1.5_r8, live_woody_loading=0.0_r8, fuel_depth=2.0_r8) + + call this%AddFuelModel(fuel_model_index=104, carrier='GR', fuel_model_name='moderate load dry climate grass', & + wind_adj_factor=0.42_r8, hr1_loading=0.3_r8, hr10_loading=0.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=1.9_r8, live_woody_loading=0.0_r8, fuel_depth=2.0_r8) + + call this%AddFuelModel(fuel_model_index=105, carrier='GR', fuel_model_name='low load humid climate grass', & + wind_adj_factor=0.39_r8, hr1_loading=0.4_r8, hr10_loading=0.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=2.5_r8, live_woody_loading=0.0_r8, fuel_depth=1.5_r8) + + call this%AddFuelModel(fuel_model_index=106, carrier='GR', fuel_model_name='moderate load humid climate grass', & + wind_adj_factor=0.39_r8, hr1_loading=0.1_r8, hr10_loading=0.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=3.4_r8, live_woody_loading=0.0_r8, fuel_depth=1.5_r8) + + call this%AddFuelModel(fuel_model_index=107, carrier='GR', fuel_model_name='high load dry climate grass', & + wind_adj_factor=0.46_r8, hr1_loading=1.0_r8, hr10_loading=0.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=5.4_r8, live_woody_loading=0.0_r8, fuel_depth=3.0_r8) + + call this%AddFuelModel(fuel_model_index=108, carrier='GR', fuel_model_name='high load humid climate grass', & + wind_adj_factor=0.49_r8, hr1_loading=0.5_r8, hr10_loading=1.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=7.3_r8, live_woody_loading=0.0_r8, fuel_depth=4.0_r8) + + call this%AddFuelModel(fuel_model_index=109, carrier='GR', fuel_model_name='very high load humid climate grass-shrub', & + wind_adj_factor=0.52_r8, hr1_loading=1.0_r8, hr10_loading=1.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=9.0_r8, live_woody_loading=0.0_r8, fuel_depth=5.0_r8) + + call this%AddFuelModel(fuel_model_index=121, carrier='GS', fuel_model_name='low load dry climate grass-shrub', & + wind_adj_factor=0.35_r8, hr1_loading=0.2_r8, hr10_loading=0.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.5_r8, live_woody_loading=0.7_r8, fuel_depth=0.9_r8) + + call this%AddFuelModel(fuel_model_index=122, carrier='GS', fuel_model_name='moderate load dry climate grass-shrub', & + wind_adj_factor=0.39_r8, hr1_loading=0.5_r8, hr10_loading=0.5_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.6_r8, live_woody_loading=1.0_r8, fuel_depth=1.5_r8) + + call this%AddFuelModel(fuel_model_index=123, carrier='GS', fuel_model_name='moderate load humid climate grass-shrub', & + wind_adj_factor=0.41_r8, hr1_loading=0.3_r8, hr10_loading=0.3_r8, hr100_loading=0.0_r8, & + live_herb_loading=1.5_r8, live_woody_loading=1.3_r8, fuel_depth=1.8_r8) + + call this%AddFuelModel(fuel_model_index=124, carrier='GS', fuel_model_name='high load humid climate grass-shrub', & + wind_adj_factor=0.42_r8, hr1_loading=1.9_r8, hr10_loading=0.3_r8, hr100_loading=0.1_r8, & + live_herb_loading=3.4_r8, live_woody_loading=7.1_r8, fuel_depth=2.1_r8) + + call this%AddFuelModel(fuel_model_index=141, carrier='SH', fuel_model_name='low load dry climate shrub', & + wind_adj_factor=0.36_r8, hr1_loading=0.3_r8, hr10_loading=0.3_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.2_r8, live_woody_loading=1.3_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=142, carrier='SH', fuel_model_name='moderate load dry climate shrub', & + wind_adj_factor=0.36_r8, hr1_loading=1.4_r8, hr10_loading=2.4_r8, hr100_loading=0.8_r8, & + live_herb_loading=0.0_r8, live_woody_loading=3.9_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=143, carrier='SH', fuel_model_name='moderate load humid climate shrub', & + wind_adj_factor=0.44_r8, hr1_loading=0.5_r8, hr10_loading=3.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=6.2_r8, fuel_depth=2.4_r8) + + call this%AddFuelModel(fuel_model_index=144, carrier='SH', fuel_model_name='low load humid climate timber-shrub', & + wind_adj_factor=0.46_r8, hr1_loading=0.9_r8, hr10_loading=1.2_r8, hr100_loading=0.2_r8, & + live_herb_loading=0.0_r8, live_woody_loading=2.6_r8, fuel_depth=3.0_r8) + + call this%AddFuelModel(fuel_model_index=145, carrier='SH', fuel_model_name='high load dry climate shrub', & + wind_adj_factor=0.55_r8, hr1_loading=3.6_r8, hr10_loading=2.1_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=2.9_r8, fuel_depth=6.0_r8) + + call this%AddFuelModel(fuel_model_index=146, carrier='SH', fuel_model_name='low load humid climate shrub', & + wind_adj_factor=0.42_r8, hr1_loading=2.9_r8, hr10_loading=1.5_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=1.4_r8, fuel_depth=2.0_r8) + + call this%AddFuelModel(fuel_model_index=147, carrier='SH', fuel_model_name='very high load dry climate shrub', & + wind_adj_factor=0.55_r8, hr1_loading=3.5_r8, hr10_loading=5.3_r8, hr100_loading=2.2_r8, & + live_herb_loading=0.0_r8, live_woody_loading=3.4_r8, fuel_depth=6.0_r8) + + call this%AddFuelModel(fuel_model_index=148, carrier='SH', fuel_model_name='high load humid climate shrub', & + wind_adj_factor=0.46_r8, hr1_loading=2.1_r8, hr10_loading=3.4_r8, hr100_loading=0.9_r8, & + live_herb_loading=0.0_r8, live_woody_loading=4.4_r8, fuel_depth=3.0_r8) + + call this%AddFuelModel(fuel_model_index=149, carrier='SH', fuel_model_name='very high load humid climate shrub', & + wind_adj_factor=0.5_r8, hr1_loading=4.5_r8, hr10_loading=2.5_r8, hr100_loading=0.0_r8, & + live_herb_loading=1.6_r8, live_woody_loading=7.0_r8, fuel_depth=4.4_r8) + + call this%AddFuelModel(fuel_model_index=161, carrier='TU', fuel_model_name='light load dry climate timber-grass-shrub', & + wind_adj_factor=0.33_r8, hr1_loading=0.2_r8, hr10_loading=0.9_r8, hr100_loading=1.5_r8, & + live_herb_loading=0.2_r8, live_woody_loading=0.9_r8, fuel_depth=0.6_r8) + + call this%AddFuelModel(fuel_model_index=162, carrier='TU', fuel_model_name='moderate load humid climate timber-shrub', & + wind_adj_factor=0.36_r8, hr1_loading=1.0_r8, hr10_loading=1.8_r8, hr100_loading=1.3_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.2_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=163, carrier='TU', fuel_model_name='moderate load humid climate timber-grass-shrub', & + wind_adj_factor=0.38_r8, hr1_loading=1.1_r8, hr10_loading=0.2_r8, hr100_loading=0.2_r8, & + live_herb_loading=0.3_r8, live_woody_loading=0.7_r8, fuel_depth=1.3_r8) + + call this%AddFuelModel(fuel_model_index=164, carrier='TU', fuel_model_name='dwarf conifer with understory', & + wind_adj_factor=0.32_r8, hr1_loading=4.5_r8, hr10_loading=0.0_r8, hr100_loading=0.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=2.0_r8, fuel_depth=0.5_r8) + + call this%AddFuelModel(fuel_model_index=165, carrier='TU', fuel_model_name='very high load dry climate timber-shrub', & + wind_adj_factor=0.33_r8, hr1_loading=4.0_r8, hr10_loading=4.0_r8, hr100_loading=3.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=3.0_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=181, carrier='TL', fuel_model_name='low load compact conifer litter', & + wind_adj_factor=0.28_r8, hr1_loading=1.0_r8, hr10_loading=2.2_r8, hr100_loading=3.6_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.2_r8) + + call this%AddFuelModel(fuel_model_index=182, carrier='TL', fuel_model_name='low load broadleaf litter', & + wind_adj_factor=0.28_r8, hr1_loading=1.4_r8, hr10_loading=2.3_r8, hr100_loading=2.2_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.2_r8) + + call this%AddFuelModel(fuel_model_index=183, carrier='TL', fuel_model_name='moderate load conifer litter', & + wind_adj_factor=0.29_r8, hr1_loading=0.5_r8, hr10_loading=2.2_r8, hr100_loading=2.8_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.3_r8) + + call this%AddFuelModel(fuel_model_index=184, carrier='TL', fuel_model_name='small downed logs', & + wind_adj_factor=0.31_r8, hr1_loading=0.5_r8, hr10_loading=1.5_r8, hr100_loading=4.2_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.4_r8) + + call this%AddFuelModel(fuel_model_index=185, carrier='TL', fuel_model_name='high load conifer litter', & + wind_adj_factor=0.33_r8, hr1_loading=1.2_r8, hr10_loading=2.5_r8, hr100_loading=4.4_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.6_r8) + + call this%AddFuelModel(fuel_model_index=186, carrier='TL', fuel_model_name='moderate load broadleaf litter', & + wind_adj_factor=0.29_r8, hr1_loading=2.4_r8, hr10_loading=1.2_r8, hr100_loading=1.2_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.3_r8) + + call this%AddFuelModel(fuel_model_index=187, carrier='TL', fuel_model_name='large downed logs', & + wind_adj_factor=0.31_r8, hr1_loading=0.3_r8, hr10_loading=1.4_r8, hr100_loading=8.1_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.4_r8) + + call this%AddFuelModel(fuel_model_index=188, carrier='TL', fuel_model_name='long-needle litter', & + wind_adj_factor=0.29_r8, hr1_loading=5.0_r8, hr10_loading=1.4_r8, hr100_loading=1.1_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.3_r8) + + call this%AddFuelModel(fuel_model_index=189, carrier='TL', fuel_model_name='very high load broadleaf litter', & + wind_adj_factor=0.33_r8, hr1_loading=6.7_r8, hr10_loading=3.3_r8, hr100_loading=4.2_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=0.6_r8) + + call this%AddFuelModel(fuel_model_index=201, carrier='SB', fuel_model_name='low load activity fuel', & + wind_adj_factor=0.36_r8, hr1_loading=1.5_r8, hr10_loading=3.0_r8, hr100_loading=11.1_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=202, carrier='SB', fuel_model_name='moderate load activity fuel or low load blowdown', & + wind_adj_factor=0.36_r8, hr1_loading=4.5_r8, hr10_loading=4.3_r8, hr100_loading=4.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=1.0_r8) + + call this%AddFuelModel(fuel_model_index=203, carrier='SB', fuel_model_name='high load activity fuel or moderate load blowdown', & + wind_adj_factor=0.38_r8, hr1_loading=5.5_r8, hr10_loading=2.8_r8, hr100_loading=3.0_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=1.2_r8) + + call this%AddFuelModel(fuel_model_index=204, carrier='SB', fuel_model_name='high load blowdown', & + wind_adj_factor=0.45_r8, hr1_loading=5.3_r8, hr10_loading=3.5_r8, hr100_loading=5.3_r8, & + live_herb_loading=0.0_r8, live_woody_loading=0.0_r8, fuel_depth=2.7_r8) + + end subroutine GetFuelModels + + ! -------------------------------------------------------------------------------------- + +end module SyntheticFuelModels diff --git a/testing/functional_testing/fire/fuel_test.py b/testing/functional_testing/fire/fuel_test.py new file mode 100644 index 0000000000..b9cd151621 --- /dev/null +++ b/testing/functional_testing/fire/fuel_test.py @@ -0,0 +1,187 @@ +""" +Concrete class for running the fuel functional test for FATES. +""" +import os +import numpy as np +import xarray as xr +import matplotlib.pyplot as plt +from functional_class import FunctionalTest + + +class FuelTest(FunctionalTest): + """Fuel test class""" + + name = "fuel" + + def __init__(self, test_dict): + super().__init__( + FuelTest.name, + test_dict["test_dir"], + test_dict["test_exe"], + test_dict["out_file"], + test_dict["use_param_file"], + test_dict["other_args"], + ) + self.plot = True + + def plot_output(self, run_dir: str, save_figs: bool, plot_dir: str): + """Plot output associated with fuel tests + + Args: + run_dir (str): run directory + out_file (str): output file + save_figs (bool): whether or not to save the figures + plot_dir (str): plot directory + """ + + fuel_dat = xr.open_dataset(os.path.join(run_dir, self.out_file)) + + self.plot_NI_dat(fuel_dat, save_figs, plot_dir) + self.plot_moisture_dat(fuel_dat, save_figs, plot_dir) + self.plot_barchart( + fuel_dat, + "fuel_loading", + "Fuel loading", + "kgC m$^{-2}$", + save_figs, + plot_dir, + ) + self.plot_barchart( + fuel_dat, + "frac_loading", + "Fractional fuel loading", + "0-1", + save_figs, + plot_dir, + ) + self.plot_barchart( + fuel_dat, + "bulk_density", + "Fuel bulk density", + "kg m$^{-3}$", + save_figs, + plot_dir, + by_litter_type=False, + ) + self.plot_barchart( + fuel_dat, + "SAV", + "Fuel surface area to volume ratio", + "cm$^{-1}$", + save_figs, + plot_dir, + by_litter_type=False, + ) + + @staticmethod + def plot_barchart( + fuel_dat: xr.Dataset, + var: str, + varname: str, + units: str, + save_figs: bool, + plot_dir: bool, + by_litter_type: bool = True, + ): + """Plots fuel data output as a bar chart + + Args: + fuel_dat (xr.Dataset): fuel data output + var (str): variable to plot + varname (str): variable name for x axis + units (str): units description + save_figs (bool): whether or not to save figure + plot_dir (bool): where to save figure + by_litter_type (bool, optional): whether the bar chart is by litter type. Defaults to True. + """ + + litter_classes = [ + "twigs", + "small branches", + "large branches", + "trunks", + "dead leaves", + "live grass", + ] + colors = [ + "darksalmon", + "peru", + "saddlebrown", + "black", + "moccasin", + "yellowgreen", + ] + fuel_models = [str(f) for f in fuel_dat.fuel_model.values] + + if by_litter_type: + data_dict = { + lc: fuel_dat.isel(litter_class=i)[var].values + for i, lc in enumerate(litter_classes) + } + else: + data_dict = fuel_dat[var].values + + _, ax = plt.subplots() + if by_litter_type: + bottom = np.zeros(len(fuel_models)) + for i, (litter_class, dat) in enumerate(data_dict.items()): + ax.bar( + fuel_models, + dat, + 0.5, + label=litter_class, + bottom=bottom, + color=colors[i], + ) + bottom += dat + plt.legend(loc="center left", bbox_to_anchor=(1, 0.5)) + else: + ax.bar(fuel_models, data_dict, color="darkcyan") + + box = ax.get_position() + ax.set_position([box.x0, box.y0, box.width * 0.75, box.height]) + plt.ylabel(f"{varname} ({units})", fontsize=11) + plt.xticks(rotation=90) + plt.xlabel("Fuel Model") + + if save_figs: + fig_name = os.path.join(plot_dir, f"{varname}_plot.png") + plt.savefig(fig_name) + + @staticmethod + def plot_NI_dat(fuel_dat: xr.Dataset, save_figs: bool, plot_dir: str): + """Plot output for Nesterov index + + Args: + fuel_dat (Xarray Dataset): output fuel data + save_figs (bool): whether or not to save the figures + plot_dir (str): plot directory + """ + + plt.figure() + fuel_dat.NI.plot() + plt.xlabel("Time", fontsize=11) + plt.ylabel("Nesterov Index", fontsize=11) + + if save_figs: + fig_name = os.path.join(plot_dir, "Nesterov_plot.png") + plt.savefig(fig_name) + + @staticmethod + def plot_moisture_dat(fuel_dat: xr.Dataset, save_figs: bool, plot_dir: str): + """Plot output for fuel moisture + + Args: + fuel_dat (Xarray Dataset): output fuel data + save_figs (bool): whether or not to save the figures + plot_dir (str): plot directory + """ + + plt.figure() + fuel_dat.fuel_moisture.plot(hue="fuel_model") + plt.xlabel("Time", fontsize=11) + plt.ylabel("Fuel Moisture", fontsize=11) + + if save_figs: + fig_name = os.path.join(plot_dir, "fuel_moisture_plot.png") + plt.savefig(fig_name) diff --git a/testing/functional_tests.cfg b/testing/functional_tests.cfg index f6c9c3f69c..46924c307c 100644 --- a/testing/functional_tests.cfg +++ b/testing/functional_tests.cfg @@ -11,3 +11,10 @@ test_exe = FATES_math_exe out_file = quad_out.nc use_param_file = False other_args = [] + +[fuel] +test_dir = fates_fuel_ftest +test_exe = FATES_fuel_exe +out_file = fuel_out.nc +use_param_file = True +other_args = ['../testing/test_data/BONA_datm.nc'] diff --git a/testing/run_functional_tests.py b/testing/run_functional_tests.py index 5f80a3f095..8a2078628a 100755 --- a/testing/run_functional_tests.py +++ b/testing/run_functional_tests.py @@ -36,12 +36,19 @@ # add testing subclasses here from functional_class import FunctionalTest -from functional_testing.allometry.allometry_test import AllometryTest # pylint: disable=unused-import -from functional_testing.math_utils.math_utils_test import QuadraticTest # pylint: disable=unused-import +from functional_testing.allometry.allometry_test import ( + AllometryTest, +) # pylint: disable=unused-import +from functional_testing.math_utils.math_utils_test import ( + QuadraticTest, +) # pylint: disable=unused-import +from functional_testing.fire.fuel_test import FuelTest # pylint: disable=unused-import add_cime_lib_to_path() -from CIME.utils import run_cmd_no_fail # pylint: disable=wrong-import-position,import-error,wrong-import-order +from CIME.utils import ( + run_cmd_no_fail, +) # pylint: disable=wrong-import-position,import-error,wrong-import-order # constants for this script _DEFAULT_CONFIG_FILE = "functional_tests.cfg" diff --git a/testing/testing_shr/CMakeLists.txt b/testing/testing_shr/CMakeLists.txt index 8cf4eb69c0..c295bdaf64 100644 --- a/testing/testing_shr/CMakeLists.txt +++ b/testing/testing_shr/CMakeLists.txt @@ -1,6 +1,7 @@ list(APPEND fates_sources FatesUnitTestParamReaderMod.F90 FatesUnitTestIOMod.F90 + FatesArgumentUtils.F90 ) sourcelist_to_parent(fates_sources) \ No newline at end of file diff --git a/testing/testing_shr/FatesArgumentUtils.F90 b/testing/testing_shr/FatesArgumentUtils.F90 new file mode 100644 index 0000000000..ed247fa157 --- /dev/null +++ b/testing/testing_shr/FatesArgumentUtils.F90 @@ -0,0 +1,36 @@ +module FatesArgumentUtils + + implicit none + private + + public :: command_line_arg + + contains + + function command_line_arg(arg_position) + ! DESCRIPTION: + ! Retrieves a single argument from the command line at a given position + + ! ARGUMENTS: + integer, intent(in) :: arg_position ! argument position + character(len=:), allocatable :: command_line_arg ! command argument + + ! LOCALS: + integer :: n_args ! number of arguments received + integer :: arglen ! argument length + + ! get total number of arguments + n_args = command_argument_count() + + if (n_args < arg_position) then + write(*, '(a, i2, a, i2)') "Incorrect number of arguments: ", n_args, ". Should be at least", arg_position, "." + stop + end if + + call get_command_argument(arg_position, length=arglen) + allocate(character(arglen) :: command_line_arg) + call get_command_argument(arg_position, value=command_line_arg) + + end function command_line_arg + +end module FatesArgumentUtils diff --git a/testing/testing_shr/FatesUnitTestIOMod.F90 b/testing/testing_shr/FatesUnitTestIOMod.F90 index 8f6ea1141a..20dd4f198e 100644 --- a/testing/testing_shr/FatesUnitTestIOMod.F90 +++ b/testing/testing_shr/FatesUnitTestIOMod.F90 @@ -10,7 +10,8 @@ module FatesUnitTestIOMod ! LOCALS integer, public, parameter :: type_double = 1 ! type integer, public, parameter :: type_int = 2 ! type - + integer, public, parameter :: type_char = 3 ! type + interface GetVar module procedure GetVarScalarReal module procedure GetVar1DReal @@ -26,6 +27,8 @@ module FatesUnitTestIOMod module procedure WriteVar2DReal module procedure WriteVar1DInt module procedure WriteVar2DInt + module procedure WriteVar1DChar + module procedure WriteVar2DChar end interface public :: OpenNCFile @@ -473,6 +476,8 @@ subroutine RegisterVar(ncid, var_name, dimID, type, att_names, atts, num_atts, v nc_type = NF90_DOUBLE else if (type == type_int) then nc_type = NF90_INT + else if (type == type_char) then + nc_type = NF90_CHAR else write(*, *) "Must pick correct type" stop @@ -570,5 +575,39 @@ subroutine WriteVar2DInt(ncid, varID, data) end subroutine WriteVar2DInt ! ===================================================================================== + + subroutine WriteVar1DChar(ncid, varID, data) + ! + ! DESCRIPTION: + ! Write 1D character data + ! + + ! ARGUMENTS: + integer, intent(in) :: ncid ! netcdf file id + integer, intent(in) :: varID ! variable ID + character(len=*), intent(in) :: data(:) ! data to write + + call Check(nf90_put_var(ncid, varID, data(:))) + + end subroutine WriteVar1DChar + + ! ===================================================================================== + + subroutine WriteVar2DChar(ncid, varID, data) + ! + ! DESCRIPTION: + ! Write 2D character data + ! + + ! ARGUMENTS: + integer, intent(in) :: ncid ! netcdf file id + integer, intent(in) :: varID ! variable ID + character(len=*), intent(in) :: data(:,:) ! data to write + + call Check(nf90_put_var(ncid, varID, data(:,:))) + + end subroutine WriteVar2DChar + + ! ===================================================================================== end module FatesUnitTestIOMod \ No newline at end of file diff --git a/testing/unit_testing/fire_fuel_test/CMakeLists.txt b/testing/unit_testing/fire_fuel_test/CMakeLists.txt new file mode 100644 index 0000000000..3ebdaef4ad --- /dev/null +++ b/testing/unit_testing/fire_fuel_test/CMakeLists.txt @@ -0,0 +1,6 @@ +set(pfunit_sources test_FireFuel.pf) + +add_pfunit_ctest(FireFuel + TEST_SOURCES "${pfunit_sources}" + LINK_LIBRARIES fates csm_share) + \ No newline at end of file diff --git a/testing/unit_testing/fire_fuel_test/test_FireFuel.pf b/testing/unit_testing/fire_fuel_test/test_FireFuel.pf new file mode 100644 index 0000000000..ffb6e9d524 --- /dev/null +++ b/testing/unit_testing/fire_fuel_test/test_FireFuel.pf @@ -0,0 +1,100 @@ +module test_FireFuel + ! + ! DESCRIPTION: + ! Test the FATES fuel portion of the SPITFIRE model + ! + use FatesConstantsMod, only : r8 => fates_r8 + use FatesFuelMod, only : fuel_type + use FatesFuelClassesMod, only : fuel_classes, num_fuel_classes + use funit + + implicit none + + @TestCase + type, extends(TestCase) :: TestFireFuel + type(fuel_type) :: fuel + contains + procedure :: setUp + end type TestFireFuel + + real(r8), parameter :: tol = 1.e-13_r8 + + contains + + subroutine setUp(this) + class(TestFireFuel), intent(inout) :: this + call this%fuel%Init() + end subroutine setUp + + @Test + subroutine UpdateLoading_CorrectInputOrder(this) + ! test that the calculate loading subroutine correctly sets the fuel values + class(TestFireFuel), intent(inout) :: this ! fuel test object + real(r8) :: leaf_litter = 5.0_r8 ! leaf litter [kgC/m2] + real(r8) :: twig_litter = 10.0_r8 ! twig litter [kgC/m2] + real(r8) :: sm_br_litter = 15.0_r8 ! small branch litter [kgC/m2] + real(r8) :: lg_br_litter = 20.0_r8 ! large branch litter [kgC/m2] + real(r8) :: trunk_litter = 25.0_r8 ! trunk branch litter [kgC/m2] + real(r8) :: live_grass = 30.0_r8 ! live grass [kgC/m2] + + call this%fuel%UpdateLoading(leaf_litter, twig_litter, sm_br_litter, & + lg_br_litter, trunk_litter, live_grass) + + @assertEqual(this%fuel%loading(fuel_classes%dead_leaves()), leaf_litter, tolerance=tol) + @assertEqual(this%fuel%loading(fuel_classes%twigs()), twig_litter, tolerance=tol) + @assertEqual(this%fuel%loading(fuel_classes%small_branches()), sm_br_litter, tolerance=tol) + @assertEqual(this%fuel%loading(fuel_classes%large_branches()), lg_br_litter, tolerance=tol) + @assertEqual(this%fuel%loading(fuel_classes%trunks()), trunk_litter, tolerance=tol) + @assertEqual(this%fuel%loading(fuel_classes%live_grass()), live_grass, tolerance=tol) + + end subroutine UpdateLoading_CorrectInputOrder + + @Test + subroutine SumLoading_CorrectValues(this) + ! test that the fuel is summed correctly (and ignores trunks) + class(TestFireFuel), intent(inout) :: this ! fuel test object + real(r8) :: dummy_litter = 5.0_r8 ! dummy litter value [kgC/m2] + real(r8) :: trunk_litter = 100.0_r8 ! trunk branch litter [kgC/m2] + real(r8) :: non_trunk_loading ! what non-trunk loading should be [kgC/m2] + + non_trunk_loading = dummy_litter*5.0_r8 + + call this%fuel%UpdateLoading(dummy_litter, dummy_litter, dummy_litter, & + dummy_litter, trunk_litter, dummy_litter) + + call this%fuel%SumLoading() + + @assertEqual(this%fuel%non_trunk_loading, non_trunk_loading, tolerance=tol) + + end subroutine SumLoading_CorrectValues + + @Test + subroutine CalculateFractionalLoading_CorrectValues(this) + ! test that the fractional loading is calculated correctly (and ignores trunks) + class(TestFireFuel), intent(inout) :: this ! fuel test object + real(r8) :: dummy_litter = 5.0_r8 ! dummy litter value [kgC/m2] + real(r8) :: trunk_litter = 100.0_r8 ! trunk branch litter [kgC/m2] + real(r8) :: non_trunk_loading ! what non-trunk loading should be [kgC/m2] + real(r8) :: frac_loading ! what the fractional loading should be [0-1] + integer :: i ! looping index + + non_trunk_loading = dummy_litter*float(num_fuel_classes - 1) + frac_loading = dummy_litter/non_trunk_loading + + call this%fuel%UpdateLoading(dummy_litter, dummy_litter, dummy_litter, & + dummy_litter, trunk_litter, dummy_litter) + + call this%fuel%SumLoading() + call this%fuel%CalculateFractionalLoading() + + do i = 1, num_fuel_classes + if (i /= fuel_classes%trunks()) then + @assertEqual(this%fuel%frac_loading(i), frac_loading, tolerance=tol) + else + @assertEqual(this%fuel%frac_loading(i), 0.0_r8, tolerance=tol) + end if + end do + + end subroutine CalculateFractionalLoading_CorrectValues + +end module test_FireFuel diff --git a/testing/unit_tests.cfg b/testing/unit_tests.cfg index 3e5ee536f4..179b924735 100644 --- a/testing/unit_tests.cfg +++ b/testing/unit_tests.cfg @@ -1,2 +1,5 @@ [fire_weather] test_dir = fates_fire_weather_utest + +[fire_fuel] +test_dir = fates_fire_fuel_utest diff --git a/testing/utils.py b/testing/utils.py index d1fafd5838..06fd74db0b 100644 --- a/testing/utils.py +++ b/testing/utils.py @@ -10,7 +10,9 @@ add_cime_lib_to_path() -from CIME.utils import run_cmd_no_fail # pylint: disable=wrong-import-position,import-error,wrong-import-order +from CIME.utils import ( + run_cmd_no_fail, +) # pylint: disable=wrong-import-position,import-error,wrong-import-order def round_up(num: float, decimals: int = 0) -> float: