| blob | f92b041a03e9d57d870f6606b337a2372d51dbc9 |
1 !WRF:MEDIATION_LAYER:FIRE_MODEL
3 !*** Jan Mandel August 2007 - February 2010
4 !*** email: Jan.Mandel@gmail.com
6 ! Routines dealing with the atmosphere
8 module module_fr_sfire_atm
10 use module_model_constants, only: cp,xlv,g
11 use module_fr_sfire_phys, only: fire_wind_height,have_fuel_cats,fcz0,fcwh,nfuelcats,no_fuel_cat,no_fuel_cat2,windrf
12 use module_fr_sfire_util
13 USE module_dm , only : wrf_dm_sum_reals
16 use module_fr_sfire_phys, only: mfuelcats, nfuelcats ! for emissions
17 USE module_state_description, only: num_tracer
18 use module_fr_sfire_phys, only: fuel_name
20 #ifdef WRF_CHEM
21 USE module_state_description, only: num_chem
22 USE module_configure, only: &
23 p_co, &
24 p_ch4, &
25 p_h2, &
26 p_no, &
27 p_no2, &
28 p_so2, &
29 p_nh3, &
30 p_p25, &
31 p_p25i, &
32 p_p25j, &
33 p_oc1, &
34 p_oc2, &
35 p_bc1, &
36 p_bc2, &
37 p_ald, &
38 p_csl, &
39 p_eth, &
40 p_hc3, &
41 p_hc5, &
42 p_hcho, &
43 p_iso, &
44 p_ket, &
45 p_mgly, &
46 p_ol2, &
47 p_olt, &
48 p_oli, &
49 p_ora2,&
50 p_tol, &
51 p_xyl, &
52 p_bigalk, &
53 p_bigene, &
54 p_c10h16, &
55 p_c2h4, &
56 p_c2h5oh, &
57 p_c2h6, &
58 p_c3h6, &
59 p_c3h8, &
60 p_ch3cooh, &
61 p_ch3oh, &
62 p_cres, &
63 p_glyald, &
64 ! p_hyac, &
65 p_isopr, &
66 p_macr, &
67 p_mek, &
68 p_mvk, &
69 p_smoke, & ! tracer smoke exists only with CHEM
70 p_sulf, &
71 p_dms, &
72 p_msa, &
73 p_dust_1, &
74 p_dust_2, &
75 p_dust_3, &
76 p_dust_4, &
77 p_dust_5, &
78 p_seas_1, &
79 p_seas_2, &
80 p_seas_3, &
81 p_seas_4, &
82 p_p10
83 #endif
85 USE module_state_description, only: &
86 p_tr17_1, &
87 p_tr17_2, &
88 p_tr17_3, &
89 p_tr17_4, &
90 p_tr17_5, &
91 p_tr17_6, &
92 p_tr17_7, &
93 p_tr17_8
95 implicit none
97 #ifndef WRF_CHEM
98 integer, parameter, private:: num_chem=0
99 #endif
101 logical, save :: have_wind_log_interpolation = .false. ! status
103 ! emission tables
104 REAL, dimension(mfuelcats), save:: &
105 co=0., &
106 ch4=0., &
107 h2=0., &
108 no=0., &
109 no2=0., &
110 so2=0., &
111 nh3=0., &
112 oc1=0., &
113 oc2=0., &
114 bc1=0., &
115 bc2=0., &
116 ald=0., &
117 csl=0., &
118 eth=0., &
119 hc3=0., &
120 p25=0., &
121 p25i=0., &
122 p25j=0., &
123 hc5=0., &
124 hcho=0., &
125 iso=0., &
126 ket=0., &
127 mgly=0., &
128 ol2=0., &
129 olt=0., &
130 oli=0., &
131 ora2=0.,&
132 tol=0., &
133 xyl=0., &
134 bigalk=0., &
135 bigene=0., &
136 c10h16=0., &
137 c2h4=0., &
138 c2h5oh=0., &
139 c2h6=0., &
140 c3h6=0., &
141 c3h8=0., &
142 ch3cooh=0., &
143 ch3oh=0., &
144 cres=0., &
145 glyald=0., &
146 ! hyac=0., &
147 isopr=0., &
148 macr=0., &
149 mek=0., &
150 mvk=0., &
151 smoke=0., &
152 sulf=0., &
153 dms=0., &
154 msa=0., &
155 dust_1=0., &
156 dust_2=0., &
157 dust_3=0., &
158 dust_4=0., &
159 dust_5=0., &
160 seas_1=0., &
161 seas_2=0., &
162 seas_3=0., &
163 seas_4=0., &
164 p10=0., &
165 tr17_1=0., &
166 tr17_2=0., &
167 tr17_3=0., &
168 tr17_4=0., &
169 tr17_5=0., &
170 tr17_6=0., &
171 tr17_7=0., &
172 tr17_8=0.
174 real, parameter:: & ! reciprocal molecular weights (mol/g)
175 imw_co = 1./28.010,&
176 imw_ch4 = 1./16.04,&
177 imw_h2 = 1./2.016,&
178 imw_no = 1./30.006,&
179 imw_no2 = 1./46.006,&
180 imw_so2 = 1./64.066,&
181 imw_nh3 = 1./17.031
183 real, parameter:: mw_air=28.97 ! molecular weight of air g/mol
185 ! should be declared in the registry and stored in the state in future because of restarts
187 real, pointer, save, dimension(:)::c_chem
188 real, pointer, save, dimension(:)::c_fuel
189 real, pointer, save, dimension(:)::c_tracer
190 logical, save:: emis_read = .false.
191 integer, save:: msglevel =1,printsums=0 ! when to print sums
192 integer, parameter:: line=5 ! number of species, emissions, etc. per line
194 contains
196 ! chem arrays are chem tracer
197 ! indices p_species are generated in inc/scalar_indices.inc and included in frame/module_configure.F
199 subroutine read_emissions_table(chem_opt,tracer_opt)
200 implicit none
201 integer, intent(in)::chem_opt,tracer_opt
202 logical, external:: wrf_dm_on_monitor
203 external::wrf_dm_bcast_integer , wrf_dm_bcast_real
204 integer, dimension(10)::compatible_chem_opt
205 integer:: iounit,ierr,i
206 character(len=128)::msg
207 namelist/emissions/ compatible_chem_opt, printsums, &
208 co, &
209 ch4, &
210 h2, &
211 no, &
212 no2, &
213 so2, &
214 nh3, &
215 p25, &
216 p25i, &
217 p25j, &
218 oc1, &
219 oc2, &
220 bc1, &
221 bc2, &
222 ald, &
223 csl, &
224 eth, &
225 hc3, &
226 hc5, &
227 hcho, &
228 iso, &
229 ket, &
230 mgly, &
231 ol2, &
232 olt, &
233 oli, &
234 ora2,&
235 tol, &
236 xyl, &
237 bigalk, &
238 bigene, &
239 c10h16, &
240 c2h4, &
241 c2h5oh, &
242 c2h6, &
243 c3h6, &
244 c3h8, &
245 ch3cooh, &
246 ch3oh, &
247 cres, &
248 glyald, &
249 ! hyac, &
250 isopr, &
251 macr, &
252 mek, &
253 mvk, &
254 smoke, &
255 sulf, &
256 dms, &
257 msa, &
258 dust_1, &
259 dust_2, &
260 dust_3, &
261 dust_4, &
262 dust_5, &
263 seas_1, &
264 seas_2, &
265 seas_3, &
266 seas_4, &
267 p10, &
268 tr17_1, &
269 tr17_2, &
270 tr17_3, &
271 tr17_4, &
272 tr17_5, &
273 tr17_6, &
274 tr17_7, &
275 tr17_8
277 !$ if (OMP_GET_THREAD_NUM() .ne. 0)then
278 !$ call crash('read_emissions_table: must be called from master thread')
279 !$ endif
281 IF ( wrf_dm_on_monitor() ) THEN
282 ! we are the master task
284 iounit=open_text_file('namelist.fire_emissions','read')
285 compatible_chem_opt=0
286 read(iounit,emissions,iostat=ierr)
287 if(ierr.ne.0)call crash('read_emissions_table: error reading namelist emissions in file namelist.fire_emissions')
288 CLOSE(iounit)
289 write(msg,'(a,i3,a)')'reading emissions table for',nfuelcats,' fuel categories'
290 call message(msg,level=0)
291 if (.not.any(compatible_chem_opt.eq.chem_opt))then
292 write(msg,'(a,i4,a)')'read_emissions_table: chem_opt=',chem_opt,' not between given compatible_chem_opt in namelist.fire_emissions'
293 call message(msg,level=0)
294 write(msg,'(a,10i4)')'compatible_chem_opt=', compatible_chem_opt
295 call message(msg,level=0)
296 call crash('chem_opt in namelist.input is not consistent with namelist.fire_emissions')
297 endif
298 ENDIF
299 call wrf_dm_bcast_integer(printsums, 1)
300 call wrf_dm_bcast_real(co, nfuelcats)
301 call wrf_dm_bcast_real(ch4, nfuelcats)
302 call wrf_dm_bcast_real(h2, nfuelcats)
303 call wrf_dm_bcast_real(no, nfuelcats)
304 call wrf_dm_bcast_real(no2, nfuelcats)
305 call wrf_dm_bcast_real(so2, nfuelcats)
306 call wrf_dm_bcast_real(nh3, nfuelcats)
307 call wrf_dm_bcast_real(p25, nfuelcats)
308 call wrf_dm_bcast_real(p25i, nfuelcats)
309 call wrf_dm_bcast_real(p25j, nfuelcats)
310 call wrf_dm_bcast_real(oc1, nfuelcats)
311 call wrf_dm_bcast_real(oc2, nfuelcats)
312 call wrf_dm_bcast_real(bc1, nfuelcats)
313 call wrf_dm_bcast_real(bc2, nfuelcats)
314 call wrf_dm_bcast_real(ald, nfuelcats)
315 call wrf_dm_bcast_real(csl, nfuelcats)
316 call wrf_dm_bcast_real(eth, nfuelcats)
317 call wrf_dm_bcast_real(hc3, nfuelcats)
318 call wrf_dm_bcast_real(hc5, nfuelcats)
319 call wrf_dm_bcast_real(hcho, nfuelcats)
320 call wrf_dm_bcast_real(iso, nfuelcats)
321 call wrf_dm_bcast_real(ket, nfuelcats)
322 call wrf_dm_bcast_real(mgly, nfuelcats)
323 call wrf_dm_bcast_real(ol2, nfuelcats)
324 call wrf_dm_bcast_real(olt, nfuelcats)
325 call wrf_dm_bcast_real(oli, nfuelcats)
326 call wrf_dm_bcast_real(ora2,nfuelcats)
327 call wrf_dm_bcast_real(tol, nfuelcats)
328 call wrf_dm_bcast_real(xyl, nfuelcats)
329 call wrf_dm_bcast_real(bigalk, nfuelcats)
330 call wrf_dm_bcast_real(bigene, nfuelcats)
331 call wrf_dm_bcast_real(c10h16, nfuelcats)
332 call wrf_dm_bcast_real(c2h4, nfuelcats)
333 call wrf_dm_bcast_real(c2h5oh, nfuelcats)
334 call wrf_dm_bcast_real(c2h6, nfuelcats)
335 call wrf_dm_bcast_real(c3h6, nfuelcats)
336 call wrf_dm_bcast_real(c3h8, nfuelcats)
337 call wrf_dm_bcast_real(ch3cooh, nfuelcats)
338 call wrf_dm_bcast_real(ch3oh, nfuelcats)
339 call wrf_dm_bcast_real(cres, nfuelcats)
340 call wrf_dm_bcast_real(glyald, nfuelcats)
341 ! call wrf_dm_bcast_real(hyac, nfuelcats)
342 call wrf_dm_bcast_real(isopr, nfuelcats)
343 call wrf_dm_bcast_real(macr, nfuelcats)
344 call wrf_dm_bcast_real(mek, nfuelcats)
345 call wrf_dm_bcast_real(mvk, nfuelcats)
346 call wrf_dm_bcast_real(smoke, nfuelcats)
347 call wrf_dm_bcast_real(sulf, nfuelcats)
348 call wrf_dm_bcast_real(dms, nfuelcats)
349 call wrf_dm_bcast_real(msa, nfuelcats)
350 call wrf_dm_bcast_real(dust_1, nfuelcats)
351 call wrf_dm_bcast_real(dust_2, nfuelcats)
352 call wrf_dm_bcast_real(dust_3, nfuelcats)
353 call wrf_dm_bcast_real(dust_4, nfuelcats)
354 call wrf_dm_bcast_real(dust_5, nfuelcats)
355 call wrf_dm_bcast_real(seas_1, nfuelcats)
356 call wrf_dm_bcast_real(seas_2, nfuelcats)
357 call wrf_dm_bcast_real(seas_3, nfuelcats)
358 call wrf_dm_bcast_real(seas_4, nfuelcats)
359 call wrf_dm_bcast_real(p10, nfuelcats)
360 call wrf_dm_bcast_real(tr17_1, nfuelcats)
361 call wrf_dm_bcast_real(tr17_2, nfuelcats)
362 call wrf_dm_bcast_real(tr17_3, nfuelcats)
363 call wrf_dm_bcast_real(tr17_4, nfuelcats)
364 call wrf_dm_bcast_real(tr17_5, nfuelcats)
365 call wrf_dm_bcast_real(tr17_6, nfuelcats)
366 call wrf_dm_bcast_real(tr17_7, nfuelcats)
367 call wrf_dm_bcast_real(tr17_8, nfuelcats)
369 if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
370 ! should be stored in the registry future because of restarts
371 write(msg,'(3(a,i3,1x))')'allocating c_chem size',num_chem,'c_tracer size',num_tracer,'c_fuel size',nfuelcats
372 call message(msg,level=2)
373 if(num_chem>0)then
374 allocate(c_chem(num_chem))
375 c_chem=0. ! cumulative burnt
376 endif
377 if(num_tracer>0)then
378 allocate(c_tracer(num_tracer))
379 c_tracer=0. ! cumulative burnt
380 endif
381 allocate(c_fuel(nfuelcats))
382 c_fuel=0. ! total per timestep, rate burnt, cumulative burnt
383 write(msg,'(a,i3,a,i3)')'allocated c_chem size',size(c_chem),' c_fuel size',size(c_fuel)
384 call message(msg,level=2)
385 endif
387 emis_read=.true.
389 end subroutine read_emissions_table
391 subroutine add_fire_emissions(chem_opt,tracer_opt,dt,dx,dy, &
392 ifms,ifme,jfms,jfme, &
393 ifts,ifte,jtfs,jfte, &
394 ids,ide,kds,kde,jds,jde, &
395 ims,ime,kms,kme,jms,jme, &
396 its,ite,kts,kte,jts,jte, &
397 rho,dz8w, & ! input on atmosphere mesh
398 fgip, fuel_frac_burnt, nfuel_cat, & ! input on fire mesh
399 chem,tracer) ! output
402 implicit none
404 !*** purpose
405 ! average fire emissions from fire mesh to coarser atmosphereic mesh and add to chemistry arrays
407 !*** arguments
408 ! the dimensions are in cells, not nodes!
410 ! input
411 integer, intent(in)::chem_opt,tracer_opt
412 real, intent(in):: dt,dx,dy ! time step & mesh spacing
413 integer, intent(in)::its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme &! atm grid dims
414 ,ids,ide,kds,kde,jds,jde
415 integer, intent(in)::ifts,ifte,jtfs,jfte,ifms,ifme,jfms,jfme ! fire grid dims
416 real, intent(in)::rho(ims:ime,kms:kme,jms:jme), & ! air density kg/m^3
417 dz8w(ims:ime,kms:kme,jms:jme) ! layer height
418 real, intent(in), dimension(ifms:ifme,jfms:jfme):: fgip, & ! initial fuel load kg/m^2
419 fuel_frac_burnt, & ! fuel fraction burned this step kg/kg
420 nfuel_cat ! fuel category (Anderson= 1 to 13)
421 ! update
422 real, intent(inout)::chem(ims:ime,kms:kme,jms:jme,num_chem),tracer(ims:ime,kms:kme,jms:jme,num_tracer)
424 !*** local
425 integer:: i,i_f,j,j_f,ir,jr,isz1,isz2,jsz1,jsz2,ioff,joff,ibase,jbase,cat,k1,areaw,m,k,k_p,errors
426 real::fuel_burnt,vol,air,conv, avgw,emis
427 character(len=128)msg
429 real, dimension(mfuelcats)::s_fuel,t_fuel,r_fuel ! total per timestep, rate burnt
430 #ifdef WRF_CHEM
431 real, dimension(num_chem) ::s_chem,t_chem,r_chem ! total per timestep, rate burnt
432 real, dimension(num_chem) ::a_chem,g_chem ! concentration in ground level 1
433 integer, parameter:: chem_np=59
434 integer:: chem_pointers(chem_np)
435 character(len=8)::chem_names(chem_np)
436 #endif
437 real, dimension(num_tracer) ::s_tracer,t_tracer,r_tracer ! total per timestep, rate burnt
439 integer, parameter:: tracer_np=8
440 integer:: tracer_pointers(tracer_np)
441 character(len=8)::tracer_names(tracer_np)
443 !*** executable
445 !check mesh dimensions and domain dimensions
446 call check_mesh_2dim(its,ite,jts,jte,ims,ime,jms,jme)
447 call check_mesh_2dim(ifts,ifte,jtfs,jfte,ifms,ifme,jfms,jfme)
449 if(.not.emis_read)call crash('add_fire_emissions: read_emissions_table must be called first')
450 write(msg,'(a,i3,a,i3,a,i3)')'add_fire_emissions: chem_opt=',chem_opt,' species ',num_chem,' tracers',num_tracer
451 call message(msg)
453 #ifdef WRF_CHEM
454 chem_pointers= (/ &
455 p_co, &
456 p_ch4, &
457 p_h2, &
458 p_no, &
459 p_no2, &
460 p_so2, &
461 p_nh3, &
462 p_p25, &
463 p_p25i, &
464 p_p25j, &
465 p_oc1, &
466 p_oc2, &
467 p_bc1, &
468 p_bc2, &
469 p_ald, &
470 p_csl, &
471 p_eth, &
472 p_hc3, &
473 p_hc5, &
474 p_hcho, &
475 p_iso, &
476 p_ket, &
477 p_mgly, &
478 p_ol2, &
479 p_olt, &
480 p_oli, &
481 p_ora2,&
482 p_tol, &
483 p_xyl, &
484 p_bigalk, &
485 p_bigene, &
486 p_c10h16, &
487 p_c2h4, &
488 p_c2h5oh, &
489 p_c2h6, &
490 p_c3h6, &
491 p_c3h8, &
492 p_ch3cooh, &
493 p_ch3oh, &
494 p_cres, &
495 p_glyald, &
496 ! p_hyac, &
497 p_isopr, &
498 p_macr, &
499 p_mek, &
500 p_mvk, &
501 p_smoke, &
502 p_sulf, &
503 p_dms, &
504 p_msa, &
505 p_dust_1, &
506 p_dust_2, &
507 p_dust_3, &
508 p_dust_4, &
509 p_dust_5, &
510 p_seas_1, &
511 p_seas_2, &
512 p_seas_3, &
513 p_seas_4, &
514 p_p10 /)
516 chem_names= (/ &
517 'co ', &
518 'ch4 ', &
519 'h2 ', &
520 'no ', &
521 'no2 ', &
522 'so2 ', &
523 'nh3 ', &
524 'p25 ', &
525 'p25i ', &
526 'p25j ', &
527 'oc1 ', &
528 'oc2 ', &
529 'bc1 ', &
530 'bc2 ', &
531 'ald ', &
532 'csl ', &
533 'eth ', &
534 'hc3 ', &
535 'hc5 ', &
536 'hcho ', &
537 'iso ', &
538 'ket ', &
539 'mgly ', &
540 'ol2 ', &
541 'olt ', &
542 'oli ', &
543 'ora2 ',&
544 'tol ', &
545 'xyl ', &
546 'bigalk ', &
547 'bigene ', &
548 'c10h16 ', &
549 'c2h4 ', &
550 'c2h5oh ', &
551 'c2h6 ', &
552 'c3h6 ', &
553 'c3h8 ', &
554 'ch3cooh ', &
555 'ch3oh ', &
556 'cres ', &
557 'glyald ', &
558 ! 'hyac ', &
559 'isopr ', &
560 'macr ', &
561 'mek ', &
562 'mvk ', &
563 'smoke ', &
564 'sulf ', &
565 'dms ', &
566 'msa ', &
567 'dust_1 ', &
568 'dust_2 ', &
569 'dust_3 ', &
570 'dust_4 ', &
571 'dust_5 ', &
572 'seas_1 ', &
573 'seas_2 ', &
574 'seas_3 ', &
575 'seas_4 ', &
576 'p10 ' /)
579 call check_pointers('chem',chem,chem_names,chem_pointers)
580 #endif
582 tracer_pointers= (/ &
583 p_tr17_1, &
584 p_tr17_2, &
585 p_tr17_3, &
586 p_tr17_4, &
587 p_tr17_5, &
588 p_tr17_6, &
589 p_tr17_7, &
590 p_tr17_8 /)
592 tracer_names= (/ &
593 'tr17_1 ', &
594 'tr17_2 ', &
595 'tr17_3 ', &
596 'tr17_4 ', &
597 'tr17_5 ', &
598 'tr17_6 ', &
599 'tr17_7 ', &
600 'tr17_8 ' /)
602 call check_pointers('tracer',tracer,tracer_names,tracer_pointers)
604 ! compute mesh sizes
605 isz1 = ite-its+1
606 jsz1 = jte-jts+1
607 isz2 = ifte-ifts+1
608 jsz2 = jfte-jtfs+1
611 ! check mesh sizes
612 if(isz1.le.0.or.jsz1.le.0.or.isz2.le.0.or.jsz2.le.0)then
613 call message('all mesh sizes must be positive',level=0)
614 goto 9
615 endif
617 ! compute mesh ratios
618 ir=isz2/isz1
619 jr=jsz2/jsz1
621 if(isz2.ne.isz1*ir .or. jsz2.ne.jsz1*jr)then
622 call message('input mesh size must be multiple of output mesh size',level=0)
623 goto 9
624 endif
626 avgw = 1.0/(ir*jr) ! averaging weight = 1/number of fire cells per atm cell
628 ! initialize emissions statistics per timestep
629 #ifdef WRF_CHEM
630 t_chem = 0.
631 #endif
632 t_fuel = 0.
633 do i=1,num_tracer
634 t_tracer(i)=0.
635 enddo
637 ! conversion fuel_burnt kg/m^2 -> chem_X 1e6*mol/mol
638 ! emis_X(g/kg)*fuel_burnt(kg/m^2)/mw_X(g/mol) = emissions (mol/m^2)
639 ! rho(kg/m^3)*dz8w(m))/mw_air(28.97e-3 kg/mol) = dry air in the 1st layer (mol/m^2)
641 k1 = kts
642 !$OMP CRITICAL(SFIRE_ATM_CRIT)
643 write(msg,'(a,i3)')'Fire tracers into atmosphere level',k1
644 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
645 call message(msg,level=msglevel)
647 #ifdef WRF_CHEM
648 if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
649 ! sum ground concentrations and check for nans
650 a_chem=0.0
651 g_chem=0.0
652 errors=0
653 do j=jts,jte
654 do k=1,chem_np
655 k_p=chem_pointers(k)
656 do i=its,ite
657 if(chem(i,k1,j,k_p ).ne.chem(i,k1,j,k_p ))errors=errors+1
658 a_chem(k_p)=a_chem(k_p)+chem(i,k1,j,k_p )
659 enddo
660 enddo
661 enddo
662 if(errors>0)call crash('NaN before chem update')
663 call wrf_dm_sum_reals(a_chem,g_chem)
664 call message('Layer1 raw sums before adding fire emissions',level=msglevel)
665 do i=1,chem_np,line
666 m=min(i+line-1,chem_np)
667 write(msg,80)'Emissions ',(trim(chem_names(j)),j=i,m)
668 call message(msg,level=msglevel)
669 write(msg,81)'Layer1 beg',(g_chem(chem_pointers(j)),j=i,m)
670 call message(msg,level=msglevel)
671 call message(' ',level=msglevel)
672 enddo
673 endif
674 #endif
676 call message("Inserting chemical species",level=msglevel)
677 do j=max(jds+1,jts),min(jte,jde-1) ! safe distance from domain boundary
678 jbase=jtfs+jr*(j-jts) ! indexing
679 do i=max(ids+1,its),min(ite,ide-1)
680 ibase=ifts+ir*(i-its) ! indexing
681 !air = 1e6*mw_air/rho(i,kds,j) ! 1e6*mw_air/air density
682 !vol = avgw/dz8w(i,kds,j) ! averaging volume factor / 1st layer depth
684 do joff=0,jr-1
685 j_f=joff+jbase
686 do ioff=0,ir-1
687 i_f=ioff+ibase
690 !*** fire cell (i_f,j_f) contributes to atmosphere cell (i,j) at ground level
692 fuel_burnt = fgip(i_f,j_f) * fuel_frac_burnt(i_f,j_f) ! kg/m^2
693 cat = nfuel_cat(i_f, j_f) ! usually 1 to 13
695 if(cat.lt.no_fuel_cat)t_fuel(cat)=t_fuel(cat) + fuel_burnt
698 !*** chem compounds emissions given in g/kg
700 ! fuel_burnt kg/m^2 * table g/kg -> ppmv = 1e6*mol/mol in 1st layer
701 !AK rho is in kg/m3 so the conversion factor must be scaled by a 1000 to match
702 !emissions in grams
703 ! conv = avgw*1e6*mw_air/(rho(i,k1,j)*dz8w(i,k1,j))
704 conv = avgw*1e3*mw_air/(rho(i,k1,j)*dz8w(i,k1,j))
705 #ifdef WRF_CHEM
707 emis=co (cat)*fuel_burnt ! emission from fire cell in g/m^2
708 t_chem(p_co) = t_chem(p_co) + emis ! add to total
709 ! if(isnan(chem(i,k1,j,p_co )))call crash('NaN before')
710 chem(i,k1,j,p_co )=chem(i,k1,j,p_co ) + emis*conv*imw_co ! add to chem
711 ! if(isnan(chem(i,k1,j,p_co )))call crash('NaN after')
713 emis=ch4 (cat)*fuel_burnt ! emission from fire cell in g/m^2
714 t_chem(p_ch4) = t_chem(p_ch4) + emis ! add to total
715 chem(i,k1,j,p_ch4 )=chem(i,k1,j,p_ch4 ) + emis*conv*imw_ch4 ! add to chem
717 emis=h2 (cat)*fuel_burnt ! emission from fire cell in g/m^2
718 t_chem(p_h2) = t_chem(p_h2) + emis ! add to total
719 chem(i,k1,j,p_h2 )=chem(i,k1,j,p_h2 ) + emis*conv*imw_h2 ! add to chem
721 emis=no (cat)*fuel_burnt ! emission from fire cell in g/m^2
722 t_chem(p_no) = t_chem(p_no) + emis ! add to total
723 chem(i,k1,j,p_no )=chem(i,k1,j,p_no ) + emis*conv*imw_no ! add to chem
725 emis=no2 (cat)*fuel_burnt ! emission from fire cell in g/m^2
726 t_chem(p_no2) = t_chem(p_no2) + emis ! add to total
727 chem(i,k1,j,p_no2 )=chem(i,k1,j,p_no2 ) + emis*conv*imw_no2 ! add to chem
729 emis=so2 (cat)*fuel_burnt ! emission from fire cell in g/m^2
730 t_chem(p_so2) = t_chem(p_so2) + emis ! ads to total
731 chem(i,k1,j,p_so2 )=chem(i,k1,j,p_so2 ) + emis*conv*imw_so2 ! add to chem
733 emis=nh3 (cat)*fuel_burnt ! emission from fire cell in g/m^2
734 t_chem(p_nh3) = t_chem(p_nh3) + emis ! add to total
735 chem(i,k1,j,p_nh3 )=chem(i,k1,j,p_nh3 ) + emis*conv*imw_nh3 ! add to chem
738 !*** other emissions already given in mol/kg
739 ! fuel_burnt kg/m^2 * table mol/kg -> ppmv = 1e6*mol/mol in 1st layer dry air in 1st layer
741 ! same conversion factor but we will not divide by the molecular weight of the compound
742 ! conv = avgw*1e6*mw_air/(rho(i,kds,j)*dz8w(i,kds,j))
744 emis=ald (cat)*fuel_burnt ! emission from fire cell
745 t_chem(p_ald) = t_chem(p_ald) + emis ! add to total
746 chem(i,k1,j,p_ald )=chem(i,k1,j,p_ald )+emis*conv
748 emis=csl (cat)*fuel_burnt ! emission from fire cell
749 t_chem(p_csl) = t_chem(p_csl) + emis ! add to total
750 chem(i,k1,j,p_csl )=chem(i,k1,j,p_csl )+emis*conv
752 emis=eth (cat)*fuel_burnt ! emission from fire cell
753 t_chem(p_eth) = t_chem(p_eth) + emis ! add to total
754 chem(i,k1,j,p_eth )=chem(i,k1,j,p_eth )+emis*conv
756 emis=hc3 (cat)*fuel_burnt ! emission from fire cell
757 t_chem(p_hc3) = t_chem(p_hc3) + emis ! add to total
758 chem(i,k1,j,p_hc3 )=chem(i,k1,j,p_hc3 )+emis*conv
760 emis=hc5 (cat)*fuel_burnt ! emission from fire cell
761 t_chem(p_hc5) = t_chem(p_hc5) + emis ! add to total
762 chem(i,k1,j,p_hc5 )=chem(i,k1,j,p_hc5 )+emis*conv
764 emis=hcho (cat)*fuel_burnt ! emission from fire cell
765 t_chem(p_hcho) = t_chem(p_hcho) + emis ! add to total
766 chem(i,k1,j,p_hcho)=chem(i,k1,j,p_hcho)+emis*conv
768 emis=iso (cat)*fuel_burnt ! emission from fire cell
769 t_chem(p_iso) = t_chem(p_iso) + emis ! add to total
770 chem(i,k1,j,p_iso )=chem(i,k1,j,p_iso )+emis*conv
772 emis=ket (cat)*fuel_burnt ! emission from fire cell
773 t_chem(p_ket) = t_chem(p_ket) + emis ! add to total
774 chem(i,k1,j,p_ket )=chem(i,k1,j,p_ket )+emis*conv
776 emis=mgly (cat)*fuel_burnt ! emission from fire cell
777 t_chem(p_mgly) = t_chem(p_mgly) + emis ! add to total
778 chem(i,k1,j,p_mgly)=chem(i,k1,j,p_mgly)+emis*conv
780 emis=ol2 (cat)*fuel_burnt ! emission from fire cell
781 t_chem(p_ol2) = t_chem(p_ol2) + emis ! add to total
782 chem(i,k1,j,p_ol2 )=chem(i,k1,j,p_ol2 )+emis*conv
784 emis=olt (cat)*fuel_burnt ! emission from fire cell
785 t_chem(p_olt) = t_chem(p_olt) + emis ! add to total
786 chem(i,k1,j,p_olt )=chem(i,k1,j,p_olt )+emis*conv
788 emis=oli (cat)*fuel_burnt ! emission from fire cell
789 t_chem(p_oli) = t_chem(p_oli) + emis ! add to total
790 chem(i,k1,j,p_oli )=chem(i,k1,j,p_oli )+emis*conv
792 emis=ora2 (cat)*fuel_burnt ! emission from fire cell
793 t_chem(p_ora2) = t_chem(p_ora2) + emis ! add to total
794 chem(i,k1,j,p_ora2)=chem(i,k1,j,p_ora2)+emis*conv
796 emis=tol (cat)*fuel_burnt ! emission from fire cell
797 t_chem(p_tol) = t_chem(p_tol) + emis ! add to total
798 chem(i,k1,j,p_tol )=chem(i,k1,j,p_tol )+emis*conv
800 emis=xyl (cat)*fuel_burnt ! emission from fire cell
801 t_chem(p_xyl) = t_chem(p_xyl) + emis ! add to total
802 chem(i,k1,j,p_xyl )=chem(i,k1,j,p_xyl )+emis*conv
804 emis=bigalk (cat)*fuel_burnt ! emission from fire cell
805 t_chem(p_bigalk) = t_chem(p_bigalk) + emis ! add to total
806 chem(i,k1,j,p_bigalk )=chem(i,k1,j,p_bigalk )+emis*conv
808 emis=bigene (cat)*fuel_burnt ! emission from fire cell
809 t_chem(p_bigene) = t_chem(p_bigene) + emis ! add to total
810 chem(i,k1,j,p_bigene )=chem(i,k1,j,p_bigene )+emis*conv
812 emis=c10h16 (cat)*fuel_burnt ! emission from fire cell
813 t_chem(p_c10h16) = t_chem(p_c10h16) + emis ! add to total
814 chem(i,k1,j,p_c10h16 )=chem(i,k1,j,p_c10h16 )+emis*conv
816 emis=c2h4 (cat)*fuel_burnt ! emission from fire cell
817 t_chem(p_c2h4) = t_chem(p_c2h4) + emis ! add to total
818 chem(i,k1,j,p_c2h4 )=chem(i,k1,j,p_c2h4 )+emis*conv
820 emis=c2h5oh (cat)*fuel_burnt ! emission from fire cell
821 t_chem(p_c2h5oh) = t_chem(p_c2h5oh) + emis ! add to total
822 chem(i,k1,j,p_c2h5oh )=chem(i,k1,j,p_c2h5oh )+emis*conv
824 emis=c2h6 (cat)*fuel_burnt ! emission from fire cell
825 t_chem(p_c2h6) = t_chem(p_c2h6) + emis ! add to total
826 chem(i,k1,j,p_c2h6 )=chem(i,k1,j,p_c2h6 )+emis*conv
828 emis=c3h6 (cat)*fuel_burnt ! emission from fire cell
829 t_chem(p_c3h6) = t_chem(p_c3h6) + emis ! add to total
830 chem(i,k1,j,p_c3h6 )=chem(i,k1,j,p_c3h6 )+emis*conv
832 emis=c3h8 (cat)*fuel_burnt ! emission from fire cell
833 t_chem(p_c3h8) = t_chem(p_c3h8) + emis ! add to total
834 chem(i,k1,j,p_c3h8 )=chem(i,k1,j,p_c3h8 )+emis*conv
836 emis=ch3cooh (cat)*fuel_burnt ! emission from fire cell
837 t_chem(p_ch3cooh) = t_chem(p_ch3cooh) + emis ! add to total
838 chem(i,k1,j,p_ch3cooh )=chem(i,k1,j,p_ch3cooh )+emis*conv
840 emis=ch3oh (cat)*fuel_burnt ! emission from fire cell
841 t_chem(p_ch3oh) = t_chem(p_ch3oh) + emis ! add to total
842 chem(i,k1,j,p_ch3oh )=chem(i,k1,j,p_ch3oh )+emis*conv
844 emis=cres (cat)*fuel_burnt ! emission from fire cell
845 t_chem(p_cres) = t_chem(p_cres) + emis ! add to total
846 chem(i,k1,j,p_cres )=chem(i,k1,j,p_cres )+emis*conv
848 emis=glyald (cat)*fuel_burnt ! emission from fire cell
849 t_chem(p_glyald) = t_chem(p_glyald) + emis ! add to total
850 chem(i,k1,j,p_glyald )=chem(i,k1,j,p_glyald )+emis*conv
852 emis=isopr (cat)*fuel_burnt ! emission from fire cell
853 t_chem(p_isopr) = t_chem(p_isopr) + emis ! add to total
854 chem(i,k1,j,p_isopr )=chem(i,k1,j,p_isopr )+emis*conv
856 emis=macr (cat)*fuel_burnt ! emission from fire cell
857 t_chem(p_macr) = t_chem(p_macr) + emis ! add to total
858 chem(i,k1,j,p_macr )=chem(i,k1,j,p_macr )+emis*conv
860 emis=mek (cat)*fuel_burnt ! emission from fire cell
861 t_chem(p_mek) = t_chem(p_mek) + emis ! add to total
862 chem(i,k1,j,p_mek )=chem(i,k1,j,p_mek )+emis*conv
864 emis=mvk (cat)*fuel_burnt ! emission from fire cell
865 t_chem(p_mvk) = t_chem(p_mvk) + emis ! add to total
866 chem(i,k1,j,p_mvk )=chem(i,k1,j,p_mvk )+emis*conv
868 ! aerosols
869 ! fuel_burnt kg/m^2 * table g/kg -> ug/kg dry air in 1st layer
871 ! see also chem/emissions_driver.F line 515
873 conv = avgw*1e6/(rho(i,k1,j)*dz8w(i,k1,j))
875 emis=p25 (cat)*fuel_burnt ! emission from fire cell
876 t_chem(p_p25) = t_chem(p_p25) + emis ! add to total
877 chem(i,k1,j,p_p25 )=chem(i,k1,j,p_p25 )+emis*conv
879 emis=p25i (cat)*fuel_burnt ! emission from fire cell
880 t_chem(p_p25i) = t_chem(p_p25i) + emis ! add to total
881 chem(i,k1,j,p_p25i )=chem(i,k1,j,p_p25i )+emis*conv
883 emis=p25j (cat)*fuel_burnt ! emission from fire cell
884 t_chem(p_p25j) = t_chem(p_p25j) + emis ! add to total
885 chem(i,k1,j,p_p25j )=chem(i,k1,j,p_p25j )+emis*conv
887 emis=oc1 (cat)*fuel_burnt ! emission from fire cell
888 t_chem(p_oc1 ) = t_chem(p_oc1 ) + emis ! add to total
889 chem(i,k1,j,p_oc1 )=chem(i,k1,j,p_oc1 )+emis*conv
891 emis=oc2 (cat)*fuel_burnt ! emission from fire cell
892 t_chem(p_oc2 ) = t_chem(p_oc2 ) + emis ! add to total
893 chem(i,k1,j,p_oc2 )=chem(i,k1,j,p_oc2 )+emis*conv
895 emis=bc1 (cat)*fuel_burnt ! emission from fire cell
896 t_chem(p_bc1 ) = t_chem(p_bc1 ) + emis ! add to total
897 chem(i,k1,j,p_bc1 )=chem(i,k1,j,p_bc1 )+emis*conv
899 emis=bc2 (cat)*fuel_burnt ! emission from fire cell
900 t_chem(p_bc2 ) = t_chem(p_bc2 ) + emis ! add to total
901 chem(i,k1,j,p_bc2 )=chem(i,k1,j,p_bc2 )+emis*conv
903 emis=sulf (cat)*fuel_burnt ! emission from fire cell
904 t_chem(p_sulf ) = t_chem(p_sulf ) + emis ! add to total
905 chem(i,k1,j,p_sulf )=chem(i,k1,j,p_sulf )+emis*conv
907 emis=dms (cat)*fuel_burnt ! emission from fire cell
908 t_chem(p_dms ) = t_chem(p_dms ) + emis ! add to total
909 chem(i,k1,j,p_dms )=chem(i,k1,j,p_dms )+emis*conv
911 emis=msa (cat)*fuel_burnt ! emission from fire cell
912 t_chem(p_msa ) = t_chem(p_msa ) + emis ! add to total
913 chem(i,k1,j,p_msa )=chem(i,k1,j,p_msa )+emis*conv
915 emis=dust_1 (cat)*fuel_burnt ! emission from fire cell
916 t_chem(p_dust_1 ) = t_chem(p_dust_1 ) + emis ! add to total
917 chem(i,k1,j,p_dust_1 )=chem(i,k1,j,p_dust_1 )+emis*conv
919 emis=dust_2 (cat)*fuel_burnt ! emission from fire cell
920 t_chem(p_dust_2 ) = t_chem(p_dust_2 ) + emis ! add to total
921 chem(i,k1,j,p_dust_2 )=chem(i,k1,j,p_dust_2 )+emis*conv
923 emis=dust_3 (cat)*fuel_burnt ! emission from fire cell
924 t_chem(p_dust_3 ) = t_chem(p_dust_3 ) + emis ! add to total
925 chem(i,k1,j,p_dust_3 )=chem(i,k1,j,p_dust_3 )+emis*conv
927 emis=dust_4 (cat)*fuel_burnt ! emission from fire cell
928 t_chem(p_dust_4 ) = t_chem(p_dust_4 ) + emis ! add to total
929 chem(i,k1,j,p_dust_4 )=chem(i,k1,j,p_dust_4 )+emis*conv
931 emis=dust_5 (cat)*fuel_burnt ! emission from fire cell
932 t_chem(p_dust_5 ) = t_chem(p_dust_5 ) + emis ! add to total
933 chem(i,k1,j,p_dust_5 )=chem(i,k1,j,p_dust_5 )+emis*conv
935 emis=seas_1 (cat)*fuel_burnt ! emission from fire cell
936 t_chem(p_seas_1 ) = t_chem(p_seas_1 ) + emis ! add to total
937 chem(i,k1,j,p_seas_1 )=chem(i,k1,j,p_seas_1 )+emis*conv
939 emis=seas_2 (cat)*fuel_burnt ! emission from fire cell
940 t_chem(p_seas_2 ) = t_chem(p_seas_2 ) + emis ! add to total
941 chem(i,k1,j,p_seas_2 )=chem(i,k1,j,p_seas_2 )+emis*conv
943 emis=seas_3 (cat)*fuel_burnt ! emission from fire cell
944 t_chem(p_seas_3 ) = t_chem(p_seas_3 ) + emis ! add to total
945 chem(i,k1,j,p_seas_3 )=chem(i,k1,j,p_seas_3 )+emis*conv
947 emis=seas_4 (cat)*fuel_burnt ! emission from fire cell
948 t_chem(p_seas_4 ) = t_chem(p_seas_4 ) + emis ! add to total
949 chem(i,k1,j,p_seas_4 )=chem(i,k1,j,p_seas_4 )+emis*conv
951 emis=p10 (cat)*fuel_burnt ! emission from fire cell
952 t_chem(p_p10 ) = t_chem(p_p10 ) + emis ! add to total
953 chem(i,k1,j,p_p10 )=chem(i,k1,j,p_p10 )+emis*conv
955 #endif
956 if (num_tracer >0)then
958 ! treat tracers exactly the same as aerosols, emissions g/kg fuel burned, tracer concentration ug/kg dry air
959 conv = avgw*1e6/(rho(i,k1,j)*dz8w(i,k1,j))
961 emis=tr17_1 (cat)*fuel_burnt ! emission from fire cell
962 t_tracer(p_tr17_1) = t_tracer(p_tr17_1) + emis ! add to total
963 tracer(i,k1,j,p_tr17_1 )=tracer(i,k1,j,p_tr17_1 )+emis*conv
965 emis=tr17_2 (cat)*fuel_burnt ! emission from fire cell
966 t_tracer(p_tr17_2) = t_tracer(p_tr17_2) + emis ! add to total
967 tracer(i,k1,j,p_tr17_2 )=tracer(i,k1,j,p_tr17_2 )+emis*conv
969 emis=tr17_3 (cat)*fuel_burnt ! emission from fire cell
970 t_tracer(p_tr17_3) = t_tracer(p_tr17_3) + emis ! add to total
971 tracer(i,k1,j,p_tr17_3 )=tracer(i,k1,j,p_tr17_3 )+emis*conv
973 emis=tr17_4 (cat)*fuel_burnt ! emission from fire cell
974 t_tracer(p_tr17_4) = t_tracer(p_tr17_4) + emis ! add to total
975 tracer(i,k1,j,p_tr17_4 )=tracer(i,k1,j,p_tr17_4 )+emis*conv
977 emis=tr17_5 (cat)*fuel_burnt ! emission from fire cell
978 t_tracer(p_tr17_5) = t_tracer(p_tr17_5) + emis ! add to total
979 tracer(i,k1,j,p_tr17_5 )=tracer(i,k1,j,p_tr17_5 )+emis*conv
981 emis=tr17_6 (cat)*fuel_burnt ! emission from fire cell
982 t_tracer(p_tr17_6) = t_tracer(p_tr17_6) + emis ! add to total
983 tracer(i,k1,j,p_tr17_6 )=tracer(i,k1,j,p_tr17_6 )+emis*conv
985 emis=tr17_7 (cat)*fuel_burnt ! emission from fire cell
986 t_tracer(p_tr17_7) = t_tracer(p_tr17_7) + emis ! add to total
987 tracer(i,k1,j,p_tr17_7 )=tracer(i,k1,j,p_tr17_7 )+emis*conv
989 emis=tr17_8 (cat)*fuel_burnt ! emission from fire cell
990 t_tracer(p_tr17_8) = t_tracer(p_tr17_8) + emis ! add to total
991 tracer(i,k1,j,p_tr17_8 )=tracer(i,k1,j,p_tr17_8 )+emis*conv
992 endif
993 enddo
994 enddo
995 enddo
996 enddo
998 #ifdef WRF_CHEM
999 if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
1000 ! sum ground concentrations and check for nans
1001 a_chem=0.0
1002 g_chem=0.0
1003 errors=0
1004 do j=jts,jte
1005 do k=1,chem_np
1006 k_p=chem_pointers(k)
1007 do i=its,ite
1008 if(chem(i,k1,j,k_p ).ne.chem(i,k1,j,k_p ))errors=errors+1
1009 a_chem(k_p)=a_chem(k_p)+chem(i,k1,j,k_p )
1010 enddo
1011 enddo
1012 enddo
1013 if(errors>0)call crash('NaN after chem update')
1014 call wrf_dm_sum_reals(a_chem,g_chem)
1015 call message('Layer1 raw sums after adding fire emissions',level=msglevel)
1016 do i=1,chem_np,line
1017 m=min(i+line-1,chem_np)
1018 write(msg,80)'Emissions ',(trim(chem_names(j)),j=i,m)
1019 call message(msg,level=msglevel)
1020 write(msg,81)'Layer1 end',(g_chem(chem_pointers(j)),j=i,m)
1021 call message(msg,level=msglevel)
1022 call message(' ',level=msglevel)
1023 enddo
1024 endif
1025 #endif
1027 if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
1028 call message("Computing totals over all processes",level=msglevel)
1029 ! sum over processes and add to cumulative sums
1031 ! fuel burned
1032 call wrf_dm_sum_reals(t_fuel,s_fuel)
1033 ! scale
1034 s_fuel = s_fuel*dx*dy
1035 ! get rates
1036 r_fuel = s_fuel/dt
1037 ! add to cumulative totals
1038 if(size(c_fuel).ne.nfuelcats)call crash('add_fire_emissions: bad size c_fuel')
1039 c_fuel = c_fuel + s_fuel
1041 #ifdef WRF_CHEM
1042 call wrf_dm_sum_reals(a_chem,g_chem)
1043 ! chem
1044 call wrf_dm_sum_reals(t_chem,s_chem)
1045 s_chem = s_chem*dx*dy
1046 r_chem = s_chem/dt
1047 if(size(c_chem).ne.num_chem)call crash('add_fire_emissions: bad size c_chem')
1048 c_chem = c_chem + s_chem
1050 call message('Total emissions in g or mol per the table',level=msglevel)
1051 do i=1,chem_np,line
1052 m=min(i+line-1,chem_np)
1053 write(msg,80)'Emissions ',(trim(chem_names(j)),j=i,m)
1054 call message(msg,level=msglevel)
1055 write(msg,81)'Cumulative',(c_chem(chem_pointers(j)),j=i,m)
1056 call message(msg,level=msglevel)
1057 write(msg,81)'Rate per s',(r_chem(chem_pointers(j)),j=i,m)
1058 call message(msg,level=msglevel)
1059 call message(' ',level=msglevel)
1060 enddo
1061 #endif
1063 if(num_tracer >0)then
1064 call message("Computing totals of tracers over all processes",level=msglevel)
1065 ! tracer
1066 call wrf_dm_sum_reals(t_tracer,s_tracer)
1067 s_tracer = s_tracer*dx*dy
1068 r_tracer = s_tracer/dt
1069 if(size(c_tracer).ne.num_tracer)call crash('add_fire_emissions: bad size c_tracer')
1070 c_tracer = c_tracer + s_tracer
1072 call message('Total emissions in g',level=msglevel)
1073 do i=1,tracer_np,line
1074 m=min(i+line-1,tracer_np)
1075 write(msg,80)'Emissions ',(trim(tracer_names(j)),j=i,m)
1076 call message(msg,level=msglevel)
1077 write(msg,81)'Cumulative',(c_tracer(tracer_pointers(j)),j=i,m)
1078 call message(msg,level=msglevel)
1079 write(msg,81)'Rate per s',(r_tracer(tracer_pointers(j)),j=i,m)
1080 call message(msg,level=msglevel)
1081 call message(' ',level=msglevel)
1082 enddo
1083 endif
1086 do cat=1,nfuelcats
1087 if(c_fuel(cat) > 0.)then
1088 write(msg,83)fuel_name(cat),' burned',c_fuel(cat),'kg, rate',r_fuel(cat),' kg/s'
1089 call message(msg,level=msglevel)
1090 endif
1091 enddo
1092 write(msg,83)'Total fuel',' burned',sum(c_fuel),'kg, rate',sum(r_fuel),' kg/s'
1093 call message(msg,level=msglevel)
1095 endif
1096 80 format(a,8a11)
1097 81 format(a,8e11.3)
1098 83 format(a30,a,g14.4,a,g14.4,a,a)
1100 return
1102 9 continue
1103 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1104 write(msg,91)ifts,ifte,jtfs,jfte,ifms,ifme,jfms,jfme
1105 call message(msg,level=0)
1106 write(msg,91)its,ite,jts,jte,ims,ime,jms,jme
1107 call message(msg,level=0)
1108 write(msg,92)'input mesh size:',isz2,jsz2
1109 call message(msg,level=0)
1110 91 format('dimensions: ',8i8)
1111 write(msg,92)'output mesh size:',isz1,jsz1
1112 call message(msg,level=0)
1113 92 format(a,2i8)
1114 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1115 call crash('add_fire_emissions: bad mesh sizes')
1117 end subroutine add_fire_emissions
1119 !
1120 !***
1121 !
1123 subroutine check_pointers(array_name,array,pointer_names,pointers)
1124 implicit none
1126 !*** arguments
1127 character(len=*)::array_name
1128 real, dimension(:,:,:,:)::array
1129 character(len=*), dimension(:)::pointer_names
1130 integer, dimension(:)::pointers
1132 !*** local
1133 integer::np,i,m,j
1134 character(len=256)::msg
1136 !** executable
1137 np=ubound(pointers,1)
1140 993 format(3a,4(1x,i3,':',i3))
1141 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1142 write(msg,993)'array ',array_name,' has dimensions ',&
1143 lbound(array,1),ubound(array,1), &
1144 lbound(array,2),ubound(array,2), &
1145 lbound(array,3),ubound(array,3), &
1146 lbound(array,4),ubound(array,4)
1147 call message(msg)
1149 do i=1,np,line
1150 m=min(i+line-1,np)
1151 write(msg,'(a,8(1x,a8))')'Species',(trim(pointer_names(j)),j=i,m)
1152 call message(msg,level=msglevel)
1153 write(msg,'(a,8i9)') 'Pointer',(pointers(j),j=i,m)
1154 call message(msg,level=msglevel)
1155 call message(' ')
1156 enddo
1158 if(maxval(pointers) > ubound(array,4) .or. minval(pointers) < lbound(array,4))then
1159 write(msg,'(3a)')'add_fire_emissions: a ',array_name,' pointer is out of bounds'
1160 call crash(msg)
1161 endif
1162 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1163 end subroutine check_pointers
1166 !
1167 !***
1168 !
1171 SUBROUTINE fire_tendency( &
1172 ids,ide, kds,kde, jds,jde, & ! dimensions
1173 ims,ime, kms,kme, jms,jme, &
1174 its,ite, kts,kte, jts,jte, &
1175 grnhfx,grnqfx,canhfx,canqfx, & ! heat fluxes summed up to atm grid
1176 alfg,alfc,z1can, & ! coeffients, properties, geometry
1177 zs,z_at_w,dz8w,mu,rho, &
1178 rthfrten,rqvfrten) ! theta and Qv tendencies
1180 ! This routine is atmospheric physics
1181 ! it does NOT go into module_fr_sfire_phys because it is not fire physics
1183 ! taken from the code by Ned Patton, only order of arguments change to the convention here
1184 ! --- this routine takes fire generated heat and moisture fluxes and
1185 ! calculates their influence on the theta and water vapor
1186 ! --- note that these tendencies are valid at the Arakawa-A location
1188 IMPLICIT NONE
1190 ! --- incoming variables
1192 INTEGER , INTENT(in) :: ids,ide, kds,kde, jds,jde, &
1193 ims,ime, kms,kme, jms,jme, &
1194 its,ite, kts,kte, jts,jte
1196 REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: grnhfx,grnqfx ! W/m^2
1197 REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: canhfx,canqfx ! W/m^2
1198 REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: zs ! topography (m abv sealvl)
1199 REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: mu ! dry air mass (Pa)
1201 REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: z_at_w ! m abv sealvl
1202 REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: dz8w ! dz across w-lvl
1203 REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: rho ! density
1205 REAL, INTENT(in) :: alfg ! extinction depth ground fire heat (m)
1206 REAL, INTENT(in) :: alfc ! extinction depth crown fire heat (m)
1207 REAL, INTENT(in) :: z1can ! height of crown fire heat release (m)
1209 ! --- outgoing variables
1211 REAL, INTENT(out), DIMENSION( ims:ime,kms:kme,jms:jme ) :: &
1212 rthfrten, & ! theta tendency from fire (in mass units)
1213 rqvfrten ! Qv tendency from fire (in mass units)
1214 ! --- local variables
1216 INTEGER :: i,j,k
1217 INTEGER :: i_st,i_en, j_st,j_en, k_st,k_en
1219 REAL :: cp_i
1220 REAL :: rho_i
1221 REAL :: xlv_i
1222 REAL :: z_w
1223 REAL :: fact_g, fact_c
1224 REAL :: alfg_i, alfc_i
1226 REAL, DIMENSION( its:ite,kts:kte,jts:jte ) :: hfx,qfx
1228 !! character(len=128)::msg
1231 do j=jts,jte
1232 do k=kts,min(kte+1,kde)
1233 do i=its,ite
1234 rthfrten(i,k,j)=0.
1235 rqvfrten(i,k,j)=0.
1236 enddo
1237 enddo
1238 enddo
1241 ! --- set some local constants
1244 cp_i = 1./cp ! inverse of specific heat
1245 xlv_i = 1./xlv ! inverse of latent heat
1246 alfg_i = 1./alfg
1247 alfc_i = 1./alfc
1249 !!write(msg,'(8e11.3)')cp,cp_i,xlv,xlv_i,alfg,alfc,z1can
1250 !!call message(msg)
1252 call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnhfx,'fire_tendency:grnhfx')
1253 call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnqfx,'fire_tendency:grnqfx')
1255 ! --- set loop indicies : note that
1257 i_st = MAX(its,ids+1)
1258 i_en = MIN(ite,ide-1)
1259 k_st = kts
1260 k_en = MIN(kte,kde-1)
1261 j_st = MAX(jts,jds+1)
1262 j_en = MIN(jte,jde-1)
1264 ! --- distribute fluxes
1266 DO j = j_st,j_en
1267 DO k = k_st,k_en
1268 DO i = i_st,i_en
1270 ! --- set z (in meters above ground)
1272 z_w = z_at_w(i,k,j) - zs(i,j) ! should be zero when k=k_st
1274 ! --- heat flux
1276 fact_g = cp_i * EXP( - alfg_i * z_w )
1277 IF ( z_w < z1can ) THEN
1278 fact_c = cp_i
1279 ELSE
1280 fact_c = cp_i * EXP( - alfc_i * (z_w - z1can) )
1281 END IF
1282 hfx(i,k,j) = fact_g * grnhfx(i,j) + fact_c * canhfx(i,j)
1284 !! write(msg,2)i,k,j,z_w,grnhfx(i,j),hfx(i,k,j)
1285 !!2 format('hfx:',3i4,6e11.3)
1286 !! call message(msg)
1288 ! --- vapor flux
1290 fact_g = xlv_i * EXP( - alfg_i * z_w )
1291 IF (z_w < z1can) THEN
1292 fact_c = xlv_i
1293 ELSE
1294 fact_c = xlv_i * EXP( - alfc_i * (z_w - z1can) )
1295 END IF
1296 qfx(i,k,j) = fact_g * grnqfx(i,j) + fact_c * canqfx(i,j)
1298 !! if(hfx(i,k,j).ne.0. .or. qfx(i,k,j) .ne. 0.)then
1299 !! write(msg,1)i,k,j,hfx(i,k,j),qfx(i,k,j)
1300 !!1 format('tend:',3i6,2e11.3)
1301 !! call message(msg)
1302 ! endif
1304 END DO
1305 END DO
1306 END DO
1308 ! --- add flux divergence to tendencies
1309 !
1310 ! multiply by dry air mass (mu) to eliminate the need to
1311 ! call sr. calculate_phy_tend (in dyn_em/module_em.F)
1313 ! print *,'fire_tendency:',i_st,i_en,j_st,j_en
1314 DO j = j_st,j_en
1315 DO k = k_st,k_en-1
1316 DO i = i_st,i_en
1318 rho_i = 1./rho(i,k,j)
1320 rthfrten(i,k,j) = - mu(i,j) * rho_i * (hfx(i,k+1,j)-hfx(i,k,j)) / dz8w(i,k,j)
1321 rqvfrten(i,k,j) = - mu(i,j) * rho_i * (qfx(i,k+1,j)-qfx(i,k,j)) / dz8w(i,k,j)
1323 ! print *,i,j,k,rthfrten(i,k,j)
1325 END DO
1326 END DO
1327 END DO
1329 call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rthfrten,'fire_tendency:rthfrten')
1330 call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rqvfrten,'fire_tendency:rqvfrten')
1332 RETURN
1334 END SUBROUTINE fire_tendency
1336 !
1337 !***
1338 !
1339 subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no output
1340 ids,ide, kds,kde, jds,jde, & ! atm grid dimensions
1341 ims,ime, kms,kme, jms,jme, &
1342 ips,ipe,jps,jpe, &
1343 its,ite,jts,jte, &
1344 ifds, ifde, jfds, jfde, & ! fire grid dimensions
1345 ifms, ifme, jfms, jfme, &
1346 ifps, ifpe, jfps, jfpe, & ! fire patch bounds
1347 ifts,ifte,jfts,jfte, &
1348 ir,jr, & ! atm/fire grid ratio
1349 u_frame, v_frame, & ! velocity frame correction
1350 u,v, & ! atm grid arrays in
1351 ph,phb, &
1352 z0,zs, &
1353 uah,vah, &
1354 uf,vf) ! fire grid arrays out
1356 implicit none
1357 ! Jan Mandel, October 2010
1358 !*** purpose: interpolate winds and height
1360 !*** arguments
1361 integer, intent(in)::id
1362 integer, intent(in):: &
1363 ids,ide, kds,kde, jds,jde, & ! atm domain bounds
1364 ims,ime, kms,kme, jms,jme, & ! atm memory bounds
1365 ips,ipe,jps,jpe, &
1366 its,ite,jts,jte, & ! atm tile bounds
1367 ifds, ifde, jfds, jfde, & ! fire domain bounds
1368 ifms, ifme, jfms, jfme, & ! fire memory bounds
1369 ifps, ifpe, jfps, jfpe, & ! fire patch bounds
1370 ifts,ifte,jfts,jfte, & ! fire tile bounds
1371 ir,jr ! atm/fire grid refinement ratio
1372 real, intent(in):: u_frame, v_frame ! velocity frame correction
1373 real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::&
1374 u,v, & ! atm wind velocity, staggered
1375 ph,phb ! geopotential
1376 real,intent(in),dimension(ims:ime,jms:jme)::&
1377 z0, & ! roughness height
1378 zs ! terrain height
1379 real,intent(out),dimension(ims:ime,jms:jme)::&
1380 uah, & ! atm wind at fire_wind_height, diagnostics
1381 vah !
1382 real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
1383 uf,vf ! wind velocity fire grid nodes
1386 !*** local
1387 character(len=256)::msg
1388 #define TDIMS its-2,ite+2,jts-2,jte+2
1389 real, dimension(its-2:ite+2,jts-2:jte+2):: ua,va ! atm winds, interpolated over height, still staggered grid
1390 real, dimension(its-2:ite+2,kds:kde,jts-2:jte+2):: altw,altub,altvb,hgtu,hgtv ! altitudes
1391 integer:: i,j,k,ifts1,ifte1,jfts1,jfte1,ite1,jte1
1392 integer::itst,itet,jtst,jtet,itsu,iteu,jtsu,jteu,itsv,itev,jtsv,jtev
1393 integer::kdmax,its1,jts1,ips1,jps1
1394 integer::itsou,iteou,jtsou,jteou,itsov,iteov,jtsov,jteov
1395 real:: ground,loght,loglast,logz0,logfwh,ht,zr
1396 real::r_nan
1397 integer::i_nan
1398 equivalence (i_nan,r_nan)
1400 !*** executable
1402 ! debug init local arrays
1403 i_nan=2147483647
1404 ua=r_nan
1405 va=r_nan
1406 altw=r_nan
1407 altub=r_nan
1408 hgtu=r_nan
1409 hgtv=r_nan
1412 if(kds.ne.1)then
1413 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1414 write(msg,*)'WARNING: bottom index kds=',kds,' should be 1?'
1415 call message(msg)
1416 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1417 endif
1419 ! ^ j
1420 ! ------------ |
1421 ! | | ----> i
1422 ! u p |
1423 ! | | nodes in cell (i,j)
1424 ! ------v----- view from top
1425 !
1426 ! v(ide,jde+1)
1427 ! -------x------
1428 ! | |
1429 ! | |
1430 ! u(ide,jde) x x x u(ide+1,jde)
1431 ! | p(ide,hde) |
1432 ! | | p=ph,phb,z0,...
1433 ! -------x------
1434 ! v(ide,jde)
1435 !
1436 ! staggered values set u(ids:ide+1,jds:jde) v(ids:ide,jds:jde+1)
1437 ! p=ph+phb set at (ids:ide,jds:jde)
1438 ! location of u(i,j) needs p(i-1,j) and p(i,j)
1439 ! location of v(i,j) needs p(i,j-1) and p(i,j)
1440 ! *** NOTE need HALO in ph, phb
1441 ! so we can compute only u(ids+1:ide,jds:jde) v(ids:ide,jds+1,jde)
1442 ! unless we extend p at the boundary
1443 ! but because we care about the fire way in the inside it does not matter
1444 ! if the fire gets close to domain boundary the simulation is over anyway
1446 ite1=snode(ite,ide,1)
1447 jte1=snode(jte,jde,1)
1448 ! do this in any case to check for nans
1449 call print_3d_stats(its,ite1,kds,kde,jts,jte,ims,ime,kms,kme,jms,jme,u,'wind U in')
1450 call print_3d_stats(its,ite,kds,kde,jts,jte1,ims,ime,kms,kme,jms,jme,v,'wind V in')
1452 if(fire_print_msg.gt.0)then
1453 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1454 write(msg,'(a,f7.2,a)')'interpolate_atm2fire: log-interpolation of wind to',fire_wind_height,' m'
1455 call message(msg)
1456 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1457 endif
1459 ! indexing
1461 ! for w
1462 itst=ifval(ids.eq.its,its,its-1)
1463 itet=ifval(ide.eq.ite,ite,ite+1)
1464 jtst=ifval(jds.ge.jts,jts,jts-1)
1465 jtet=ifval(jde.eq.jte,jte,jte+1)
1467 ! for u
1468 itsu=ifval(ids.eq.its,its+1,its) ! staggered direction
1469 iteu=ifval(ide.eq.ite,ite,ite+1) ! staggered direction
1470 jtsu=ifval(jds.ge.jts,jts,jts-1)
1471 jteu=ifval(jde.eq.jte,jte,jte+1)
1473 ! for v
1474 jtsv=ifval(jds.eq.jts,jts+1,jts) ! staggered direction
1475 jtev=ifval(jde.eq.jte,jte,jte+1) ! staggered direction
1476 itsv=ifval(ids.ge.its,its,its-1)
1477 itev=ifval(ide.eq.ite,ite,ite+1)
1480 if(fire_print_msg.ge.1)then
1481 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1482 write(msg,7001)'atm input ','tile',its,ite,jts,jte
1483 call message(msg)
1484 write(msg,7001)'altw ','tile',itst,itet,jtst,jtet
1485 call message(msg)
1486 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1487 endif
1488 7001 format(a,' dimensions ',a4,':',i6,' to ',i6,' by ',i6,' to ',i6)
1490 !**********************************************************
1491 !* *
1492 !* find the altitude of the w points *
1493 !* *
1494 !**********************************************************
1495 !! in future, replace by z8w & test if the same
1497 kdmax=kde-1 ! max layer to interpolate from, can be less
1499 do j = jtst,jtet
1500 do k=kds,kdmax+1
1501 do i = itst,itet
1502 altw(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g ! altitude of the bottom w-point
1503 enddo
1504 enddo
1505 enddo
1507 ! values at u points
1508 if(fire_print_msg.ge.1)then
1509 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1510 write(msg,7001)'u interp at','tile',itsu,iteu,jtsu,jteu
1511 call message(msg)
1512 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1513 endif
1515 !**********************************************************
1516 !* *
1517 !* interpolate the altitude from w points to u points *
1518 !* *
1519 !**********************************************************
1521 do j = jtsu,jteu
1522 do k=kds,kdmax+1
1523 do i = itsu,iteu
1524 altub(i,k,j)= 0.5*(altw(i-1,k,j)+altw(i,k,j)) ! altitude of the bottom point under u-point
1525 enddo
1526 enddo
1527 do k=kds,kdmax
1528 do i = itsu,iteu
1529 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
1530 enddo
1531 enddo
1532 enddo
1534 ! values at v points
1535 if(fire_print_msg.ge.1)then
1536 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1537 write(msg,7001)'v interp at','tile',itsv,itev,jtsv,jtev
1538 call message(msg)
1539 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1540 endif
1542 !**********************************************************
1543 !* *
1544 !* interpolate the altitude from w points to v points *
1545 !* *
1546 !**********************************************************
1548 do j = jtsv,jtev
1549 do k=kds,kdmax+1
1550 do i = itsv,itev
1551 altvb(i,k,j)= 0.5*(altw(i,k,j-1)+altw(i,k,j)) ! altitude of the bottom point under v-point
1552 enddo
1553 enddo
1554 do k=kds,kdmax
1555 do i = itsv,itev
1556 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
1557 enddo
1558 enddo
1559 enddo
1561 #ifdef DEBUG_OUT
1562 call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,altub,'altub',id)
1563 call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,altvb,'altvb',id)
1564 call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,hgtu,'hgtu',id)
1565 call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,hgtv,'hgtv',id)
1566 #endif
1568 logfwh = log(fire_wind_height)
1570 !**********************************************************
1571 !* *
1572 !* interpolate u vertically to fire_wind_height *
1573 !* *
1574 !**********************************************************
1576 ! interpolate u, staggered in X
1578 do j = jtsu,jteu ! compute on domain by 1 smaller, otherwise z0 and ph not available
1579 do i = itsu,iteu ! compute with halo 2
1580 zr = 0.5*(z0(i,j)+z0(i-1,j)) ! interpolated roughness length under this u point
1581 if(fire_wind_height > zr)then !
1582 do k=kds,kdmax
1583 ht = hgtu(i,k,j) ! height of this u point above the ground
1584 if( .not. ht < fire_wind_height) then ! found layer k this point is in
1585 loght = log(ht)
1586 if(k.eq.kds)then ! first layer, log linear interpolation from 0 at zr
1587 logz0 = log(zr)
1588 ua(i,j)= u(i,k,j)*(logfwh-logz0)/(loght-logz0)
1589 else ! log linear interpolation
1590 loglast=log(hgtu(i,k-1,j))
1591 ua(i,j)= u(i,k-1,j) + (u(i,k,j) - u(i,k-1,j)) * ( logfwh - loglast) / (loght - loglast)
1592 endif
1593 goto 10
1594 endif
1595 if(k.eq.kdmax)then ! last layer, still not high enough
1596 ua(i,j)=u(i,k,j)
1597 endif
1598 enddo
1599 10 continue
1600 else ! roughness higher than the fire wind height
1601 ua(i,j)=0.
1602 endif
1603 enddo
1604 enddo
1607 !**********************************************************
1608 !* *
1609 !* interpolate v vertically to fire_wind_height *
1610 !* *
1611 !**********************************************************
1613 ! interpolate v, staggered in Y
1615 do j = jtsv,jtev
1616 do i = itsv,itev
1617 zr = 0.5*(z0(i,j-1)+z0(i,j)) ! roughness length under this v point
1618 if(fire_wind_height > zr)then !
1619 do k=kds,kdmax
1620 ht = hgtv(i,k,j) ! height of this u point above the ground
1621 if( .not. ht < fire_wind_height) then ! found layer k this point is in
1622 loght = log(ht)
1623 if(k.eq.kds)then ! first layer, log linear interpolation from 0 at zr
1624 logz0 = log(zr)
1625 va(i,j)= v(i,k,j)*(logfwh-logz0)/(loght-logz0)
1626 else ! log linear interpolation
1627 loglast=log(hgtv(i,k-1,j))
1628 va(i,j)= v(i,k-1,j) + (v(i,k,j) - v(i,k-1,j)) * ( logfwh - loglast) / (loght - loglast)
1629 endif
1630 goto 11
1631 endif
1632 if(k.eq.kdmax)then ! last layer, still not high enough
1633 va(i,j)=v(i,k,j)
1634 endif
1635 enddo
1636 11 continue
1637 else ! roughness higher than the fire wind height
1638 va(i,j)=0.
1639 endif
1640 enddo
1641 enddo
1643 #ifdef DEBUG_OUT
1644 ! store the output for diagnostics
1645 do j = jts,jte1
1646 do i = its,ite1
1647 uah(i,j)=ua(i,j)
1648 vah(i,j)=va(i,j)
1649 enddo
1650 enddo
1652 call write_array_m(its,ite1,jts,jte,ims,ime,jms,jme,uah,'uah_n',id) ! no reflection
1653 call write_array_m(its,ite,jts,jte1,ims,ime,jms,jme,vah,'vah_n',id)
1654 #endif
1656 !**********************************************************
1657 !* *
1658 !* interpolate ua,va vertically to the fire mesh *
1659 !* *
1660 !**********************************************************
1663 ips1 = ifval(ips.eq.ids,ips+1,ips)
1664 call continue_at_boundary(1,1,0., & ! x direction
1665 TDIMS, &! memory dims atm grid tile
1666 ids+1,ide,jds,jde, & ! domain dims - where u defined
1667 ips1,ipe,jps,jpe, & ! patch dims
1668 itsu,iteu,jtsu,jteu, & ! tile dims - in nonextended direction one beyond if at patch boundary but not domain
1669 itsou,iteou,jtsou,jteou, & ! out, where set
1670 ua) ! array
1672 jps1 = ifval(jps.eq.jds,jps+1,jps)
1673 call continue_at_boundary(1,1,0., & ! y direction
1674 TDIMS, & ! memory dims atm grid tile
1675 ids,ide,jds+1,jde, & ! domain dims - where v wind defined
1676 ips,ipe,jps1,jpe, & ! patch dims
1677 itsv,itev,jtsv,jtev, & ! tile dims
1678 itsov,iteov,jtsov,jteov, & ! where set
1679 va) ! array
1681 ! store the output for diagnostics
1682 do j = jts,jte1
1683 do i = its,ite1
1684 uah(i,j)=ua(i,j)
1685 vah(i,j)=va(i,j)
1686 enddo
1687 enddo
1689 #ifdef DEBUG_OUT
1690 call write_array_m(itsou,iteou,jtsou,jteou,TDIMS,ua,'ua',id)
1691 call write_array_m(itsov,iteov,jtsov,jteov,TDIMS,va,'va',id)
1692 #endif
1694 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1695 ! don't have all values valid, don't check
1696 write(msg,12)'atm mesh wind U at',fire_wind_height,' m'
1697 call print_2d_stats(itsou,iteou,jtsou,jteou,TDIMS,ua,msg)
1698 write(msg,12)'atm mesh wind V at',fire_wind_height,' m'
1699 call print_2d_stats(itsov,iteov,jtsov,jteov,TDIMS,va,msg)
1700 12 format(a,f6.2,a)
1701 call print_2d_stats(its,ite1,jts,jte,ims,ime,jms,jme,uah,'UAH')
1702 call print_2d_stats(its,ite,jts,jte1,ims,ime,jms,jme,vah,'VAH')
1703 !call write_array_m(its,ite1,jts,jte,ims,ime,jms,jme,uah,'uah',id)
1704 !call write_array_m(its,ite,jts,jte1,ims,ime,jms,jme,vah,'vah',id)
1705 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1707 ! ---------------
1708 ! | F | F | F | F | Example of atmospheric and fire grid with
1709 ! |-------|-------| ir=jr=4.
1710 ! | F | F | F | F | Winds are given at the midpoints of the sides of the atmosphere grid,
1711 ! ua------z-------| interpolated to midpoints of the cells of the fine fire grid F.
1712 ! | F | F | F | F | This is (1,1) cell of atmosphere grid, and [*] is the (1,1) cell of the fire grid.
1713 ! |---------------| ua(1,1) <--> uf(0.5,2.5)
1714 ! | * | F | F | F | va(1,1) <--> vf(2.5,0.5)
1715 ! -------va------ za(1,1) <--> zf(2.5,2.5)
1716 !
1717 ! ^ x2
1718 ! | --------va(1,2)---------
1719 ! | | | | Example of atmospheric and fire grid with
1720 ! | | | | ir=jr=1.
1721 ! | | za,zf | Winds are given at the midpoints of the sides of the atmosphere grid,
1722 ! | ua(1,1)----uf,vf-----ua(2,1) interpolated to midpoints of the cells of the (the same) fire grid
1723 ! | | (1,1) | ua(1,1) <--> uf(0.5,1)
1724 ! | | | | va(1,1) <--> vf(1,0.5)
1725 ! | | | | za(1,1) <--> zf(1,1)
1726 ! | --------va(1,1)---------
1727 ! |--------------------> x1
1728 !
1729 ! Meshes are aligned by the lower left cell of the domain. Then in the above figure
1730 ! u = node with the ua component of the wind at (ids,jds), midpoint of side
1731 ! v = node with the va component of the wind at (ids,jds), midpoint of side
1732 ! * = fire grid node at (ifds,jfds)
1733 ! z = node with height, midpoint of cell
1734 !
1735 ! ua(ids,jds)=uf(ifds-0.5,jfds+jr*0.5-0.5) = uf(ifds-0.5,jfds+(jr-1)*0.5)
1736 ! va(ids,jds)=vf(ifds+ir*0.5-0.5,jfds-0.5) = vf(ifds+(ir-1)*0.5,jfds-0.5)
1737 ! 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)
1739 ! ifts1=max(snode(ifts,ifps,-1),ifds) ! go 1 beyond patch boundary but not at domain boundary
1740 ! ifte1=min(snode(ifte,ifpe,+1),ifde)
1741 ! jfts1=max(snode(jfts,jfps,-1),jfds)
1742 ! jfte1=min(snode(jfte,jfpe,+1),jfde)
1744 call interpolate_2d( &
1745 TDIMS, & ! memory dims atm grid tile
1746 itsou,iteou,jtsou,jteou,& ! where set
1747 ifms,ifme,jfms,jfme, & ! array dims fire grid
1748 ifts,ifte,jfts,jfte,& ! dimensions on the fire grid to interpolate to
1749 ir,jr, & ! refinement ratio
1750 real(ids),real(jds),ifds-0.5,jfds+(jr-1)*0.5, & ! line up by lower left corner of domain
1751 ua, & ! in atm grid
1752 uf) ! out fire grid
1754 call interpolate_2d( &
1755 TDIMS, & ! memory dims atm grid tile
1756 itsov,iteov,jtsov,jteov,& ! where set
1757 ifms,ifme,jfms,jfme, & ! array dims fire grid
1758 ifts,ifte,jfts,jfte,& ! dimensions on the fire grid to interpolate to
1759 ir,jr, & ! refinement ratio
1760 real(ids),real(jds),ifds+(ir-1)*0.5,jfds-0.5, & ! line up by lower left corner of domain
1761 va, & ! in atm grid
1762 vf) ! out fire grid
1764 call print_2d_stats_vec(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,vf,'fire wind (m/s)')
1765 ! call print_2d_stats_vec(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,vf,'fire wind extended')
1766 #ifdef DEBUG_OUT
1767 call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,'uf1',id)
1768 call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,vf,'vf1',id)
1769 ! call write_array_m(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,'uf1',id)
1770 ! call write_array_m(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,vf,'vf1',id)
1771 #endif
1773 return
1775 end subroutine interpolate_atm2fire
1777 subroutine apply_windrf( &
1778 ifms,ifme,jfms,jfme, &
1779 ifts,ifte,jfts,jfte, &
1780 nfuel_cat,uf,vf)
1781 integer:: &
1782 ifms, ifme, jfms, jfme, & ! fire memory bounds
1783 ifts, ifte, jfts, jfte ! fire tile bounds
1784 real,intent(in),dimension(ifms:ifme,jfms:jfme)::nfuel_cat
1785 real,intent(inout),dimension(ifms:ifme,jfms:jfme)::uf,vf
1786 !*** local
1787 integer::i,j,k
1788 !*** executable
1790 do j = jfts,jfte
1791 do i = ifts,ifte
1792 k=int( nfuel_cat(i,j) )
1793 if(k.lt.no_fuel_cat)then
1794 uf(i,j)=uf(i,j)*windrf(k)
1795 vf(i,j)=vf(i,j)*windrf(k)
1796 else
1797 uf(i,j)=0.
1798 vf(i,j)=0.
1799 endif
1800 enddo
1801 enddo
1803 end subroutine apply_windrf
1805 !
1806 !***
1807 !
1809 subroutine setup_wind_log_interpolation( &
1810 ids,ide, jds,jde, & ! atm grid dimensions
1811 ims,ime, jms,jme, &
1812 ips,ipe,jps,jpe, &
1813 its,ite,jts,jte, &
1814 ifds, ifde, jfds, jfde, & ! fire grid dimensions
1815 ifms, ifme, jfms, jfme, &
1816 ifts, ifte, jfts, jfte, &
1817 ir,jr, & ! atm/fire grid ratio
1818 z0, & ! atm grid arrays in
1819 nfuel_cat, & ! fire array in
1820 fz0,fwh) ! fire arrays out
1821 !*** arguments
1822 integer, intent(in):: &
1823 ids,ide, jds,jde, & ! atm domain bounds
1824 ims,ime, jms,jme, & ! atm memory bounds
1825 ips,ipe,jps,jpe, &
1826 its,ite,jts,jte, & ! atm tile bounds
1827 ifds, ifde, jfds, jfde, & ! fire domain bounds
1828 ifms, ifme, jfms, jfme, & ! fire memory bounds
1829 ifts, ifte, jfts, jfte, & ! fire tile bounds
1830 ir,jr ! atm/fire grid refinement ratio
1831 real,intent(in),dimension(ims:ime,jms:jme)::z0 ! landuse roughness length
1832 real,intent(in),dimension(ifms:ifme,jfms:jfme)::nfuel_cat ! fuel category
1833 real,intent(out),dimension(ifms:ifme,jfms:jfme)::&
1834 fz0, & ! roughness height
1835 fwh ! height to read the wind from
1836 !*** local
1837 integer::i,j,ii,jj,k,id=0
1838 character(len=128)::msg
1839 real::r
1840 !*** executable
1842 if(.not.have_fuel_cats)call crash('setup_wind_log_interpolation: fuel categories not yet set')
1844 select case(fire_wind_log_interp)
1846 case(1)
1847 call message('fire_wind_log_interp=1: log interpolation on fire mesh, roughness and wind height from fuel categories')
1848 do j=jfts,jfte
1849 do i=ifts,ifte
1850 k = int(nfuel_cat(i,j))
1851 if(k.ge.no_fuel_cat.and.k.le.no_fuel_cat2)then ! no fuel
1852 fz0(i,j)=-1.
1853 fwh(i,j)=-1.
1854 elseif(k < 1 .or. k > nfuelcats) then
1855 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1856 write(msg,*)'i,j,nfuel_cat,nfuelcats=',i,j,k,nfuelcats
1857 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1858 call message(msg)
1859 call crash('setup_wind_log_interpolation: fuel category out of bounds')
1860 else
1861 fz0(i,j)=fcz0(k)
1862 fwh(i,j)=fcwh(k)
1863 endif
1864 enddo
1865 enddo
1867 case(2)
1868 call message('fire_wind_log_interp=2: log interpolation on fire mesh' // &
1869 'piecewise constant roughness from landuse, constant fire_wind_height')
1870 do j=jts,jte
1871 do i=its,ite
1872 do jj=(j-1)*jr+1,(j-1)*jr+jr
1873 do ii=(i-1)*ir+1,(i-1)*ir+ir
1874 fz0(ii,jj)=z0(i,j)
1875 enddo
1876 enddo
1877 enddo
1878 enddo
1879 do j=jfts,jfte
1880 do i=ifts,ifte
1881 k = int(nfuel_cat(i,j))
1882 if(k.lt.no_fuel_cat)then ! no fuel, interpolation does not matter
1883 fwh(i,j)=fcwh(k)
1884 else
1885 fz0(i,j)=-1.
1886 fwh(i,j)=-1.
1887 endif
1888 enddo
1889 enddo
1891 case(3)
1892 call message('fire_wind_log_interp=3: log interpolation on fire mesh,' // &
1893 ' interpolated roughness from landuse, constant fire_wind_height')
1894 call interpolate_z2fire(id,1, & ! for debug output, <= 0 no output
1895 ids,ide, jds,jde, & ! atm grid dimensions
1896 ims,ime, jms,jme, &
1897 ips,ipe,jps,jpe, &
1898 its,ite,jts,jte, &
1899 ifds, ifde, jfds, jfde, & ! fire grid dimensions
1900 ifms, ifme, jfms, jfme, &
1901 ifts,ifte,jfts,jfte, &
1902 ir,jr, & ! atm/fire grid ratio
1903 z0, & ! atm grid arrays in
1904 fz0) ! fire grid arrays out
1905 do j=jfts,jfte
1906 do i=ifts,ifte
1907 k = int(nfuel_cat(i,j))
1908 if(k.ne.no_fuel_cat)then ! no fuel, interpolation does not matter
1909 fwh(i,j)=fcwh(k)
1910 else
1911 fz0(i,j)=-1.
1912 fwh(i,j)=-1.
1913 endif
1914 enddo
1915 enddo
1917 case(4)
1918 call message('fire_wind_log_interp=4: log interpolation on atmospheric' // &
1919 ' mesh, roughness from landuse, constant fire_wind_height')
1920 return
1922 case default
1923 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1924 write(msg,*)'setup_wind_log_interpolation: invalid fire_wind_log_interp=',fire_wind_log_interp
1925 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1926 call crash('msg')
1928 end select
1930 select case(fire_use_windrf)
1932 case(0)
1933 call message('setup_wind_log_interpolation: not using wind reduction factors')
1935 case(1)
1936 call message('setup_wind_log_interpolation: multiplying wind by reduction factors')
1938 case(2)
1939 call message('setup_wind_log_interpolation: resetting wind interpolation height from wind reduction factors')
1940 do j=jfts,jfte
1941 do i=ifts,ifte
1942 k = int(nfuel_cat(i,j))
1943 if(k.ne.no_fuel_cat)then
1944 fwh(i,j) = fz0(i,j) ** (1.-windrf(k)) * fire_wind_height ** windrf(k) ! GMD paper eq. (26)
1946 if (.not. fz0(i,j) > 0. .or. .not. fwh(i,j) > fz0(i,j))then
1947 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1948 write(msg,*)'category ',k,'windrf=',windrf(k),' fire_wind_height=',fire_wind_height
1949 call message(msg,level=-1)
1950 write(msg,*)'i=',i,' j=',j,' fz0(i,j)=',fz0(i,j),' fwh(i,j)=',fwh(i,j)
1951 call message(msg,level=-1)
1952 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1953 call crash('setup_wind_log_interpolation: must have fwh > fz0 > 0')
1954 endif
1956 endif
1957 enddo
1958 enddo
1960 case(3)
1961 if(fire_wind_log_interp.eq.2.or.fire_wind_log_interp.eq.3)then
1962 call message('setup_wind_log_interpolation: adjusting wind interpolation height for LANDUSE roughness height')
1963 do j=jfts,jfte
1964 do i=ifts,ifte
1965 k = int(nfuel_cat(i,j))
1966 if(k.lt.no_fuel_cat)then
1967 r = log(fcwh(k)/fcz0(k))/log(fire_wind_height/fcz0(k))! GMD paper eq. (25)
1968 fwh(i,j) = fz0(i,j) ** (1.-r) * fire_wind_height ** r ! GMD paper eq. (26)
1970 if (.not. fz0(i,j) > 0. .or. .not. fwh(i,j) > fz0(i,j))then
1971 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1972 write(msg,*)'category ',k, 'roughness ',fcz0(k),' midflame height ',fcwh(k),' fire_wind_height=',fire_wind_height
1973 call message(msg,level=-1)
1974 write(msg,*)'computed wind reduction factor ',r
1975 call message(msg,level=-1)
1976 write(msg,*)'i=',i,' j=',j,' fz0(i,j)=',fz0(i,j),' fwh(i,j)=',fwh(i,j)
1977 call message(msg,level=-1)
1978 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1979 call crash('setup_wind_log_interpolation: must have fwh > fz0 > 0')
1980 endif
1982 endif
1983 enddo
1984 enddo
1985 else
1986 call message('setup_wind_log_interpolation: not using wind reduction factors')
1987 endif
1989 case default
1990 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1991 write(msg,*)'setup_wind_log_interpolation: invalid fire_use_windrf=',fire_use_windrf
1992 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1993 call crash('msg')
1995 end select
1997 ! consistency check
1998 do j=jfts,jfte
1999 do i=ifts,ifte
2000 k = int(nfuel_cat(i,j))
2001 if(k.lt.no_fuel_cat)then
2002 if (.not. fz0(i,j) > 0. .or. .not. fwh(i,j) > fz0(i,j))then
2003 !$OMP CRITICAL(SFIRE_ATM_CRIT)
2004 write(msg,*)'i=',i,' j=',j,' fz0(i,j)=',fz0(i,j),' fwh(i,j)=',fwh(i,j)
2005 call message(msg,level=-1)
2006 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
2007 call crash('setup_wind_log_interpolation: must have fwh > fz0 > 0')
2008 endif
2009 else
2010 if(.not.fwh(i,j)<0.)then
2011 !$OMP CRITICAL(SFIRE_ATM_CRIT)
2012 write(msg,*)'i=',i,' j=',j,' fz0(i,j)=',fz0(i,j),' fwh(i,j)=',fwh(i,j)
2013 call message(msg,level=-1)
2014 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
2015 call crash('setup_wind_log_interpolation: no fuel must be signalled by fwh<0')
2016 endif
2017 endif
2018 enddo
2019 enddo
2021 have_wind_log_interpolation = .true.
2023 call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fz0,'setup_wind_log:fz0')
2024 call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fz0,'setup_wind_log:fwh')
2026 end subroutine setup_wind_log_interpolation
2028 !
2029 !***
2030 !
2032 subroutine interpolate_wind2fire_height(id, & ! to identify debugging prints and files if needed
2033 ids,ide, kds,kde, jds,jde, & ! atm grid dimensions
2034 ims,ime, kms,kme, jms,jme, &
2035 ips,ipe,jps,jpe, &
2036 its,ite,jts,jte, &
2037 ifds, ifde, jfds, jfde, & ! fire grid dimensions
2038 ifms, ifme, jfms, jfme, &
2039 ifps, ifpe, jfps, jfpe, & ! fire patch bounds
2040 ifts,ifte,jfts,jfte, &
2041 ir,jr, & ! atm/fire grid ratio
2042 u_frame, v_frame, & ! velocity frame correction
2043 u,v,ph,phb, & ! input atmospheric arrays
2044 fz0,fwh, & ! input fire arrays
2045 uf,vf) ! output fire arrays
2047 implicit none
2048 !*** arguments
2049 integer, intent(in):: id, & ! debug identification
2050 ids,ide, kds,kde, jds,jde, & ! atm domain bounds
2051 ims,ime, kms,kme, jms,jme, & ! atm memory bounds
2052 ips,ipe,jps,jpe, &
2053 its,ite,jts,jte, & ! atm tile bounds
2054 ifds, ifde, jfds, jfde, & ! fire domain bounds
2055 ifms, ifme, jfms, jfme, & ! fire memory bounds
2056 ifps, ifpe, jfps, jfpe, & ! fire patch bounds
2057 ifts,ifte,jfts,jfte, & ! fire tile bounds
2058 ir,jr ! atm/fire grid refinement ratio
2059 real, intent(in):: u_frame, v_frame ! velocity frame correction
2060 real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::&
2061 u,v, & ! atm wind velocity, staggered
2062 ph,phb ! geopotential
2063 real,intent(in),dimension(ifms:ifme,jfms:jfme)::&
2064 fz0, & ! roughness height
2065 fwh ! height to read the wind from
2066 real,intent(out),dimension(ifms:ifme,jfms:jfme)::&
2067 uf, & ! atm wind at fire_wind_height, diagnostics
2068 vf !
2070 !*** local
2071 integer:: i,j,k,jcb,jcm,icb,icm,kdmax,kmin,kmax
2072 integer::itst,itet,jtst,jtet
2073 integer::iftst,iftet,jftst,jftet
2074 real:: wjcb,wjcm,wicb,wicm,ht,i_g,loght,zr,ht_last,logwh,wh,loght_last,uk,vk,uk1,vk1,z0,logz0
2075 real, dimension (its-1:ite+1,kds:kde,jts-1:jte+1):: z
2076 character(len=128)::msg
2078 !*** executable
2080 if(.not. have_wind_log_interpolation) call crash('interpolate_wind2fire_height: wind_log_interpolation must be set up first')
2082 ! print *,'interpolate_wind2fire_height start, id=',id
2084 kdmax=kde-1 ! max layer to use
2086 ! find the altitude of atm cell centers
2088 ! index bounds for cell centers - need to go one beyond if at end of tile but not domain
2089 itst=ifval(ids.eq.its,its,its-1)
2090 itet=ifval(ide.eq.ite,ite,ite+1)
2091 jtst=ifval(jds.ge.jts,jts,jts-1)
2092 jtet=ifval(jde.eq.jte,jte,jte+1)
2094 ! print *,'its, ite, jts, jte =',its, ite, jts, jte
2095 ! print *,'itst, itet, jtst, jtet=',itst, itet, jtst, jtet
2097 ! get altitudes
2098 i_g = 1./g
2099 do j = jtst,jtet
2100 do k=kds,kdmax+1
2101 do i = itst,itet
2102 z(i,k,j) = (ph(i,k,j)+phb(i,k,j))*i_g ! altitude of the bottom w-point
2104 ! print *,'i,k,j=',i,k,j,'ph, phb, z=',ph(i,k,j),phb(i,k,j),z(i,k,j)
2106 enddo
2107 enddo
2108 do k=kds,kdmax
2109 do i = itst,itet
2110 z(i,k,j) = (z(i,k,j)+z(i,k+1,j))*0.5 - z(i,kds,j) ! height of the cell center
2111 enddo
2112 enddo
2113 enddo
2115 ! index bounds for fire mesh cell centers
2116 ! to prevent setting values from uninitialized memory
2117 iftst=ifval(ifds.eq.ifts,ifts+ir/2,ifts)
2118 iftet=ifval(ifde.eq.ifte,ifte-ir/2,ifte)
2119 jftst=ifval(jfds.ge.jfts,jfts+jr/2,jfts)
2120 jftet=ifval(jfde.eq.jfte,jfte-jr/2,jfte)
2122 ! print *,'iftst, iftet, jftst, jftet=',iftst, iftet, jftst, jftet
2124 ! zero out first, to prevent unitialized values on strips along domain boundaries
2125 ! it would be faster but longer code to do just cleanup loop on the strips
2126 do j = jfts,jfte
2127 do i = ifts,ifte
2128 uf(i,j)=0.
2129 vf(i,j)=0.
2130 enddo
2131 enddo
2133 ! vertical and horizontal interpolation
2135 kmin=kde ! init stats
2136 kmax=kds
2138 loop_j: do j = jftst,jftet
2139 call coarse(j,jr,-2,jcb,wjcb) ! get interpolation coefficients from the boundary
2140 call coarse(j,jr,ir,jcm,wjcm) ! get interpolation coefficients from the midpoint
2141 loop_i: do i = iftst,iftet
2142 call coarse(i,ir,-2,icb,wicb) ! get interpolation coefficients from the boundary
2143 call coarse(i,ir,ir,icm,wicm) ! get interpolation coefficients from the midpoint
2144 z0 = fz0(i,j) ! roughness length
2145 wh = fwh(i,j) ! wind height
2147 ! print *,'i=',i,' j=',j,' icb=',icb,' jcb=',jcb,' z0=',z0,' wh=',wh
2150 if( wh > z0 .and. z0 > 0)then
2152 ht_last=z0 ! initialize starting height of this layer
2153 loop_k: do k=kds,kdmax ! search for layer k such that ht(k-1)<=wh<ht(k), ht(0)=z0
2154 ! interpolate height from atmospheric cell midpoints
2155 ht=interpolate_h(its-1,ite+1,kds,kde,jts-1,jts+1,icm,k,jcm,wicm,wjcm,z)
2157 ! print *,'i=',i,' j=',j,'k=',k,' ht=',ht
2159 if(.not. ht < wh) exit loop_k ! found layer k this point is in
2160 ht_last = ht
2161 enddo loop_k
2162 if(k .gt. kdmax) then
2163 goto 91 ! run out of vertical levels, this must be wrong
2164 endif
2165 kmin=min(k,kmin)
2166 kmax=max(k,kmax)
2167 ! found layer k, ht_last < wh <= ht
2168 logz0 = log(z0)
2169 logwh= log(wh)
2170 loght_last = log(ht_last)
2171 loght = log(ht)
2173 ! interpolate u at level k from staggered coarse grid: boundary in i, midpoints in j
2174 uk=interpolate_h(ims,ime,kms,kme,jms,jme,icb,k,jcm,wicb,wjcm,u) - u_frame
2176 ! print *,'i=',i,' j=',j,' uk=',uk
2178 ! interpolate v at level k from staggered coarse grid: midpoints in i, boundary in j
2179 vk=interpolate_h(ims,ime,kms,kme,jms,jme,icm,k,jcb,wicm,wjcb,v) - v_frame
2181 ! print *,'i=',i,' j=',j,' vk=',vk
2183 if(k.gt.kds)then ! interpolate u,v horizontaly at the previous layer, k-1
2184 ! interpolate u at level k-1 from staggered coarse grid: boundary in i, midpoints in j
2185 uk1=interpolate_h(ims,ime,kms,kme,jms,jme,icb,k-1,jcm,wicb,wjcm,u)
2187 ! print *,'i=',i,' j=',j,' uk1=',uk1
2189 ! interpolate v at level k-1 from staggered coarse grid: midpoints in i, boundary in j
2190 vk1=interpolate_h(ims,ime,kms,kme,jms,jme,icm,k-1,jcb,wicm,wjcb,v)
2192 ! print *,'i=',i,' j=',j,' uk1=',uk1
2194 else ! there is no previous layer, wind at roughness is assumed 0
2195 uk1=0.
2196 vk1=0.
2197 endif
2199 ! log interpolate the wind to wh
2200 uf(i,j)= uk1 + (uk - uk1) * ( logwh - loght_last) / (loght - loght_last)
2201 vf(i,j)= vk1 + (vk - vk1) * ( logwh - loght_last) / (loght - loght_last)
2203 ! print *,'i=',i,' j=',j,' uk=',uk,' vk=',vk,' uk1=',uk1,' vk1=',vk1,' uf=',uf(i,j),' vf=',vf(i,j)
2205 else
2206 ! no fuel, no wind
2207 uf(i,j) = 0.
2208 vf(i,j) = 0.
2210 endif
2212 enddo loop_i
2213 enddo loop_j
2215 ! print *,'interpolate_wind2fire_height complete, id=',id
2217 !$OMP CRITICAL(SFIRE_ATM_CRIT)
2218 write(msg,*)'wind interpolated from layers',kmin,' to ',kmax
2219 call message(msg)
2220 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
2222 return
2224 91 call crash('interpolate_wind2fire_height: fire wind height too large, increase kdmax or atm height')
2225 92 continue
2226 !$OMP CRITICAL(SFIRE_ATM_CRIT)
2227 write(msg,*)'fz0(',i,j,')=',fz0(i,j),'fwh(',i,j,')=',fwh(i,j)
2228 call message(msg)
2229 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
2230 call crash('interpolate_wind2fire_height: must have fire wind height > roughness height > 0')
2232 contains
2234 real function interpolate_h(ims,ime,kms,kme,jms,jme,ic,kc,jc,wic,wjc,a)
2235 !*** purpose: bilinear interpolation from a(ic:ic+1,k,jc:jc+1) with weights wicm wjcm
2236 implicit none
2237 !*** arguments
2238 integer, intent(in)::ims,ime,jms,kms,kme,jme,ic,kc,jc
2239 real, intent(in)::wic,wjc,a(ims:ime,kms:kme,jms:jme)
2240 !*** executable
2241 interpolate_h = &
2242 a(ic,kc,jc) *wic *wjc + &
2243 a(ic,kc,jc+1) *wic *(1.-wjc) + &
2244 a(ic+1,kc,jc) *(1.-wic)*wjc + &
2245 a(ic+1,kc,jc+1)*(1.-wic)*(1.-wjc)
2246 end function interpolate_h
2249 subroutine coarse(ix,nr,ia,ic,w)
2250 !*** find coarse mesh index and interpolation weight
2251 !*** arguments
2252 implicit none
2253 integer, intent(in)::ix,nr,ia
2254 integer, intent(out)::ic
2255 real, intent(out)::w
2256 !*** description
2257 ! given fine mesh with nr cells for each coarse mesh cell and such that
2258 ! coarse mesh node 1 is aligned with the fine mesh at (na+1)/2
2259 ! for fine mesh node ix find its coarse mesh coordinate c and return
2260 ! ic=floor(c), the index of the nearest coarse mesh node below
2261 ! w =1-(c-ic), the interpolation weight from coarse mesh node ic to fine mesh node ix
2262 !
2263 ! Intended use:
2264 ! fine mesh nodes are always at the middle of their cells
2265 !
2266 ! the alignment when the coarse nodes are on the boundary of coarse cells:
2267 ! |---1---|---2---|.......|--nr---| fine mesh
2268 ! 1-------------------------------2 coarse mesh
2269 ! ia = -2 because coarse node 1 is aligned with the fine mesh at -1/2 = (-2 + 1)/2
2270 !
2271 ! the alignment when the coarse node is at the midpoint of coarse cell:
2272 ! |---1---|---2---|---3---|---4---| fine mesh, here nr=4
2273 ! |---------------1---------------| coarse mesh
2274 ! ia = nr because coarse node 1 is aligned with the fine mesh at (nr + 1)/2
2275 ! here, (4 + 1)/2 = 2.5
2276 !
2278 !*** local
2279 real:: c,a
2281 !*** executable
2283 a = (ia + 1)*0.5 ! the location on the fine grid where coarse node 1 is aligned
2284 c = 1 + (ix - a)/nr ! coarse mesh coordinate of ix
2285 ic= floor(c) ! nearest coarse node to the left
2286 w = (1 + ic) - c ! interpolation weight, 1-(c-ic)
2288 end subroutine coarse
2290 end subroutine interpolate_wind2fire_height
2292 end module module_fr_sfire_atm