Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_fr_sfire_driver.F @ 1419

Last change on this file since 1419 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 62.0 KB

Line
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	!
15
16	module module_fr_sfire_driver
17	! use this module for standalone call, you only need to provide some mock-up wrf modules
18
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
23
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
29
30	implicit none
31
32
33	private
34	public sfire_driver_em,use_atm_vars
35
36	logical:: use_atm_vars=.true. ! interpolate wind from atm mesh and average output fluxes back
37
38	contains
39
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 )
49
50	!*** purpose: driver from grid structure
51
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
59
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
72
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
84
85
86	!*** executable
87
88	! populate our structures from wrf
89
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
102
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)
106
107	! copy configuration flags to sfire internal structures
108	call set_flags(config_flags)
109
110	! refinement r
111	ir=grid%sr_x ! refinement ratio
112	jr=grid%sr_y
113	itimestep=grid%itimestep
114	restart=config_flags%restart
115
116
117
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)
123
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
126
127	need_lfn_update=.false.
128	do fire_ifun=fire_ifun_start,fire_ifun_end
129
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
136
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
142
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"
152
153	elseif(fire_ifun.eq.2)then
154	! halo exchange on zsf width 2
155	#include "HALO_FIRE_ZSF.inc"
156
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"
161
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"
165
166	elseif(fire_ifun.eq.6)then
167	! computing fuel_left needs ignition time from neighbors
168	#include "HALO_FIRE_TIGN.inc"
169
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	)
212
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
219
220
221
222	enddo
223	enddo
224	if(tsteps>0)call crash('sfire_driver_em: test run of uncoupled fire model completed')
225
226	end subroutine sfire_driver_em
227
228	!
229	!*******************
230	!
231
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
268
269	implicit none
270
271	!*** arguments
272
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
289
290	logical, intent(in)::restart
291
292
293	logical, intent(out)::need_lfn_update
294
295	integer,dimension(num_tiles),intent(in) :: i_start,i_end,j_start,j_end ! atm grid tiling
296
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
306
307
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
311
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
319
320	real, dimension(ims:ime, jms:jme), intent(inout)::xlong, xlat ! inout because of extension at bdry
321
322	real, intent(inout), dimension(ifms:ifme,jfms:jfme):: &
323	nfuel_cat ! fuel data; can be also set inside (cell based, fire grid)
324
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
328
329	real, intent(out), dimension(ifms:ifme, jfms:jfme):: &
330	fire_area ! fraction of each cell burning
331
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)
338
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)
345
346	! *** data (constant in time) ***
347
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
350
351	type(fire_params),intent(inout)::fp
352
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
362
363
364	!*** executable
365
366	! time - assume dt does not change
367	time_start = itimestep * dt
368
369	! fire mesh step
370	dxf=dx/ir
371	dyf=dy/jr
372
373
374
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)
393
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
403
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
410
411
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
416
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
422
423	! fake atm tile bounds
424	kts=kps
425	kte=kpe
426
427	! staggered atm patch bounds
428	ipe1=ifval(ipe.eq.ide,ipe+1,ipe)
429	jpe1=ifval(jpe.eq.jde,jpe+1,jpe)
430
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
435
436	id = ifval(pid.ne.0,pid+ij*10000,0) ! for print
437
438	ignitions_done_tile(ij)=0
439
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
449
450	! staggered atm tile bounds
451	ite1=ifval(ite.eq.ide,ite+1,ite)
452	jte1=ifval(jte.eq.jde,jte+1,jte)
453
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)
462
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
467
468
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)
475
476	if(ifun.eq.1)then ! set terrain
477
478	if(restart)then
479
480	call message('restart - topo initialization skipped')
481
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
504
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
521
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
537
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
549
550	endif
551
552	endif
553
554	elseif(ifun.eq.2)then
555
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')
558
559	elseif(ifun.eq.3)then ! interpolate winds to the fire grid
560
561	if(use_atm_vars)then
562
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
586
587	endif
588
589	! the following executes in any case
590
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	)
611
612	if(ifun.eq.6)then ! heat fluxes into the atmosphere
613
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')
616
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)
641
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_feedbackgrnhfx(i,j)s
652	grnqfx(i,j)=fire_atm_feedbackgrnqfx(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
658
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
662
663	endif ! ifun=6
664
665	enddo ! tiles
666	!$OMP END PARALLEL DO
667
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
676
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
697
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
715
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
724
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
740
741	end subroutine sfire_driver_phys
742	!
743	!*******************
744	!
745
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')
770
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
841
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
852
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_ctrpi2/360.)unit_fxlat ! latitude
866	endif
867
868	end subroutine fire_ignition_convert
869
870	!
871	!*****************************
872	!
873
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
885
886	implicit none
887	!*** purpose: interpolate height
888
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
902
903
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
907
908	! terrain height
909
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
920
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
928
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)
934
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
944
945	end subroutine interpolate_z2fire
946	!
947	!*****************************
948	!
949
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
967
968	implicit none
969	!*** purpose: interpolate winds and height
970
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
996
997
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)
1011
1012	!*** executable
1013
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
1022
1023
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
1030
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
1057
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')
1063
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
1070
1071	! indexing
1072
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)
1078
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)
1088
1089	kdmax=kde-1 ! max layer to interpolate from, can be less
1090
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
1098
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)
1104
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
1111
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
1124
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)
1130
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
1149
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)
1157
1158	! interpolate u, staggered in X
1159
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
1187
1188	! interpolate v, staggered in Y
1189
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
1217
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
1226
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
1230
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
1239
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
1248
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
1256
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
1261
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)
1274
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+jr0.5-0.5) = uf(ifds-0.5,jfds+(jr-1)0.5)
1304	! va(ids,jds)=vf(ifds+ir0.5-0.5,jfds-0.5) = vf(ifds+(ir-1)0.5,jfds-0.5)
1305	! za(ids,jds)=zf(ifds+ir0.5-0.5,jfds+jr0.5-0.5) = zf(ifds+(ir-1)0.5,jfds+(jr-1)0.5)
1306
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)
1311
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
1321
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
1331
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
1341
1342	end subroutine interpolate_atm2fire
1343
1344	!
1345	!***
1346	!
1347
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
1365
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
1374
1375
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
1395
1396
1397
1398
1399	end subroutine set_flags
1400
1401	end module module_fr_sfire_driver

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format