diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index d75dbaa211..3dedd6e98c 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -3289,7 +3289,6 @@ subroutine fuse_2_patches(csite, dp, rp) 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%effect_wspeed = (dp%effect_wspeed*dp%area + rp%effect_wspeed*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%tfc_ros = (dp%tfc_ros*dp%area + rp%tfc_ros*rp%area) * inv_sum_area diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index cc6c77ded1..999bd73228 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -6,6 +6,7 @@ module FatesPatchMod use FatesConstantsMod, only : primaryland, secondaryland use FatesConstantsMod, only : n_landuse_cats use FatesConstantsMod, only : TRS_regeneration + use FatesConstantsMod, only : itrue use FatesGlobals, only : fates_log use FatesGlobals, only : endrun => fates_endrun use FatesUtilsMod, only : check_hlm_list @@ -16,6 +17,8 @@ module FatesPatchMod use FatesLitterMod, only : litter_type use PRTGenericMod, only : num_elements use PRTGenericMod, only : element_list + use PRTParametersMod, only : prt_params + use FatesConstantsMod, only : nocomp_bareground use EDParamsMod, only : nlevleaf, nclmax, maxpft use FatesConstantsMod, only : n_dbh_bins, n_dist_types use FatesConstantsMod, only : t_water_freeze_k_1atm @@ -100,6 +103,7 @@ module FatesPatchMod ! used to determine attenuation of parameters during photosynthesis real(r8) :: total_canopy_area ! area that is covered by vegetation [m2] real(r8) :: total_tree_area ! area that is covered by woody vegetation [m2] + real(r8) :: total_grass_area ! area that is covered by non-woody vegetation [m2] real(r8) :: zstar ! height of smallest canopy tree, only meaningful in "strict PPA" mode [m] ! exposed leaf area in each canopy layer, pft, and leaf layer [m2 leaf/m2 contributing crown area] @@ -213,7 +217,6 @@ module FatesPatchMod ! fire spread real(r8) :: ros_front ! rate of forward spread of fire [m/min] real(r8) :: ros_back ! rate of backward spread of fire [m/min] - real(r8) :: effect_wspeed ! windspeed modified by fraction of relative grass and tree cover [m/min] real(r8) :: tau_l ! duration of lethal heating [min] 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] @@ -241,6 +244,7 @@ module FatesPatchMod procedure :: InitRunningMeans procedure :: InitLitter procedure :: Create + procedure :: UpdateTreeGrassArea procedure :: FreeMemory procedure :: Dump procedure :: CheckVars @@ -451,7 +455,8 @@ subroutine NanValues(this) this%pft_agb_profile(:,:) = nan this%canopy_layer_tlai(:) = nan this%total_canopy_area = nan - this%total_tree_area = nan + this%total_tree_area = nan + this%total_grass_area = nan this%zstar = nan @@ -511,7 +516,6 @@ subroutine NanValues(this) this%litter_moisture(:) = nan this%ros_front = nan this%ros_back = nan - this%effect_wspeed = nan this%tau_l = nan this%fi = nan this%fire = fates_unset_int @@ -565,7 +569,8 @@ subroutine ZeroValues(this) ! LEAF ORGANIZATION this%canopy_layer_tlai(:) = 0.0_r8 - this%total_tree_area = 0.0_r8 + this%total_tree_area = 0.0_r8 + this%total_grass_area = 0.0_r8 this%zstar = 0.0_r8 this%c_stomata = 0.0_r8 @@ -574,7 +579,6 @@ subroutine ZeroValues(this) ! RADIATION this%rad_error(:) = 0.0_r8 - this%tr_soil_dir_dif(:) = 0.0_r8 this%fab(:) = 0.0_r8 this%fabi(:) = 0.0_r8 @@ -606,7 +610,6 @@ subroutine ZeroValues(this) this%litter_moisture(:) = 0.0_r8 this%ros_front = 0.0_r8 this%ros_back = 0.0_r8 - this%effect_wspeed = 0.0_r8 this%tau_l = 0.0_r8 this%fi = 0.0_r8 this%fd = 0.0_r8 @@ -763,6 +766,42 @@ end subroutine Create !=========================================================================== + subroutine UpdateTreeGrassArea(this) + ! + ! DESCRIPTION: + ! calculate and update the total tree area and grass area (by canopy) on patch + ! + + ! ARGUMENTS: + class(fates_patch_type), intent(inout) :: this ! patch object + + ! LOCALS: + 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] + + if (this%nocomp_pft_label /= nocomp_bareground) then + tree_area = 0.0_r8 + grass_area = 0.0_r8 + + currentCohort => this%tallest + do while(associated(currentCohort)) + if (prt_params%woody(currentCohort%pft) == itrue) then + tree_area = tree_area + currentCohort%c_area + else + grass_area = grass_area + currentCohort%c_area + end if + currentCohort => currentCohort%shorter + end do + + this%total_tree_area = min(tree_area, this%area) + this%total_grass_area = min(grass_area, this%area) + end if + + end subroutine UpdateTreeGrassArea + + !=========================================================================== + subroutine FreeMemory(this, regeneration_model, numpft) ! ! DESCRIPTION: @@ -913,6 +952,7 @@ subroutine Dump(this) write(fates_log(),*) 'pa%ncl_p = ',this%ncl_p write(fates_log(),*) 'pa%total_canopy_area = ',this%total_canopy_area write(fates_log(),*) 'pa%total_tree_area = ',this%total_tree_area + write(fates_log(),*) 'pa%total_grass_area = ',this%total_grass_area write(fates_log(),*) 'pa%zstar = ',this%zstar write(fates_log(),*) 'pa%solar_zenith_flag = ',this%solar_zenith_flag write(fates_log(),*) 'pa%solar_zenith_angle = ',this%solar_zenith_angle diff --git a/fire/SFFireWeatherMod.F90 b/fire/SFFireWeatherMod.F90 index e39834ad1e..11b8602fa4 100644 --- a/fire/SFFireWeatherMod.F90 +++ b/fire/SFFireWeatherMod.F90 @@ -6,27 +6,64 @@ module SFFireWeatherMod private type, abstract, public :: fire_weather + real(r8) :: fire_weather_index ! fire weather index + real(r8) :: effective_windspeed ! effective wind speed, corrected for by tree/grass cover [m/min] + contains + procedure(initialize_fire_weather), public, deferred :: Init - procedure(update_fire_weather), public, deferred :: Update + procedure(update_fire_weather), public, deferred :: UpdateIndex + procedure, public :: UpdateEffectiveWindSpeed + end type fire_weather abstract interface subroutine initialize_fire_weather(this) + import :: fire_weather + class(fire_weather), intent(inout) :: this + end subroutine initialize_fire_weather subroutine update_fire_weather(this, temp_C, precip, rh, wind) + use FatesConstantsMod, only : r8 => fates_r8 + import :: fire_weather + class(fire_weather), intent(inout) :: this real(r8), intent(in) :: temp_C real(r8), intent(in) :: precip real(r8), intent(in) :: rh real(r8), intent(in) :: wind + end subroutine update_fire_weather end interface + contains + + subroutine UpdateEffectiveWindSpeed(this, wind_speed, tree_fraction, grass_fraction, & + bare_fraction) + ! + ! DESCRIPTION: + ! Calculates effective wind speed + + ! CONSTANTS: + real(r8), parameter :: wind_atten_treed = 0.4_r8 ! wind attenuation factor for tree fraction + real(r8), parameter :: wind_atten_grass = 0.6_r8 ! wind attenuation factor for grass fraction + + ! ARGUMENTS + class(fire_weather), intent(inout) :: this ! fire weather class + real(r8), intent(in) :: wind_speed ! wind speed [m/min] + real(r8), intent(in) :: tree_fraction ! tree fraction [0-1] + real(r8), intent(in) :: grass_fraction ! grass fraction [0-1] + real(r8), intent(in) :: bare_fraction ! bare ground fraction [0-1] + + this%effective_windspeed = wind_speed*(tree_fraction*wind_atten_treed + & + (grass_fraction + bare_fraction)*wind_atten_grass) + + end subroutine UpdateEffectiveWindSpeed + end module SFFireWeatherMod \ No newline at end of file diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index f1e7de14ea..eeb37e79f7 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -9,8 +9,7 @@ module SFMainMod 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 EDTypesMod , only : numWaterMem + use FatesInterfaceTypesMod, only : hlm_masterproc ! 1= master process, 0=not master process use FatesGlobals , only : fates_log use FatesInterfaceTypesMod, only : hlm_spitfire_mode use FatesInterfaceTypesMod, only : hlm_sf_nofire_def @@ -39,13 +38,9 @@ module SFMainMod use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : carbon12_element use PRTGenericMod, only : leaf_organ - use PRTGenericMod, only : fnrt_organ use PRTGenericMod, only : sapw_organ - use PRTGenericMod, only : store_organ - use PRTGenericMod, only : repro_organ use PRTGenericMod, only : struct_organ - use PRTGenericMod, only : SetState - use FatesInterfaceTypesMod , only : numpft + use FatesInterfaceTypesMod, only : numpft use FatesAllometryMod, only : CrownDepth use FatesConstantsMod, only : nearzero @@ -56,7 +51,6 @@ module SFMainMod public :: charecteristics_of_fuel public :: rate_of_spread public :: ground_fuel_consumption - public :: wind_effect public :: area_burnt_intensity public :: crown_scorching public :: crown_damage @@ -95,7 +89,6 @@ subroutine fire_model(currentSite, bc_in) if (hlm_spitfire_mode > hlm_sf_nofire_def) then call UpdateFireWeather(currentSite, bc_in) - call wind_effect(currentSite, bc_in) call charecteristics_of_fuel(currentSite) call rate_of_spread(currentSite) call ground_fuel_consumption(currentSite) @@ -113,33 +106,41 @@ end subroutine fire_model subroutine UpdateFireWeather(currentSite, bc_in) ! ! DESCRIPTION: - ! Updates the site's fire weather index + ! 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 use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - use FatesConstantsMod, only : sec_per_day + 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 + type(bc_in_type), intent(in) :: bc_in ! LOCALS: - type(fates_patch_type), pointer :: currentPatch ! patch object - real(r8) :: temp_C ! daily averaged temperature [deg C] - real(r8) :: precip ! daily precip [mm/day] - real(r8) :: rh ! daily relative humidity [%] - real(r8) :: wind ! wind speed [m/s] - integer :: iofp ! index of oldest the fates patch + type(fates_patch_type), pointer :: currentPatch ! patch object + real(r8) :: temp_C ! daily averaged temperature [deg C] + real(r8) :: precip ! daily precip [mm/day] + real(r8) :: rh ! daily relative humidity [%] + real(r8) :: wind ! wind speed [m/s] + real(r8) :: tree_fraction ! site-level tree fraction [0-1] + 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 inpact, unless we decide to ever calculated fire weather for each 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 ! for the iofp index (i.e. the next younger patch) - if(currentPatch%nocomp_pft_label .eq. nocomp_bareground)then + if (currentPatch%nocomp_pft_label == nocomp_bareground) then currentPatch => currentPatch%younger endif @@ -149,8 +150,18 @@ subroutine UpdateFireWeather(currentSite, bc_in) rh = bc_in%relhumid24_pa(iofp) wind = bc_in%wind24_pa(iofp) + ! convert to m/min + currentSite%wind = wind*sec_per_min + ! update fire weather index - call currentSite%fireWeather%Update(temp_C, precip, rh, wind) + call currentSite%fireWeather%UpdateIndex(temp_C, precip, rh, wind) + + ! calculate site-level tree, grass, and bare fraction + call CalculateTreeGrassAreaSite(currentSite, tree_fraction, grass_fraction, bare_fraction) + + ! update effective wind speed + call currentSite%fireWeather%UpdateEffectiveWindSpeed(wind*sec_per_min, tree_fraction, & + grass_fraction, bare_fraction) end subroutine UpdateFireWeather @@ -339,108 +350,7 @@ subroutine charecteristics_of_fuel ( currentSite ) end subroutine charecteristics_of_fuel - !***************************************************************** - subroutine wind_effect ( currentSite, bc_in) - !*****************************************************************. - - ! Routine called daily from within ED within a site loop. - ! Calculates the effective windspeed based on vegetation charecteristics. - ! currentSite%wind is daily wind converted to m/min for Spitfire units - - use FatesConstantsMod, only : sec_per_min - - type(ed_site_type) , intent(inout), target :: currentSite - type(bc_in_type) , intent(in) :: bc_in - - type(fates_patch_type) , pointer :: currentPatch - type(fates_cohort_type), pointer :: currentCohort - - real(r8) :: total_grass_area ! per patch,in m2 - real(r8) :: tree_fraction ! site level. no units - real(r8) :: grass_fraction ! site level. no units - real(r8) :: bare_fraction ! site level. no units - integer :: iofp ! index of oldest fates patch - - - currentPatch => currentSite%oldest_patch - - ! If the oldest patch is a bareground patch (i.e. nocomp mode is on) use the first vegetated patch - ! for the iofp index (i.e. the next younger patch) - if(currentPatch%nocomp_pft_label .eq. nocomp_bareground)then - currentPatch => currentPatch%younger - endif - - ! note - this is a patch level temperature, which probably won't have much inpact, - ! unless we decide to ever calculated the NI for each patch. - iofp = currentPatch%patchno - currentSite%wind = bc_in%wind24_pa(iofp) * sec_per_min !Convert to m/min for SPITFIRE - - if(write_SF == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'wind24', currentSite%wind - endif - ! --- influence of wind speed, corrected for surface roughness---- - ! --- averaged over the whole grid cell to prevent extreme divergence - ! average_wspeed = 0.0_r8 - tree_fraction = 0.0_r8 - grass_fraction = 0.0_r8 - currentPatch=>currentSite%oldest_patch; - do while(associated(currentPatch)) - - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - - currentPatch%total_tree_area = 0.0_r8 - total_grass_area = 0.0_r8 - currentCohort => currentPatch%tallest - - do while(associated(currentCohort)) - if (debug) write(fates_log(),*) 'SF currentCohort%c_area ',currentCohort%c_area - if( prt_params%woody(currentCohort%pft) == itrue)then - currentPatch%total_tree_area = currentPatch%total_tree_area + currentCohort%c_area - else - total_grass_area = total_grass_area + currentCohort%c_area - endif - currentCohort => currentCohort%shorter - enddo - tree_fraction = tree_fraction + min(currentPatch%area,currentPatch%total_tree_area)/AREA - grass_fraction = grass_fraction + min(currentPatch%area,total_grass_area)/AREA - - if(debug)then - write(fates_log(),*) 'SF currentPatch%area ',currentPatch%area - write(fates_log(),*) 'SF currentPatch%total_area ',currentPatch%total_tree_area - write(fates_log(),*) 'SF total_grass_area ',tree_fraction,grass_fraction - write(fates_log(),*) 'SF AREA ',AREA - endif - - endif !nocomp_pft_label check - - currentPatch => currentPatch%younger - enddo !currentPatch loop - - !if there is a cover of more than one, then the grasses are under the trees - grass_fraction = min(grass_fraction,1.0_r8-tree_fraction) - bare_fraction = 1.0_r8 - tree_fraction - grass_fraction - if(write_sf == itrue)then - if ( hlm_masterproc == itrue ) write(fates_log(),*) 'grass, trees, bare', & - grass_fraction, tree_fraction, bare_fraction - endif - - currentPatch=>currentSite%oldest_patch; - - do while(associated(currentPatch)) - if(currentPatch%nocomp_pft_label .ne. nocomp_bareground)then - - currentPatch%total_tree_area = min(currentPatch%total_tree_area,currentPatch%area) - ! effect_wspeed in units m/min - currentPatch%effect_wspeed = currentSite%wind * (tree_fraction*0.4_r8+(grass_fraction+bare_fraction)*0.6_r8) - endif ! nocomp_pft_label check - - currentPatch => currentPatch%younger - enddo !end patch loop - - end subroutine wind_effect - - !***************************************************************** subroutine rate_of_spread ( currentSite ) !*****************************************************************. !Routine called daily from within ED within a site loop. @@ -520,8 +430,6 @@ subroutine rate_of_spread ( currentSite ) if (debug) then if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) 'SF - c ',c - if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) 'SF - currentPatch%effect_wspeed ', & - currentPatch%effect_wspeed if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) 'SF - b ',b if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) 'SF - beta_ratio ',beta_ratio if ( hlm_masterproc == itrue .and.debug) write(fates_log(),*) 'SF - e ',e @@ -530,7 +438,7 @@ subroutine rate_of_spread ( currentSite ) ! Equation A5 in Thonicke et al. 2010 ! phi_wind (unitless) ! convert current_wspeed (wind at elev relevant to fire) from m/min to ft/min for Rothermel ROS eqn - phi_wind = c * ((3.281_r8*currentPatch%effect_wspeed)**b)*(beta_ratio**(-e)) + phi_wind = c * ((3.281_r8*currentSite%fireWeather%effective_windspeed)**b)*(beta_ratio**(-e)) ! ---propagating flux---- @@ -571,7 +479,6 @@ subroutine rate_of_spread ( currentSite ) 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',currentPatch%ROS_front,phi_wind,currentPatch%effect_wspeed ! write(fates_log(),*) 'ros calcs',currentPatch%fuel_bulkd,ir,xi,eps,q_ig endif ! Equation 10 in Thonicke et al. 2010 @@ -795,16 +702,16 @@ subroutine area_burnt_intensity ( currentSite, bc_in ) write(fates_log(),*) 'SF AREA ',AREA endif - if ((currentPatch%effect_wspeed*m_per_min__to__km_per_hour) < 1._r8) then !16.67m/min = 1km/hr + if ((currentSite%fireWeather%effective_windspeed*m_per_min__to__km_per_hour) < 1._r8) then !16.67m/min = 1km/hr lb = 1.0_r8 else if (tree_fraction_patch > forest_grassland_lengthtobreadth_threshold) then !benchmark forest cover, Staver 2010 ! EQ 79 forest fuels (Canadian Forest Fire Behavior Prediction System Ont.Inf.Rep. ST-X-3, 1992) lb = (1.0_r8 + (8.729_r8 * & - ((1.0_r8 -(exp(-0.03_r8 * m_per_min__to__km_per_hour * currentPatch%effect_wspeed)))**2.155_r8))) + ((1.0_r8 -(exp(-0.03_r8 * m_per_min__to__km_per_hour * currentSite%fireWeather%effective_windspeed)))**2.155_r8))) else ! EQ 80 grass fuels (CFFBPS Ont.Inf.Rep. ST-X-3, 1992, but with a correction from an errata published within ! Information Report GLC-X-10 by Wotton et al., 2009 for a typo in CFFBPS Ont.Inf.Rep. ST-X-3, 1992) - lb = (1.1_r8*((m_per_min__to__km_per_hour * currentPatch%effect_wspeed)**0.464_r8)) + lb = (1.1_r8*((m_per_min__to__km_per_hour * currentSite%fireWeather%effective_windspeed)**0.464_r8)) endif endif diff --git a/fire/SFNesterovMod.F90 b/fire/SFNesterovMod.F90 index 85d1d94396..f8001c2552 100644 --- a/fire/SFNesterovMod.F90 +++ b/fire/SFNesterovMod.F90 @@ -1,7 +1,10 @@ module SFNesterovMod use FatesConstantsMod, only : r8 => fates_r8 + use FatesGlobals, only : endrun => fates_endrun + use FatesGlobals, only : fates_log use SFFireWeatherMod, only : fire_weather + use shr_log_mod, only : errMsg => shr_log_errMsg implicit none private @@ -11,7 +14,7 @@ module SFNesterovMod contains procedure, public :: Init => init_nesterov_fire_weather - procedure, public :: Update => update_nesterov_index + procedure, public :: UpdateIndex => update_nesterov_index end type nesterov_index @@ -29,6 +32,7 @@ subroutine init_nesterov_fire_weather(this) ! initialize values to 0.0 this%fire_weather_index = 0.0_r8 + this%effective_windspeed = 0.0_r8 end subroutine init_nesterov_fire_weather diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 95da8cbfa8..15bbb3934c 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -980,7 +980,6 @@ subroutine init_patches( nsites, sites, bc_in) currentPatch%fuel_sav = 0._r8 currentPatch%fuel_mef = 0._r8 currentPatch%ros_front = 0._r8 - currentPatch%effect_wspeed = 0._r8 currentPatch%tau_l = 0._r8 currentPatch%fuel_frac(:) = 0._r8 currentPatch%tfc_ros = 0._r8 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 14865a56bd..10ef14734c 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -481,6 +481,7 @@ module EDTypesMod ! Make public necessary subroutines and functions public :: dump_site + public :: CalculateTreeGrassAreaSite contains @@ -537,6 +538,43 @@ end subroutine ZeroMassBalFlux ! ===================================================================================== + subroutine CalculateTreeGrassAreaSite(csite, tree_fraction, grass_fraction, bare_fraction) + ! + ! DESCRIPTION: + ! Calculates total grass, tree, and bare fractions for a site + + use FatesConstantsMod, only : nocomp_bareground + + ! ARGUMENTS: + type(ed_site_type), intent(inout) :: csite ! site object + real(r8), intent(out) :: tree_fraction ! total site tree fraction + real(r8), intent(out) :: grass_fraction ! total site grass fraction + real(r8), intent(out) :: bare_fraction ! total site bare fraction + + ! LOCALS: + type(fates_patch_type), pointer :: currentPatch ! patch object + + tree_fraction = 0.0_r8 + grass_fraction = 0.0_r8 + + currentPatch => csite%oldest_patch + do while(associated(currentPatch)) + if (currentPatch%nocomp_pft_label /= nocomp_bareground) then + call currentPatch%UpdateTreeGrassArea() + tree_fraction = tree_fraction + currentPatch%total_tree_area/AREA + grass_fraction = grass_fraction + currentPatch%total_grass_area/AREA + end if + currentPatch => currentPatch%younger + end do + + ! if cover > 1.0, grasses are under the trees + grass_fraction = min(grass_fraction, 1.0_r8 - tree_fraction) + bare_fraction = 1.0_r8 - tree_fraction - grass_fraction + + end subroutine CalculateTreeGrassAreaSite + + !--------------------------------------------------------------------------------------- + subroutine dump_site(csite) type(ed_site_type),intent(in),target :: csite diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index fd3c09dfd5..1f23c18951 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -2530,6 +2530,8 @@ subroutine update_history_dyn1(this,nc,nsites,sites,bc_in) ! Nesterov index (unitless) hio_nesterov_fire_danger_si(io_si) = sites(s)%fireWeather%fire_weather_index + + hio_effect_wspeed_si(io_si) = sites(s)%fireWeather%effective_windspeed/sec_per_min ! number of ignitions [#/km2/day -> #/m2/s] hio_fire_nignitions_si(io_si) = sites(s)%NF_successful / m2_per_km2 / & @@ -2606,7 +2608,7 @@ subroutine update_history_dyn1(this,nc,nsites,sites,bc_in) sites(s)%fmort_crownarea_ustory + & sites(s)%term_crownarea_ustory * days_per_year + & sites(s)%imort_crownarea - + flux_diags_c => sites(s)%flux_diags(element_pos(carbon12_element)) hio_litter_in_si(io_si) = (sum(flux_diags_c%cwd_ag_input(:)) + & @@ -2670,7 +2672,6 @@ subroutine update_history_dyn1(this,nc,nsites,sites,bc_in) ! Update Fire Variables hio_spitfire_ros_si(io_si) = hio_spitfire_ros_si(io_si) + cpatch%ROS_front * cpatch%area * AREA_INV / sec_per_min - hio_effect_wspeed_si(io_si) = hio_effect_wspeed_si(io_si) + cpatch%effect_wspeed * cpatch%area * AREA_INV / sec_per_min 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 diff --git a/testing/unit_testing/fire_weather_test/test_FireWeather.pf b/testing/unit_testing/fire_weather_test/test_FireWeather.pf index c5111394bc..65149a46e3 100644 --- a/testing/unit_testing/fire_weather_test/test_FireWeather.pf +++ b/testing/unit_testing/fire_weather_test/test_FireWeather.pf @@ -12,9 +12,7 @@ module test_FireWeather @TestCase type, extends(TestCase) :: TestFireWeather - class(fire_weather), allocatable :: fireWeatherNesterov - contains procedure :: setUp procedure :: tearDown @@ -36,32 +34,81 @@ module test_FireWeather end subroutine tearDown @Test - subroutine zero_NI_rain(this) + subroutine UpdateIndex_OverPrecipThreshold_ZerosNI(this) ! test that over 3 mm of rain is 0.0 - class(TestFireWeather), intent(inout) :: this ! fire weather object + class(TestFireWeather), intent(inout) :: this ! fire weather test object + real(r8) :: tempC = 25.0_r8 ! temperature [degC] + real(r8) :: precip = 3.1_r8 ! precipitation [mm] + real(r8) :: rh = 10.0_r8 ! relative humidity [%] + real(r8) :: wind = 0.0_r8 ! wind speed [m/s] + + + call this%fireWeatherNesterov%UpdateIndex(tempC, precip, rh, wind) - call this%fireWeatherNesterov%Update(25.0_r8, 3.1_r8, 10.0_r8, 0.0_r8) @assertEqual(this%fireWeatherNesterov%fire_weather_index, 0.0_r8, tolerance=tol) - this%fireWeatherNesterov%fire_weather_index = 0.0_r8 - end subroutine zero_NI_rain + + end subroutine UpdateIndex_OverPrecipThreshold_ZerosNI @Test - subroutine NI_rain_min(this) - ! test that at 3 mm is over zero - class(TestFireWeather), intent(inout) :: this ! fire weather object + subroutine UpdateIndex_AtPrecipThreshold_AccumulateNI(this) + ! test that at 3 mm is over 0.0 + class(TestFireWeather), intent(inout) :: this ! fire weather test object + real(r8) :: tempC = 25.0_r8 ! temperature [degC] + real(r8) :: precip = 3.0_r8 ! precipitation [mm] + real(r8) :: rh = 10.0_r8 ! relative humidity [%] + real(r8) :: wind = 0.0_r8 ! wind speed [m/s] + + call this%fireWeatherNesterov%UpdateIndex(tempC, precip, rh, wind) - call this%fireWeatherNesterov%Update(25.0_r8, 3.0_r8, 10.0_r8, 0.0_r8) @assertGreaterThan(this%fireWeatherNesterov%fire_weather_index, 0.0_r8, tolerance=tol) - this%fireWeatherNesterov%fire_weather_index = 0.0_r8 - end subroutine NI_rain_min + + end subroutine UpdateIndex_AtPrecipThreshold_AccumulateNI @Test - subroutine NI_not_negative(this) + subroutine UpdateIndex_NegativeCalculation_ZerosNI(this) ! test that NI is not negative - class(TestFireWeather), intent(inout) :: this ! fire weather object + class(TestFireWeather), intent(inout) :: this ! fire weather test object + real(r8) :: tempC = -30.0_r8 ! temperature [degC] + real(r8) :: precip = 0.0_r8 ! precipitation [mm] + real(r8) :: rh = 99.0_r8 ! relative humidity [%] + real(r8) :: wind = 0.0_r8 ! wind speed [m/s] + + + call this%fireWeatherNesterov%UpdateIndex(tempC, precip, rh, wind) - call this%fireWeatherNesterov%Update(-30.0_r8, 0.0_r8, 99.0_r8, 0.0_r8) @assertEqual(this%fireWeatherNesterov%fire_weather_index, 0.0_r8, tolerance=tol) - end subroutine NI_not_negative + + end subroutine UpdateIndex_NegativeCalculation_ZerosNI + + @Test + subroutine UpdateEffectiveWindSpeed_ZeroWindSpeed_ZeroEffectiveWindSpeed(this) + ! test that effective wind speed is zero when wind speed is zero + class(TestFireWeather), intent(inout) :: this ! fire weather test object + real(r8) :: wind_speed = 0.0_r8 ! wind speed [m/s] + real(r8) :: tree_fraction = 0.5_r8 ! tree fraction [0-1] + real(r8) :: grass_fraction = 0.5_r8 ! grass fraction [0-1] + real(r8) :: bare_fraction = 0.0_r8 ! bare fraction [0-1] + + call this%fireWeatherNesterov%UpdateEffectiveWindSpeed(wind_speed, tree_fraction, & + grass_fraction, bare_fraction) + @assertEqual(this%fireWeatherNesterov%effective_windspeed, 0.0_r8, tolerance=tol) + + end subroutine UpdateEffectiveWindSpeed_ZeroWindSpeed_ZeroEffectiveWindSpeed + + @Test + subroutine UpdateEffectiveWindSpeed_WindSpeed_AttenuatesWindSpeed(this) + ! test that effective wind speed is less than input wind speed + class(TestFireWeather), intent(inout) :: this ! fire weather test object + real(r8) :: wind_speed = 50.0_r8 ! wind speed [m/s] + real(r8) :: tree_fraction = 0.5_r8 ! tree fraction [0-1] + real(r8) :: grass_fraction = 0.5_r8 ! grass fraction [0-1] + real(r8) :: bare_fraction = 0.0_r8 ! bare fraction [0-1] + + call this%fireWeatherNesterov%UpdateEffectiveWindSpeed(wind_speed, tree_fraction, & + grass_fraction, bare_fraction) + + @assertLessThan(this%fireWeatherNesterov%effective_windspeed, wind_speed, tolerance=tol) + + end subroutine UpdateEffectiveWindSpeed_WindSpeed_AttenuatesWindSpeed end module test_FireWeather \ No newline at end of file