diff --git a/Registry/registry.fire b/Registry/registry.fire
index a092b7bff2..aa503ac283 100644
--- a/Registry/registry.fire
+++ b/Registry/registry.fire
@@ -12,7 +12,7 @@ package fire_sfire ifire==2 - state:nfuel_ca
# level function history support
dimspec ign 2 constant=1 z i_lfn_history
-state real lfn_hist *i{ign}*j fire 1 Z i012hr "LFN_HIST" "level function history" "1"
+state real lfn_hist *i*j fire 1 Z i012hr "LFN_HIST" "level function history" "1"
state real lfn_time {ign} fire 1 - i012hr "LFN_TIME" "level function history time" "s"
@@ -32,18 +32,28 @@ state real rthfrten ikj fire 1 z hr "RTHFRTEN"
state real rqvfrten ikj fire 1 z hr "RQVFRTEN" "humidity tendency"
# diagnostics only
-state real avg_fuel_frac ij fire 1 z hr "AVG_FUEL_FRAC" "fuel remaining averaged to atmospheric grid" "1"
-state real grnhfx ij fire 1 z hr "GRNHFX" "heat flux from ground fire" "W/m^2"
-state real grnqfx ij fire 1 z hr "GRNQFX" "moisture flux from ground fire" "W/m^2"
-state real canhfx ij fire 1 z hr "CANHFX" "heat flux from crown fire" "W/m^2"
-state real canqfx ij fire 1 z hr "CANQFX" "moisture flux from crown fire" "W/m^2"
-state real uah ij fire 1 X hr "UAH" "wind at fire_wind_height" "m/s"
-state real vah ij fire 1 Y hr "VAH" "wind at fire_wind_height" "m/s"
+state real avg_fuel_frac ij fire 1 z hr "AVG_FUEL_FRAC" "fuel remaining averaged to atmospheric grid" "1"
+state real grnhfx ij fire 1 z hr "GRNHFX" "heat flux from ground fire" "W/m^2"
+state real grnqfx ij fire 1 z hr "GRNQFX" "moisture flux from ground fire" "W/m^2"
+state real canhfx ij fire 1 z hr "CANHFX" "heat flux from crown fire" "W/m^2"
+state real canqfx ij fire 1 z hr "CANQFX" "moisture flux from crown fire" "W/m^2"
+state real uah ij fire 1 X hr "UAH" "wind at fire_wind_height" "m/s"
+state real vah ij fire 1 Y hr "VAH" "wind at fire_wind_height" "m/s"
+state real grnhfx_fu ij fire 1 z r "GRNHFX_FU" "heat flux from ground fire (feedback unsensitive)" "W/m^2"
+state real grnqfx_fu ij fire 1 z r "GRNQFX_FU" "moisture flux from ground fire (feedback unsensitive)" "W/m^2"
# sfire variables on fire grid
# (also using tign_g,zs,z_at_w,dz8w,nfuel_cat,fluxes,zsf)
#
state real lfn *i*j fire 1 z hr "LFN" "level function" "1"
+state real lfn_0 *i*j fire 1 z r "LFN_0" "level function for time integration, step 0" "1"
+state real lfn_1 *i*j fire 1 z r "LFN_1" "level function for time integration, step 1" "1"
+state real lfn_2 *i*j fire 1 z r "LFN_2" "level function for time integration, step 2" "1"
+state real lfn_s0 *i*j fire 1 z r "LFN_S0" "level set sign function from LSM integration" "1"
+state real lfn_s1 *i*j fire 1 z r "LFN_S1" "level set function for reinitialization integration" "1"
+state real lfn_s2 *i*j fire 1 z r "LFN_S2" "level set function for reinitialization integration" "1"
+state real lfn_s3 *i*j fire 1 z r "LFN_S3" "level set function for reinitialization integration" "1"
+#
state real fuel_frac *i*j fire 1 z hr "FUEL_FRAC" "fuel remaining" "1"
state real fire_area *i*j fire 1 z hr "FIRE_AREA" "fraction of cell area on fire" "1"
state real uf *i*j fire 1 z hr "UF" "fire wind" "m/s"
@@ -52,7 +62,10 @@ state real fgrnhfx *i*j fire 1 z hr "FGRNHFX"
state real fgrnqfx *i*j fire 1 z hr "FGRNQFX" "moisture flux from ground fire" "W/m^2"
state real fcanhfx *i*j fire 1 z hr "FCANHFX" "heat flux from crown fire" "W/m^2"
state real fcanqfx *i*j fire 1 z hr "FCANQFX" "moisture flux from crown fire" "W/m^2"
-state real ros *i*j fire 1 z hr "ROS" "rate of spread" "m/s"
+state real ros *i*j fire 1 z r "ROS" "rate of spread" "m/s"
+state real burnt_area_dt *i*j fire 1 z hr "BURNT_AREA_DT" "fraction of cell area burnt on current dt" "-"
+state real flame_length *i*j fire 1 z hr "FLAME_LENGTH" "fire flame length" "m"
+state real ros_front *i*j fire 1 z hr "ROS_FRONT" "rate of spread at fire front" "m/s"
# constant data arrays
state real fxlong *i*j fire 1 z ihr "FXLONG" "longitude of midpoints of fire cells" "degrees"
@@ -64,14 +77,16 @@ state real phiwc *i*j fire 1 z hr "PHIWC"
state real r_0 *i*j fire 1 z hr "R_0" "fuel"
state real fgip *i*j fire 1 z hr "FGIP" "fuel"
state real ischap *i*j fire 1 z hr "ISCHAP" "fuel"
+state real fz0 *i*j fire 1 z ihr "FZ0" "roughness length of fire cells" "m"
+state real iboros *i*j fire 1 z hr "IBOROS" "fire intensity over rate of spread" "kJ/m^2"
#
# fire configure namelist variables
#
#
rconfig integer ifire namelist,fire max_domains 0
-rconfig integer fire_boundary_guard namelist,fire max_domains 2 - "fire_boundary_guard" "cells to stop when fire close to domain boundary"
-# ignition for sfire
+rconfig integer fire_boundary_guard namelist,fire max_domains 8 - "fire_boundary_guard" "cells to stop when fire close to domain boundary"
+# ignition for WRF-Fire
rconfig integer fire_num_ignitions namelist,fire max_domains 0 - "fire_num_ignitions" "number of ignition lines"
rconfig real fire_ignition_ros1 namelist,fire max_domains 0.01 - "fire_ignition_ros1" "rate of spread during ignition" "m/s"
rconfig real fire_ignition_start_lon1 namelist,fire max_domains 0. - "fire_ignition_start_long1" "long coord of start of ignition line" "deg"
@@ -148,13 +163,12 @@ rconfig integer fire_fuel_cat namelist,fire max_domains
# sfire switches
rconfig integer fire_print_msg namelist,fire max_domains 0 - "fire_write_msg" "write fire statistics, 0 no writes, 1+ for more" ""
rconfig integer fire_print_file namelist,fire max_domains 0 - "fire_write_file" "write fire output text files, 0 no writes, 1+ for more" ""
-# method selection`
+# method selection
rconfig integer fire_fuel_left_method namelist,fire max_domains 1 - "fire_fuel_left_method" "1 or 2, compute fuel_left" ""
rconfig integer fire_fuel_left_irl namelist,fire max_domains 2 - "fire_fuel_left_irl" "submesh to compute fuel lwft, even, at least 2" ""
rconfig integer fire_fuel_left_jrl namelist,fire max_domains 2 - "fire_fuel_left_jrl" "submesh to compute fuel lwft, even, at least 2" ""
-rconfig real fire_back_weight namelist,fire max_domains 0.5 - "fire_back_weight" "RK timestepping coefficient, 0=forward, 0.5=Heun" "1"
rconfig integer fire_grows_only namelist,fire max_domains 1 - "fire_grows_only" "if >0 level set function cannot increase = fire can only grow" "1"
-rconfig integer fire_upwinding namelist,fire max_domains 3 - "fire_upwinding" "upwind normal spread: 1=standard, 2=godunov, 3=eno, 4=sethian" "1"
+rconfig integer fire_upwinding namelist,fire max_domains 9 - "fire_upwinding" "upwind normal spread: 1=standard, 2=godunov, 3=eno, 4=sethian, 5=2nd-order, 6=WENO3, 7=WENO5, 8=hybrid WENO3/ENO1, 9=hybrid WENO5/ENO1" "1"
rconfig integer fire_upwind_split namelist,fire max_domains 0 - "fire_upwind_split" "1=upwind advection separately from normal direction spread" "1"
rconfig real fire_viscosity namelist,fire max_domains 0.4 - "fire_viscosity" "artificial viscosity in level set method" "1"
rconfig real fire_lfn_ext_up namelist,fire max_domains 1.0 - "fire_lfn_ext_up" "0.=extend level set function at boundary by reflection, 1.=always up" "1"
@@ -200,10 +214,32 @@ rconfig real sfc_vegfra namelist,fire max_domains
rconfig real sfc_canwat namelist,fire max_domains 0 - "sfc_canwat" "canopy water" ""
rconfig integer sfc_ivgtyp namelist,fire max_domains 18 - "sfc_ivgtyp" "dominant vegetation category in the LSM scheme" ""
rconfig integer sfc_isltyp namelist,fire max_domains 7 - "sfc_isltyp" "dominant soil category in the LSM scheme" ""
+# namelist options for new WRF-Fire capabilities DME
+rconfig logical fire_lsm_reinit namelist,fire max_domains .true. - "flag to activate reinitialization of level set method"
+rconfig integer fire_lsm_reinit_iter namelist,fire max_domains 1 - "number of iterations for the reinitialization PDE"
+rconfig integer fire_upwinding_reinit namelist,fire max_domains 4 - "numerical scheme (space) for reinitialization PDE: 1=WENO3, 2=WENO5, 3=hybrid WENO3-ENO1, 4=hybrid WENO5-ENO1"
+rconfig logical fire_is_real_perim namelist,fire max_domains .false. - ".false. = point/line ignition, .true. = observed perimeter"
+rconfig integer fire_lsm_band_ngp namelist,fire max_domains 4 - "number of grid points around lfn=0 that WENO5/3 is used (ENO1 elsewhere), for fire_upwinding_reinit=4,5 and fire_upwinding=8,9 options"
+rconfig logical fire_lsm_zcoupling namelist,fire max_domains .false. - "flag to activate reference velocity at a different height from fire_wind_height"
+rconfig real fire_lsm_zcoupling_ref namelist,fire max_domains 50. - "reference height from wich u at fire_wind_hegiht is calculated using a logarithmic profile" "m"
+rconfig real fire_tracer_smoke namelist,fire max_domains 0.02 - "parts per unit of burned fuel becoming smoke (tracer_opt=3)" "g_smoke/kg_air"
+rconfig real fire_viscosity_bg namelist,fire max_domains 0.4 - "fire_viscosity_bg" "artificial viscosity in the near-front region" "1"
+rconfig real fire_viscosity_band namelist,fire max_domains 0.5 - "fire_viscosity_band" "number of times the hybrid advection band to transition from fire_viscosity_bg to fire_viscosity" "1"
+rconfig integer fire_viscosity_ngp namelist,fire max_domains 2 - "number of grid points around lfn=0 where low artificial viscosity is used = fire_viscosity_bg"
+rconfig real fire_slope_factor namelist,fire max_domains 1.0 - "slope correction factor" "-"
+# tracers for WRF-Fire
+state real fire_smoke ikjftb tracer 1 - irhusdf=(bdy_interp:dt) "fire_smoke" "fire_smoke" "g_smoke/kg_air"
+package tracer_fire tracer_opt==3 - tracer:fire_smoke
#
# Fire halo descriptions
#
-halo HALO_FIRE_LFN dyn_em 24:lfn
+halo HALO_FIRE_LFN dyn_em 48:lfn
+halo HALO_FIRE_LFN_0 dyn_em 48:lfn_0
+halo HALO_FIRE_LFN_1 dyn_em 48:lfn_1
+halo HALO_FIRE_LFN_2 dyn_em 48:lfn_2
+halo HALO_FIRE_LFN_S1 dyn_em 48:lfn_s1
+halo HALO_FIRE_LFN_S2 dyn_em 48:lfn_s2
+halo HALO_FIRE_LFN_S3 dyn_em 48:lfn_s3
halo HALO_FIRE_TIGN dyn_em 8:tign_g
halo HALO_FIRE_HT dyn_em 8:ht
halo HALO_FIRE_PHB dyn_em 8:phb
@@ -218,5 +254,4 @@ halo HALO_FIRE_FUEL dyn_em 8:fuel_frac,fuel_time,bbb,betafl,phiwc,r_0,fgip,
# ----------------------------------------
# end fire variables and configuration
# ----------------------------------------
-
##
diff --git a/dyn_em/module_initialize_fire.F b/dyn_em/module_initialize_fire.F
index 7fe293caf9..a3518ee294 100644
--- a/dyn_em/module_initialize_fire.F
+++ b/dyn_em/module_initialize_fire.F
@@ -445,11 +445,6 @@ SUBROUTINE init_domain_rk ( grid &
ifms,ifme,jfms,jfme, &
ifts,ifte,jfts,jfte, &
grid%fxlong,grid%fxlat )
- call set_ideal_coord( grid%dx,grid%dy, &
- ids,ide,jds,jde, &
- ims,ime,jms,jme, &
- its,ite,jts,jte, &
- grid%xlong,grid%xlat )
! set terrain height
diff --git a/phys/module_fr_fire_atm.F b/phys/module_fr_fire_atm.F
index 47cba70b59..9980343fd9 100644
--- a/phys/module_fr_fire_atm.F
+++ b/phys/module_fr_fire_atm.F
@@ -10,9 +10,69 @@ module module_fr_fire_atm
use module_model_constants, only: cp,xlv
use module_fr_fire_util
+use module_state_description, only: num_tracer
+use module_state_description, only: p_fire_smoke
contains
+! DME subroutine for passive tracers
+subroutine add_fire_tracer_emissions( &
+ tracer_opt,dt,dx,dy, &
+ ifms,ifme,jfms,jfme, &
+ ifts,ifte,jtfs,jfte, &
+ ids,ide,kds,kde,jds,jde, &
+ ims,ime,kms,kme,jms,jme, &
+ its,ite,kts,kte,jts,jte, &
+ rho,dz8w, &
+ burnt_area_dt,fgip, &
+ tracer,fire_tracer_smoke &
+)
+
+implicit none
+! arguments
+integer,intent(in)::tracer_opt
+real,intent(in)::fire_tracer_smoke
+real,intent(in)::dt,dx,dy
+integer,intent(in)::ifms,ifme,jfms,jfme,ifts,ifte,jtfs,jfte,ids,ide,kds,kde,jds,jde,ims,ime,kms,kme,jms,jme,its,ite,kts,kte,jts,jte
+real,intent(in)::rho(ims:ime,kms:kme,jms:jme),dz8w(ims:ime,kms:kme,jms:jme)
+real,intent(in),dimension(ifms:ifme,jfms:jfme)::burnt_area_dt,fgip
+real,intent(inout)::tracer(ims:ime,kms:kme,jms:jme,num_tracer)
+! local
+integer::isz1,jsz1,isz2,jsz2,ir,jr
+integer::i,j,ibase,jbase,i_f,ioff,j_f,joff
+real::avgw,emis,conv
+
+isz1 = ite-its+1
+jsz1 = jte-jts+1
+isz2 = ifte-ifts+1
+jsz2 = jfte-jtfs+1
+ir=isz2/isz1
+jr=jsz2/jsz1
+avgw = 1.0/(ir*jr)
+
+do j=max(jds+1,jts),min(jte,jde-2)
+ jbase=jtfs+jr*(j-jts)
+ do i=max(ids+1,its),min(ite,ide-2)
+ ibase=ifts+ir*(i-its)
+ do joff=0,jr-1
+ j_f=joff+jbase
+ do ioff=0,ir-1
+ i_f=ioff+ibase
+ if (num_tracer >0)then
+ emis=avgw*fire_tracer_smoke*burnt_area_dt(i_f,j_f)*fgip(i_f,j_f)*1000/(rho(i,kts,j)*dz8w(i,kts,j)) ! g_smoke/kg_air
+ tracer(i,kts,j,p_fire_smoke)=tracer(i,kts,j,p_fire_smoke)+emis
+ endif
+ enddo
+ enddo
+ enddo
+enddo
+
+end subroutine add_fire_tracer_emissions
+
+!
+!***
+!
+
SUBROUTINE fire_tendency( &
ids,ide, kds,kde, jds,jde, & ! dimensions
ims,ime, kms,kme, jms,jme, &
diff --git a/phys/module_fr_fire_core.F b/phys/module_fr_fire_core.F
index 44da2d5c77..e6c9f38363 100644
--- a/phys/module_fr_fire_core.F
+++ b/phys/module_fr_fire_core.F
@@ -92,7 +92,7 @@ subroutine ignite_fire( ifds,ifde,jfds,jfde, & ! fire domain
lfn,tign,ignited)
implicit none
-!*** purpose: ignite a circular/line fire
+!*** purpose: ignite a circular/line/prescribed fire
!*** description
! ignite fire in the region within radius r from the line (sx,sy) to (ex,ey).
@@ -142,21 +142,19 @@ subroutine ignite_fire( ifds,ifde,jfds,jfde, & ! fire domain
end_time = ignition_line%end_time! ignition time for the end poin from simulation start (s)
radius = ignition_line%radius ! all within this radius ignites immediately
ros = ignition_line%ros ! rate of spread
-
-
tos = radius/ros ! time of spread to the given radius
st = start_time ! the start time of ignition considered in this time step
et = min(end_ts,end_time) ! the end time of the ignition segment in this time step
! this should be called whenever (start_ts, end_ts) \subset (start_time, end_time + tos)
if(start_ts>et+tos .or. end_ts0.) then
ignited=ignited+1 ! count
endif
- if(lfn(i,j)>0. .and. .not. lfn_new > 0.) then
+ if(lfn(i,j)>0. .and. .not. lfn_new > 0.) then
tign(i,j)=time_ign + d/ros ! newly ignited now
if(fire_print_msg.ge.3)then
!$OMP CRITICAL(FIRE_CORE_CRIT)
@@ -382,9 +380,6 @@ subroutine fuel_left(&
integer::m,omp_get_thread_num
-call check_mesh_2dim(its-1,ite+1,jts-1,jte+1,ims,ime,jms,jme)
-call check_mesh_2dim(its,ite,jts,jte,ifs,ife,jfs,jfe)
-
! refinement
ir=fuel_left_irl
jr=fuel_left_jrl
@@ -1273,20 +1268,44 @@ subroutine matmatp(A,mA,nA,B,mB,nB,C,mC,nC,nrow,ncolumn,nP)
!
#endif
-subroutine prop_ls( id, & ! for debug
- ids,ide,jds,jde, & ! domain dims
- ims,ime,jms,jme, & ! memory dims
- ips,ipe,jps,jpe, & ! patch - nodes owned by this process
- its,ite,jts,jte, & ! tile dims
- ts,dt,dx,dy, & ! scalars in
- tbound, & ! scalars out
- lfn_in,lfn_out,tign,ros, & ! arrays inout
- fp &
+subroutine prop_ls_rk3(id, &
+ ifds,ifde,jfds,jfde, &
+ ifms,ifme,jfms,jfme, &
+ ifps,ifpe,jfps,jfpe, &
+ ifts,ifte,jfts,jfte, &
+ ts,dt,dx,dy, &
+ tbound, &
+ lfn_in, &
+ lfn_0,lfn_1,lfn_2, &
+ lfn_out,tign,ros, &
+ fp, &
+ grid, &
+ ids,ide,jds,jde,kds,kde, &
+ ims,ime,jms,jme,kms,kme, &
+ ips,ipe,jps,jpe,kps,kpe &
)
+
+
+#ifdef DM_PARALLEL
+ USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks
+ USE module_comm_dm , ONLY : halo_fire_lfn_1_sub, halo_fire_lfn_2_sub, halo_fire_lfn_0_sub
+#endif
+USE module_domain , only: domain
+
implicit none
!*** purpose: advance level function in time
+!***************************************************************************************!
+!*** Level-set method recoded for true mpi communications and added ***!
+!*** 3rd-order Runge-Kutta time integration and 3rd/5th-order WENO advection schemes ***!
+!*** Implemented by: Domingo Munoz-Esparza (NCAR/RAL, April 2016) ***!
+!*** Reference: D. Munoz-Esparza, B. Kosovic, P. Jimenez, J. Coen: "An accurate ***!
+!*** fire-spread algorithm in the Weather Research and Forecasting model using the ***!
+!*** level-set method", Journal of Advances in Modeling Earth Systems, 2018 ***!
+!*** doi:AAA
+!***************************************************************************************!
+
!*** description
!
! Propagation of closed curve by a level function method. The level function
@@ -1329,47 +1348,28 @@ subroutine prop_ls( id, & ! for debug
! lfn_in,lfn_out inout,out the level set function at nodes
! tign inout the ignition time at nodes
-! The dimensions are cell-based, the nodal value is associated with the south-west corner.
-! The whole computation is on domain indices ids:ide+1,jds:jde+1.
-!
-! The region where new lfn and tign are computed is the tile its:ite,jts:jte
-! except when the tile is at domain upper boundary, an extra band of points is added:
-! if ite=ide then region goes up to ite+1, if jte=jde then region goes up to jte+1.
-
-! The time step requires values from 2 rows of nodes beyond the region except when at the
-! domain boundary one-sided derivatives are used. This is implemented by extending the input
-! beyond the domain boundary so sufficient memory bounds must be allocated.
-! The update on all tiles can be done in parallel. To avoid the race condition (different regions
-! of the same array updated by different threads), the in and out versions of the
-! arrays lft and tign are distinct. If the time step dt is larger
-! that the returned tbound, the routine should be called again with timestep td<=tbound, and then
-! having distinct inputs and outputs comes handy.
-
!*** calls
!
-! tend_ls
+! DME added for halo update
+type(domain) , target :: grid
+integer, intent(in):: ids,ide,jds,jde,kds,kde, &
+ ims,ime,jms,jme,kms,kme, &
+ ips,ipe,jps,jpe,kps,kpe
!
-integer,intent(in)::id,ims,ime,jms,jme,ids,ide,jds,jde,its,ite,jts,jte,ips,ipe,jps,jpe
-real,dimension(ims:ime,jms:jme),intent(inout)::lfn_in,tign
-real,dimension(ims:ime,jms:jme),intent(out)::lfn_out,ros
+integer,intent(in)::id,ifms,ifme,jfms,jfme,ifds,ifde,jfds,jfde,ifts,ifte,jfts,jfte,ifps,ifpe,jfps,jfpe
+real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_in,tign
+real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_1,lfn_2,lfn_0
+real,dimension(ifms:ifme,jfms:jfme),intent(out)::lfn_out,ros
real,intent(in)::dx,dy,ts,dt
real,intent(out)::tbound
type(fire_params),intent(in)::fp
-!*** local
-! arrays
-#define IMTS its-1
-#define IMTE ite+1
-#define JMTS jts-1
-#define JMTE jte+1
-real,dimension(IMTS:IMTE,JMTS:JMTE):: tend, lfn1 ! region-sized with halo
+real,dimension(ifms:ifme,jfms:jfme):: tend
! scalars
-real::grad2,rr,tbound2,a,a1 ! a=0 euler, a=0.5 heun
+real::grad2,rr,tbound2,tbound3
real::gradx,grady,aspeed,err,aerr,time_now
-integer::ihs,ihe,jhs,jhe
-integer::ihs2,ihe2,jhs2,jhe2
integer::i,j,its1,ite1,jts1,jte1,k,kk,id1
character(len=128)::msg
integer::nfirenodes,nfireline
@@ -1383,98 +1383,110 @@ subroutine prop_ls( id, & ! for debug
! f90 intrinsic function
intrinsic max,min,sqrt,nint,epsilon,tiny,huge
-
+
!*** executable
!$OMP CRITICAL(FIRE_CORE_CRIT)
-write(msg,'(a8,i5,a6,i5,3(a1,i5))')'prop_ls:',id,' tile ',its,':',ite,',',jts,':',jte
+write(msg,'(a8,i5,a6,i5,3(a1,i5))')'prop_ls:',id,' tile ',ifts,':',ifte,',',jfts,':',jfte
!$OMP END CRITICAL(FIRE_CORE_CRIT)
call message(msg)
- a=fire_back_weight ! from module_fr_fire_util
- a1=1. - a
-
- ! tend = F(lfn)
-
- ihs2=max(its-2,ids) ! need lfn two beyond the tile but not outside the domain
- ihe2=min(ite+2,ide)
- jhs2=max(jts-2,jds)
- jhe2=min(jte+2,jde)
+!!!!!!!!!!!!!!!!!!!!!!
+! Runge-Kutta step 1 !
+!!!!!!!!!!!!!!!!!!!!!!
- ihs=max(its-1,ids) ! compute tend one beyond the tile but not outside the domain
- ihe=min(ite+1,ide)
- jhs=max(jts-1,jds)
- jhe=min(jte+1,jde)
+ do j=jfts,jfte
+ do i=ifts,ifte
+ lfn_0(i,j) = lfn_in(i,j)
+ enddo
+ enddo
-#ifdef DEBUG_OUT
- call write_array_m(ihs,ihe,jhs,jhe,ims,ime,jms,jme,lfn_in,'lfn_in',id)
+#ifdef DM_PARALLEL
+# include "HALO_FIRE_LFN_0.inc"
#endif
- ! check array dimensions
- call check_mesh_2dim(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme)
- call print_2d_stats(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme, &
- lfn_in,'prop_ls: lfn in')
-
- ! NOTE: tend_ls will extrapolate to one node strip at domain boundaries
- ! so that it can compute gradient at domain boundaries.
- ! To avoid copying, lfn_in is declared inout.
- ! At tile boundaries that are not domain boundaries values of lfn_in two nodes
- ! outside of the tile are needed.
id1 = id ! for debug prints
if(id1.ne.0)id1=id1+1000
- call tend_ls( id1, &
- ims,ime,jms,jme, & ! memory dims for lfn_in
- IMTS,IMTE,JMTS,JMTE, & ! memory dims for tend
- ids,ide,jds,jde, & ! domain dims - where lfn exists
- ips,ipe,jps,jpe, & ! patch - nodes owned by this process
- ihs,ihe,jhs,jhe, & ! where tend computed
- ims,ime,jms,jme, & ! memory dims for ros
- its,ite,jts,jte, & ! tile dims - where to set ros
- ts,dt,dx,dy, & ! scalars in
- lfn_in, & ! arrays in
- tbound, & ! scalars out
- tend, ros, & ! arrays out
- fp & ! params
+ call tend_ls(id1, &
+ ifds,ifde,jfds,jfde, & ! domain dims
+ ifps,ifpe,jfps,jfpe, & ! patch dims
+ ifts,ifte,jfts,jfte, & ! tile dims
+ ifms,ifme,jfms,jfme, & ! memory dims
+ ts,dt,dx,dy, & ! scalars in
+ lfn_0, & ! arrays in
+ tbound, & ! scalars out
+ tend, ros, & ! arrays out
+ fp & ! params
)
-#ifdef DEBUG_OUT
- call write_array_m(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE,tend,'tend1',id)
+ do j=jfts,jfte
+ do i=ifts,ifte
+ lfn_1(i,j) = lfn_0(i,j) + (dt/3.0)*tend(i,j)
+ enddo
+ enddo
+
+#ifdef DM_PARALLEL
+# include "HALO_FIRE_LFN_1.inc"
#endif
- ! Euler method, the half-step, same region as ted
- do j=jhs,jhe
- do i=ihs,ihe
- lfn1(i,j) = lfn_in(i,j) + dt*tend(i,j)
+!!!!!!!!!!!!!!!!!!!!!!
+! Runge-Kutta step 2 !
+!!!!!!!!!!!!!!!!!!!!!!
+
+ if(id1.ne.0)id1=id1+1000
+ call tend_ls(id1, &
+ ifds,ifde,jfds,jfde, &
+ ifps,ifpe,jfps,jfpe, &
+ ifts,ifte,jfts,jfte, &
+ ifms,ifme,jfms,jfme, &
+ ts+dt,dt,dx,dy, &
+ lfn_1, &
+ tbound2, &
+ tend,ros, &
+ fp &
+)
+
+ do j=jfts,jfte
+ do i=ifts,ifte
+ lfn_2(i,j) = lfn_0(i,j) + (dt/2.0)*tend(i,j)
enddo
- enddo
-
- call print_2d_stats(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE, &
- lfn1,'prop_ls: lfn1')
- ! tend = F(lfn1) on the tile (not beyond)
+ enddo
+
+#ifdef DM_PARALLEL
+# include "HALO_FIRE_LFN_2.inc"
+#endif
+
+!!!!!!!!!!!!!!!!!!!!!!
+! Runge-Kutta step 3 !
+!!!!!!!!!!!!!!!!!!!!!!
if(id1.ne.0)id1=id1+1000
- call tend_ls( id1,&
- IMTS,IMTE,JMTS,JMTE, & ! memory dims for lfn
- IMTS,IMTE,JMTS,JMTE, & ! memory dims for tend
- ids,ide,jds,jde, & ! domain dims - where lfn exists
- ips,ipe,jps,jpe, & ! patch - nodes owned by this process
- its,ite,jts,jte, & ! tile dims - where is tend computed
- ims,ime,jms,jme, & ! memory dims for ros
- its,ite,jts,jte, & ! tile dims - where is ros set
- ts+dt,dt,dx,dy, & ! scalars in
- lfn1, & ! arrays in
- tbound2, & ! scalars out
- tend,ros, & ! arrays out
- fp &
+ call tend_ls(id1, &
+ ifds,ifde,jfds,jfde, &
+ ifps,ifpe,jfps,jfpe, &
+ ifts,ifte,jfts,jfte, &
+ ifms,ifme,jfms,jfme, &
+ ts+dt,dt,dx,dy, &
+ lfn_2, &
+ tbound3, &
+ tend,ros, &
+ fp &
)
-#ifdef DEBUG_OUT
- call write_array_m(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'tend2',id)
+ do j=jfts,jfte
+ do i=ifts,ifte
+ lfn_out(i,j) = lfn_0(i,j) + dt*tend(i,j)
+ lfn_2(i,j) = lfn_out(i,j) ! lfn_2=lfn_out (needed for reinitialization purposes)
+ enddo
+ enddo
+
+#ifdef DM_PARALLEL
+# include "HALO_FIRE_LFN_2.inc"
#endif
- call print_2d_stats(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'prop_ls: tend2')
-
- tbound=min(tbound,tbound2)
+! CFL check
+
+ tbound=min(tbound,tbound2,tbound3)
!$OMP CRITICAL(FIRE_CORE_CRIT)
write(msg,'(a,f10.2,4(a,f7.2))')'prop_ls: time',ts,' dt=',dt,' bound',min(tbound,999.99), &
@@ -1489,146 +1501,398 @@ subroutine prop_ls( id, & ! for debug
call message(msg)
endif
- ! combine lfn1 and lfn_in + dt*tend -> lfn_out
-
- do j=jts,jte
- do i=its,ite
- lfn_out(i,j) = a1*lfn1(i,j) + a*(lfn_in(i,j) + dt*tend(i,j))
- lfn_out(i,j) = min(lfn_out(i,j),lfn_in(i,j)) ! fire area can only increase
- enddo
- enddo
-
- ! compute ignition time by interpolation
- ! the node was not burning at start but it is burning at end
- ! interpolate from the level functions at start and at end
- ! lfn_in is the level set function value at time ts
- ! lfn_out is the level set function value at time ts+dt
- ! 0 is the level set function value at time tign(i,j)
- ! thus assuming the level function is approximately linear =>
- ! tign(i,j)= ts + ((ts + td) - ts) * lfn_in / (lfn_in - lfn_out)
- ! = ts + dt * lfn_in / (lfn_in - lfn_out)
-
- time_now=ts+dt
- time_now = time_now + abs(time_now)*epsilon(time_now)*2.
- do j=jts,jte
- do i=its,ite
- ! interpolate the cross-over time
- if (.not. lfn_out(i,j)>0 .and. lfn_in(i,j)>0)then
- tign(i,j) = ts + dt * lfn_in(i,j) / (lfn_in(i,j) - lfn_out(i,j))
- endif
- ! set the ignition time outside of burning region
- if(lfn_out(i,j)>0.)tign(i,j)=time_now
+end subroutine prop_ls_rk3
+
+!
+!*****************************
+!
+
+subroutine reinit_ls_rk3(id, &
+ ifts,ifte,jfts,jfte, &
+ ifms,ifme,jfms,jfme, &
+ ifds,ifde,jfds,jfde, &
+ ifps,ifpe,jfps,jfpe, &
+ ts,dt,dx,dy, &
+ lfn_in, &
+ lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3, &
+ lfn_out,tign, &
+ grid, &
+ ids,ide,jds,jde,kds,kde, &
+ ims,ime,jms,jme,kms,kme, &
+ ips,ipe,jps,jpe,kps,kpe &
+ )
+
+
+#ifdef DM_PARALLEL
+ USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks
+ USE module_comm_dm , ONLY : halo_fire_lfn_s1_sub,halo_fire_lfn_s2_sub,halo_fire_lfn_s3_sub
+#endif
+USE module_domain , only: domain
+USE module_configure, only: grid_config_rec_type
+
+implicit none
+
+!*************************************************************************************!
+!*** Level-set function reinitialization following: ***!
+!*** Sussman, Smereka, Osher. Journal of Computational Physics 114, 146-159 (1994) ***!
+!*** implemented by Domingo Munoz-Esparza (NCAR/RAL, April 2016) ***!
+!*** Reference: D. Munoz-Esparza, B. Kosovic, P. Jimenez, J. Coen: "An accurate ***!
+!*** fire-spread algorithm in the Weather Research and Forecasting model using the ***!
+!*** level-set method", Journal of Advances in Modeling Earth Systems, 2018 ***!
+!*** doi:AAA
+!*************************************************************************************!
+
+!*** purpose: reinitialize the level-set function
+
+type(domain) , target :: grid
+integer, intent(in):: ids,ide,jds,jde,kds,kde, &
+ ims,ime,jms,jme,kms,kme, &
+ ips,ipe,jps,jpe,kps,kpe
+
+integer,intent(in)::ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,id
+integer,intent(in)::ifds,ifde,jfds,jfde,ifps,ifpe,jfps,jfpe
+real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_in,tign
+real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3
+real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_out
+real,intent(in)::dx,dy,ts,dt
+real::dt_s
+real,dimension(ifts:ifte,jfts:jfte):: tend_1,tend_2,tend_3
+real::diffLx,diffLy,diffRx,diffRy,diff2x,diff2y,grad,time_now
+integer::nts,i,j,k,kk
+intrinsic epsilon
+character(len=128)::msg
+integer::itso,iteo,jtso,jteo
+real::threshold_HLl,threshold_HLu
+integer,parameter::bdy_eno1=10
+!*** executable
+
+threshold_HLl=-fire_lsm_band_ngp*dx
+threshold_HLu=fire_lsm_band_ngp*dx
+
+! Define S0 based on current lfn values
+
+ do j=jfts,jfte
+ do i=ifts,ifte
+ lfn_s0(i,j) = lfn_2(i,j)/sqrt(lfn_2(i,j)**2.0+dx**2.0)
+ lfn_s3(i,j) = lfn_2(i,j)
enddo
enddo
+#ifdef DM_PARALLEL
+# include "HALO_FIRE_LFN_S3.inc"
+#endif
+
+ call continue_at_boundary(1,1,fire_lfn_ext_up, & !extend by extrapolation but never down
+ ifms,ifme,jfms,jfme, & ! memory dims
+ ifds,ifde,jfds,jfde, & ! domain
+ ifps,ifpe,jfps,jfpe, & ! patch
+ ifts,ifte,jfts,jfte, & ! tile
+ itso,iteo,jtso,jteo, & ! where set now
+ lfn_s3)
+
- ! check local speed error and stats
- ! may not work correctly in parallel
- ! init stats
- nfirenodes=0
- nfireline=0
- sum_err=0.
- min_err=rmax
- max_err=rmin
- sum_aerr=0.
- min_aerr=rmax
- max_aerr=rmin
- its1=its+1
- jts1=jts+1
- ite1=ite-1
- jte1=jte-1
- ! loop over right inside of the domain
- ! cannot use values outside of the domain, would have to reflect and that
- ! would change lfn_out
- ! cannot use values outside of tile, not synchronized yet
- ! so in parallel mode the statistics is not accurate
- do j=jts1,jte1
- do i=its1,ite1
- if(lfn_out(i,j)>0.0)then ! a point out of burning region
- if(lfn_out(i+1,j)<=0.or.lfn_out(i,j+1)<=0.or. & ! neighbor in burning region
- lfn_out(i-1,j)<=0.or.lfn_out(i,j-1)<=0)then ! point next to fireline
- gradx=(lfn_out(i+1,j)-lfn_out(i-1,j))/(2.0*dx) ! central differences
- grady=(lfn_out(i,j+1)-lfn_out(i,j-1))/(2.0*dy)
- grad2=sqrt(gradx*gradx+grady*grady)
- aspeed = (lfn_in(i,j)-lfn_out(i,j))/(dt*max(grad2,rmin))
- rr = speed_func(gradx,grady,dx,dy,i,j,fp)
- err=aspeed-rr
- sum_err=sum_err+err
- min_err=min(min_err,err)
- max_err=max(max_err,err)
- aerr=abs(err)
- sum_aerr=sum_aerr+aerr
- min_aerr=min(min_aerr,aerr)
- max_aerr=max(max_aerr,aerr)
- nfireline=nfireline+1
- endif
- else
- nfirenodes=nfirenodes+1
- endif
+dt_s=0.01*dx
+do nts=1,fire_lsm_reinit_iter ! iterate to solve to steady state (1 iter each time step is enoguh)
+ ! Runge-Kutta step 1 !
+
+ call advance_ls_reinit( &
+ ifms,ifme,jfms,jfme, &
+ ifds,ifde,jfds,jfde, &
+ ifts,ifte,jfts,jfte, &
+ dx,dy,dt_s,bdy_eno1,threshold_HLu, &
+ lfn_s0,lfn_s3,lfn_s3,lfn_s1,1.0/3.0) ! sign funcition, initial ls, current stage ls, next stage advanced ls, RK coefficient
+
+#ifdef DM_PARALLEL
+# include "HALO_FIRE_LFN_S1.inc"
+#endif
+
+ call continue_at_boundary(1,1,fire_lfn_ext_up, &
+ ifms,ifme,jfms,jfme, &
+ ifds,ifde,jfds,jfde, &
+ ifps,ifpe,jfps,jfpe, &
+ ifts,ifte,jfts,jfte, &
+ itso,iteo,jtso,jteo, &
+ lfn_s1)
+
+ ! Runge-Kutta step 2 !
+
+ call advance_ls_reinit( &
+ ifms,ifme,jfms,jfme, &
+ ifds,ifde,jfds,jfde, &
+ ifts,ifte,jfts,jfte, &
+ dx,dy,dt_s,bdy_eno1,threshold_HLu, &
+ lfn_s0,lfn_s3,lfn_s1,lfn_s2,1.0/2.0)
+
+#ifdef DM_PARALLEL
+# include "HALO_FIRE_LFN_S2.inc"
+#endif
+
+ call continue_at_boundary(1,1,fire_lfn_ext_up, &
+ ifms,ifme,jfms,jfme, &
+ ifds,ifde,jfds,jfde, &
+ ifps,ifpe,jfps,jfpe, &
+ ifts,ifte,jfts,jfte, &
+ itso,iteo,jtso,jteo, &
+ lfn_s2)
+
+ ! Runge-Kutta step 3 !
+
+ call advance_ls_reinit( &
+ ifms,ifme,jfms,jfme, &
+ ifds,ifde,jfds,jfde, &
+ ifts,ifte,jfts,jfte, &
+ dx,dy,dt_s,bdy_eno1,threshold_HLu, &
+ lfn_s0,lfn_s3,lfn_s2,lfn_s3,1.0)
+
+#ifdef DM_PARALLEL
+# include "HALO_FIRE_LFN_S3.inc"
+#endif
+
+ call continue_at_boundary(1,1,fire_lfn_ext_up, &
+ ifms,ifme,jfms,jfme, &
+ ifds,ifde,jfds,jfde, &
+ ifps,ifpe,jfps,jfpe, &
+ ifts,ifte,jfts,jfte, &
+ itso,iteo,jtso,jteo, &
+ lfn_s3)
+
+enddo ! end iterations for steady state solution of reinitialization PDE
+
+ do j=jfts,jfte
+ do i=ifts,ifte
+ lfn_out(i,j)=lfn_s3(i,j) ! assing to lfn_out the reinitialized level-set function
+ lfn_out(i,j)=min(lfn_out(i,j),lfn_in(i,j)) ! fire area can only increase
enddo
enddo
-!$OMP CRITICAL(FIRE_CORE_CRIT)
- write(msg,'(2(a,i6,f8.4))')'prop_ls: nodes burning',nfirenodes, &
- (100.*nfirenodes)/((ite1-its1+1)*(jte1-jts1+1)),'% next to fireline',nfireline
-!$OMP END CRITICAL(FIRE_CORE_CRIT)
- call message(msg)
- if(nfireline>0)then
- call print_stat_line('speed error',its1,ite1,jts1,jte1,min_err,max_err,sum_err/nfireline)
- call print_stat_line('abs(speed error)',its1,ite1,jts1,jte1,min_aerr,max_aerr,sum_aerr/nfireline)
- endif
- ! check if the fire did not get to the domain boundary
- do k=-1,1,2
- ! do kk=1,(boundary_guard*(ide-ids+1))/100 ! in %
- do kk=1,boundary_guard ! measured in cells
- i=ids+k*kk
- if(i.ge.its.and.i.le.ite)then
- do j=jts,jte
- if(lfn_out(i,j)<=0.)goto 9
- enddo
- endif
- enddo
- ! do kk=1,(boundary_guard*(jde-jds+1))/100
- do kk=1,boundary_guard ! measured in cells
- j=jds+k*kk
- if(j.ge.jts.and.j.le.jte)then
- do i=its,ite
- if(lfn_out(i,j)<=0.)goto 9
- enddo
- endif
+end subroutine reinit_ls_rk3
+
+!
+!*****************************
+!
+
+subroutine advance_ls_reinit(ifms,ifme,jfms,jfme, &
+ ifds,ifde,jfds,jfde, &
+ ifts,ifte,jfts,jfte, &
+ dx,dy,dt_s,bdy_eno1,threshold_HLu, &
+ lfn_s0,lfn_ini,lfn_curr,lfn_fin,rk_coeff)
+
+! Calculates right-hand-side forcing and advances a RK-stage the level-set reinitialization PDE
+! Domingo Munoz-Esparza, NCAR/RAL
+! Reference: D. Munoz-Esparza, B. Kosovic, P. Jimenez, J. Coen: "An accurate
+! fire-spread algorithm in the Weather Research and Forecasting model using the
+! level-set method", Journal of Advances in Modeling Earth Systems, 2018
+! doi:AAA
+
+
+implicit none
+
+integer,intent(in)::ifms,ifme,jfms,jfme,ifts,ifte,jfts,jfte,ifds,ifde,jfds,jfde,bdy_eno1
+real,dimension(ifms:ifme,jfms:jfme),intent(in)::lfn_s0,lfn_ini,lfn_curr
+real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_fin
+real,intent(in)::dx,dy,dt_s,threshold_HLu,rk_coeff
+!
+integer::i,j
+real::diffLx,diffLy,diffRx,diffRy,diff2x,diff2y,grad,tend_r
+!
+
+ do j=jfts,jfte
+ do i=ifts,ifte
+ if (i.lt.ifds+bdy_eno1 .OR. i.gt.ifde-bdy_eno1 .OR. j.lt.jfds+bdy_eno1 .OR. j.gt.jfde-bdy_eno1) then !
+ diffLx=(lfn_curr(i,j)-lfn_curr(i-1,j))/dx
+ diffLy=(lfn_curr(i,j)-lfn_curr(i,j-1))/dy
+ diffRx=(lfn_curr(i+1,j)-lfn_curr(i,j))/dx
+ diffRy=(lfn_curr(i,j+1)-lfn_curr(i,j))/dy
+ diff2x=select_eno(diffLx,diffRx)
+ diff2y=select_eno(diffLy,diffRy)
+ else
+ select case(fire_upwinding_reinit)
+ case(1)
+ diff2x=select_4th(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j))
+ diff2y=select_4th(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2))
+ diff2x=select_weno3(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j),lfn_s0(i,j)*diff2x)
+ diff2y=select_weno3(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2),lfn_s0(i,j)*diff2y)
+ case(2)
+ diff2x=select_4th(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j))
+ diff2y=select_4th(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2))
+ diff2x=select_weno5(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i-3,j),lfn_curr(i+1,j),lfn_curr(i+2,j),lfn_curr(i+3,j),lfn_s0(i,j)*diff2x)
+ diff2y=select_weno5(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j-3),lfn_curr(i,j+1),lfn_curr(i,j+2),lfn_curr(i,j+3),lfn_s0(i,j)*diff2y)
+ case(3)
+ if (lfn_curr(i,j).lt.threshold_HLu) then
+ diff2x=select_4th(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j))
+ diff2y=select_4th(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2))
+ diff2x=select_weno3(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j),lfn_s0(i,j)*diff2x)
+ diff2y=select_weno3(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2),lfn_s0(i,j)*diff2y)
+ else
+ diffLx=(lfn_curr(i,j)-lfn_curr(i-1,j))/dx
+ diffLy=(lfn_curr(i,j)-lfn_curr(i,j-1))/dy
+ diffRx=(lfn_curr(i+1,j)-lfn_curr(i,j))/dx
+ diffRy=(lfn_curr(i,j+1)-lfn_curr(i,j))/dy
+ diff2x=select_eno(diffLx,diffRx)
+ diff2y=select_eno(diffLy,diffRy)
+ endif
+ case(4)
+ if (lfn_curr(i,j).lt.threshold_HLu) then
+ diff2x=select_4th(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j))
+ diff2y=select_4th(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2))
+ diff2x=select_weno5(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i-3,j),lfn_curr(i+1,j),lfn_curr(i+2,j),lfn_curr(i+3,j),lfn_s0(i,j)*diff2x)
+ diff2y=select_weno5(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j-3),lfn_curr(i,j+1),lfn_curr(i,j+2),lfn_curr(i,j+3),lfn_s0(i,j)*diff2y)
+ else
+ diffLx=(lfn_curr(i,j)-lfn_curr(i-1,j))/dx
+ diffLy=(lfn_curr(i,j)-lfn_curr(i,j-1))/dy
+ diffRx=(lfn_curr(i+1,j)-lfn_curr(i,j))/dx
+ diffRy=(lfn_curr(i,j+1)-lfn_curr(i,j))/dy
+ diff2x=select_eno(diffLx,diffRx)
+ diff2y=select_eno(diffLy,diffRy)
+ endif
+ case default
+ if (lfn_curr(i,j).lt.threshold_HLu) then
+ diff2x=select_4th(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j))
+ diff2y=select_4th(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2))
+ diff2x=select_weno5(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i-3,j),lfn_curr(i+1,j),lfn_curr(i+2,j),lfn_curr(i+3,j),lfn_s0(i,j)*diff2x)
+ diff2y=select_weno5(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j-3),lfn_curr(i,j+1),lfn_curr(i,j+2),lfn_curr(i,j+3),lfn_s0(i,j)*diff2y)
+ else
+ diffLx=(lfn_curr(i,j)-lfn_curr(i-1,j))/dx
+ diffLy=(lfn_curr(i,j)-lfn_curr(i,j-1))/dy
+ diffRx=(lfn_curr(i+1,j)-lfn_curr(i,j))/dx
+ diffRy=(lfn_curr(i,j+1)-lfn_curr(i,j))/dy
+ diff2x=select_eno(diffLx,diffRx)
+ diff2y=select_eno(diffLy,diffRy)
+ endif
+ end select
+ endif
+ grad=sqrt(diff2x*diff2x+diff2y*diff2y)
+ tend_r=lfn_s0(i,j)*(1.0-grad)
+ lfn_fin(i,j)=lfn_ini(i,j)+(dt_s*rk_coeff)*tend_r
enddo
enddo
- goto 10
-9 continue
-!$OMP CRITICAL(FIRE_CORE_CRIT)
- write(msg,'(a,i2,a,2i8)')'prop_ls: fire',boundary_guard, &
- ' cells from domain boundary at node ',i,j
-!$OMP END CRITICAL(FIRE_CORE_CRIT)
- call message(msg)
- call crash('prop_ls: increase the fire region')
-10 continue
- call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme, &
- lfn_out,'prop_ls: lfn out')
+end subroutine advance_ls_reinit
+
+!
+!*****************************
+!
+
+subroutine tign_update(ifts,ifte,jfts,jfte, & ! tile dims
+ ifms,ifme,jfms,jfme, & ! memory dims
+ ifds,jfds,ifde,jfde, & ! domain dims
+ ts,dt, & ! scalars in
+ lfn_in,lfn_out,tign & ! arrays inout
+ )
+
+implicit none
+
+integer,intent(in)::ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,ifds,jfds,ifde,jfde
+real,dimension(ifms:ifme,jfms:jfme),intent(inout)::tign
+real,dimension(ifms:ifme,jfms:jfme),intent(in)::lfn_in,lfn_out
+real,intent(in)::ts,dt
+real::time_now
+integer::i,j,k,kk
+intrinsic epsilon
+character(len=128)::msg
+!*** executable
+
+! compute ignition time by interpolation
+! the node was not burning at start but it is burning at end
+! interpolate from the level functions at start and at end
+! lfn_in==lfn is the level set function value at time ts
+! lfn_out is the level set function value at time ts+dt (after reinitialization)
+! 0 is the level set function value at time tign(i,j)
+! thus assuming the level function is approximately linear =>
+! tign(i,j)= ts + ((ts + td) - ts) * lfn_in / (lfn_in - lfn_out)
+! = ts + dt * lfn_in / (lfn_in - lfn_out)
+
+time_now=ts+dt
+time_now = time_now + abs(time_now)*epsilon(time_now)*2.
+do j=jfts,jfte
+ do i=ifts,ifte
+ ! interpolate the cross-over time
+ if (.not. lfn_out(i,j)>0 .and. lfn_in(i,j)>0)then
+ tign(i,j) = ts + dt * lfn_in(i,j) / (lfn_in(i,j) - lfn_out(i,j))
+ endif
+ ! set the ignition time outside of burning region
+ if(lfn_out(i,j)>0.)tign(i,j)=time_now
+ enddo
+enddo
+! stop simulation if fire is within boundary_guard grid points from the domain boundaries
+do j=jfts,jfte
+ if (j.le.boundary_guard .or. j.gt.(jfde-boundary_guard)) then
+ do i=ifts,ifte
+ if (lfn_out(i,j).lt.0.) then
+ WRITE(msg,*)'j-boundary reached'
+ CALL wrf_debug(0,msg)
+ WRITE(msg,*)'i,j,lfn_out=',i,j,lfn_out(i,j)
+ CALL wrf_debug(0,msg)
+ call crash('wrf: SUCCESS COMPLETE WRF. Fire has reached domain boundary.')
+ endif
+ enddo
+ endif
+enddo
+do i=ifts,ifte
+ if (i.le.boundary_guard .or. i.gt.(ifde-boundary_guard)) then
+ do j=jfts,jfte
+ if (lfn_out(i,j).lt.0.) then
+ WRITE(msg,*)'i-boundary reached'
+ CALL wrf_debug(0,msg)
+ WRITE(msg,*)'i,j,lfn_out=',i,j,lfn_out(i,j)
+ CALL wrf_debug(0,msg)
+ call crash('wrf: SUCCESS COMPLETE WRF. Fire has reached domain boundary.')
+ endif
+ enddo
+ endif
+enddo
+
+call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme, &
+ lfn_out,'prop_ls: lfn out')
-end subroutine prop_ls
+end subroutine tign_update
+
+!
+!*****************************
+!
+
+subroutine calc_flame_length(ifts,ifte,jfts,jfte, &
+ ifms,ifme,jfms,jfme, &
+ ros,iboros,flame_length,ros_fl,fire_area)
+
+! Diagnoses flame length
+! Domingo Munoz-Esparza, NCAR/RAL
+
+implicit none
+
+integer,intent(in)::ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme
+real,dimension(ifms:ifme,jfms:jfme),intent(in)::ros,iboros,fire_area
+real,dimension(ifms:ifme,jfms:jfme),intent(out)::flame_length,ros_fl
+! local
+integer::i,j
+
+do j=jfts,jfte
+ do i=ifts,ifte
+ if (fire_area(i,j).gt.0.0 .AND. fire_area(i,j).lt.1.0) then
+ flame_length(i,j)=0.0775*(iboros(i,j)*ros(i,j))**0.46 ! flame length according to Byram (1959)
+ ros_fl(i,j)=ros(i,j)
+ else
+ flame_length(i,j)=0.0
+ ros_fl(i,j)=0.0
+ endif
+ enddo
+enddo
+
+end subroutine calc_flame_length
!
!*****************************
!
subroutine tend_ls( id, &
- lims,lime,ljms,ljme, & ! memory dims for lfn
- tims,time,tjms,tjme, & ! memory dims for tend
- ids,ide,jds,jde, & ! domain - nodes where lfn defined
- ips,ipe,jps,jpe, & ! patch - nodes owned by this process
- ints,inte,jnts,jnte, & ! region - nodes where tend computed
- ims,ime,jms,jme, & ! memory dims for ros
- its,ite,jts,jte, & ! tile dims - where is ros set
- t,dt,dx,dy, & ! scalars in
- lfn, & ! arrays in
- tbound, & ! scalars out
- tend, ros, & ! arrays out
- fp &
+ ids,ide,jds,jde, & ! domain - nodes where lfn defined
+ ips,ipe,jps,jpe, & ! patch - nodes owned by this process
+ its,ite,jts,jte, & ! region - nodes where tend computed
+ ims,ime,jms,jme, & ! memory dims for ros
+ t,dt,dx,dy, & ! scalars in
+ lfn, & ! arrays in
+ tbound, & ! scalars out
+ tend, ros, & ! arrays out
+ fp &
)
implicit none
@@ -1636,15 +1900,15 @@ subroutine tend_ls( id, &
! compute the right hand side of the level set equation
!*** arguments
-integer,intent(in)::id,lims,lime,ljms,ljme,tims,time,tjms,tjme
+integer,intent(in)::id
integer,intent(in)::ims,ime,jms,jme,its,ite,jts,jte
-integer, intent(in)::ids,ide,jds,jde,ints,inte,jnts,jnte,ips,ipe,jps,jpe
-real,intent(in)::t ! time
-real,intent(in)::dt,dx,dy ! mesh step
-real,dimension(lims:lime,ljms:ljme),intent(inout)::lfn ! level set function
-real,intent(out)::tbound ! max allowed time step
-real,dimension(tims:time,tjms:tjme),intent(out)::tend ! tendency (rhs of the level set pde)
-real,dimension(ims:ime,jms:jme),intent(out)::ros ! rate of spread
+integer, intent(in)::ids,ide,jds,jde,ips,ipe,jps,jpe
+real,intent(in)::t ! time
+real,intent(in)::dt,dx,dy ! mesh step
+real,dimension(ims:ime,jms:jme),intent(inout)::lfn ! level set function
+real,dimension(ims:ime,jms:jme),intent(out)::tend ! tendency (rhs of the level set pde)
+real,dimension(ims:ime,jms:jme),intent(out)::ros ! rate of spread
+real,intent(out)::tbound ! max allowed time step
type(fire_params),intent(in)::fp
!*** local
@@ -1660,41 +1924,35 @@ subroutine tend_ls( id, &
!intrinsic epsilon
real, parameter:: zero=0.,one=1.,tol=100*eps, &
safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
-
+!
+real :: diff2xn, diff2yn
+real :: a_valor, signo_x, signo_y
+real :: threshold_HLl,threshold_HLu
+real :: threshold_AV,fire_viscosity_var
+integer,parameter :: bdy_eno1=10
! f90 intrinsic function
intrinsic max,min,sqrt,nint,tiny,huge
-#ifdef DEBUG_OUT
-real,dimension(tims:time,tjms:tjme)::rra,grada,speeda,tanphia
-#endif
-
!*** executable
-
- ! check array dimensions
- call check_mesh_2dim(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme)
- call check_mesh_2dim(ints,inte,jnts,jnte,tims,time,tjms,tjme)
+
+threshold_HLl=-fire_lsm_band_ngp*dx
+threshold_HLu=fire_lsm_band_ngp*dx
+threshold_AV=fire_viscosity_ngp*dx
call continue_at_boundary(1,1,fire_lfn_ext_up, & !extend by extrapolation but never down
- lims,lime,ljms,ljme, & ! memory dims
- ids,ide,jds,jde, & ! domain - nodes where lfn defined
- ips,ipe,jps,jpe, & ! patch - nodes owned by this process
- ints,inte,jnts,jnte, & ! tile - nodes update by this thread
- itso,iteo,jtso,jteo, & ! where set now
- lfn) ! array
-
- call print_2d_stats(itso,iteo,jtso,jteo,lims,lime,ljms,ljme, &
- lfn,'tend_ls: lfn cont')
-
-#ifdef DEBUG_OUT
- call write_array_m(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme,lfn,'tend_lfn_in',id)
-#endif
-
+ ims,ime,jms,jme, & ! memory dims
+ ids,ide,jds,jde, & ! domain
+ ips,ipe,jps,jpe, & ! patch
+ its,ite,jts,jte, & ! tile
+ itso,iteo,jtso,jteo, & ! where set now
+ lfn) ! field
+
tbound=0
- do j=jnts,jnte
- do i=ints,inte
+ do j=jts,jte
+ do i=its,ite
! one sided differences
diffRx = (lfn(i+1,j)-lfn(i,j))/dx
diffLx = (lfn(i,j)-lfn(i-1,j))/dx
@@ -1702,47 +1960,95 @@ subroutine tend_ls( id, &
diffLy = (lfn(i,j)-lfn(i,j-1))/dy
diffCx = diffLx+diffRx ! TWICE CENTRAL DIFFERENCE
diffCy = diffLy+diffRy
-
- !upwinding - select right or left derivative
- select case(fire_upwinding)
- case(0) ! none
- grad=sqrt(diffCx**2 + diffCy**2)
- case(1) ! standard
- diff2x=select_upwind(diffLx,diffRx)
- diff2y=select_upwind(diffLy,diffRy)
- grad=sqrt(diff2x*diff2x + diff2y*diff2y)
- case(2) ! godunov per osher/fedkiw
- diff2x=select_godunov(diffLx,diffRx)
- diff2y=select_godunov(diffLy,diffRy)
- grad=sqrt(diff2x*diff2x + diff2y*diff2y)
- case(3) ! eno
- diff2x=select_eno(diffLx,diffRx)
- diff2y=select_eno(diffLy,diffRy)
- grad=sqrt(diff2x*diff2x + diff2y*diff2y)
- case(4) ! Sethian - twice stronger pushdown of bumps
- grad=sqrt(max(diffLx,0.)**2+min(diffRx,0.)**2 &
- + max(diffLy,0.)**2+min(diffRy,0.)**2)
- case default
- grad=0.
- end select
-
- ! normal direction, from central differences
- scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps)
- nvx=diffCx/scale
- nvy=diffCy/scale
-
+
+ if (i.lt.ids+bdy_eno1 .OR. i.gt.ide-bdy_eno1 .OR. j.lt.jds+bdy_eno1 .OR. j.gt.jde-bdy_eno1) then ! DME: use eno1 near domain boundaries
+ diff2x=select_eno(diffLx,diffRx)
+ diff2y=select_eno(diffLy,diffRy)
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ else
+ select case(fire_upwinding)
+ case(0) ! none
+ grad=sqrt(diffCx**2 + diffCy**2)
+ case(1) ! standard
+ diff2x=select_upwind(diffLx,diffRx)
+ diff2y=select_upwind(diffLy,diffRy)
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ case(2) ! godunov per osher/fedkiw
+ diff2x=select_godunov(diffLx,diffRx)
+ diff2y=select_godunov(diffLy,diffRy)
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ case(3) ! eno
+ diff2x=select_eno(diffLx,diffRx)
+ diff2y=select_eno(diffLy,diffRy)
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ case(4) ! Sethian - twice stronger pushdown of bumps
+ grad=sqrt(max(diffLx,0.)**2+min(diffRx,0.)**2 &
+ + max(diffLy,0.)**2+min(diffRy,0.)**2)
+ case(5) ! 2nd order (DME)
+ diff2x=select_2nd(rr,dx,lfn(i,j),lfn(i-1,j),lfn(i+1,j))
+ diff2y=select_2nd(rr,dy,lfn(i,j),lfn(i,j-1),lfn(i,j+1))
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ case(6) ! 3rd-order WENO (DME)
+ a_valor=select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))*fp%vx(i,j)+ &
+ select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))*fp%vy(i,j)
+ signo_x=a_valor*select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))
+ signo_y=a_valor*select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))
+ diff2x=select_weno3(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j),signo_x)
+ diff2y=select_weno3(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2),signo_y)
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ case(7) ! 5th-order WENO (DME)
+ a_valor=select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))*fp%vx(i,j)+ &
+ select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))*fp%vy(i,j)
+ signo_x=a_valor*select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))
+ signo_y=a_valor*select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))
+ diff2x=select_weno5(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i-3,j),lfn(i+1,j),lfn(i+2,j),lfn(i+3,j),signo_x)
+ diff2y=select_weno5(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j-3),lfn(i,j+1),lfn(i,j+2),lfn(i,j+3),signo_y)
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ case(8) ! high-order local band method WENO3/ENO1 (DME)
+ if (abs(lfn(i,j)).lt.threshold_HLu) then
+ a_valor=select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))*fp%vx(i,j)+ &
+ select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))*fp%vy(i,j)
+ signo_x=a_valor*select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))
+ signo_y=a_valor*select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))
+ diff2x=select_weno3(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j),signo_x)
+ diff2y=select_weno3(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2),signo_y)
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ else
+ diff2x=select_eno(diffLx,diffRx)
+ diff2y=select_eno(diffLy,diffRy)
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ endif
+ case(9) ! high-order local band method WENO5/ENO1 (DME)
+ if (abs(lfn(i,j)).lt.threshold_HLu) then
+ a_valor=select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))*fp%vx(i,j)+ &
+ select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))*fp%vy(i,j)
+ signo_x=a_valor*select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))
+ signo_y=a_valor*select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))
+ diff2x=select_weno5(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i-3,j),lfn(i+1,j),lfn(i+2,j),lfn(i+3,j),signo_x)
+ diff2y=select_weno5(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j-3),lfn(i,j+1),lfn(i,j+2),lfn(i,j+3),signo_y)
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ else
+ diff2x=select_eno(diffLx,diffRx)
+ diff2y=select_eno(diffLy,diffRy)
+ grad=sqrt(diff2x*diff2x + diff2y*diff2y)
+ endif
+ case default
+ grad=0.
+ end select
+ endif
+! DME use the same discretization scheme from fire_upwinding option to calculate the normals
+ scale=sqrt(grad**2.0+eps)
+ nvx=diff2x/scale
+ nvy=diff2y/scale
+
! wind speed in direction of spread
speed = fp%vx(i,j)*nvx + fp%vy(i,j)*nvy
-
-
! get rate of spread from wind speed and slope
-
call fire_ros(ros_base,ros_wind,ros_slope, &
nvx,nvy,i,j,fp)
+ rr=ros_base + ros_wind + fire_slope_factor*ros_slope
- rr=ros_base + ros_wind + ros_slope
if(fire_grows_only.gt.0)rr=max(rr,0.)
-
! set ros for output
if(i.ge.its.and.i.le.ite.and.j.ge.jts.and.j.le.jte)ros(i,j)=rr
@@ -1796,35 +2102,20 @@ subroutine tend_ls( id, &
tbound = max(tbound,rr*(abs(diff2x)/dx+abs(diff2y)/dy)/grad)
endif
- ! add numerical viscosity
- te=te + fire_viscosity*abs(rr)*((diffRx-diffLx)+(diffRy-diffLy))
-
+ ! DME: adds spatially variable artificial viscosity
+ if (abs(lfn(i,j)).lt.threshold_AV .AND. (i.gt.ids+bdy_eno1 .and. i.lt.ide-bdy_eno1) .AND. (j.gt.jds+bdy_eno1 .and. j.lt.jde-bdy_eno1)) then
+ fire_viscosity_var=fire_viscosity_bg
+ elseif (abs(lfn(i,j)).ge.threshold_AV .AND. abs(lfn(i,j)).lt.threshold_AV*(1.0+fire_viscosity_band) .AND. (i.gt.ids+10 .and. i.lt.ide-10) .AND. (j.gt.jds+10 .and. j.lt.jde-10)) then
+ fire_viscosity_var=min(fire_viscosity_bg+(fire_viscosity-fire_viscosity_bg)*(abs(lfn(i,j))-threshold_AV)/(fire_viscosity_band*threshold_AV),fire_viscosity)
+ else
+ fire_viscosity_var=fire_viscosity
+ endif
+ te=te + fire_viscosity_var*abs(rr)*((diffRx-diffLx)+(diffRy-diffLy))
tend(i,j)=te
-#ifdef DEBUG_OUT
- rra(i,j)=rr
- grada(i,j)=grad
- speeda(i,j)=speed
- tanphia(i,j)=tanphi
-#endif
- !write(msg,*)i,j,grad,dzdx,dzdy
- !call message(msg)
-
- !if(abs(te)>1e4)then
- ! write(msg,'(a,2i5,e12.3)')'tend_ls: tend out of bounds at ',i,j,te
- ! call crash(msg)
- !endif
enddo
enddo
-#ifdef DEBUG_OUT
- call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,rra,'rr',id)
- call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,grada,'grad',id)
- call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,speeda,'speed',id)
- call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,tanphia,'tanphi',id)
- call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,tend,'tend',id)
-#endif
-
- call print_2d_stats(ints,inte,jnts,jnte,tims,time,tjms,tjme, &
+ call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme, &
tend,'tend_ls: tend out')
! the final CFL bound
@@ -1835,6 +2126,7 @@ end subroutine tend_ls
!
!**************************
!
+!
real function select_upwind(diffLx,diffRx)
implicit none
@@ -1855,7 +2147,6 @@ end function select_upwind
!**************************
!
-
real function select_godunov(diffLx,diffRx)
implicit none
real, intent(in):: diffLx, diffRx
@@ -1906,6 +2197,178 @@ end function select_eno
!**************************
!
+
+real function select_2nd(ros,dx,lfn_i,lfn_im1,lfn_ip1)
+implicit none
+real, intent(in):: lfn_i, lfn_im1, lfn_ip1
+real, intent(in):: ros, dx
+real:: diff2x_p, diff2x_m
+
+! 2nd-order advection scheme in the x,y-direction (DME)
+
+diff2x_p=0.
+diff2x_m=0.
+diff2x_p=(lfn_ip1+lfn_i)/(2.*dx)
+diff2x_m=(lfn_i+lfn_im1)/(2.*dx)
+select_2nd=diff2x_p-diff2x_m
+end function select_2nd
+
+!
+!**************************
+!
+
+
+real function select_4th(dx,lfn_i,lfn_im1,lfn_im2,lfn_ip1,lfn_ip2)
+implicit none
+real, intent(in):: lfn_i,lfn_im1,lfn_im2,lfn_ip1,lfn_ip2
+real, intent(in):: dx
+real:: diff2x_p,diff2x_m
+
+! 4th-order advection scheme in the x,y-direction (DME)
+
+diff2x_p=0.
+diff2x_m=0.
+diff2x_p=(7.0*lfn_ip1+7.0*lfn_i-lfn_ip2-lfn_im1)/(12.*dx)
+diff2x_m=(7.0*lfn_i+7.0*lfn_im1-lfn_ip1-lfn_im2)/(12.*dx)
+select_4th=diff2x_p-diff2x_m
+end function select_4th
+
+!
+!**************************
+!
+
+
+real function select_weno3(dx,lfn_it,lfn_im1t,lfn_im2t,lfn_ip1t,lfn_ip2t,uf)
+implicit none
+real, intent(in)::lfn_it,lfn_im1t,lfn_im2t,lfn_ip1t,lfn_ip2t,uf
+real:: lfn_i,lfn_im1,lfn_im2,lfn_ip1,lfn_ip2
+real, intent(in):: dx
+real:: flux_p,flux_m
+real, parameter:: gamma1=1./3.,gamma2=2./3.,tol=1e-6
+real:: w1,w2,w1t,w2t,beta1,beta2
+real:: fh_1,fh_2
+
+! 3rd-order advection WENO scheme in the x,y-direction (DME)
+
+flux_p=0.
+flux_m=0.
+
+if (uf .ge. 0.0 ) then
+ lfn_i=lfn_it
+ lfn_im1=lfn_im1t
+ lfn_im2=lfn_im2t
+ lfn_ip1=lfn_ip1t
+else
+ lfn_i=lfn_it
+ lfn_im1=lfn_ip1t
+ lfn_im2=lfn_ip2t
+ lfn_ip1=lfn_im1t
+endif
+
+! numerical flux at i,j+1/2 face
+fh_1=-0.5*lfn_im1+1.5*lfn_i
+fh_2=0.5*lfn_i+0.5*lfn_ip1
+beta1=(lfn_i-lfn_im1)**2
+beta2=(lfn_ip1-lfn_i)**2
+w1t=gamma1/(beta1+tol)**2
+w2t=gamma2/(beta2+tol)**2
+w1=w1t/(w1t+w2t)
+w2=w2t/(w1t+w2t)
+flux_p=w1*fh_1+w2*fh_2
+! numerical flux at i,j-1/2 face
+fh_1=-0.5*lfn_im2+1.5*lfn_im1
+fh_2=0.5*lfn_im1+0.5*lfn_i
+beta1=(lfn_im1-lfn_im2)**2
+beta2=(lfn_i-lfn_im1)**2
+w1t=gamma1/(beta1+tol)**2
+w2t=gamma2/(beta2+tol)**2
+w1=w1t/(w1t+w2t)
+w2=w2t/(w1t+w2t)
+flux_m=w1*fh_1+w2*fh_2
+
+if (uf .ge. 0.0 ) then
+ select_weno3=(flux_p-flux_m)/dx
+else
+ select_weno3=(flux_m-flux_p)/dx
+endif
+end function select_weno3
+
+!
+!**************************
+!
+
+
+real function select_weno5(dx,lfn_it,lfn_im1t,lfn_im2t,lfn_im3t,lfn_ip1t,lfn_ip2t,lfn_ip3t,uf)
+implicit none
+real, intent(in):: lfn_it,lfn_im1t,lfn_im2t,lfn_im3t,lfn_ip1t,lfn_ip2t,lfn_ip3t,uf
+real:: lfn_i,lfn_im1,lfn_im2,lfn_im3,lfn_ip1,lfn_ip2,lfn_ip3
+real, intent(in):: dx
+real:: flux_p,flux_m
+real, parameter:: gamma1=1./10.,gamma2=3./5.,gamma3=3./10.,tol=1e-6
+real:: w1,w2,w3,w1t,w2t,w3t,beta1,beta2,beta3
+real:: fh_1,fh_2,fh_3
+
+! 5th-order advection WENO scheme in the x,y-direction (DME)
+
+flux_p=0.
+flux_m=0.
+
+if (uf .ge. 0.0 ) then
+ lfn_i=lfn_it
+ lfn_im1=lfn_im1t
+ lfn_im2=lfn_im2t
+ lfn_im3=lfn_im3t
+ lfn_ip1=lfn_ip1t
+ lfn_ip2=lfn_ip2t
+else
+ lfn_i=lfn_it
+ lfn_im1=lfn_ip1t
+ lfn_im2=lfn_ip2t
+ lfn_im3=lfn_ip3t
+ lfn_ip1=lfn_im1t
+ lfn_ip2=lfn_im2t
+endif
+
+! numerical flux at i,j+1/2 face
+fh_1=(2.0*lfn_im2-7.0*lfn_im1+11.0*lfn_i)/6.0
+fh_2=(-1.0*lfn_im1+5.0*lfn_i+2.0*lfn_ip1)/6.0
+fh_3=(2.0*lfn_i+5.0*lfn_ip1-1.0*lfn_ip2)/6.0
+beta1=(13.0/12.0)*(lfn_im2-2.0*lfn_im1+lfn_i)**2+0.25*(lfn_im2-4.0*lfn_im1+3.0*lfn_i)**2
+beta2=(13.0/12.0)*(lfn_im1-2.0*lfn_i+lfn_ip1)**2+0.25*(lfn_im1-lfn_ip1)**2
+beta3=(13.0/12.0)*(lfn_i-2.0*lfn_ip1+lfn_ip2)**2+0.25*(3.0*lfn_i-4.0*lfn_ip1+lfn_ip2)**2
+w1t=gamma1/(beta1+tol)**2
+w2t=gamma2/(beta2+tol)**2
+w3t=gamma3/(beta3+tol)**2
+w1=w1t/(w1t+w2t+w3t)
+w2=w2t/(w1t+w2t+w3t)
+w3=w3t/(w1t+w2t+w3t)
+flux_p=w1*fh_1+w2*fh_2+w3*fh_3
+! numerical flux at i,j-1/2 face
+fh_1=(2.0*lfn_im3-7.0*lfn_im2+11.0*lfn_im1)/6.0
+fh_2=(-1.0*lfn_im2+5.0*lfn_im1+2.0*lfn_i)/6.0
+fh_3=(2.0*lfn_im1+5.0*lfn_i-1.0*lfn_ip1)/6.0
+beta1=(13.0/12.0)*(lfn_im3-2.0*lfn_im2+lfn_im1)**2+0.25*(lfn_im3-4.0*lfn_im2+3.0*lfn_im1)**2
+beta2=(13.0/12.0)*(lfn_im2-2.0*lfn_im1+lfn_i)**2+0.25*(lfn_im2-lfn_i)**2
+beta3=(13.0/12.0)*(lfn_im1-2.0*lfn_i+lfn_ip1)**2+0.25*(3.0*lfn_im1-4.0*lfn_i+lfn_ip1)**2
+w1t=gamma1/(beta1+tol)**2
+w2t=gamma2/(beta2+tol)**2
+w3t=gamma3/(beta3+tol)**2
+w1=w1t/(w1t+w2t+w3t)
+w2=w2t/(w1t+w2t+w3t)
+w3=w3t/(w1t+w2t+w3t)
+flux_m=w1*fh_1+w2*fh_2+w3*fh_3
+
+if (uf .ge. 0.0 ) then
+ select_weno5=(flux_p-flux_m)/dx
+else
+ select_weno5=(flux_m-flux_p)/dx
+endif
+end function select_weno5
+
+!
+!**************************
+!
+
real function speed_func(diffCx,diffCy,dx,dy,i,j,fp)
!*** purpose
! the level set method speed function
diff --git a/phys/module_fr_fire_driver.F b/phys/module_fr_fire_driver.F
index 26d7a401f0..e464c948ac 100644
--- a/phys/module_fr_fire_driver.F
+++ b/phys/module_fr_fire_driver.F
@@ -20,6 +20,7 @@ module module_fr_fire_driver
use module_fr_fire_phys, only : fire_params , init_fuel_cats
use module_fr_fire_util
use module_fr_fire_core, only: ignition_line_type
+use module_fr_fire_atm, only: add_fire_tracer_emissions
USE module_domain, only: domain
USE module_configure, only: grid_config_rec_type
@@ -45,7 +46,9 @@ subroutine fire_driver_em ( grid , config_flags &
,ips,ipe, kps,kpe, jps,jpe &
,ifds,ifde, jfds,jfde &
,ifms,ifme, jfms,jfme &
- ,ifps,ifpe, jfps,jfpe )
+ ,ifps,ifpe, jfps,jfpe &
+ ,rho,z_at_w,dz8w &
+ )
!*** purpose: driver from grid structure
@@ -69,6 +72,7 @@ subroutine fire_driver_em ( grid , config_flags &
ifds,ifde, jfds,jfde, &
ifms,ifme, jfms,jfme, &
ifps,ifpe, jfps,jfpe
+ real,dimension(ims:ime, kms:kme, jms:jme),intent(in), optional::rho,z_at_w,dz8w
!*** local
INTEGER:: fire_num_ignitions
@@ -99,6 +103,7 @@ subroutine fire_driver_em ( grid , config_flags &
fp%r_0 => grid%r_0 ! a rate of spread formula variable
fp%fgip => grid%fgip ! a rate of spread formula coeff
fp%ischap => grid%ischap ! a rate of spread formula switch
+ fp%iboros => grid%iboros ! Ib divided by ROS
! get ignition data
call fire_ignition_convert (config_flags,fire_max_ignitions,fire_ignition_longlat, &
@@ -136,7 +141,7 @@ subroutine fire_driver_em ( grid , config_flags &
#ifdef DM_PARALLEL
if(need_lfn_update)then
-! halo exchange on lfn width 2
+! halo exchange on lfn width 3
#include "HALO_FIRE_LFN.inc"
endif
@@ -148,6 +153,7 @@ subroutine fire_driver_em ( grid , config_flags &
!! endif
! base geopotential and roughness
#include "HALO_FIRE_PHB.inc"
+
#include "HALO_FIRE_Z0.inc"
elseif(fire_ifun.eq.2)then
@@ -169,6 +175,8 @@ subroutine fire_driver_em ( grid , config_flags &
endif
#endif
+
+
! need domain by 1 smaller, in last row.col winds are not set properly
call fire_driver_phys ( &
fire_ifun,need_lfn_update, &
@@ -195,21 +203,33 @@ subroutine fire_driver_em ( grid , config_flags &
grid%ph_2,grid%phb, & ! geopotential
grid%z0, & ! roughness height
grid%ht, & ! terrain height
- grid%xlong,grid%xlat, & ! coordinates of atm grid centers, for ignition location
- grid%lfn,grid%tign_g,grid%fuel_frac, & ! state arrays, fire grid
+ grid%xlong,grid%xlat, & ! coordinates of atm grid centers, for ignition location
+ grid%lfn, &
+ grid%lfn_hist, & ! PAJ: to initialize fire perimeter from obs
+ config_flags%fire_is_real_perim, & ! PAJ: use obs fire perimeter.
+ grid%lfn_0,grid%lfn_1,grid%lfn_2,grid%lfn_s0,grid%lfn_s1,grid%lfn_s2,grid%lfn_s3,grid%flame_length,grid%ros_front, &
+ grid%tign_g,grid%fuel_frac, & ! state arrays, fire grid
grid%fire_area, & ! redundant, for display, fire grid
+ grid%burnt_area_dt, & ! redundant, for display
lfn_out, & ! work - one timestep
grid%avg_fuel_frac, & ! out redundant arrays, atm grid
grid%grnhfx,grid%grnqfx,grid%canhfx,grid%canqfx, & ! out redundant arrays, atm grid
+ grid%grnhfx_fu,grid%grnqfx_fu, & ! out redundant arrays, atm grid
grid%uah,grid%vah, &
grid%fgrnhfx,grid%fgrnqfx,grid%fcanhfx,grid%fcanqfx, & ! out redundant arrays, atm grid
grid%ros, & ! rate of spread
- grid%fxlong,grid%fxlat, &
+ grid%fxlong,grid%fxlat,grid%fz0, &
grid%nfuel_cat, & ! input, or internal for safekeeping
grid%fuel_time, &
- fp &
+ fp, &
+ grid, & ! DME added for halo update
+ ids,ide,jds,jde,kds,kde, &
+ ims,ime,jms,jme,kms,kme, &
+ ips,ipe,jps,jpe,kps,kpe, &
+ config_flags &
)
+
#ifdef DM_PARALLEL
if(fire_ifun.eq.2)then
! halo exchange on all fuel data width 2
@@ -217,8 +237,19 @@ subroutine fire_driver_em ( grid , config_flags &
endif
#endif
-
-
+! DME update tracer emissions
+ if(fire_ifun.eq.6 .AND. config_flags%tracer_opt.eq.3)then
+ call add_fire_tracer_emissions(config_flags%tracer_opt,grid%dt,grid%dx,grid%dy, &
+ ifms,ifme,jfms,jfme, &
+ ifps,ifpe,jfps,jfpe, &
+ ids,ide,kds,kde,jds,jde, &
+ ims,ime,kms,kme,jms,jme, &
+ ips,ipe,kps,kpe,jps,jpe, &
+ rho,dz8w, &
+ grid%burnt_area_dt,grid%fgip, &
+ grid%tracer,config_flags%fire_tracer_smoke)
+ endif
+! DME
enddo
enddo
if(tsteps>0)call crash('fire_driver_em: test run of uncoupled fire model completed')
@@ -251,23 +282,38 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
ph,phb, &
z0,zs, &
xlong,xlat, &
- lfn,tign,fuel_frac, & ! state arrays, fire grid
+ lfn, &
+ lfn_hist, & ! PAJ to init perimeter from obs
+ fire_is_real_perim, & ! PAJ to init perimeter from obs
+ lfn_0,lfn_1,lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3,flame_length,ros_front, &
+ tign,fuel_frac, & ! state arrays, fire grid
fire_area, & ! redundant state, fire grid
+ burnt_area_dt, & ! redundant state, fire grid
lfn_out, & ! out level set function
avg_fuel_frac, &
grnhfx,grnqfx,canhfx,canqfx, & ! out redundant arrays, atm grid
+ grnhfx_fu,grnqfx_fu, & ! out redundant arrays, atm grid
uah,vah, & ! out atm grid
fgrnhfx,fgrnqfx,fcanhfx,fcanqfx, & ! out redundant arrays, fire grid
ros, &
- fxlong,fxlat, & !
+ fxlong,fxlat,fz0, &
nfuel_cat, & ! in array, data, fire grid, or constant internal
fuel_time, & ! save constant internal data, fire grid
- fp & ! fire params
+ fp, & ! fire params
+ grid, &
+ ids_hu,ide_hu,jds_hu,jde_hu,kds_hu,kde_hu, &
+ ims_hu,ime_hu,jms_hu,jme_hu,kms_hu,kme_hu, &
+ ips_hu,ipe_hu,jps_hu,jpe_hu,kps_hu,kpe_hu, &
+ config_flags &
)
+
USE module_dm, only:wrf_dm_maxval
+USE module_domain, only: domain
implicit none
+ type(domain) , target :: grid ! state
+
!*** arguments
integer, intent(in)::ifun, &
@@ -287,6 +333,12 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
nfuel_cat0, & ! fuel category to initialize everything to
num_tiles ! number of tiles
+integer, intent(in)::ids_hu,ide_hu,jds_hu,jde_hu,kds_hu,kde_hu, &
+ ims_hu,ime_hu,jms_hu,jme_hu,kms_hu,kme_hu, &
+ ips_hu,ipe_hu,jps_hu,jpe_hu,kps_hu,kpe_hu
+
+type (grid_config_rec_type), intent(in) :: config_flags
+
logical, intent(in)::restart
@@ -324,11 +376,20 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
real, intent(inout), dimension(ifms:ifme, jfms:jfme):: &
lfn,tign,fuel_frac, & ! state: level function, ign time, fuel left
+ lfn_hist, & ! PAJ: to init obs. fira perimeter
lfn_out ! fire wind velocities
+real, intent(inout), dimension(ifms:ifme, jfms:jfme):: &
+ lfn_0,lfn_1,lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3,flame_length,ros_front ! level set function and reinitialization steps
+
real, intent(out), dimension(ifms:ifme, jfms:jfme):: &
fire_area ! fraction of each cell burning
+real, intent(out), dimension(ifms:ifme, jfms:jfme):: &
+ burnt_area_dt
+real, intent(out), dimension(ims:ime, jms:jme):: &
+ grnhfx_fu,grnqfx_fu
+
real, intent(out), dimension(ims:ime, jms:jme):: & ! redundant arrays, for display purposes only (atm grid)
avg_fuel_frac, & ! average fuel fraction
grnhfx, & ! heat flux from surface fire (W/m^2)
@@ -345,7 +406,7 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
! ***** data (constant in time) *****
-real, dimension(ifms:ifme, jfms:jfme), intent(inout)::fxlong,fxlat ! fire mesh coordinates
+real, dimension(ifms:ifme, jfms:jfme), intent(inout)::fxlong,fxlat,fz0 ! fire mesh coordinates
real, intent(out), dimension(ifms:ifme, jfms:jfme)::fuel_time ! fire params arrays
type(fire_params),intent(inout)::fp
@@ -360,6 +421,8 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
integer::ignitions_done_tile(num_tiles),ignited_tile(num_ignitions,num_tiles)
integer::ignitions_done,ignited_patch(num_ignitions),idex,jdex
+ ! PAJ:
+ logical :: fire_is_real_perim
!*** executable
@@ -505,11 +568,26 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
if(ignition_longlat .eq.0)then
! set ideal fire mesh coordinates - used for ignition only
! do not forget to set unit_fxlong, unit_fxlat outside of parallel loop
- !call set_ideal_coord( dxf,dyf, &
- ! ifds,ifde,jfds,jfde, &
- ! ifms,ifme,jfms,jfme, &
- ! ifts,ifte,jfts,jfte, &
- ! fxlong,fxlat )
+
+ ! DME: Next call added to fixe a bug when
+ ! initializing a nested domain (before, fxlong & fxlat where not
+ ! assigned, therefore were zero and fire does not ignite)
+ call set_ideal_coord( dxf,dyf, &
+ ifds,ifde,jfds,jfde, &
+ ifms,ifme,jfms,jfme, &
+ ifts,ifte,jfts,jfte, &
+ fxlong,fxlat )
+ call interpolate_z2fire(id, &
+ ids,ide,jds,jde, &
+ ims,ime,jms,jme, &
+ ips,ipe,jps,jpe, &
+ its,ite,jts,jte, &
+ ifds,ifde,jfds,jfde, &
+ ifms,ifme,jfms,jfme, &
+ ifts,ifte,jfts,jfte, &
+ ir,jr, &
+ z0, &
+ fz0,1)
!call set_ideal_coord( dx,dy, &
! ids,ide,jds,jde, &
! ims,ime,jms,jme, &
@@ -533,7 +611,7 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
ifts,ifte,jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
xlat, & ! atm grid arrays in
- fxlat) ! fire grid arrays out
+ fxlat,0) ! fire grid arrays out
call interpolate_z2fire(id, & ! for debug output, <= 0 no output
ids,ide, jds,jde, & ! atm grid dimensions
@@ -545,7 +623,19 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
ifts,ifte,jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
xlong, & ! atm grid arrays in
- fxlong) ! fire grid arrays out
+ fxlong,0) ! fire grid arrays out
+
+ call interpolate_z2fire(id, &
+ ids,ide, jds,jde, &
+ ims,ime, jms,jme, &
+ ips,ipe,jps,jpe, &
+ its,ite,jts,jte, &
+ ifds, ifde, jfds, jfde, &
+ ifms, ifme, jfms, jfme, &
+ ifts,ifte,jfts,jfte, &
+ ir,jr, &
+ z0, &
+ fz0,1)
endif
@@ -557,9 +647,10 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fp%zsf,'driver_phys:zsf')
elseif(ifun.eq.3)then ! interpolate winds to the fire grid
-
+
+
if(use_atm_vars)then
-
+
call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,z0,'z0',id)
call write_array_m3(its,ite1,kts,kde-1,jts,jte,ims,ime,kms,kme,jms,jme,u,'u_2',id)
call write_array_m3(its,ite,kts,kde-1,jts,jte1,ims,ime,kms,kme,jms,jme,v,'v_2',id)
@@ -581,9 +672,10 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
ph,phb, &
z0,zs, & ! 2D atm grid arrays in
uah,vah, & ! 2D atm grid out
- fp%vx,fp%vy) ! fire grid arrays out
+ fp%vx,fp%vy,fz0) ! fire grid arrays out
+
endif
-
+
endif
! the following executes in any case
@@ -600,15 +692,24 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
ignition_line, & ! description of ignition lines
ignitions_done_tile(ij),ignited_tile(1,ij), &
fxlong,fxlat,unit_fxlong,unit_fxlat, & ! fire mesh coordinates
- lfn,lfn_out,tign,fuel_frac, & ! state: level function, ign time, fuel left
+ lfn, & ! state: level function
+ lfn_hist, & ! PAJ: to init obs fire perimeter
+ fire_is_real_perim, & ! PAJ: to init obs fire perimeter
+ lfn_0,lfn_1,lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3,flame_length,ros_front, & ! state
+ lfn_out,tign,fuel_frac, & ! state: ign time, fuel left
fire_area, & ! output: fraction of cell burning
+ burnt_area_dt, &
fgrnhfx,fgrnqfx, & ! output: heat fluxes
ros, & ! output: rate of spread for display
nfuel_cat, & ! fuel data per point
fuel_time, & ! save derived internal data
- fp & ! fire coefficients
+ fp, & ! fire coefficients
+ grid, & ! DME added for halo update
+ ids_hu,ide_hu,jds_hu,jde_hu,kds_hu,kde_hu, &
+ ims_hu,ime_hu,jms_hu,jme_hu,kms_hu,kme_hu, &
+ ips_hu,ipe_hu,jps_hu,jpe_hu,kps_hu,kpe_hu &
)
-
+
if(ifun.eq.6)then ! heat fluxes into the atmosphere
call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fgrnhfx,'fire_driver:fgrnhfx')
@@ -646,6 +747,9 @@ subroutine fire_driver_phys (ifun,need_lfn_update, &
s = 1./(ir*jr)
do j=jts,jte
do i=its,ite
+ ! DME heat fluxes contribution for the case wiythout feedback
+ grnhfx_fu(i,j)=grnhfx(i,j)*s
+ grnqfx_fu(i,j)=grnqfx(i,j)*s
! scale surface fluxes to get the averages
avg_fuel_frac(i,j)=avg_fuel_frac(i,j)*s
grnhfx(i,j)=fire_atm_feedback*grnhfx(i,j)*s
@@ -881,7 +985,7 @@ subroutine interpolate_z2fire(id, & ! for debug output, <= 0 no
ifts,ifte,jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
zs, & ! atm grid arrays in
- zsf) ! fire grid arrays out
+ zsf,flag_z0) ! fire grid arrays out
implicit none
!*** purpose: interpolate height
@@ -899,6 +1003,7 @@ subroutine interpolate_z2fire(id, & ! for debug output, <= 0 no
real, intent(in), dimension(ims:ime, jms:jme):: zs ! terrain height at atm cell centers & ! terrain height
real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
zsf ! terrain height fire grid nodes
+integer,intent(in)::flag_z0
!*** local
@@ -942,6 +1047,14 @@ subroutine interpolate_z2fire(id, & ! for debug output, <= 0 no
za, & ! in atm grid
zsf) ! out fire grid
+if (flag_z0 .eq. 1 ) then
+ do j=jfts1,jfte1
+ do i=ifts1,ifte1
+ zsf(i,j)=max(zsf(i,j),0.001)
+ enddo
+ enddo
+endif
+
end subroutine interpolate_z2fire
!
!*****************************
@@ -963,7 +1076,7 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no
ph,phb, &
z0,zs, &
uah,vah, &
- uf,vf) ! fire grid arrays out
+ uf,vf,z0f) ! fire grid arrays out
implicit none
!*** purpose: interpolate winds and height
@@ -993,6 +1106,7 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no
vah !
real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
uf,vf ! wind velocity fire grid nodes
+real,intent(in),dimension(ifms:ifme,jfms:jfme)::z0f ! roughness length in fire grid
!*** local
@@ -1008,6 +1122,9 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no
real::r_nan
integer::i_nan
equivalence (i_nan,r_nan)
+real::fire_wind_height_local,z0fc
+real::ust_d,wsf,wsf1,uf_temp,vf_temp
+real,parameter::vk_kappa=0.4
!*** executable
@@ -1153,17 +1270,25 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no
call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,hgtu,'hgtu',id)
call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,hgtv,'hgtv',id)
#endif
- logfwh = log(fire_wind_height)
+! DME
+if (fire_lsm_zcoupling) then
+ logfwh = log(fire_lsm_zcoupling_ref)
+ fire_wind_height_local = fire_lsm_zcoupling_ref
+else
+ logfwh = log(fire_wind_height)
+ fire_wind_height_local = fire_wind_height
+endif
! interpolate u, staggered in X
do j = jtsu,jteu ! compute on domain by 1 smaller, otherwise z0 and ph not available
- do i = itsu,iteu ! compute with halo 2
+! do i = itsu,iteu ! compute with halo 2
+ do i = itsu,iteu
zr = 0.5*(z0(i,j)+z0(i-1,j)) ! interpolated roughness length under this u point
- if(fire_wind_height > zr)then !
+ if(fire_wind_height_local > zr)then !
do k=kds,kdmax
ht = hgtu(i,k,j) ! height of this u point above the ground
- if( .not. ht < fire_wind_height) then ! found layer k this point is in
+ if( .not. ht < fire_wind_height_local) then ! found layer k this point is in
loght = log(ht)
if(k.eq.kds)then ! first layer, log linear interpolation from 0 at zr
logz0 = log(zr)
@@ -1185,15 +1310,20 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no
enddo
enddo
+ do j = jtsu,jteu
+ ua(itsu-1,j)=ua(itsu,j)
+ enddo
+
+
! interpolate v, staggered in Y
do j = jtsv,jtev
do i = itsv,itev
zr = 0.5*(z0(i,j-1)+z0(i,j)) ! roughness length under this v point
- if(fire_wind_height > zr)then !
+ if(fire_wind_height_local > zr)then !
do k=kds,kdmax
ht = hgtv(i,k,j) ! height of this u point above the ground
- if( .not. ht < fire_wind_height) then ! found layer k this point is in
+ if( .not. ht < fire_wind_height_local) then ! found layer k this point is in
loght = log(ht)
if(k.eq.kds)then ! first layer, log linear interpolation from 0 at zr
logz0 = log(zr)
@@ -1215,6 +1345,10 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no
enddo
enddo
+ do i = itsv,itev
+ va(i,jtsv-1)=va(i,jtsv)
+ enddo
+
#ifdef DEBUG_OUT
! store the output for diagnostics
do j = jts,jte1
@@ -1229,6 +1363,7 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no
#endif
ips1 = ifval(ips.eq.ids,ips+1,ips)
+
call continue_at_boundary(1,1,0., & ! x direction
TDIMS, &! memory dims atm grid tile
ids+1,ide,jds,jde, & ! domain dims - where u defined
@@ -1238,6 +1373,7 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no
ua) ! array
jps1 = ifval(jps.eq.jds,jps+1,jps)
+
call continue_at_boundary(1,1,0., & ! y direction
TDIMS, & ! memory dims atm grid tile
ids,ide,jds+1,jde, & ! domain dims - where v wind defined
@@ -1329,6 +1465,23 @@ subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no
va, & ! in atm grid
vf) ! out fire grid
+! DME here code to extrapolate mid-flame height velocity -> fire_lsm_zcoupling = .true.
+if (fire_lsm_zcoupling) then
+ do j = jfts,jfte
+ do i = ifts,ifte
+ uf_temp=uf(i,j)
+ vf_temp=vf(i,j)
+ wsf=max(sqrt(uf_temp**2.+vf_temp**2.),0.1)
+ z0fc=z0f(i,j)
+ ust_d=wsf*vk_kappa/log(fire_lsm_zcoupling_ref/z0fc)
+ wsf1=(ust_d/vk_kappa)*log((fire_wind_height+z0fc)/z0fc)
+ uf(i,j)=wsf1*uf_temp/wsf
+ vf(i,j)=wsf1*vf_temp/wsf
+ enddo
+ enddo
+endif
+
+
call print_2d_stats_vec(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,vf,'fire wind (m/s)')
! call print_2d_stats_vec(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,vf,'fire wind extended')
#ifdef DEBUG_OUT
@@ -1383,15 +1536,23 @@ subroutine set_flags(config_flags)
fire_const_grnqfx = config_flags%fire_const_grnqfx
fire_atm_feedback = config_flags%fire_atm_feedback
boundary_guard = config_flags%fire_boundary_guard
-fire_back_weight = config_flags%fire_back_weight
fire_grows_only = config_flags%fire_grows_only
fire_upwinding = config_flags%fire_upwinding
fire_upwind_split = config_flags%fire_upwind_split
fire_viscosity = config_flags%fire_viscosity
fire_lfn_ext_up = config_flags%fire_lfn_ext_up
fire_test_steps = config_flags%fire_test_steps
-!fire_topo_from_atm = config_flags%fire_topo_from_atm
fire_advection = config_flags%fire_advection
+fire_lsm_reinit = config_flags%fire_lsm_reinit
+fire_lsm_reinit_iter = config_flags%fire_lsm_reinit_iter
+fire_upwinding_reinit = config_flags%fire_upwinding_reinit
+fire_lsm_band_ngp = config_flags%fire_lsm_band_ngp
+fire_lsm_zcoupling = config_flags%fire_lsm_zcoupling
+fire_lsm_zcoupling_ref = config_flags%fire_lsm_zcoupling_ref
+fire_viscosity_bg = config_flags%fire_viscosity_bg
+fire_viscosity_band = config_flags%fire_viscosity_band
+fire_viscosity_ngp = config_flags%fire_viscosity_ngp
+fire_slope_factor = config_flags%fire_slope_factor
diff --git a/phys/module_fr_fire_driver_wrf.F b/phys/module_fr_fire_driver_wrf.F
index ff5471efb3..9905d64d36 100644
--- a/phys/module_fr_fire_driver_wrf.F
+++ b/phys/module_fr_fire_driver_wrf.F
@@ -38,7 +38,7 @@ subroutine fire_driver_em_init (grid , config_flags &
! dummies
real,dimension(1,1,1)::rho,z_at_w,dz8w
-
+
call message('fire_driver_em_init: FIRE initialization start')
! get fire mesh dimensions
@@ -110,7 +110,9 @@ subroutine fire_driver_em_step (grid , config_flags &
,ifds,ifde, jfds,jfde &
,ifms,ifme, jfms,jfme &
,ifps,ifpe, jfps,jfpe &
- )
+ ,rho,z_at_w,dz8w &
+ )
+
! --- add heat and moisture fluxes to tendency variables by postulated decay
do ij=1,grid%num_tiles
diff --git a/phys/module_fr_fire_model.F b/phys/module_fr_fire_model.F
index 6741d7f379..7a0f03e45f 100644
--- a/phys/module_fr_fire_model.F
+++ b/phys/module_fr_fire_model.F
@@ -25,12 +25,21 @@ subroutine fire_model ( &
ignition_line, & ! small array of ignition line descriptions
ignitions_done,ignited_tile, &
coord_xf,coord_yf,unit_xf,unit_yf, & ! fire mesh coordinates
- lfn,lfn_out,tign,fuel_frac,fire_area, & ! state: level function, ign time, fuel left, area burning
+ lfn, & ! state: level function
+ lfn_hist, & ! PAJ: to init obs fire perimeter.
+ fire_is_real_perim, & ! PAJ: to init obs fire perimeter.
+ lfn_0,lfn_1,lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3,flame_length,ros_front, & ! state
+ lfn_out,tign,fuel_frac,fire_area, & ! state: level function, ign time, fuel left, area burning
+ burnt_area_dt, &
grnhfx,grnqfx, & ! output: heat fluxes
ros, & ! output: rate of spread
nfuel_cat, & ! fuel data per point
fuel_time, & ! save derived internal data
- fp &
+ fp, &
+ grid, &
+ ids_hu,ide_hu,jds_hu,jde_hu,kds_hu,kde_hu, &
+ ims_hu,ime_hu,jms_hu,jme_hu,kms_hu,kme_hu, &
+ ips_hu,ipe_hu,jps_hu,jpe_hu,kps_hu,kpe_hu &
)
! This subroutine implements the fire spread model.
@@ -69,10 +78,24 @@ subroutine fire_model ( &
!
! enddo
+USE module_domain , only: domain
+#ifdef DM_PARALLEL
+ USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks
+ USE module_comm_dm , ONLY : halo_fire_lfn_sub
+#endif
+USE module_configure, only: grid_config_rec_type
+
implicit none
!*** arguments
+! DME added for halo update
+ type(domain) , target :: grid
+ integer, intent(in):: ids_hu,ide_hu,jds_hu,jde_hu,kds_hu,kde_hu, &
+ ims_hu,ime_hu,jms_hu,jme_hu,kms_hu,kme_hu, &
+ ips_hu,ipe_hu,jps_hu,jpe_hu,kps_hu,kpe_hu
+!
+
! control switches
integer, intent(in) :: id
integer, intent(in) :: ifun ! 1 = initialize run pass 1
@@ -88,7 +111,7 @@ subroutine fire_model ( &
integer, intent(in) :: ifuelread,nfuel_cat0 ! for set_fire_params
integer, intent(in) :: ifds,ifde,jfds,jfde,& ! fire domain bounds
ifps,ifpe,jfps,jfpe ! patch - nodes owned by this process
-integer, intent(in) :: ifts,ifte,jfts,jfte ! fire tile bounds
+integer, intent(in) :: ifts,ifte,jfts,jfte ! fire tile bounds
integer, intent(in) :: ifms,ifme,jfms,jfme ! fire memory array bounds
REAL,INTENT(in) :: time_start,dt ! starting time, time step
REAL,INTENT(in) :: fdx,fdy ! spacing of the fire mesh
@@ -100,14 +123,24 @@ subroutine fire_model ( &
real, intent(in):: unit_xf,unit_yf ! coordinate units in m
! state
+ ! PAJ: to init obs fire perimeter
+ real, intent(in), dimension(ifms:ifme,jfms:jfme):: lfn_hist
+ logical, intent(in) :: fire_is_real_perim
+
REAL, INTENT(inout), dimension(ifms:ifme,jfms:jfme):: &
lfn , & ! level function: fire is where lfn<0 (node)
tign , & ! absolute time of ignition (node)
fuel_frac ! fuel fraction (node), currently redundant
+REAL, INTENT(inout), dimension(ifms:ifme,jfms:jfme):: &
+ lfn_0,lfn_1,lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3,flame_length,ros_front ! level function stages
+
REAL, INTENT(out), dimension(ifms:ifme,jfms:jfme):: &
fire_area ! fraction of each cell burning
+REAL, INTENT(out), dimension(ifms:ifme,jfms:jfme):: &
+ burnt_area_dt
+
! output
REAL, INTENT(out), dimension(ifms:ifme,jfms:jfme):: &
lfn_out, & !
@@ -129,6 +162,10 @@ subroutine fire_model ( &
logical:: freeze_fire
integer:: stat_lev=1
+ ! PAJ:
+ real :: start_time_ig, end_time_ig
+ real, parameter :: EPSILON = 0.00001
+
!*** executable
call check_mesh_2dim(ifts-1,ifte+1,jfts-1,jfte+1,ifms,ifme,jfms,jfme)
@@ -258,14 +295,50 @@ subroutine fire_model ( &
if(.not. freeze_fire)then
- call prop_ls(id, &
- ifds,ifde,jfds,jfde, & ! fire domain dims - the whole domain
- ifms,ifme,jfms,jfme, &
- ifps,ifpe,jfps,jfpe, & ! patch - nodes owned by this process
- ifts,ifte,jfts,jfte, &
- time_start,dt,fdx,fdy,tbound, &
- lfn,lfn_out,tign,ros, fp &
- )
+ call prop_ls_rk3(id, &
+ ifds,ifde,jfds,jfde, &
+ ifms,ifme,jfms,jfme, &
+ ifps,ifpe,jfps,jfpe, &
+ ifts,ifte,jfts,jfte, &
+ time_start,dt,fdx,fdy,tbound, &
+ lfn, &
+ lfn_0,lfn_1,lfn_2, &
+ lfn_out,tign,ros, fp, &
+ grid, &
+ ids_hu,ide_hu,jds_hu,jde_hu,kds_hu,kde_hu, &
+ ims_hu,ime_hu,jms_hu,jme_hu,kms_hu,kme_hu, &
+ ips_hu,ipe_hu,jps_hu,jpe_hu,kps_hu,kpe_hu &
+ )
+
+ call tign_update(ifts,ifte,jfts,jfte, &
+ ifms,ifme,jfms,jfme, &
+ ifds,jfds,ifde,jfde, &
+ time_start,dt, &
+ lfn,lfn_out,tign &
+ )
+
+ call calc_flame_length(ifts,ifte,jfts,jfte, &
+ ifms,ifme,jfms,jfme, &
+ ros,fp%iboros,flame_length,ros_front,fire_area)
+
+ if (fire_lsm_reinit) then ! DME added call to reinitialize level-set function
+
+ call reinit_ls_rk3(id, &
+ ifts,ifte,jfts,jfte, &
+ ifms,ifme,jfms,jfme, &
+ ifds,ifde,jfds,jfde, &
+ ifps,ifpe,jfps,jfpe, &
+ time_start,dt,fdx,fdy, &
+ lfn, &
+ lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3, &
+ lfn_out,tign, &
+ grid, &
+ ids_hu,ide_hu,jds_hu,jde_hu,kds_hu,kde_hu, &
+ ims_hu,ime_hu,jms_hu,jme_hu,kms_hu,kme_hu, &
+ ips_hu,ipe_hu,jps_hu,jpe_hu,kps_hu,kpe_hu &
+ )
+ endif ! fire_lsm_reinit
+
else
call message('fire_model: EXPERIMENTAL: skipping fireline propagation')
@@ -275,19 +348,36 @@ subroutine fire_model ( &
! this cannot be done in the time step itself because of race condition
! some thread may still be using lfn as input in their tile halo
- if(.not. freeze_fire)then
-
- do j=jfts,jfte
- do i=ifts,ifte
- lfn(i,j)=lfn_out(i,j)
- ! if want to try timestep again treat tign the same way here
- ! even if tign does not need a halo
- enddo
- enddo
+ if (.not. freeze_fire) then
+
+ do j=jfts,jfte
+ do i=ifts,ifte
+ lfn(i,j)=lfn_out(i,j)
+ ! if want to try timestep again treat tign the same way here
+ ! even if tign does not need a halo
+ enddo
+ enddo
endif
! check for ignitions
+!paj
+ ig = 1
+ start_time_ig = ignition_line(ig)%start_time
+ end_time_ig = ignition_line(ig)%end_time
+
+ if ( fire_is_real_perim .and. time_start >= start_time_ig .and. time_start < start_time_ig + dt) then
+ ignited = 0
+ do j = jfts, jfte
+ do i = ifts, ifte
+ lfn(i, j) = lfn_hist(i, j)
+ if (abs(lfn(i, j)) < EPSILON) then
+ tign(i, j)= time_start
+ ignited = ignited + 1
+ end if
+ enddo
+ enddo
+ elseif (.not. fire_is_real_perim) then
do ig = 1,num_ignitions
! for now, check for ignition every time step...
@@ -313,7 +403,9 @@ subroutine fire_model ( &
! endif
enddo
+ end if
+
call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme, &
lfn,'fire_model: lfn out')
@@ -339,6 +431,7 @@ subroutine fire_model ( &
do j=jfts,jfte
do i=ifts,ifte
fuel_frac_burnt(i,j)=fuel_frac(i,j)-fuel_frac_end(i,j) ! fuel lost this timestep
+ burnt_area_dt(i,j)=fuel_frac_burnt(i,j)
fuel_frac(i,j)=fuel_frac_end(i,j) ! copy new value to state array
enddo
enddo
diff --git a/phys/module_fr_fire_phys.F b/phys/module_fr_fire_phys.F
index 856536d169..af28829203 100644
--- a/phys/module_fr_fire_phys.F
+++ b/phys/module_fr_fire_phys.F
@@ -19,6 +19,7 @@ module module_fr_fire_phys
real,pointer,dimension(:,:):: bbb,betafl,phiwc,r_0 ! spread formula coefficients
real,pointer,dimension(:,:):: fgip ! init mass of surface fuel (kg/m^2)
real,pointer,dimension(:,:):: ischap ! if fuel is chaparral and want rate of spread treated differently
+real,pointer,dimension(:,:):: iboros ! fire instensity over rate of spread -> avialable fuel x heat of combustion
end type fire_params
! use as
@@ -296,6 +297,7 @@ subroutine write_fuels_m(nsteps,maxwind,maxslope)
!real,pointer,dimension(:,:):: ischap ! 1 if chapparal
!end type fire_params
real, dimension(1:2,1:nsteps), target::vx,vy,zsf,dzdxf,dzdyf,bbb,betafl,phiwc,r_0,fgip,ischap
+real, dimension(1:2,1:nsteps), target::iboros
real, dimension(1:2,1:nsteps)::nfuel_cat,fuel_time,ros
real::ros_base,ros_wind,ros_slope,propx,propy,r
@@ -309,6 +311,7 @@ subroutine write_fuels_m(nsteps,maxwind,maxslope)
fp%r_0=>r_0
fp%fgip=>fgip
fp%ischap=>ischap
+fp%iboros=>iboros
iounit = open_text_file('fuels.m','write')
@@ -434,21 +437,102 @@ subroutine set_fire_params( &
integer::nerr
character(len=128)::msg
+integer :: kk
+integer,parameter :: nf_sb = 204 ! maximum category on
+integer,dimension(1:nf_sb) :: ksb ! Anderson82 + S&B2005 fuel categories array
+
!*** executable
+!*** Scott & Burgan ***!
+! assign no fuel by default to all the categories
+do kk=1,nf_sb
+ ksb(kk)=14
+enddo
+! Anderson 1982
+ksb(1)=1
+ksb(2)=2
+ksb(3)=3
+ksb(4)=4
+ksb(5)=5
+ksb(6)=6
+ksb(7)=7
+ksb(8)=8
+ksb(9)=9
+ksb(10)=10
+ksb(11)=11
+ksb(12)=12
+ksb(13)=13
+! Scott & Burgan crosswalks
+! Short grass -- 1
+ksb(101)=1
+ksb(104)=1
+ksb(107)=1
+! Timber grass and understory -- 2
+ksb(102)=2
+ksb(121)=2
+ksb(122)=2
+ksb(123)=2
+ksb(124)=2
+! Tall grass -- 3
+ksb(103)=3
+ksb(105)=3
+ksb(106)=3
+ksb(108)=3
+ksb(109)=3
+! Chaparral -- 4
+ksb(145)=4
+ksb(147)=4
+! Brush -- 5
+ksb(142)=5
+! Dormant Brushi -- 6
+ksb(141)=6
+ksb(146)=6
+! Southern Rough -- 7
+ksb(143)=7
+ksb(144)=7
+ksb(148)=7
+ksb(149)=7
+! Compact Timber Litter -- 8
+ksb(181)=8
+ksb(183)=8
+ksb(184)=8
+ksb(187)=8
+! Hardwood Litter -- 9
+ksb(182)=9
+ksb(186)=9
+ksb(188)=9
+ksb(189)=9
+! Timber (understory) -- 10
+ksb(161)=10
+ksb(162)=10
+ksb(163)=10
+ksb(164)=10
+ksb(165)=10
+! Light Logging Slash -- 11
+ksb(185)=11
+ksb(201)=11
+! Medium Logging Slash -- 12
+ksb(202)=12
+! Heavy Logging Slash -- 13
+ksb(203)=13
+ksb(204)=13
+
+! ****** !
+
nerr=0
do j=jfts,jfte
do i=ifts,ifte
! fuel category
- k=int( nfuel_cat(i,j) )
+ k=ksb(int(nfuel_cat(i,j))) ! DME S&B05
if(k.eq.no_fuel_cat)then ! no fuel
fp%fgip(i,j)=0. ! no mass
fp%ischap(i,j)=0.
- fp%betafl(i,j)=0. ! to prevent division by zero
+ fp%betafl(i,j)=1. ! DME: set to 1.0 to prevent fp%betafl(i,j)**(-0.3) to be Inf in fire_ros
fp%bbb(i,j)=1. !
fuel_time(i,j)=7./0.85 ! does not matter, just what was there before
fp%phiwc(i,j)=0.
fp%r_0(i,j)=0. ! no fuel, no spread.
+ fp%iboros(i,j)=0. ! DME Ib/ROS zero for no fuel
else
if(k.eq.0.and.nfuel_cat0.ge.1.and.nfuel_cat0.le.nfuelcats)then
! replace k=0 by default
@@ -502,6 +586,7 @@ subroutine set_fire_params( &
etas = 0.174* se(k)**(-0.19) ! mineral damping coef
ir = gamma * wn * fuelheat * etam * etas !rxn intensity,btu/ft^2 min
! irm = ir * 1055./( 0.3048**2 * 60.) * 1.e-6 !for mw/m^2
+ fp%iboros(i,j) = ir * 1055./( 0.3048**2 * 60.) * 1.e-3 * (60.*12.6/savr(k)) ! I_R x t_r (kJ m^-2)
xifr = exp( (0.792 + 0.681*savr(k)**0.5) &
* (fp%betafl(i,j)+0.1)) /(192. + 0.2595*savr(k)) ! propagating flux ratio
@@ -759,7 +844,7 @@ subroutine fire_ros(ros_base,ros_wind,ros_slope, &
! if advection, multiply by the cosines
ros_wind=ros_wind*cor_wind
ros_slope=ros_slope*cor_slope
-!
+
! ----------note! put an 6 m/s cap on max spread rate -----------
! rosm= min(rosm, 6.) ! no faster than this cap ! param !
@@ -771,6 +856,7 @@ subroutine fire_ros(ros_base,ros_wind,ros_slope, &
ros_slope = ros_slope - excess*ros_slope/(ros_wind+ros_slope)
endif
+
!write(msg,*)i,j,' speed=',speed,' tanphi',tanphi,' ros=',ros_base,ros_wind,ros_slope
!call message(msg)
return
diff --git a/phys/module_fr_fire_util.F b/phys/module_fr_fire_util.F
index 322250c4b7..b47e9d4ca8 100644
--- a/phys/module_fr_fire_util.F
+++ b/phys/module_fr_fire_util.F
@@ -23,21 +23,31 @@ module module_fr_fire_util
fuel_left_jrl=2, &
boundary_guard=-1, & ! crash if fire gets this many cells to domain boundary, -1=off
fire_grows_only=1, & ! fire can spread out only (level set functions may not increase)
- fire_upwinding=3, & ! upwind normal spread: 1=standard, 2=godunov, 3=eno, 4=sethian
+ fire_upwinding=9, & ! upwind normal spread: 1=standard, 2=godunov, 3=eno, 4=sethian, 5=2nd-order, 6=WENO3, 7=WENO5, 8=hybrid WENO3/ENO1, 9=hybrid WENO5/ENO1
fire_upwind_split=0, & ! 1=upwind advection separately from normal direction spread
fire_test_steps=0, & ! 0=no fire, 1=normal, >1 = do specified number of steps and terminate (testing only)
fire_topo_from_atm=1, & ! 0 = expect ZSF set correctly on entry, 1 = populate by interploating from atmosphere
- fire_advection=0 ! 0 = fire spread from normal wind/slope (CAWFE), 1 = full speed projected
-
+ fire_advection=0, & ! 0 = fire spread from normal wind/slope (CAWFE), 1 = full speed projected
+ fire_upwinding_reinit=4, & ! spatial discretization for reinitialization PDE: 1=WENO3, 2=WENO5, 3=hybrid WENO3-ENO1, 4=hybrid WENO5-ENO1
+ fire_lsm_reinit_iter=1, & ! number of iterations for the reinitialization PDE
+ fire_lsm_band_ngp=4, & ! number of grid points around lfn=0 that the higher-order scheme WENO5/WENO3 is used (ENO1 elsewhere), for ire_upwinding=8,9 and fire_upwinding_reinit=3,4 and
+ fire_viscosity_ngp=4 ! number of grid points around lfn=0 where low artificial viscosity is used = fire_viscosity_bg
real, save:: &
fire_const_time=-1, & ! time from ignition to start constant heat output <0=never
fire_const_grnhfx=-1., & ! if both >=0, the constant heat flux to output then
fire_const_grnqfx=-1., & ! if both >=0, the constant heat flux to output then
fire_atm_feedback=1. , & ! 1 = normal, 0. = one way coupling atmosphere -> fire only
- fire_back_weight=0.5, & ! RK parameter, 0 = Euler method, 0.5 = Heun, 1 = fake backward Euler
fire_viscosity=0.4, & ! artificial viscosity
- fire_lfn_ext_up=1. ! 0.=extend level set function at boundary by reflection, 1.=always up
+ fire_lfn_ext_up=1., & ! 0.=extend level set function at boundary by reflection, 1.=always up
+ fire_lsm_zcoupling_ref=50., & ! Reference height from wich u at fire_wind_hegiht is calculated using a logarithmic profile
+ fire_viscosity_bg=0.4, & ! artificial viscosity in the near-front region
+ fire_viscosity_band=0.5, & ! number of times the hybrid advection band to transition from fire_viscosity_bg to fire_viscosity
+ fire_slope_factor=1.0 ! slope correction factor
+
+logical,save:: &
+ fire_lsm_reinit=.true., & ! flag to activate reinitialization of level set method
+ fire_lsm_zcoupling=.false. ! flag to activate reference velocity at a different height from fire_wind_height
integer, parameter:: REAL_SUM=10, REAL_MAX=20, RNRM_SUM=30, RNRM_MAX=40
@@ -204,10 +214,9 @@ subroutine continue_at_boundary(ix,iy,bias, & ! do x direction or y direction
integer i,j
character(len=128)::msg
integer::its1,ite1,jts1,jte1
-integer,parameter::halo=1 ! width of halo region to update
+integer,parameter::halo=1 ! only 1 domain halo is needed since ENO1 is used near domain boundaries
!*** executable
-! check if there is space for the extension
-call check_mesh_2dim(its-1,ite+1,jts-1,jte+1,ims,ime,jms,jme)
+
! for dislay only
itso=its
jtso=jts
diff --git a/test/em_fire/namelist.input_hill_simple b/test/em_fire/namelist.input_hill_simple
index 00c0d7490c..bcae2c61ac 100644
--- a/test/em_fire/namelist.input_hill_simple
+++ b/test/em_fire/namelist.input_hill_simple
@@ -93,6 +93,7 @@
time_step_sound = 20, 20, 20,
moist_adv_opt = 1, 1, 1,
scalar_adv_opt = 1, 1, 1,
+ tracer_opt = 3, 3, 3,
/
&bdy_control
@@ -166,15 +167,12 @@
!
! method switches for developers only, do not change!
!
- fire_boundary_guard = -1, ! integer, number of cells to stop when fire close to the domain boundary, -1 turn off
fire_fuel_left_irl=2, ! refinement to integrate fuel_left, must be even
fire_fuel_left_jrl=2, ! refinement to integrate fuel_left, must be even
fire_atm_feedback=1.0, ! real, multiplier for heat fluxes, 1.=normal, 0.=turn off two-way coupling
- fire_back_weight=0.5, ! RK timestepping coefficient, 0=forward, 0.5=Heun
fire_grows_only=1, ! if >0 level set function cannot increase = fire can only grow
fire_viscosity=0.4, ! artificial viscosity in level set method (max 1, needed with fire_upwinding=0)
- fire_upwinding=3, ! 0=none, 1=standard, 2=godunov, 3=eno, 4=sethian
+ fire_upwinding=9, ! 0=none, 1=standard, 2=godunov, 3=eno, 4=sethian
fire_fuel_left_method=1, ! for now, use 1 only
fire_lfn_ext_up=1.0, ! 0.=extend level set function at boundary by reflection, 1.=always up
- fire_advection=0, ! 0 = cawfe, 1 = use abs speed/slope in spread rate, then project on normal to fireline
/
diff --git a/test/em_fire/namelist.input_two_fires b/test/em_fire/namelist.input_two_fires
index 94e6c6e36b..08899134ef 100644
--- a/test/em_fire/namelist.input_two_fires
+++ b/test/em_fire/namelist.input_two_fires
@@ -94,6 +94,7 @@
time_step_sound = 20, 20, 20,
moist_adv_opt = 1, 1, 1,
scalar_adv_opt = 1, 1, 1,
+ tracer_opt = 3, 3, 3,
/
&bdy_control
@@ -173,15 +174,12 @@
!
! method switches for developers only, do not change!
!
- fire_boundary_guard = -1, ! integer, number of cells to stop when fire close to the domain boundary, -1 turn off
fire_fuel_left_irl=2, ! refinement to integrate fuel_left, must be even
fire_fuel_left_jrl=2, ! refinement to integrate fuel_left, must be even
fire_atm_feedback=1.0, ! real, multiplier for heat fluxes, 1.=normal, 0.=turn off two-way coupling
- fire_back_weight=0.5, ! RK timestepping coefficient, 0=forward, 0.5=Heun
fire_grows_only=1, ! if >0 level set function cannot increase = fire can only grow
fire_viscosity=0.4, ! artificial viscosity in level set method (max 1, needed with fire_upwinding=0)
- fire_upwinding=3, ! 0=none, 1=standard, 2=godunov, 3=eno, 4=sethian
+ fire_upwinding=9, ! 0=none, 1=standard, 2=godunov, 3=eno, 4=sethian
fire_fuel_left_method=1, ! for now, use 1 only
fire_lfn_ext_up=1.0, ! 0.=extend level set function at boundary by reflection, 1.=always up
- fire_advection=0, ! 0 = cawfe, 1 = use abs speed/slope in spread rate, then project on normal to fireline
/