| blob | c22ae6db6314601f5f23cf9b36dc4ffb4eb19ff2 |
1 !
2 ! This module is the entry point from WRF ARW to the wildland
3 ! fire module. The call to sfire_driver advances the fire module by
4 ! one timestep. The fire module inputs the wind and outputs
5 ! temperature and humidity tendencies. The fire module also inputs a
6 ! number of constant arrays (fuel data, topography). Additional
7 ! arguments are model state (for data assimilation) and constant arrays
8 ! the model gives to WRF for safekeeping because it is not allowed
9 ! to save anything.
10 !
11 ! Contributions to this wildland fire module have come from individuals at
12 ! NCAR, the U.S.D.A. Forest Service, the Australian Bureau of Meteorology,
13 ! and the University of Colorado at Denver.
14 !
16 module module_fr_sfire_driver
17 ! use this module for standalone call, you only need to provide some mock-up wrf modules
19 use module_fr_sfire_model
20 use module_fr_sfire_phys, only : fire_params , init_fuel_cats
21 use module_fr_sfire_util
22 use module_fr_sfire_core, only: ignition_line_type
24 USE module_domain, only: domain
25 USE module_configure, only: grid_config_rec_type
26 use module_model_constants, only: reradius, & ! 1/earth radiusw
27 g, & ! gravitation acceleration
28 pi2 ! 2*pi
30 implicit none
33 private
34 public sfire_driver_em,use_atm_vars
36 logical:: use_atm_vars=.true. ! interpolate wind from atm mesh and average output fluxes back
38 contains
40 #define DEBUG_OUT
41 subroutine sfire_driver_em ( grid , config_flags &
42 ,fire_ifun_start,fire_ifun_end,tsteps &
43 ,ids,ide, kds,kde, jds,jde &
44 ,ims,ime, kms,kme, jms,jme &
45 ,ips,ipe, kps,kpe, jps,jpe &
46 ,ifds,ifde, jfds,jfde &
47 ,ifms,ifme, jfms,jfme &
48 ,ifps,ifpe, jfps,jfpe )
50 !*** purpose: driver from grid structure
52 ! Driver layer modules
53 #ifdef DM_PARALLEL
54 USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks
55 USE module_comm_dm , ONLY : halo_fire_fuel_sub, halo_fire_tign_sub, halo_fire_wind_f_sub, &
56 halo_fire_wind_a_sub, halo_fire_ph_sub, halo_fire_zsf_sub, halo_fire_longlat_sub, &
57 halo_fire_phb_sub, halo_fire_z0_sub, halo_fire_lfn_sub
58 #endif
60 implicit none
61 !*** arguments
62 TYPE(domain) , TARGET :: grid ! state
63 TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags ! namelist
64 integer, intent(in):: fire_ifun_start,fire_ifun_end,tsteps ! driver cycle controls
65 integer, intent(in):: &
66 ids,ide, kds,kde, jds,jde, &
67 ims,ime, kms,kme, jms,jme, &
68 ips,ipe, kps,kpe, jps,jpe, &
69 ifds,ifde, jfds,jfde, &
70 ifms,ifme, jfms,jfme, &
71 ifps,ifpe, jfps,jfpe
73 !*** local
74 INTEGER:: fire_num_ignitions
75 integer, parameter::fire_max_ignitions=5
76 TYPE(ignition_line_type), DIMENSION(fire_max_ignitions):: ignition_line
77 integer::fire_ifun,ir,jr,fire_ignition_longlat,istep,itimestep
78 logical::need_lfn_update,restart
79 real, dimension(ifms:ifme, jfms:jfme)::lfn_out
80 real:: corner_ll,corner_ul,corner_ur,corner_lr
81 character(len=128)msg
82 real:: unit_fxlong ,unit_fxlat
83 type(fire_params)::fp
86 !*** executable
88 ! populate our structures from wrf
90 ! pointers to be passed to fire rate of spread formulas
91 fp%vx => grid%uf ! W-E winds used in fire module
92 fp%vy => grid%vf ! S-N winds used in fire module
93 fp%zsf => grid%zsf ! terrain height
94 fp%dzdxf => grid%dzdxf ! terrain grad
95 fp%dzdyf => grid%dzdyf ! terrain grad
96 fp%bbb => grid%bbb ! a rate of spread formula coeff
97 fp%betafl => grid%betafl ! a rate of spread formula variable
98 fp%phiwc => grid%phiwc ! a rate of spread formula coeff
99 fp%r_0 => grid%r_0 ! a rate of spread formula variable
100 fp%fgip => grid%fgip ! a rate of spread formula coeff
101 fp%ischap => grid%ischap ! a rate of spread formula switch
103 ! get ignition data
104 call fire_ignition_convert (config_flags,fire_max_ignitions,fire_ignition_longlat, &
105 ignition_line,fire_num_ignitions,unit_fxlong,unit_fxlat)
107 ! copy configuration flags to sfire internal structures
108 call set_flags(config_flags)
110 ! refinement r
111 ir=grid%sr_x ! refinement ratio
112 jr=grid%sr_y
113 itimestep=grid%itimestep
114 restart=config_flags%restart
118 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
119 write(msg,'(a,i1,a,i1,a,i4)') &
120 'sfire_driver_em: ifun from ',fire_ifun_start,' to ',fire_ifun_end,' test steps',tsteps
121 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
122 call message(msg)
124 do istep=0,tsteps ! istep >0 is for testing only, exit after the first call
125 itimestep = grid%itimestep + istep ! in the first call, do fire_test_steps steps of the fire model
127 need_lfn_update=.false.
128 do fire_ifun=fire_ifun_start,fire_ifun_end
130 ! 1 = initialize run pass 1: interpolate height to zsf=terrain
131 ! 2 = initialize run pass 2: set fuel data, terrain gradient
132 ! 3 = initialize timestep: interpolate winds, check for ignition
133 ! 4 = do one timestep
134 ! 5 = copy timestep output to input
135 ! 6 = compute output fluxes
137 #ifdef DM_PARALLEL
138 if(need_lfn_update)then
139 ! halo exchange on lfn width 2
140 #include "HALO_FIRE_LFN.inc"
141 endif
143 if(fire_ifun.eq.1)then
144 ! halo exchange on topography
145 #include "HALO_FIRE_LONGLAT.inc"
146 !! if(fire_topo_from_atm.eq.1)then
147 !!#include "HALO_FIRE_HT.inc"
148 !! endif
149 ! base geopotential and roughness
150 #include "HALO_FIRE_PHB.inc"
151 #include "HALO_FIRE_Z0.inc"
153 elseif(fire_ifun.eq.2)then
154 ! halo exchange on zsf width 2
155 #include "HALO_FIRE_ZSF.inc"
157 elseif(fire_ifun.eq.3)then
158 ! halo exchange on atm winds and geopotential, width 1 for interpolation
159 #include "HALO_FIRE_WIND_A.inc"
160 #include "HALO_FIRE_PH.inc"
162 elseif(fire_ifun.eq.4)then
163 ! halo exchange on fire winds width 2 for a 2-step RK method
164 #include "HALO_FIRE_WIND_F.inc"
166 elseif(fire_ifun.eq.6)then
167 ! computing fuel_left needs ignition time from neighbors
168 #include "HALO_FIRE_TIGN.inc"
170 endif
171 #endif
172 ! need domain by 1 smaller, in last row.col winds are not set properly
173 call sfire_driver_phys ( &
174 fire_ifun,need_lfn_update, &
175 ids,ide-1, kds,kde, jds,jde-1, &
176 ims,ime, kms,kme, jms,jme, &
177 ips,min(ipe,ide-1), kps,kpe, jps,min(jpe,jde-1), &
178 ifds,ifde-ir, jfds,jfde-jr, &
179 ifms,ifme, jfms,jfme, &
180 ifps,min(ifpe,ifde-ir), jfps,min(jfpe,jfde-jr), &
181 ir,jr, & ! atm/fire grid ratio
182 grid%num_tiles, & ! atm grid tiling
183 grid%i_start,min(grid%i_end,ide-1), &
184 grid%j_start,min(grid%j_end,jde-1), &
185 itimestep,restart,config_flags%fire_fuel_read,config_flags%fire_fuel_cat, & ! in scalars
186 grid%dt,grid%dx,grid%dy, &
187 grid%u_frame,grid%v_frame, &
188 unit_fxlong,unit_fxlat, & ! coordinates of grid center
189 config_flags%fire_ext_grnd,config_flags%fire_ext_crwn,config_flags%fire_crwn_hgt, &
190 config_flags%fire_wind_height, & ! height of wind input to fire spread formula
191 fire_num_ignitions, &
192 fire_ignition_longlat, &
193 ignition_line, &
194 grid%u_2,grid%v_2, & ! atm arrays in
195 grid%ph_2,grid%phb, & ! geopotential
196 grid%z0, & ! roughness height
197 grid%ht, & ! terrain height
198 grid%xlong,grid%xlat, & ! coordinates of atm grid centers, for ignition location
199 grid%lfn,grid%tign_g,grid%fuel_frac, & ! state arrays, fire grid
200 grid%fire_area, & ! redundant, for display, fire grid
201 lfn_out, & ! work - one timestep
202 grid%avg_fuel_frac, & ! out redundant arrays, atm grid
203 grid%grnhfx,grid%grnqfx,grid%canhfx,grid%canqfx, & ! out redundant arrays, atm grid
204 grid%uah,grid%vah, &
205 grid%fgrnhfx,grid%fgrnqfx,grid%fcanhfx,grid%fcanqfx, & ! out redundant arrays, atm grid
206 grid%ros, & ! rate of spread
207 grid%fxlong,grid%fxlat, &
208 grid%nfuel_cat, & ! input, or internal for safekeeping
209 grid%fuel_time, &
210 fp &
211 )
213 #ifdef DM_PARALLEL
214 if(fire_ifun.eq.2)then
215 ! halo exchange on all fuel data width 2
216 #include "HALO_FIRE_FUEL.inc"
217 endif
218 #endif
222 enddo
223 enddo
224 if(tsteps>0)call crash('sfire_driver_em: test run of uncoupled fire model completed')
226 end subroutine sfire_driver_em
228 !
229 !*******************
230 !
232 ! module_fr_sfire_driver%%sfire_driver
233 subroutine sfire_driver_phys (ifun,need_lfn_update, &
234 ids,ide, kds,kde, jds,jde, & ! atm grid dimensions
235 ims,ime, kms,kme, jms,jme, &
236 ips,ipe, kps,kpe, jps,jpe, &
237 ifds, ifde, jfds, jfde, & ! fire grid dimensions
238 ifms, ifme, jfms, jfme, &
239 ifps, ifpe, jfps, jfpe, & ! fire patch in - will use smaller
240 ir,jr, & ! atm/fire grid ratio
241 num_tiles,i_start,i_end,j_start,j_end, & ! atm grid tiling
242 itimestep,restart,ifuelread,nfuel_cat0,dt,dx,dy, & ! in scalars
243 u_frame,v_frame, &
244 unit_fxlong,unit_fxlat, & ! fxlong, fxlat units in m
245 fire_ext_grnd,fire_ext_crwn,fire_crwn_hgt, &
246 fire_wind_height, & ! for vertical wind interpolation
247 num_ignitions, &
248 ignition_longlat, &
249 ignition_line, &
250 u,v, & ! in arrays, atm grid
251 ph,phb, &
252 z0,zs, &
253 xlong,xlat, &
254 lfn,tign,fuel_frac, & ! state arrays, fire grid
255 fire_area, & ! redundant state, fire grid
256 lfn_out, & ! out level set function
257 avg_fuel_frac, &
258 grnhfx,grnqfx,canhfx,canqfx, & ! out redundant arrays, atm grid
259 uah,vah, & ! out atm grid
260 fgrnhfx,fgrnqfx,fcanhfx,fcanqfx, & ! out redundant arrays, fire grid
261 ros, &
262 fxlong,fxlat, & !
263 nfuel_cat, & ! in array, data, fire grid, or constant internal
264 fuel_time, & ! save constant internal data, fire grid
265 fp & ! fire params
266 )
267 USE module_dm, only:wrf_dm_maxval
269 implicit none
271 !*** arguments
273 integer, intent(in)::ifun, &
274 ids,ide, kds,kde, jds,jde, & ! atm domain bounds
275 ims,ime, kms,kme, jms,jme, & ! atm memory bounds
276 ips,ipe, kps,kpe, jps,jpe, & ! atm patch bounds
277 ifds, ifde, jfds, jfde, & ! fire domain bounds
278 ifms, ifme, jfms, jfme, & ! fire memory bounds
279 ifps, ifpe, jfps, jfpe, & ! fire patch bounds
280 ir,jr, & ! atm/fire grid refinement ratio
281 itimestep, & ! number of this timestep
282 ifuelread, & ! how to initialize nfuel_cat:
283 ! -1=not at all, done outside
284 ! 0=from nfuel_cat0
285 ! 1=from altitude
286 ! 2=from file
287 nfuel_cat0, & ! fuel category to initialize everything to
288 num_tiles ! number of tiles
290 logical, intent(in)::restart
293 logical, intent(out)::need_lfn_update
295 integer,dimension(num_tiles),intent(in) :: i_start,i_end,j_start,j_end ! atm grid tiling
297 real, intent(in):: &
298 dt, & ! time step
299 dx,dy, & ! atm grid step
300 u_frame,v_frame, & ! velocity offset
301 unit_fxlong,unit_fxlat, & ! fxlong, fxlat units in m
302 fire_crwn_hgt, & ! lowest height crown fire heat is released (m)
303 fire_ext_grnd, & ! extinction depth of surface fire heat flux (m)
304 fire_ext_crwn, & ! wind height for vertical interploation to fire spread
305 fire_wind_height
308 integer, intent(in):: num_ignitions ! number of ignitions, can be 0
309 integer, intent(in):: ignition_longlat ! if 1 ignition give as long/lat, otherwise as m from lower left corner
310 TYPE (ignition_line_type), DIMENSION(num_ignitions), intent(out):: ignition_line
312 real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::u,v, & ! wind velocity (m/s) (staggered atm grid)
313 ph, phb ! geopotential (w-points atm grid)
314 real,intent(in),dimension(ims:ime, jms:jme):: z0, & ! roughness height
315 zs ! terrain height
316 real,intent(out),dimension(ims:ime,jms:jme)::&
317 uah, & ! atm wind at fire_wind_height, diagnostics
318 vah ! atm wind at fire_wind_height, diagnostics
320 real, dimension(ims:ime, jms:jme), intent(inout)::xlong, xlat ! inout because of extension at bdry
322 real, intent(inout), dimension(ifms:ifme,jfms:jfme):: &
323 nfuel_cat ! fuel data; can be also set inside (cell based, fire grid)
325 real, intent(inout), dimension(ifms:ifme, jfms:jfme):: &
326 lfn,tign,fuel_frac, & ! state: level function, ign time, fuel left
327 lfn_out ! fire wind velocities
329 real, intent(out), dimension(ifms:ifme, jfms:jfme):: &
330 fire_area ! fraction of each cell burning
332 real, intent(out), dimension(ims:ime, jms:jme):: & ! redundant arrays, for display purposes only (atm grid)
333 avg_fuel_frac, & ! average fuel fraction
334 grnhfx, & ! heat flux from surface fire (W/m^2)
335 grnqfx, & ! moisture flux from surface fire (W/m^2)
336 canhfx, & ! heat flux from crown fire (W/m^2)
337 canqfx ! moisture flux from crown fire (W/m^2)
339 real, intent(out), dimension(ifms:ifme, jfms:jfme):: & ! redundant arrays, for display only, fire grid
340 fgrnhfx, & ! heat flux from surface fire (W/m^2)
341 fgrnqfx, & ! moisture flux from surface fire (W/m^2)
342 fcanhfx, & ! heat flux from crown fire (W/m^2)
343 fcanqfx, & ! moisture flux from crown fire (W/m^2)
344 ros ! fire rate of spread (m/s)
346 ! ***** data (constant in time) *****
348 real, dimension(ifms:ifme, jfms:jfme), intent(inout)::fxlong,fxlat ! fire mesh coordinates
349 real, intent(out), dimension(ifms:ifme, jfms:jfme)::fuel_time ! fire params arrays
351 type(fire_params),intent(inout)::fp
353 !*** local
354 real :: dxf,dyf,time_start,latm, s
355 integer :: its,ite,jts,jte,kts,kte, & ! tile
356 ij,i,j,k,id,pid,ipe1,jpe1,ite1,jte1, &
357 ifts,ifte,jfts,jfte ! fire tile
358 character(len=128)::msg
359 character(len=3)::kk
360 integer::ignitions_done_tile(num_tiles),ignited_tile(num_ignitions,num_tiles)
361 integer::ignitions_done,ignited_patch(num_ignitions),idex,jdex
364 !*** executable
366 ! time - assume dt does not change
367 time_start = itimestep * dt
369 ! fire mesh step
370 dxf=dx/ir
371 dyf=dy/jr
375 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
376 write(msg,'(a,2f15.6)')'atmosphere mesh step:',dx,dy
377 call message(msg)
378 write(msg,'(a,2f15.6)')'fire mesh step: ',dxf,dyf
379 call message(msg)
380 write(msg,7001)'atm domain ','ids',ids,ide,jds,jde
381 call message(msg)
382 write(msg,7001)'atm memory ','ims',ims,ime,jms,jme
383 call message(msg)
384 write(msg,7001)'atm patch ','ips',ips,ipe,jps,jpe
385 call message(msg)
386 write(msg,7001)'fire domain ','ifds',ifds,ifde,jfds,jfde
387 call message(msg)
388 write(msg,7001)'fire memory ','ifms',ifms,ifme,jfms,jfme
389 call message(msg)
390 write(msg,7001)'fire patch ','ifps',ifps,ifpe,jfps,jfpe
391 call message(msg)
392 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
394 ! check mesh dimensions
395 call check_fmesh(ids,ide,ifds,ifde,ir,'id') ! check if atm and fire grids line up
396 call check_fmesh(jds,jde,jfds,jfde,jr,'jd')
397 call check_fmesh(ips,ipe,ifps,ifpe,ir,'ip')
398 call check_fmesh(jps,jpe,jfps,jfpe,jr,'jp')
399 call check_mesh_2dim(ips,ipe,jps,jpe,ims,ime,jms,jme) ! check if atm patch fits in atm array
400 call check_mesh_2dim(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme) ! check if fire patch fits in fire array
401 call check_mesh_2dim(ips,ipe,jps,jpe,ids,ide,jds,jde) ! check if atm patch fits in atm domain
402 call check_mesh_2dim(ifps,ifpe,jfps,jfpe,ifds,ifde,jfds,jfde) ! check if fire patch fits in fire domain
404 !$OMP SINGLE
405 ! init rest of fuel tables with derived constants
406 if(ifun.eq.1) then
407 call init_fuel_cats ! common for all threads
408 endif
409 !$OMP END SINGLE
412 pid=0
413 if(fire_print_file.gt.0)then
414 if(itimestep.le.fire_print_file.or.mod(itimestep,fire_print_file).eq.0)pid=itimestep ! print 1-fire_print_file then every fire_print_file-th
415 endif
417 if(ifun.eq.3)then
418 call print_chsum(itimestep,ims,ime,kms,kme,jms,jme,ids,ide,kds,kde,jds,jde,ips,ipe,kps,kpe,jps,jpe,1,0,0,u,'u')
419 call print_chsum(itimestep,ims,ime,kms,kme,jms,jme,ids,ide,kds,kde,jds,jde,ips,ipe,kps,kpe,jps,jpe,0,0,1,v,'v')
420 call print_chsum(itimestep,ims,ime,kms,kme,jms,jme,ids,ide,kds,kde,jds,jde,ips,ipe,kps,kpe,jps,jpe,0,1,0,ph,'ph')
421 endif
423 ! fake atm tile bounds
424 kts=kps
425 kte=kpe
427 ! staggered atm patch bounds
428 ipe1=ifval(ipe.eq.ide,ipe+1,ipe)
429 jpe1=ifval(jpe.eq.jde,jpe+1,jpe)
431 ! set up fire tiles & interpolate to fire grid
432 !$OMP PARALLEL DO PRIVATE(ij,its,ite,jts,jte,ite1,jte1,ifts,ifte,jfts,jfte,msg,id) &
433 !$OMP SCHEDULE(STATIC)
434 do ij=1,num_tiles
436 id = ifval(pid.ne.0,pid+ij*10000,0) ! for print
438 ignitions_done_tile(ij)=0
440 ! set up tile bounds
441 its = i_start(ij) ! start atmospheric tile in i
442 ite = i_end(ij) ! end atmospheric tile in i
443 jts = j_start(ij) ! start atmospheric tile in j
444 jte = j_end(ij) ! end atmospheric tile in j
445 ifts= (its-ids)*ir+ifds ! start fire tile in i
446 ifte= (ite-ids+1)*ir+ifds-1 ! end fire tile in i
447 jfts= (jts-jds)*jr+jfds ! start fire tile in j
448 jfte= (jte-jds+1)*jr+jfds-1 ! end fire tile in j
450 ! staggered atm tile bounds
451 ite1=ifval(ite.eq.ide,ite+1,ite)
452 jte1=ifval(jte.eq.jde,jte+1,jte)
454 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
455 write(msg,'(a,i3,1x,a,i7,1x,a,i3)')'tile=',ij,' id=',id,' ifun=',ifun
456 call message(msg)
457 write(msg,7001)'atm tile ','its',its,ite,jts,jte
458 call message(msg)
459 write(msg,7001)'fire tile ','ifts',ifts,ifte,jfts,jfte
460 call message(msg)
461 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
463 ! check the tiles
464 call check_mesh_2dim(its,ite,jts,jte,ips,ipe,jps,jpe) ! check if atm tile fits in atm patch
465 call check_mesh_2dim(ifts,ifte,jfts,jfte,ifps,ifpe,jfps,jfpe) ! check if fire tile fits in fire patch
466 call check_mesh_2dim(ifts-2,ifte+2,jfts-2,jfte+2,ifms,ifme,jfms,jfme)! check if fire node tile fits in memory
469 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
470 write(msg,'(a,i6,a,2(f15.6,a))')'time step',itimestep,' at',time_start,' duration',dt,'s'
471 call message(msg)
472 7001 format(a,' dimensions ',a4,':',i6,' to ',i6,' by ',i6,' to ',i6)
473 write(msg,'(a,2i9)')'refinement ratio:',ir,jr
474 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
476 if(ifun.eq.1)then ! set terrain
478 if(restart)then
480 call message('restart - topo initialization skipped')
482 else
483 !
484 ! call print_2d_stats(ips,ipe,jps,jpe,ims,ime,jms,jme,zs,'driver:zs')
485 !
486 ! ! interpolate terrain height
487 ! if(fire_topo_from_atm.eq.1)then
488 ! call interpolate_z2fire(id, & ! for debug output, <= 0 no output
489 ! ids,ide, jds,jde, & ! atm grid dimensions
490 ! ims,ime, jms,jme, &
491 ! ips,ipe,jps,jpe, &
492 ! its,ite,jts,jte, &
493 ! ifds, ifde, jfds, jfde, & ! fire grid dimensions
494 ! ifms, ifme, jfms, jfme, &
495 ! ifts,ifte,jfts,jfte, &
496 ! ir,jr, & ! atm/fire grid ratio
497 ! zs, & ! atm grid arrays in
498 ! fp%zsf) ! fire grid arrays out
499 ! else
500 !!$OMP CRITICAL(SFIRE_DRIVER_CRIT)
501 ! write(msg,'(a,i3,a)')'fire_topo_from_atm=',fire_topo_from_atm,' assuming ZSF set, interpolation skipped'
502 !!$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
503 ! endif
505 if(ignition_longlat .eq.0)then
506 ! set ideal fire mesh coordinates - used for ignition only
507 ! do not forget to set unit_fxlong, unit_fxlat outside of parallel loop
508 !call set_ideal_coord( dxf,dyf, &
509 ! ifds,ifde,jfds,jfde, &
510 ! ifms,ifme,jfms,jfme, &
511 ! ifts,ifte,jfts,jfte, &
512 ! fxlong,fxlat )
513 !call set_ideal_coord( dx,dy, &
514 ! ids,ide,jds,jde, &
515 ! ims,ime,jms,jme, &
516 ! its,ite,jts,jte, &
517 ! xlong,xlat )
518 elseif(use_atm_vars)then
519 ! assume halo xlong xlat
520 ! interpolate nodal coordinates
522 #ifdef DEBUG_OUT
523 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,xlat,'xlat',id)
524 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,xlong,'xlong',id)
525 #endif
526 call interpolate_z2fire(id, & ! for debug output, <= 0 no output
527 ids,ide, jds,jde, & ! atm grid dimensions
528 ims,ime, jms,jme, &
529 ips,ipe,jps,jpe, &
530 its,ite,jts,jte, &
531 ifds, ifde, jfds, jfde, & ! fire grid dimensions
532 ifms, ifme, jfms, jfme, &
533 ifts,ifte,jfts,jfte, &
534 ir,jr, & ! atm/fire grid ratio
535 xlat, & ! atm grid arrays in
536 fxlat) ! fire grid arrays out
538 call interpolate_z2fire(id, & ! for debug output, <= 0 no output
539 ids,ide, jds,jde, & ! atm grid dimensions
540 ims,ime, jms,jme, &
541 ips,ipe,jps,jpe, &
542 its,ite,jts,jte, &
543 ifds, ifde, jfds, jfde, & ! fire grid dimensions
544 ifms, ifme, jfms, jfme, &
545 ifts,ifte,jfts,jfte, &
546 ir,jr, & ! atm/fire grid ratio
547 xlong, & ! atm grid arrays in
548 fxlong) ! fire grid arrays out
550 endif
552 endif
554 elseif(ifun.eq.2)then
556 ! after the loop where zsf created exited and all synced
557 call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fp%zsf,'driver_phys:zsf')
559 elseif(ifun.eq.3)then ! interpolate winds to the fire grid
561 if(use_atm_vars)then
563 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,z0,'z0',id)
564 call write_array_m3(its,ite1,kts,kde-1,jts,jte,ims,ime,kms,kme,jms,jme,u,'u_2',id)
565 call write_array_m3(its,ite,kts,kde-1,jts,jte1,ims,ime,kms,kme,jms,jme,v,'v_2',id)
566 call write_array_m3(its,ite,kts,kde,jts,jte,ims,ime,kms,kme,jms,jme,ph,'ph_2',id)
567 call write_array_m3(its,ite,kts,kde,jts,jte,ims,ime,kms,kme,jms,jme,phb,'phb',id)
568 call interpolate_atm2fire(id, & ! flag for debug output
569 fire_wind_height, & ! height to interpolate to
570 ids,ide, kds,kde, jds,jde, & ! atm grid dimensions
571 ims,ime, kms,kme, jms,jme, &
572 ips,ipe, jps,jpe, &
573 its,ite,jts,jte, &
574 ifds, ifde, jfds, jfde, & ! fire grid dimensions
575 ifms, ifme, jfms, jfme, &
576 ifps, ifpe, jfps, jfpe, & ! fire patch bounds
577 ifts, ifte, jfts, jfte, &
578 ir,jr, & ! atm/fire grid ratio
579 u_frame, v_frame, & ! velocity frame correction
580 u,v, & ! 3D atm grid arrays in
581 ph,phb, &
582 z0,zs, & ! 2D atm grid arrays in
583 uah,vah, & ! 2D atm grid out
584 fp%vx,fp%vy) ! fire grid arrays out
585 endif
587 endif
589 ! the following executes in any case
591 call sfire_model (id,ifun,restart,need_lfn_update, &
592 num_ignitions, & ! switches
593 ifuelread,nfuel_cat0, & ! initialize fuel categories
594 ifds,ifde,jfds,jfde, & ! fire domain dims
595 ifms,ifme,jfms,jfme, & ! fire memory dims
596 ifps,ifpe,jfps,jfpe, &
597 ifts,ifte,jfts,jfte, & ! fire patch dims
598 time_start,dt, & ! time and increment
599 dxf,dyf, & ! fire mesh spacing
600 ignition_line, & ! description of ignition lines
601 ignitions_done_tile(ij),ignited_tile(1,ij), &
602 fxlong,fxlat,unit_fxlong,unit_fxlat, & ! fire mesh coordinates
603 lfn,lfn_out,tign,fuel_frac, & ! state: level function, ign time, fuel left
604 fire_area, & ! output: fraction of cell burning
605 fgrnhfx,fgrnqfx, & ! output: heat fluxes
606 ros, & ! output: rate of spread for display
607 nfuel_cat, & ! fuel data per point
608 fuel_time, & ! save derived internal data
609 fp & ! fire coefficients
610 )
612 if(ifun.eq.6)then ! heat fluxes into the atmosphere
614 call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fgrnhfx,'fire_driver:fgrnhfx')
615 call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fgrnqfx,'fire_driver:fgrnqfx')
617 ! sum the fluxes over atm cells
618 if(use_atm_vars)then
619 call sum_2d_cells( &
620 ifms,ifme,jfms,jfme, &
621 ifts,ifte,jfts,jfte, &
622 fuel_frac, &
623 ims, ime, jms, jme, &
624 its,ite,jts,jte, &
625 avg_fuel_frac)
626 call sum_2d_cells( &
627 ifms,ifme,jfms,jfme, &
628 ifts,ifte,jfts,jfte, &
629 fgrnhfx, &
630 ims, ime, jms, jme, &
631 its,ite,jts,jte, &
632 grnhfx)
633 !comment out the next call to get results as before commit 55fd92051196b796891b60cb7ec1c4bdb8800078
634 call sum_2d_cells( &
635 ifms,ifme,jfms,jfme, &
636 ifts,ifte,jfts,jfte, &
637 fgrnqfx, &
638 ims, ime, jms, jme, &
639 its,ite,jts,jte, &
640 grnqfx)
642 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
643 write(msg,'(a,f6.3)')'fire-atmosphere feedback scaling ',fire_atm_feedback
644 !$OMP end CRITICAL(SFIRE_DRIVER_CRIT)
645 call message(msg)
646 s = 1./(ir*jr)
647 do j=jts,jte
648 do i=its,ite
649 ! scale surface fluxes to get the averages
650 avg_fuel_frac(i,j)=avg_fuel_frac(i,j)*s
651 grnhfx(i,j)=fire_atm_feedback*grnhfx(i,j)*s
652 grnqfx(i,j)=fire_atm_feedback*grnqfx(i,j)*s
653 ! we do not have canopy fluxes yet...
654 canhfx(i,j)=0
655 canqfx(i,j)=0
656 enddo
657 enddo
659 call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnhfx,'fire_driver:grnhfx')
660 call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnqfx,'fire_driver:grnqfx')
661 endif
663 endif ! ifun=6
665 enddo ! tiles
666 !$OMP END PARALLEL DO
668 #ifdef DEBUG_OUT
669 if(ifun.eq.1)then
670 if(pid.ne.0)then
671 call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,zs,'zs',pid)
672 call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fp%zsf,'zsf',pid)
673 endif
674 endif
675 #endif
677 if (ifun.eq.3)then
678 ignitions_done=ignitions_done_tile(1) ! all should be the same
679 do i=1,ignitions_done
680 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
681 write(msg,'(2(a,i4,1x))')'sfire_driver_phys: checking ignition ',i,' of ',ignitions_done
682 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
683 call message(msg)
684 ignited_patch(i)=0
685 do ij=1,num_tiles
686 ignited_patch(i)=ignited_patch(i)+ignited_tile(i,ij)
687 enddo
688 #ifdef DM_PARALLEL
689 call message('sfire_driver_phys: checking ignitions, collect counts')
690 call wrf_dm_maxval(ignited_patch(i),idex,jdex)
691 call message('sfire_driver_phys: collected')
692 #endif
693 if(ignited_patch(i).eq.0)then
694 call crash('sfire_driver_phys: Ignition failed, no nodes ignited. Bad coordinates?')
695 endif
696 enddo
698 call print_chsum(itimestep,ims,ime,1,1,jms,jme,ids,ide,1,1,jds,jde,ips,ipe,1,1,jps,jpe,1,0,0,uah,'uah')
699 call print_chsum(itimestep,ims,ime,1,1,jms,jme,ids,ide,1,1,jds,jde,ips,ipe,1,1,jps,jpe,0,0,1,vah,'vah')
700 call print_chsum(itimestep,ifms,ifme,1,1,jfms,jfme,ifds,ifde,1,1,jfds,jfde,ifps,ifpe,1,1,jfps,jfpe,0,0,0,fp%vx,'uf')
701 call print_chsum(itimestep,ifms,ifme,1,1,jfms,jfme,ifds,ifde,1,1,jfds,jfde,ifps,ifpe,1,1,jfps,jfpe,0,0,0,fp%vy,'vf')
702 #ifdef DEBUG_OUT
703 if(pid.gt.0)then
704 call write_array_m(ips,ipe1,jps,jpe,ims,ime,jms,jme,uah,'uah',pid)
705 call write_array_m(ips,ipe,jps,jpe1,ims,ime,jms,jme,vah,'vah',pid)
706 call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnhfx,'grnhfx',pid)
707 call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnqfx,'grnqfx',pid)
708 call write_array_m3(ips,ipe1,kds,kde+1,jps,jpe,ims,ime,kms,kme,jms,jme,u,'u',pid)
709 call write_array_m3(ips,ipe,kds,kde+1,jps,jpe1,ims,ime,kms,kme,jms,jme,v,'v',pid)
710 call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fp%vx,'uf',pid)
711 call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fp%vy,'vf',pid)
712 endif
713 #endif
714 endif
716 if(ifun.eq.5)then
717 #ifdef DEBUG_OUT
718 if(pid.gt.0)then
719 call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,lfn,'lfn',pid)
720 call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,tign,'tign',pid)
721 endif
722 #endif
723 endif
725 if(ifun.eq.6)then
726 call print_chsum(itimestep,ifms,ifme,1,1,jfms,jfme,ifds,ifde,1,1,jfds,jfde,ifps,ifpe,1,1,jfps,jfpe,0,0,0,fgrnhfx,'fgrnhfx')
727 call print_chsum(itimestep,ifms,ifme,1,1,jfms,jfme,ifds,ifde,1,1,jfds,jfde,ifps,ifpe,1,1,jfps,jfpe,0,0,0,fgrnqfx,'fgrnqfx')
728 call print_chsum(itimestep,ims,ime,1,1,jms,jme,ids,ide,1,1,jds,jde,ips,ipe,1,1,jps,jpe,0,0,0,grnhfx,'grnhfx')
729 call print_chsum(itimestep,ims,ime,1,1,jms,jme,ids,ide,1,1,jds,jde,ips,ipe,1,1,jps,jpe,0,0,0,grnqfx,'grnqfx')
730 #ifdef DEBUG_OUT
731 if(pid.gt.0)then
732 call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnhfx,'grnhfx',pid)
733 call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnqfx,'grnqfx',pid)
734 call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fuel_frac,'fuel_frac',pid)
735 call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fgrnhfx,'fgrnhfx',pid)
736 call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fgrnqfx,'fgrnqfx',pid)
737 endif
738 #endif
739 endif
741 end subroutine sfire_driver_phys
742 !
743 !*******************
744 !
746 subroutine fire_ignition_convert (config_flags,fire_max_ignitions,fire_ignition_longlat, &
747 ignition_line,fire_num_ignitions,unit_fxlong,unit_fxlat)
748 USE module_configure, only : grid_config_rec_type
749 implicit none
750 ! create ignition arrays from scalar flags
751 !*** arguments
752 TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags
753 integer, intent(in)::fire_max_ignitions
754 TYPE (ignition_line_type), DIMENSION(fire_max_ignitions), intent(out):: ignition_line ! any values from input discarded
755 integer, intent(out)::fire_num_ignitions,fire_ignition_longlat
756 real, intent(out)::unit_fxlong,unit_fxlat
757 !*** local
758 integer::i
759 logical:: real,ideal
760 real::lat_ctr,lon_ctr
761 !*** executable
762 ! this is only until I figure out how to input arrays through the namelist...
763 if(fire_max_ignitions.lt.5)call crash('fire_max_ignitions too small')
764 ! figure out which kind of coordinates from the first given
765 ideal=config_flags%fire_ignition_start_x1 .ne.0. .or. config_flags%fire_ignition_start_y1 .ne. 0.
766 real=config_flags%fire_ignition_start_lon1 .ne. 0. .or. config_flags%fire_ignition_start_lat1 .ne. 0.
767 if(ideal)call message('Using ideal ignition coordinates, m from the lower left domain corner')
768 if(real)call message('Using real ignition coordinates, longitude and latitude')
769 if(ideal.and.real)call crash('Only one of the ideal or real coordinates may be given')
771 fire_ignition_longlat=0 ! default, if no ignition
772 if(ideal)then
773 ! use values from _x and _y variables
774 fire_ignition_longlat=0
775 ignition_line(1)%start_x=config_flags%fire_ignition_start_x1
776 ignition_line(1)%start_y=config_flags%fire_ignition_start_y1
777 ignition_line(1)%end_x=config_flags%fire_ignition_end_x1
778 ignition_line(1)%end_y=config_flags%fire_ignition_end_y1
779 ignition_line(2)%start_x=config_flags%fire_ignition_start_x2
780 ignition_line(2)%start_y=config_flags%fire_ignition_start_y2
781 ignition_line(2)%end_x=config_flags%fire_ignition_end_x2
782 ignition_line(2)%end_y=config_flags%fire_ignition_end_y2
783 ignition_line(3)%start_x=config_flags%fire_ignition_start_x3
784 ignition_line(3)%start_y=config_flags%fire_ignition_start_y3
785 ignition_line(3)%end_x=config_flags%fire_ignition_end_x3
786 ignition_line(3)%end_y=config_flags%fire_ignition_end_y3
787 ignition_line(4)%start_x=config_flags%fire_ignition_start_x4
788 ignition_line(4)%start_y=config_flags%fire_ignition_start_y4
789 ignition_line(4)%end_x=config_flags%fire_ignition_end_x4
790 ignition_line(4)%end_y=config_flags%fire_ignition_end_y4
791 ignition_line(5)%start_x=config_flags%fire_ignition_start_x5
792 ignition_line(5)%start_y=config_flags%fire_ignition_start_y5
793 ignition_line(5)%end_x=config_flags%fire_ignition_end_x5
794 ignition_line(5)%end_y=config_flags%fire_ignition_end_y5
795 endif
796 if(real)then
797 ! use values from _long and _lat
798 fire_ignition_longlat=1
799 ignition_line(1)%start_x=config_flags%fire_ignition_start_lon1
800 ignition_line(1)%start_y=config_flags%fire_ignition_start_lat1
801 ignition_line(1)%end_x=config_flags%fire_ignition_end_lon1
802 ignition_line(1)%end_y=config_flags%fire_ignition_end_lat1
803 ignition_line(2)%start_x=config_flags%fire_ignition_start_lon2
804 ignition_line(2)%start_y=config_flags%fire_ignition_start_lat2
805 ignition_line(2)%end_x=config_flags%fire_ignition_end_lon2
806 ignition_line(2)%end_y=config_flags%fire_ignition_end_lat2
807 ignition_line(3)%start_x=config_flags%fire_ignition_start_lon3
808 ignition_line(3)%start_y=config_flags%fire_ignition_start_lat3
809 ignition_line(3)%end_x=config_flags%fire_ignition_end_lon3
810 ignition_line(3)%end_y=config_flags%fire_ignition_end_lat3
811 ignition_line(4)%start_x=config_flags%fire_ignition_start_lon4
812 ignition_line(4)%start_y=config_flags%fire_ignition_start_lat4
813 ignition_line(4)%end_x=config_flags%fire_ignition_end_lon4
814 ignition_line(4)%end_y=config_flags%fire_ignition_end_lat4
815 ignition_line(5)%start_x=config_flags%fire_ignition_start_lon5
816 ignition_line(5)%start_y=config_flags%fire_ignition_start_lat5
817 ignition_line(5)%end_x=config_flags%fire_ignition_end_lon5
818 ignition_line(5)%end_y=config_flags%fire_ignition_end_lat5
819 endif
820 ! common to both cases
821 ignition_line(1)%ros=config_flags%fire_ignition_ros1
822 ignition_line(1)%radius=config_flags%fire_ignition_radius1
823 ignition_line(1)%start_time=config_flags%fire_ignition_start_time1
824 ignition_line(1)%end_time=config_flags%fire_ignition_end_time1
825 ignition_line(2)%ros=config_flags%fire_ignition_ros2
826 ignition_line(2)%radius=config_flags%fire_ignition_radius2
827 ignition_line(2)%start_time=config_flags%fire_ignition_start_time2
828 ignition_line(2)%end_time=config_flags%fire_ignition_end_time2
829 ignition_line(3)%ros=config_flags%fire_ignition_ros3
830 ignition_line(3)%radius=config_flags%fire_ignition_radius3
831 ignition_line(3)%start_time=config_flags%fire_ignition_start_time3
832 ignition_line(3)%end_time=config_flags%fire_ignition_end_time3
833 ignition_line(4)%ros=config_flags%fire_ignition_ros4
834 ignition_line(4)%radius=config_flags%fire_ignition_radius4
835 ignition_line(4)%start_time=config_flags%fire_ignition_start_time4
836 ignition_line(4)%end_time=config_flags%fire_ignition_end_time4
837 ignition_line(5)%ros=config_flags%fire_ignition_ros5
838 ignition_line(5)%radius=config_flags%fire_ignition_radius5
839 ignition_line(5)%start_time=config_flags%fire_ignition_start_time5
840 ignition_line(5)%end_time=config_flags%fire_ignition_end_time5
842 !
843 fire_num_ignitions=0
844 do i=1,min(5,config_flags%fire_num_ignitions)
845 ! count the ignitions
846 if(ignition_line(i)%radius.gt.0.)fire_num_ignitions=i
847 ! expand ignition data given as zero
848 if(ignition_line(i)%end_x.eq.0.)ignition_line(i)%end_x=ignition_line(i)%start_x
849 if(ignition_line(i)%end_y.eq.0.)ignition_line(i)%end_y=ignition_line(i)%start_y
850 if(ignition_line(i)%end_time.eq.0.)ignition_line(i)%end_time=ignition_line(i)%start_time
851 enddo
853 if(fire_ignition_longlat .eq. 0)then
854 ! ideal
855 ! ignition is in m
856 unit_fxlong=1.
857 unit_fxlat=1.
858 ! will set fire mesh coordinates to uniform mesh below
859 else
860 ! real
861 lat_ctr=config_flags%cen_lat
862 lon_ctr=config_flags%cen_lon
863 ! 1 degree in m (approximate OK)
864 unit_fxlat=pi2/(360.*reradius) ! earth circumference in m / 360 degrees
865 unit_fxlong=cos(lat_ctr*pi2/360.)*unit_fxlat ! latitude
866 endif
868 end subroutine fire_ignition_convert
870 !
871 !*****************************
872 !
874 subroutine interpolate_z2fire(id, & ! for debug output, <= 0 no output
875 ids,ide, jds,jde, & ! atm grid dimensions
876 ims,ime, jms,jme, &
877 ips,ipe,jps,jpe, &
878 its,ite,jts,jte, &
879 ifds, ifde, jfds, jfde, & ! fire grid dimensions
880 ifms, ifme, jfms, jfme, &
881 ifts,ifte,jfts,jfte, &
882 ir,jr, & ! atm/fire grid ratio
883 zs, & ! atm grid arrays in
884 zsf) ! fire grid arrays out
886 implicit none
887 !*** purpose: interpolate height
889 !*** arguments
890 integer, intent(in)::id, &
891 ids,ide, jds,jde, & ! atm domain bounds
892 ims,ime,jms,jme, & ! atm memory bounds
893 ips,ipe,jps,jpe, &
894 its,ite,jts,jte, & ! atm tile bounds
895 ifds, ifde, jfds, jfde, & ! fire domain bounds
896 ifms, ifme, jfms, jfme, & ! fire memory bounds
897 ifts,ifte,jfts,jfte, & ! fire tile bounds
898 ir,jr ! atm/fire grid refinement ratio
899 real, intent(in), dimension(ims:ime, jms:jme):: zs ! terrain height at atm cell centers & ! terrain height
900 real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
901 zsf ! terrain height fire grid nodes
904 !*** local
905 real, dimension(its-2:ite+2,jts-2:jte+2):: za ! terrain height
906 integer:: i,j,jts1,jte1,its1,ite1,jfts1,jfte1,ifts1,ifte1,itso,jtso,iteo,jteo
908 ! terrain height
910 jts1=max(jts-1,jds) ! lower loop limit by one less when at end of domain
911 its1=max(its-1,ids) ! ASSUMES THE HALO IS THERE if patch != domain
912 jte1=min(jte+1,jde)
913 ite1=min(ite+1,ide)
914 do j = jts1,jte1
915 do i = its1,ite1
916 ! copy to local array
917 za(i,j)=zs(i,j)
918 enddo
919 enddo
921 call continue_at_boundary(1,1,0., & ! do x direction or y direction
922 its-2,ite+2,jts-2,jte+2, & ! memory dims
923 ids,ide,jds,jde, & ! domain dims - winds defined up to +1
924 ips,ipe,jps,jpe, & ! patch dims - winds defined up to +1
925 its1,ite1,jts1,jte1, & ! tile dims
926 itso,jtso,iteo,jteo, &
927 za) ! array
929 ! interpolate to tile plus strip along domain boundary if at boundary
930 jfts1=snode(jfts,jfds,-1) ! lower loop limit by one less when at end of domain
931 ifts1=snode(ifts,ifds,-1)
932 jfte1=snode(jfte,jfde,+1)
933 ifte1=snode(ifte,ifde,+1)
935 call interpolate_2d( &
936 its-2,ite+2,jts-2,jte+2, & ! memory dims atm grid tile
937 its1-1,ite1+1,jts1-1,jte1+1, & ! where atm grid values set
938 ifms,ifme,jfms,jfme, & ! array dims fire grid
939 ifts1,ifte1,jfts1,jfte1, & ! dimensions fire grid tile
940 ir,jr, & ! refinement ratio
941 real(ids),real(jds),ifds+(ir-1)*0.5,jfds+(jr-1)*0.5, & ! line up by lower left corner of domain
942 za, & ! in atm grid
943 zsf) ! out fire grid
945 end subroutine interpolate_z2fire
946 !
947 !*****************************
948 !
950 subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no output
951 fire_wind_height, & ! interpolation height
952 ids,ide, kds,kde, jds,jde, & ! atm grid dimensions
953 ims,ime, kms,kme, jms,jme, &
954 ips,ipe,jps,jpe, &
955 its,ite,jts,jte, &
956 ifds, ifde, jfds, jfde, & ! fire grid dimensions
957 ifms, ifme, jfms, jfme, &
958 ifps, ifpe, jfps, jfpe, & ! fire patch bounds
959 ifts,ifte,jfts,jfte, &
960 ir,jr, & ! atm/fire grid ratio
961 u_frame, v_frame, & ! velocity frame correction
962 u,v, & ! atm grid arrays in
963 ph,phb, &
964 z0,zs, &
965 uah,vah, &
966 uf,vf) ! fire grid arrays out
968 implicit none
969 !*** purpose: interpolate winds and height
971 !*** arguments
972 integer, intent(in)::id
973 real, intent(in):: fire_wind_height ! height above the terrain for vertical interpolation
974 integer, intent(in):: &
975 ids,ide, kds,kde, jds,jde, & ! atm domain bounds
976 ims,ime, kms,kme, jms,jme, & ! atm memory bounds
977 ips,ipe,jps,jpe, &
978 its,ite,jts,jte, & ! atm tile bounds
979 ifds, ifde, jfds, jfde, & ! fire domain bounds
980 ifms, ifme, jfms, jfme, & ! fire memory bounds
981 ifps, ifpe, jfps, jfpe, & ! fire patch bounds
982 ifts,ifte,jfts,jfte, & ! fire tile bounds
983 ir,jr ! atm/fire grid refinement ratio
984 real, intent(in):: u_frame, v_frame ! velocity frame correction
985 real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::&
986 u,v, & ! atm wind velocity, staggered
987 ph,phb ! geopotential
988 real,intent(in),dimension(ims:ime,jms:jme)::&
989 z0, & ! roughness height
990 zs ! terrain height
991 real,intent(out),dimension(ims:ime,jms:jme)::&
992 uah, & ! atm wind at fire_wind_height, diagnostics
993 vah !
994 real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
995 uf,vf ! wind velocity fire grid nodes
998 !*** local
999 character(len=256)::msg
1000 #define TDIMS its-2,ite+2,jts-2,jte+2
1001 real, dimension(its-2:ite+2,jts-2:jte+2):: ua,va ! atm winds, interpolated over height, still staggered grid
1002 real, dimension(its-2:ite+2,kds:kde,jts-2:jte+2):: altw,altub,altvb,hgtu,hgtv ! altitudes
1003 integer:: i,j,k,ifts1,ifte1,jfts1,jfte1,ite1,jte1
1004 integer::itst,itet,jtst,jtet,itsu,iteu,jtsu,jteu,itsv,itev,jtsv,jtev
1005 integer::kdmax,its1,jts1,ips1,jps1
1006 integer::itsou,iteou,jtsou,jteou,itsov,iteov,jtsov,jteov
1007 real:: ground,loght,loglast,logz0,logfwh,ht,zr
1008 real::r_nan
1009 integer::i_nan
1010 equivalence (i_nan,r_nan)
1012 !*** executable
1014 ! debug init local arrays
1015 i_nan=2147483647
1016 ua=r_nan
1017 va=r_nan
1018 altw=r_nan
1019 altub=r_nan
1020 hgtu=r_nan
1021 hgtv=r_nan
1024 if(kds.ne.1)then
1025 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1026 write(msg,*)'WARNING: bottom index kds=',kds,' should be 1?'
1027 call message(msg)
1028 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1029 endif
1031 ! ^ j
1032 ! ------------ |
1033 ! | | ----> i
1034 ! u p |
1035 ! | | nodes in cell (i,j)
1036 ! ------v----- view from top
1037 !
1038 ! v(ide,jde+1)
1039 ! -------x------
1040 ! | |
1041 ! | |
1042 ! u(ide,jde) x x x u(ide+1,jde)
1043 ! | p(ide,hde) |
1044 ! | | p=ph,phb,z0,...
1045 ! -------x------
1046 ! v(ide,jde)
1047 !
1048 ! staggered values set u(ids:ide+1,jds:jde) v(ids:ide,jds:jde+1)
1049 ! p=ph+phb set at (ids:ide,jds:jde)
1050 ! location of u(i,j) needs p(i-1,j) and p(i,j)
1051 ! location of v(i,j) needs p(i,j-1) and p(i,j)
1052 ! *** NOTE need HALO in ph, phb
1053 ! so we can compute only u(ids+1:ide,jds:jde) v(ids:ide,jds+1,jde)
1054 ! unless we extend p at the boundary
1055 ! but because we care about the fire way in the inside it does not matter
1056 ! if the fire gets close to domain boundary the simulation is over anyway
1058 ite1=snode(ite,ide,1)
1059 jte1=snode(jte,jde,1)
1060 ! do this in any case to check for nans
1061 call print_3d_stats(its,ite1,kds,kde,jts,jte,ims,ime,kms,kme,jms,jme,u,'wind U in')
1062 call print_3d_stats(its,ite,kds,kde,jts,jte1,ims,ime,kms,kme,jms,jme,v,'wind V in')
1064 if(fire_print_msg.gt.0)then
1065 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1066 write(msg,'(a,f7.2,a)')'interpolate_atm2fire: log-interpolation of wind to',fire_wind_height,' m'
1067 call message(msg)
1068 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1069 endif
1071 ! indexing
1073 ! file for u
1074 itst=ifval(ids.eq.its,its,its-1)
1075 itet=ifval(ide.eq.ite,ite,ite+1)
1076 jtst=ifval(jds.ge.jts,jts,jts-1)
1077 jtet=ifval(jde.eq.jte,jte,jte+1)
1079 if(fire_print_msg.ge.1)then
1080 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1081 write(msg,7001)'atm input ','tile',its,ite,jts,jte
1082 call message(msg)
1083 write(msg,7001)'altw ','tile',itst,itet,jtst,jtet
1084 call message(msg)
1085 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1086 endif
1087 7001 format(a,' dimensions ',a4,':',i6,' to ',i6,' by ',i6,' to ',i6)
1089 kdmax=kde-1 ! max layer to interpolate from, can be less
1091 do j = jtst,jtet
1092 do k=kds,kdmax+1
1093 do i = itst,itet
1094 altw(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g ! altitude of the bottom w-point
1095 enddo
1096 enddo
1097 enddo
1099 ! values at u points
1100 itsu=ifval(ids.eq.its,its+1,its) ! staggered direction
1101 iteu=ifval(ide.eq.ite,ite,ite+1) ! staggered direction
1102 jtsu=ifval(jds.ge.jts,jts,jts-1)
1103 jteu=ifval(jde.eq.jte,jte,jte+1)
1105 if(fire_print_msg.ge.1)then
1106 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1107 write(msg,7001)'u interp at','tile',itsu,iteu,jtsu,jteu
1108 call message(msg)
1109 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1110 endif
1112 do j = jtsu,jteu
1113 do k=kds,kdmax+1
1114 do i = itsu,iteu
1115 altub(i,k,j)= 0.5*(altw(i-1,k,j)+altw(i,k,j)) ! altitude of the bottom point under u-point
1116 enddo
1117 enddo
1118 do k=kds,kdmax
1119 do i = itsu,iteu
1120 hgtu(i,k,j) = 0.5*(altub(i,k,j)+altub(i,k+1,j)) - altub(i,kds,j) ! height of the u-point above the ground
1121 enddo
1122 enddo
1123 enddo
1125 ! values at v points
1126 jtsv=ifval(jds.eq.jts,jts+1,jts) ! staggered direction
1127 jtev=ifval(jde.eq.jte,jte,jte+1) ! staggered direction
1128 itsv=ifval(ids.ge.its,its,its-1)
1129 itev=ifval(ide.eq.ite,ite,ite+1)
1131 if(fire_print_msg.ge.1)then
1132 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1133 write(msg,7001)'v interp at','tile',itsv,itev,jtsv,jtev
1134 call message(msg)
1135 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1136 endif
1137 do j = jtsv,jtev
1138 do k=kds,kdmax+1
1139 do i = itsv,itev
1140 altvb(i,k,j)= 0.5*(altw(i,k,j-1)+altw(i,k,j)) ! altitude of the bottom point under v-point
1141 enddo
1142 enddo
1143 do k=kds,kdmax
1144 do i = itsv,itev
1145 hgtv(i,k,j) = 0.5*(altvb(i,k,j)+altvb(i,k+1,j)) - altvb(i,kds,j) ! height of the v-point above the ground
1146 enddo
1147 enddo
1148 enddo
1150 #ifdef DEBUG_OUT
1151 call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,altub,'altub',id)
1152 call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,altvb,'altvb',id)
1153 call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,hgtu,'hgtu',id)
1154 call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,hgtv,'hgtv',id)
1155 #endif
1156 logfwh = log(fire_wind_height)
1158 ! interpolate u, staggered in X
1160 do j = jtsu,jteu ! compute on domain by 1 smaller, otherwise z0 and ph not available
1161 do i = itsu,iteu ! compute with halo 2
1162 zr = 0.5*(z0(i,j)+z0(i-1,j)) ! interpolated roughness length under this u point
1163 if(fire_wind_height > zr)then !
1164 do k=kds,kdmax
1165 ht = hgtu(i,k,j) ! height of this u point above the ground
1166 if( .not. ht < fire_wind_height) then ! found layer k this point is in
1167 loght = log(ht)
1168 if(k.eq.kds)then ! first layer, log linear interpolation from 0 at zr
1169 logz0 = log(zr)
1170 ua(i,j)= u(i,k,j)*(logfwh-logz0)/(loght-logz0)
1171 else ! log linear interpolation
1172 loglast=log(hgtu(i,k-1,j))
1173 ua(i,j)= u(i,k-1,j) + (u(i,k,j) - u(i,k-1,j)) * ( logfwh - loglast) / (loght - loglast)
1174 endif
1175 goto 10
1176 endif
1177 if(k.eq.kdmax)then ! last layer, still not high enough
1178 ua(i,j)=u(i,k,j)
1179 endif
1180 enddo
1181 10 continue
1182 else ! roughness higher than the fire wind height
1183 ua(i,j)=0.
1184 endif
1185 enddo
1186 enddo
1188 ! interpolate v, staggered in Y
1190 do j = jtsv,jtev
1191 do i = itsv,itev
1192 zr = 0.5*(z0(i,j-1)+z0(i,j)) ! roughness length under this v point
1193 if(fire_wind_height > zr)then !
1194 do k=kds,kdmax
1195 ht = hgtv(i,k,j) ! height of this u point above the ground
1196 if( .not. ht < fire_wind_height) then ! found layer k this point is in
1197 loght = log(ht)
1198 if(k.eq.kds)then ! first layer, log linear interpolation from 0 at zr
1199 logz0 = log(zr)
1200 va(i,j)= v(i,k,j)*(logfwh-logz0)/(loght-logz0)
1201 else ! log linear interpolation
1202 loglast=log(hgtv(i,k-1,j))
1203 va(i,j)= v(i,k-1,j) + (v(i,k,j) - v(i,k-1,j)) * ( logfwh - loglast) / (loght - loglast)
1204 endif
1205 goto 11
1206 endif
1207 if(k.eq.kdmax)then ! last layer, still not high enough
1208 va(i,j)=v(i,k,j)
1209 endif
1210 enddo
1211 11 continue
1212 else ! roughness higher than the fire wind height
1213 va(i,j)=0.
1214 endif
1215 enddo
1216 enddo
1218 #ifdef DEBUG_OUT
1219 ! store the output for diagnostics
1220 do j = jts,jte1
1221 do i = its,ite1
1222 uah(i,j)=ua(i,j)
1223 vah(i,j)=va(i,j)
1224 enddo
1225 enddo
1227 call write_array_m(its,ite1,jts,jte,ims,ime,jms,jme,uah,'uah_n',id) ! no reflection
1228 call write_array_m(its,ite,jts,jte1,ims,ime,jms,jme,vah,'vah_n',id)
1229 #endif
1231 ips1 = ifval(ips.eq.ids,ips+1,ips)
1232 call continue_at_boundary(1,1,0., & ! x direction
1233 TDIMS, &! memory dims atm grid tile
1234 ids+1,ide,jds,jde, & ! domain dims - where u defined
1235 ips1,ipe,jps,jpe, & ! patch dims
1236 itsu,iteu,jtsu,jteu, & ! tile dims - in nonextended direction one beyond if at patch boundary but not domain
1237 itsou,iteou,jtsou,jteou, & ! out, where set
1238 ua) ! array
1240 jps1 = ifval(jps.eq.jds,jps+1,jps)
1241 call continue_at_boundary(1,1,0., & ! y direction
1242 TDIMS, & ! memory dims atm grid tile
1243 ids,ide,jds+1,jde, & ! domain dims - where v wind defined
1244 ips,ipe,jps1,jpe, & ! patch dims
1245 itsv,itev,jtsv,jtev, & ! tile dims
1246 itsov,iteov,jtsov,jteov, & ! where set
1247 va) ! array
1249 ! store the output for diagnostics
1250 do j = jts,jte1
1251 do i = its,ite1
1252 uah(i,j)=ua(i,j)
1253 vah(i,j)=va(i,j)
1254 enddo
1255 enddo
1257 #ifdef DEBUG_OUT
1258 call write_array_m(itsou,iteou,jtsou,jteou,TDIMS,ua,'ua',id)
1259 call write_array_m(itsov,iteov,jtsov,jteov,TDIMS,va,'va',id)
1260 #endif
1262 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1263 ! don't have all values valid, don't check
1264 write(msg,12)'atm mesh wind U at',fire_wind_height,' m'
1265 call print_2d_stats(itsou,iteou,jtsou,jteou,TDIMS,ua,msg)
1266 write(msg,12)'atm mesh wind V at',fire_wind_height,' m'
1267 call print_2d_stats(itsov,iteov,jtsov,jteov,TDIMS,va,msg)
1268 12 format(a,f6.2,a)
1269 call print_2d_stats(its,ite1,jts,jte,ims,ime,jms,jme,uah,'UAH')
1270 call print_2d_stats(its,ite,jts,jte1,ims,ime,jms,jme,vah,'VAH')
1271 !call write_array_m(its,ite1,jts,jte,ims,ime,jms,jme,uah,'uah',id)
1272 !call write_array_m(its,ite,jts,jte1,ims,ime,jms,jme,vah,'vah',id)
1273 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1275 ! ---------------
1276 ! | F | F | F | F | Example of atmospheric and fire grid with
1277 ! |-------|-------| ir=jr=4.
1278 ! | F | F | F | F | Winds are given at the midpoints of the sides of the atmosphere grid,
1279 ! ua------z-------| interpolated to midpoints of the cells of the fine fire grid F.
1280 ! | F | F | F | F | This is (1,1) cell of atmosphere grid, and [*] is the (1,1) cell of the fire grid.
1281 ! |---------------| ua(1,1) <--> uf(0.5,2.5)
1282 ! | * | F | F | F | va(1,1) <--> vf(2.5,0.5)
1283 ! -------va------ za(1,1) <--> zf(2.5,2.5)
1284 !
1285 ! ^ x2
1286 ! | --------va(1,2)---------
1287 ! | | | | Example of atmospheric and fire grid with
1288 ! | | | | ir=jr=1.
1289 ! | | za,zf | Winds are given at the midpoints of the sides of the atmosphere grid,
1290 ! | ua(1,1)----uf,vf-----ua(2,1) interpolated to midpoints of the cells of the (the same) fire grid
1291 ! | | (1,1) | ua(1,1) <--> uf(0.5,1)
1292 ! | | | | va(1,1) <--> vf(1,0.5)
1293 ! | | | | za(1,1) <--> zf(1,1)
1294 ! | --------va(1,1)---------
1295 ! |--------------------> x1
1296 !
1297 ! Meshes are aligned by the lower left cell of the domain. Then in the above figure
1298 ! u = node with the ua component of the wind at (ids,jds), midpoint of side
1299 ! v = node with the va component of the wind at (ids,jds), midpoint of side
1300 ! * = fire grid node at (ifds,jfds)
1301 ! z = node with height, midpoint of cell
1302 !
1303 ! ua(ids,jds)=uf(ifds-0.5,jfds+jr*0.5-0.5) = uf(ifds-0.5,jfds+(jr-1)*0.5)
1304 ! va(ids,jds)=vf(ifds+ir*0.5-0.5,jfds-0.5) = vf(ifds+(ir-1)*0.5,jfds-0.5)
1305 ! za(ids,jds)=zf(ifds+ir*0.5-0.5,jfds+jr*0.5-0.5) = zf(ifds+(ir-1)*0.5,jfds+(jr-1)*0.5)
1307 ! ifts1=max(snode(ifts,ifps,-1),ifds) ! go 1 beyond patch boundary but not at domain boundary
1308 ! ifte1=min(snode(ifte,ifpe,+1),ifde)
1309 ! jfts1=max(snode(jfts,jfps,-1),jfds)
1310 ! jfte1=min(snode(jfte,jfpe,+1),jfde)
1312 call interpolate_2d( &
1313 TDIMS, & ! memory dims atm grid tile
1314 itsou,iteou,jtsou,jteou,& ! where set
1315 ifms,ifme,jfms,jfme, & ! array dims fire grid
1316 ifts,ifte,jfts,jfte,& ! dimensions on the fire grid to interpolate to
1317 ir,jr, & ! refinement ratio
1318 real(ids),real(jds),ifds-0.5,jfds+(jr-1)*0.5, & ! line up by lower left corner of domain
1319 ua, & ! in atm grid
1320 uf) ! out fire grid
1322 call interpolate_2d( &
1323 TDIMS, & ! memory dims atm grid tile
1324 itsov,iteov,jtsov,jteov,& ! where set
1325 ifms,ifme,jfms,jfme, & ! array dims fire grid
1326 ifts,ifte,jfts,jfte,& ! dimensions on the fire grid to interpolate to
1327 ir,jr, & ! refinement ratio
1328 real(ids),real(jds),ifds+(ir-1)*0.5,jfds-0.5, & ! line up by lower left corner of domain
1329 va, & ! in atm grid
1330 vf) ! out fire grid
1332 call print_2d_stats_vec(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,vf,'fire wind (m/s)')
1333 ! call print_2d_stats_vec(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,vf,'fire wind extended')
1334 #ifdef DEBUG_OUT
1335 call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,'uf1',id)
1336 call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,vf,'vf1',id)
1337 ! call write_array_m(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,'uf1',id)
1338 ! call write_array_m(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,vf,'vf1',id)
1339 #endif
1340 return
1342 end subroutine interpolate_atm2fire
1344 !
1345 !***
1346 !
1348 subroutine check_fmesh(ids,ide,ifds,ifde,ir,s)
1349 !*** purpose: check if fire and atm meshes line up
1350 implicit none
1351 !*** arguments
1352 integer, intent(in)::ids,ide,ifds,ifde,ir
1353 character(len=*),intent(in)::s
1354 !*** local
1355 character(len=128)msg
1356 !*** executable
1357 if ((ide-ids+1)*ir.ne.(ifde-ifds+1))then
1358 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1359 write(msg,1)s,ids,ide,ifds,ifde,ir
1360 1 format('module_fr_sfire_driver: incompatible bounds ',a,' atm ',i5,':',i5,' fire ',i5,':',i5,' ratio ',i3)
1361 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1362 call crash(msg)
1363 endif
1364 end subroutine check_fmesh
1366 !
1367 !*****************************
1368 !
1369 subroutine set_flags(config_flags)
1370 implicit none
1371 TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags
1372 ! copy flags from wrf to module_fr_sfire_util
1373 ! for instructions how to add a flag see the top of module_fr_sfire_util.F
1376 fire_print_msg = config_flags%fire_print_msg
1377 fire_print_file = config_flags%fire_print_file
1378 fuel_left_method = config_flags%fire_fuel_left_method
1379 fuel_left_irl = config_flags%fire_fuel_left_irl
1380 fuel_left_jrl = config_flags%fire_fuel_left_jrl
1381 fire_const_time = config_flags%fire_const_time
1382 fire_const_grnhfx = config_flags%fire_const_grnhfx
1383 fire_const_grnqfx = config_flags%fire_const_grnqfx
1384 fire_atm_feedback = config_flags%fire_atm_feedback
1385 boundary_guard = config_flags%fire_boundary_guard
1386 fire_back_weight = config_flags%fire_back_weight
1387 fire_grows_only = config_flags%fire_grows_only
1388 fire_upwinding = config_flags%fire_upwinding
1389 fire_upwind_split = config_flags%fire_upwind_split
1390 fire_viscosity = config_flags%fire_viscosity
1391 fire_lfn_ext_up = config_flags%fire_lfn_ext_up
1392 fire_test_steps = config_flags%fire_test_steps
1393 !fire_topo_from_atm = config_flags%fire_topo_from_atm
1394 fire_advection = config_flags%fire_advection
1399 end subroutine set_flags
1401 end module module_fr_sfire_driver