wrfv2_fire/phys/module_fr_sfire_atm.F

   1 !WRF:MEDIATION_LAYER:FIRE_MODEL
   2
   3 !*** Jan Mandel August 2007 - February 2010
   4 !*** email: Jan.Mandel@gmail.com
   5
   6 ! Routines dealing with the atmosphere
   7
   8 module module_fr_sfire_atm
   9
  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
  14
  15
  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
  19
  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
  84
  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
  94
  95 implicit none
  96
  97 #ifndef WRF_CHEM
  98 integer, parameter, private:: num_chem=0
  99 #endif
 100
 101 logical, save :: have_wind_log_interpolation = .false.  ! status
 102
 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.
 173
 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
 182
 183 real, parameter:: mw_air=28.97   ! molecular weight of air g/mol
 184
 185 ! should be declared in the registry and stored in the state in future because of restarts
 186
 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
 193
 194 contains
 195
 196 ! chem arrays are chem tracer
 197 ! indices p_species are generated in inc/scalar_indices.inc and included in frame/module_configure.F
 198
 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
 276
 277 !$  if (OMP_GET_THREAD_NUM() .ne. 0)then
 278 !$     call crash('read_emissions_table: must be called from master thread')
 279 !$  endif
 280
 281 IF ( wrf_dm_on_monitor() ) THEN
 282     ! we are the master task
 283
 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)
 368
 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
 386
 387 emis_read=.true.
 388
 389 end subroutine read_emissions_table
 390
 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
 400
 401
 402 implicit none
 403
 404 !*** purpose
 405 ! average fire emissions from fire mesh to coarser atmosphereic mesh and add to chemistry arrays
 406
 407 !*** arguments
 408 ! the dimensions are in cells, not nodes!
 409
 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)
 423
 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
 428
 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
 438
 439 integer, parameter:: tracer_np=8
 440 integer:: tracer_pointers(tracer_np)
 441 character(len=8)::tracer_names(tracer_np)
 442
 443 !*** executable
 444
 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)
 448
 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)
 452
 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 /)
 515
 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     ' /)
 577
 578
 579 call check_pointers('chem',chem,chem_names,chem_pointers)
 580 #endif
 581
 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  /)
 591
 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  ' /)
 601
 602 call check_pointers('tracer',tracer,tracer_names,tracer_pointers)
 603
 604 ! compute mesh sizes
 605 isz1 = ite-its+1
 606 jsz1 = jte-jts+1
 607 isz2 = ifte-ifts+1
 608 jsz2 = jfte-jtfs+1
 609
 610
 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
 616
 617 ! compute mesh ratios
 618 ir=isz2/isz1
 619 jr=jsz2/jsz1
 620
 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
 625
 626 avgw = 1.0/(ir*jr)     ! averaging weight = 1/number of fire cells per atm cell
 627
 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
 636
 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)
 640
 641 k1 = kts
 642 write(msg,'(a,i3)')'Fire emissions inserted into atmosphere level',k1
 643 call message(msg,level=msglevel)
 644
 645 #ifdef WRF_CHEM
 646 if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
 647     ! sum ground concentrations and check for nans
 648     a_chem=0.0
 649     g_chem=0.0
 650     errors=0
 651     do j=jts,jte
 652        do k=1,chem_np
 653           k_p=chem_pointers(k)
 654           do i=its,ite
 655              if(chem(i,k1,j,k_p ).ne.chem(i,k1,j,k_p ))errors=errors+1
 656              a_chem(k_p)=a_chem(k_p)+chem(i,k1,j,k_p )
 657           enddo
 658        enddo
 659     enddo
 660     if(errors>0)call crash('NaN before chem update')
 661     call wrf_dm_sum_reals(a_chem,g_chem)
 662     call message('Layer1 raw sums before adding fire emissions',level=msglevel)
 663     do i=1,chem_np,line
 664         m=min(i+line-1,chem_np)
 665         write(msg,80)'Emissions ',(trim(chem_names(j)),j=i,m)
 666         call message(msg,level=msglevel)
 667         write(msg,81)'Layer1 beg',(g_chem(chem_pointers(j)),j=i,m)
 668         call message(msg,level=msglevel)
 669         call message(' ',level=msglevel)
 670     enddo
 671 endif
 672 #endif
 673
 674 do j=max(jds+1,jts),min(jte,jde-1)          ! safe distance from domain boundary
 675     jbase=jtfs+jr*(j-jts)                   ! indexing
 676     do i=max(ids+1,its),min(ite,ide-1)
 677        ibase=ifts+ir*(i-its)                ! indexing
 678        !air = 1e6*mw_air/rho(i,kds,j)    ! 1e6*mw_air/air density
 679        !vol = avgw/dz8w(i,kds,j)             ! averaging volume factor / 1st layer depth
 680
 681        do joff=0,jr-1
 682            j_f=joff+jbase
 683            do ioff=0,ir-1
 684                i_f=ioff+ibase
 685
 686
 687                !*** fire cell (i_f,j_f) contributes to atmosphere cell (i,j) at ground level
 688
 689                fuel_burnt = fgip(i_f,j_f) * fuel_frac_burnt(i_f,j_f) ! kg/m^2
 690                cat = nfuel_cat(i_f, j_f)                       ! usually 1 to 13
 691
 692                if(cat.lt.no_fuel_cat)t_fuel(cat)=t_fuel(cat) + fuel_burnt
 693
 694
 695                !*** chem compounds emissions given in g/kg
 696
 697                ! fuel_burnt kg/m^2 *  table g/kg ->  ppmv = 1e6*mol/mol in 1st layer
 698 !AK rho is in kg/m3 so the conversion factor must be scaled by a 1000 to match
 699 !emissions in grams
 700 !                conv = avgw*1e6*mw_air/(rho(i,k1,j)*dz8w(i,k1,j))
 701                  conv = avgw*1e3*mw_air/(rho(i,k1,j)*dz8w(i,k1,j))
 702 #ifdef WRF_CHEM
 703
 704                emis=co  (cat)*fuel_burnt                                    ! emission from fire cell in g/m^2
 705                t_chem(p_co) = t_chem(p_co) + emis                           ! add to total
 706                ! if(isnan(chem(i,k1,j,p_co  )))call crash('NaN before')
 707                chem(i,k1,j,p_co  )=chem(i,k1,j,p_co  ) + emis*conv*imw_co   ! add to chem
 708                ! if(isnan(chem(i,k1,j,p_co  )))call crash('NaN after')
 709
 710                emis=ch4  (cat)*fuel_burnt                                    ! emission from fire cell in g/m^2
 711                t_chem(p_ch4) = t_chem(p_ch4) + emis                           ! add to total
 712                chem(i,k1,j,p_ch4  )=chem(i,k1,j,p_ch4  ) + emis*conv*imw_ch4   ! add to chem
 713
 714                emis=h2  (cat)*fuel_burnt                                    ! emission from fire cell in g/m^2
 715                t_chem(p_h2) = t_chem(p_h2) + emis                           ! add to total
 716                chem(i,k1,j,p_h2  )=chem(i,k1,j,p_h2  ) + emis*conv*imw_h2   ! add to chem
 717
 718                emis=no  (cat)*fuel_burnt                                    ! emission from fire cell in g/m^2
 719                t_chem(p_no) = t_chem(p_no) + emis                           ! add to total
 720                chem(i,k1,j,p_no  )=chem(i,k1,j,p_no  ) + emis*conv*imw_no   ! add to chem
 721
 722                emis=no2 (cat)*fuel_burnt                                    ! emission from fire cell in g/m^2
 723                t_chem(p_no2) = t_chem(p_no2) + emis                         ! add to total
 724                chem(i,k1,j,p_no2 )=chem(i,k1,j,p_no2 ) + emis*conv*imw_no2  ! add to chem
 725
 726                emis=so2 (cat)*fuel_burnt                                    ! emission from fire cell in g/m^2
 727                t_chem(p_so2) = t_chem(p_so2) + emis                         ! ads to total
 728                chem(i,k1,j,p_so2 )=chem(i,k1,j,p_so2 ) + emis*conv*imw_so2  ! add to chem
 729
 730                emis=nh3 (cat)*fuel_burnt                                    ! emission from fire cell in g/m^2
 731                t_chem(p_nh3) = t_chem(p_nh3) + emis                         ! add to total
 732                chem(i,k1,j,p_nh3 )=chem(i,k1,j,p_nh3 ) + emis*conv*imw_nh3  ! add to chem
 733
 734
 735                !*** other emissions already given in mol/kg
 736                ! fuel_burnt kg/m^2 *  table mol/kg -> ppmv = 1e6*mol/mol in 1st layer dry air in 1st layer
 737
 738                ! same conversion factor but we will not divide by the molecular weight of the compound
 739                ! conv = avgw*1e6*mw_air/(rho(i,kds,j)*dz8w(i,kds,j))
 740
 741                emis=ald (cat)*fuel_burnt                                    ! emission from fire cell
 742                t_chem(p_ald) = t_chem(p_ald) + emis                         ! add to total
 743                chem(i,k1,j,p_ald )=chem(i,k1,j,p_ald )+emis*conv
 744
 745                emis=csl (cat)*fuel_burnt                                    ! emission from fire cell
 746                t_chem(p_csl) = t_chem(p_csl) + emis                         ! add to total
 747                chem(i,k1,j,p_csl )=chem(i,k1,j,p_csl )+emis*conv
 748
 749                emis=eth (cat)*fuel_burnt                                    ! emission from fire cell
 750                t_chem(p_eth) = t_chem(p_eth) + emis                         ! add to total
 751                chem(i,k1,j,p_eth )=chem(i,k1,j,p_eth )+emis*conv
 752
 753                emis=hc3 (cat)*fuel_burnt                                    ! emission from fire cell
 754                t_chem(p_hc3) = t_chem(p_hc3) + emis                         ! add to total
 755                chem(i,k1,j,p_hc3 )=chem(i,k1,j,p_hc3 )+emis*conv
 756
 757                emis=hc5 (cat)*fuel_burnt                                    ! emission from fire cell
 758                t_chem(p_hc5) = t_chem(p_hc5) + emis                         ! add to total
 759                chem(i,k1,j,p_hc5 )=chem(i,k1,j,p_hc5 )+emis*conv
 760
 761                emis=hcho (cat)*fuel_burnt                                    ! emission from fire cell
 762                t_chem(p_hcho) = t_chem(p_hcho) + emis                         ! add to total
 763                chem(i,k1,j,p_hcho)=chem(i,k1,j,p_hcho)+emis*conv
 764
 765                emis=iso (cat)*fuel_burnt                                    ! emission from fire cell
 766                t_chem(p_iso) = t_chem(p_iso) + emis                         ! add to total
 767                chem(i,k1,j,p_iso )=chem(i,k1,j,p_iso )+emis*conv
 768
 769                emis=ket (cat)*fuel_burnt                                    ! emission from fire cell
 770                t_chem(p_ket) = t_chem(p_ket) + emis                         ! add to total
 771                chem(i,k1,j,p_ket )=chem(i,k1,j,p_ket )+emis*conv
 772
 773                emis=mgly (cat)*fuel_burnt                                    ! emission from fire cell
 774                t_chem(p_mgly) = t_chem(p_mgly) + emis                         ! add to total
 775                chem(i,k1,j,p_mgly)=chem(i,k1,j,p_mgly)+emis*conv
 776
 777                emis=ol2 (cat)*fuel_burnt                                    ! emission from fire cell
 778                t_chem(p_ol2) = t_chem(p_ol2) + emis                         ! add to total
 779                chem(i,k1,j,p_ol2 )=chem(i,k1,j,p_ol2 )+emis*conv
 780
 781                emis=olt (cat)*fuel_burnt                                    ! emission from fire cell
 782                t_chem(p_olt) = t_chem(p_olt) + emis                         ! add to total
 783                chem(i,k1,j,p_olt )=chem(i,k1,j,p_olt )+emis*conv
 784
 785                emis=oli (cat)*fuel_burnt                                    ! emission from fire cell
 786                t_chem(p_oli) = t_chem(p_oli) + emis                         ! add to total
 787                chem(i,k1,j,p_oli )=chem(i,k1,j,p_oli )+emis*conv
 788
 789                emis=ora2 (cat)*fuel_burnt                                    ! emission from fire cell
 790                t_chem(p_ora2) = t_chem(p_ora2) + emis                         ! add to total
 791                chem(i,k1,j,p_ora2)=chem(i,k1,j,p_ora2)+emis*conv
 792
 793                emis=tol (cat)*fuel_burnt                                    ! emission from fire cell
 794                t_chem(p_tol) = t_chem(p_tol) + emis                         ! add to total
 795                chem(i,k1,j,p_tol )=chem(i,k1,j,p_tol )+emis*conv
 796
 797                emis=xyl (cat)*fuel_burnt                                    ! emission from fire cell
 798                t_chem(p_xyl) = t_chem(p_xyl) + emis                         ! add to total
 799                chem(i,k1,j,p_xyl )=chem(i,k1,j,p_xyl )+emis*conv
 800
 801                emis=bigalk (cat)*fuel_burnt                                    ! emission from fire cell
 802                t_chem(p_bigalk) = t_chem(p_bigalk) + emis                         ! add to total
 803                chem(i,k1,j,p_bigalk )=chem(i,k1,j,p_bigalk )+emis*conv
 804
 805                emis=bigene (cat)*fuel_burnt                                    ! emission from fire cell
 806                t_chem(p_bigene) = t_chem(p_bigene) + emis                         ! add to total
 807                chem(i,k1,j,p_bigene )=chem(i,k1,j,p_bigene )+emis*conv
 808
 809                emis=c10h16 (cat)*fuel_burnt                                    ! emission from fire cell
 810                t_chem(p_c10h16) = t_chem(p_c10h16) + emis                         ! add to total
 811                chem(i,k1,j,p_c10h16 )=chem(i,k1,j,p_c10h16 )+emis*conv
 812
 813                emis=c2h4 (cat)*fuel_burnt                                    ! emission from fire cell
 814                t_chem(p_c2h4) = t_chem(p_c2h4) + emis                         ! add to total
 815                chem(i,k1,j,p_c2h4 )=chem(i,k1,j,p_c2h4 )+emis*conv
 816
 817                emis=c2h5oh (cat)*fuel_burnt                                    ! emission from fire cell
 818                t_chem(p_c2h5oh) = t_chem(p_c2h5oh) + emis                         ! add to total
 819                chem(i,k1,j,p_c2h5oh )=chem(i,k1,j,p_c2h5oh )+emis*conv
 820
 821                emis=c2h6 (cat)*fuel_burnt                                    ! emission from fire cell
 822                t_chem(p_c2h6) = t_chem(p_c2h6) + emis                         ! add to total
 823                chem(i,k1,j,p_c2h6 )=chem(i,k1,j,p_c2h6 )+emis*conv
 824
 825                emis=c3h6 (cat)*fuel_burnt                                    ! emission from fire cell
 826                t_chem(p_c3h6) = t_chem(p_c3h6) + emis                         ! add to total
 827                chem(i,k1,j,p_c3h6 )=chem(i,k1,j,p_c3h6 )+emis*conv
 828
 829                emis=c3h8 (cat)*fuel_burnt                                    ! emission from fire cell
 830                t_chem(p_c3h8) = t_chem(p_c3h8) + emis                         ! add to total
 831                chem(i,k1,j,p_c3h8 )=chem(i,k1,j,p_c3h8 )+emis*conv
 832
 833                emis=ch3cooh (cat)*fuel_burnt                                    ! emission from fire cell
 834                t_chem(p_ch3cooh) = t_chem(p_ch3cooh) + emis                         ! add to total
 835                chem(i,k1,j,p_ch3cooh )=chem(i,k1,j,p_ch3cooh )+emis*conv
 836
 837                emis=ch3oh (cat)*fuel_burnt                                    ! emission from fire cell
 838                t_chem(p_ch3oh) = t_chem(p_ch3oh) + emis                         ! add to total
 839                chem(i,k1,j,p_ch3oh )=chem(i,k1,j,p_ch3oh )+emis*conv
 840
 841                emis=cres (cat)*fuel_burnt                                    ! emission from fire cell
 842                t_chem(p_cres) = t_chem(p_cres) + emis                         ! add to total
 843                chem(i,k1,j,p_cres )=chem(i,k1,j,p_cres )+emis*conv
 844
 845                emis=glyald (cat)*fuel_burnt                                    ! emission from fire cell
 846                t_chem(p_glyald) = t_chem(p_glyald) + emis                         ! add to total
 847                chem(i,k1,j,p_glyald )=chem(i,k1,j,p_glyald )+emis*conv
 848
 849                emis=isopr (cat)*fuel_burnt                                    ! emission from fire cell
 850                t_chem(p_isopr) = t_chem(p_isopr) + emis                         ! add to total
 851                chem(i,k1,j,p_isopr )=chem(i,k1,j,p_isopr )+emis*conv
 852
 853                emis=macr (cat)*fuel_burnt                                    ! emission from fire cell
 854                t_chem(p_macr) = t_chem(p_macr) + emis                         ! add to total
 855                chem(i,k1,j,p_macr )=chem(i,k1,j,p_macr )+emis*conv
 856
 857                emis=mek (cat)*fuel_burnt                                    ! emission from fire cell
 858                t_chem(p_mek) = t_chem(p_mek) + emis                         ! add to total
 859                chem(i,k1,j,p_mek )=chem(i,k1,j,p_mek )+emis*conv
 860
 861                emis=mvk (cat)*fuel_burnt                                    ! emission from fire cell
 862                t_chem(p_mvk) = t_chem(p_mvk) + emis                         ! add to total
 863                chem(i,k1,j,p_mvk )=chem(i,k1,j,p_mvk )+emis*conv
 864
 865                ! aerosols
 866                ! fuel_burnt kg/m^2 *  table g/kg -> ug/kg dry air in 1st layer
 867
 868                ! see also chem/emissions_driver.F line 515
 869
 870                conv = avgw*1e6/(rho(i,k1,j)*dz8w(i,k1,j))
 871
 872                emis=p25 (cat)*fuel_burnt                                    ! emission from fire cell
 873                t_chem(p_p25) = t_chem(p_p25) + emis                         ! add to total
 874                chem(i,k1,j,p_p25 )=chem(i,k1,j,p_p25 )+emis*conv
 875
 876                emis=p25i (cat)*fuel_burnt                                    ! emission from fire cell
 877                t_chem(p_p25i) = t_chem(p_p25i) + emis                         ! add to total
 878                chem(i,k1,j,p_p25i )=chem(i,k1,j,p_p25i )+emis*conv
 879
 880                emis=p25j (cat)*fuel_burnt                                    ! emission from fire cell
 881                t_chem(p_p25j) = t_chem(p_p25j) + emis                         ! add to total
 882                chem(i,k1,j,p_p25j )=chem(i,k1,j,p_p25j )+emis*conv
 883
 884                emis=oc1  (cat)*fuel_burnt                                    ! emission from fire cell
 885                t_chem(p_oc1 ) = t_chem(p_oc1 ) + emis                         ! add to total
 886                chem(i,k1,j,p_oc1  )=chem(i,k1,j,p_oc1  )+emis*conv
 887
 888                emis=oc2  (cat)*fuel_burnt                                    ! emission from fire cell
 889                t_chem(p_oc2 ) = t_chem(p_oc2 ) + emis                         ! add to total
 890                chem(i,k1,j,p_oc2  )=chem(i,k1,j,p_oc2  )+emis*conv
 891
 892                emis=bc1  (cat)*fuel_burnt                                    ! emission from fire cell
 893                t_chem(p_bc1 ) = t_chem(p_bc1 ) + emis                         ! add to total
 894                chem(i,k1,j,p_bc1  )=chem(i,k1,j,p_bc1  )+emis*conv
 895
 896                emis=bc2  (cat)*fuel_burnt                                    ! emission from fire cell
 897                t_chem(p_bc2 ) = t_chem(p_bc2 ) + emis                         ! add to total
 898                chem(i,k1,j,p_bc2  )=chem(i,k1,j,p_bc2  )+emis*conv
 899 #endif
 900                if (num_tracer >0)then
 901
 902                      ! treat tracers exactly the same as aerosols, emissions g/kg fuel burned, tracer concentration ug/kg dry air
 903                      conv = avgw*1e6/(rho(i,k1,j)*dz8w(i,k1,j))
 904
 905                      emis=tr17_1 (cat)*fuel_burnt                                    ! emission from fire cell
 906                      t_tracer(p_tr17_1) = t_tracer(p_tr17_1) + emis                         ! add to total
 907                      tracer(i,k1,j,p_tr17_1 )=tracer(i,k1,j,p_tr17_1 )+emis*conv
 908
 909                      emis=tr17_2 (cat)*fuel_burnt                                    ! emission from fire cell
 910                      t_tracer(p_tr17_2) = t_tracer(p_tr17_2) + emis                         ! add to total
 911                      tracer(i,k1,j,p_tr17_2 )=tracer(i,k1,j,p_tr17_2 )+emis*conv
 912
 913                      emis=tr17_3 (cat)*fuel_burnt                                    ! emission from fire cell
 914                      t_tracer(p_tr17_3) = t_tracer(p_tr17_3) + emis                         ! add to total
 915                      tracer(i,k1,j,p_tr17_3 )=tracer(i,k1,j,p_tr17_3 )+emis*conv
 916
 917                      emis=tr17_4 (cat)*fuel_burnt                                    ! emission from fire cell
 918                      t_tracer(p_tr17_4) = t_tracer(p_tr17_4) + emis                         ! add to total
 919                      tracer(i,k1,j,p_tr17_4 )=tracer(i,k1,j,p_tr17_4 )+emis*conv
 920
 921                      emis=tr17_5 (cat)*fuel_burnt                                    ! emission from fire cell
 922                      t_tracer(p_tr17_5) = t_tracer(p_tr17_5) + emis                         ! add to total
 923                      tracer(i,k1,j,p_tr17_5 )=tracer(i,k1,j,p_tr17_5 )+emis*conv
 924
 925                      emis=tr17_6 (cat)*fuel_burnt                                    ! emission from fire cell
 926                      t_tracer(p_tr17_6) = t_tracer(p_tr17_6) + emis                         ! add to total
 927                      tracer(i,k1,j,p_tr17_6 )=tracer(i,k1,j,p_tr17_6 )+emis*conv
 928
 929                      emis=tr17_7 (cat)*fuel_burnt                                    ! emission from fire cell
 930                      t_tracer(p_tr17_7) = t_tracer(p_tr17_7) + emis                         ! add to total
 931                      tracer(i,k1,j,p_tr17_7 )=tracer(i,k1,j,p_tr17_7 )+emis*conv
 932
 933                      emis=tr17_8 (cat)*fuel_burnt                                    ! emission from fire cell
 934                      t_tracer(p_tr17_8) = t_tracer(p_tr17_8) + emis                         ! add to total
 935                      tracer(i,k1,j,p_tr17_8 )=tracer(i,k1,j,p_tr17_8 )+emis*conv
 936                 endif
 937            enddo
 938        enddo
 939     enddo
 940 enddo
 941
 942 #ifdef WRF_CHEM
 943 if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
 944     ! sum ground concentrations and check for nans
 945     a_chem=0.0
 946     g_chem=0.0
 947     errors=0
 948     do j=jts,jte
 949        do k=1,chem_np
 950           k_p=chem_pointers(k)
 951           do i=its,ite
 952              if(chem(i,k1,j,k_p ).ne.chem(i,k1,j,k_p ))errors=errors+1
 953              a_chem(k_p)=a_chem(k_p)+chem(i,k1,j,k_p )
 954           enddo
 955        enddo
 956     enddo
 957     if(errors>0)call crash('NaN after chem update')
 958     call wrf_dm_sum_reals(a_chem,g_chem)
 959     call message('Layer1 raw sums after adding fire emissions',level=msglevel)
 960     do i=1,chem_np,line
 961         m=min(i+line-1,chem_np)
 962         write(msg,80)'Emissions ',(trim(chem_names(j)),j=i,m)
 963         call message(msg,level=msglevel)
 964         write(msg,81)'Layer1 end',(g_chem(chem_pointers(j)),j=i,m)
 965         call message(msg,level=msglevel)
 966         call message(' ',level=msglevel)
 967     enddo
 968 endif
 969 #endif
 970
 971 if(fire_print_msg .ge. msglevel .and.printsums .gt. 0)then
 972     ! sum over processes and add to cumulative sums
 973
 974     ! fuel burned
 975     call wrf_dm_sum_reals(t_fuel,s_fuel)
 976     ! scale
 977     s_fuel = s_fuel*dx*dy
 978     ! get rates
 979     r_fuel = s_fuel/dt
 980     ! add to cumulative totals
 981     if(size(c_fuel).ne.nfuelcats)call crash('add_fire_emissions: bad size c_fuel')
 982     c_fuel = c_fuel + s_fuel
 983
 984 #ifdef WRF_CHEM
 985     call wrf_dm_sum_reals(a_chem,g_chem)
 986     ! chem
 987     call wrf_dm_sum_reals(t_chem,s_chem)
 988     s_chem = s_chem*dx*dy
 989     r_chem = s_chem/dt
 990     if(size(c_chem).ne.num_chem)call crash('add_fire_emissions: bad size c_chem')
 991     c_chem = c_chem + s_chem
 992
 993     call message('Total emissions in g or mol per the table',level=msglevel)
 994     do i=1,chem_np,line
 995         m=min(i+line-1,chem_np)
 996         write(msg,80)'Emissions ',(trim(chem_names(j)),j=i,m)
 997         call message(msg,level=msglevel)
 998         write(msg,81)'Cumulative',(c_chem(chem_pointers(j)),j=i,m)
 999         call message(msg,level=msglevel)
1000         write(msg,81)'Rate per s',(r_chem(chem_pointers(j)),j=i,m)
1001         call message(msg,level=msglevel)
1002         call message(' ',level=msglevel)
1003     enddo
1004 #endif
1005
1006     if(num_tracer >0)then
1007         ! tracer
1008         call wrf_dm_sum_reals(t_tracer,s_tracer)
1009         s_tracer = s_tracer*dx*dy
1010         r_tracer = s_tracer/dt
1011         if(size(c_tracer).ne.num_tracer)call crash('add_fire_emissions: bad size c_tracer')
1012         c_tracer = c_tracer + s_tracer
1013
1014         call message('Total emissions in g',level=msglevel)
1015         do i=1,tracer_np,line
1016             m=min(i+line-1,tracer_np)
1017             write(msg,80)'Emissions ',(trim(tracer_names(j)),j=i,m)
1018             call message(msg,level=msglevel)
1019             write(msg,81)'Cumulative',(c_tracer(tracer_pointers(j)),j=i,m)
1020             call message(msg,level=msglevel)
1021             write(msg,81)'Rate per s',(r_tracer(tracer_pointers(j)),j=i,m)
1022             call message(msg,level=msglevel)
1023             call message(' ',level=msglevel)
1024         enddo
1025     endif
1026
1027
1028     do cat=1,nfuelcats
1029         if(c_fuel(cat) > 0.)then
1030             write(msg,83)fuel_name(cat),' burned',c_fuel(cat),'kg, rate',r_fuel(cat),' kg/s'
1031             call message(msg,level=msglevel)
1032         endif
1033     enddo
1034     write(msg,83)'Total fuel',' burned',sum(c_fuel),'kg, rate',sum(r_fuel),' kg/s'
1035     call message(msg,level=msglevel)
1036
1037 endif
1038 80 format(a,8a11)
1039 81 format(a,8e11.3)
1040 83 format(a30,a,g14.4,a,g14.4,a,a)
1041
1042 return
1043
1044 9 continue
1045 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1046 write(msg,91)ifts,ifte,jtfs,jfte,ifms,ifme,jfms,jfme
1047 call message(msg,level=0)
1048 write(msg,91)its,ite,jts,jte,ims,ime,jms,jme
1049 call message(msg,level=0)
1050 write(msg,92)'input  mesh size:',isz2,jsz2
1051 call message(msg,level=0)
1052 91 format('dimensions: ',8i8)
1053 write(msg,92)'output mesh size:',isz1,jsz1
1054 call message(msg,level=0)
1055 92 format(a,2i8)
1056 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1057 call crash('add_fire_emissions: bad mesh sizes')
1058
1059 end subroutine add_fire_emissions
1060
1061 !
1062 !***
1063 !
1064
1065 subroutine check_pointers(array_name,array,pointer_names,pointers)
1066 implicit none
1067
1068 !*** arguments
1069 character(len=*)::array_name
1070 real, dimension(:,:,:,:)::array
1071 character(len=*), dimension(:)::pointer_names
1072 integer, dimension(:)::pointers
1073
1074 !*** local
1075 integer::np,i,m,j
1076 character(len=256)::msg
1077
1078 !** executable
1079 np=ubound(pointers,1)
1080
1081
1082 993 format(3a,4(1x,i3,':',i3))
1083 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1084 write(msg,993)'array ',array_name,' has dimensions ',&
1085          lbound(array,1),ubound(array,1), &
1086          lbound(array,2),ubound(array,2), &
1087          lbound(array,3),ubound(array,3), &
1088          lbound(array,4),ubound(array,4)
1089 call message(msg)
1090
1091 do i=1,np,line
1092    m=min(i+line-1,np)
1093    write(msg,'(a,8(1x,a8))')'Species',(trim(pointer_names(j)),j=i,m)
1094    call message(msg,level=msglevel)
1095    write(msg,'(a,8i9)')     'Pointer',(pointers(j),j=i,m)
1096    call message(msg,level=msglevel)
1097    call message(' ')
1098 enddo
1099
1100 if(maxval(pointers) > ubound(array,4) .or. minval(pointers) < lbound(array,4))then
1101     write(msg,'(3a)')'add_fire_emissions: a ',array_name,' pointer is out of bounds'
1102     call crash(msg)
1103 endif
1104 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1105 end subroutine check_pointers
1106
1107
1108 !
1109 !***
1110 !
1111
1112
1113 SUBROUTINE fire_tendency( &
1114     ids,ide, kds,kde, jds,jde,   & ! dimensions
1115     ims,ime, kms,kme, jms,jme,   &
1116     its,ite, kts,kte, jts,jte,   &
1117     grnhfx,grnqfx,canhfx,canqfx, & ! heat fluxes summed up to  atm grid
1118     alfg,alfc,z1can,             & ! coeffients, properties, geometry
1119     zs,z_at_w,dz8w,mu,rho,       &
1120     rthfrten,rqvfrten)             ! theta and Qv tendencies
1121
1122 ! This routine is atmospheric physics
1123 ! it does NOT go into module_fr_sfire_phys because it is not fire physics
1124
1125 ! taken from the code by Ned Patton, only order of arguments change to the convention here
1126 ! --- this routine takes fire generated heat and moisture fluxes and
1127 !     calculates their influence on the theta and water vapor
1128 ! --- note that these tendencies are valid at the Arakawa-A location
1129
1130    IMPLICIT NONE
1131
1132 ! --- incoming variables
1133
1134    INTEGER , INTENT(in) :: ids,ide, kds,kde, jds,jde, &
1135                            ims,ime, kms,kme, jms,jme, &
1136                            its,ite, kts,kte, jts,jte
1137
1138    REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: grnhfx,grnqfx  ! W/m^2
1139    REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: canhfx,canqfx  ! W/m^2
1140    REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: zs  ! topography (m abv sealvl)
1141    REAL, INTENT(in), DIMENSION( ims:ime,jms:jme ) :: mu  ! dry air mass (Pa)
1142
1143    REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: z_at_w ! m abv sealvl
1144    REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: dz8w   ! dz across w-lvl
1145    REAL, INTENT(in), DIMENSION( ims:ime,kms:kme,jms:jme ) :: rho    ! density
1146
1147    REAL, INTENT(in) :: alfg ! extinction depth ground fire heat (m)
1148    REAL, INTENT(in) :: alfc ! extinction depth crown  fire heat (m)
1149    REAL, INTENT(in) :: z1can    ! height of crown fire heat release (m)
1150
1151 ! --- outgoing variables
1152
1153    REAL, INTENT(out), DIMENSION( ims:ime,kms:kme,jms:jme ) ::   &
1154        rthfrten, & ! theta tendency from fire (in mass units)
1155        rqvfrten    ! Qv tendency from fire (in mass units)
1156 ! --- local variables
1157
1158    INTEGER :: i,j,k
1159    INTEGER :: i_st,i_en, j_st,j_en, k_st,k_en
1160
1161    REAL :: cp_i
1162    REAL :: rho_i
1163    REAL :: xlv_i
1164    REAL :: z_w
1165    REAL :: fact_g, fact_c
1166    REAL :: alfg_i, alfc_i
1167
1168    REAL, DIMENSION( its:ite,kts:kte,jts:jte ) :: hfx,qfx
1169
1170 !!   character(len=128)::msg
1171
1172
1173         do j=jts,jte
1174             do k=kts,min(kte+1,kde)
1175                do i=its,ite
1176                    rthfrten(i,k,j)=0.
1177                    rqvfrten(i,k,j)=0.
1178                enddo
1179             enddo
1180         enddo
1181
1182
1183 ! --- set some local constants
1184
1185
1186    cp_i = 1./cp     ! inverse of specific heat
1187    xlv_i = 1./xlv   ! inverse of latent heat
1188    alfg_i = 1./alfg
1189    alfc_i = 1./alfc
1190
1191 !!write(msg,'(8e11.3)')cp,cp_i,xlv,xlv_i,alfg,alfc,z1can
1192 !!call message(msg)
1193
1194    call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnhfx,'fire_tendency:grnhfx')
1195    call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnqfx,'fire_tendency:grnqfx')
1196
1197 ! --- set loop indicies : note that
1198
1199    i_st = MAX(its,ids+1)
1200    i_en = MIN(ite,ide-1)
1201    k_st = kts
1202    k_en = MIN(kte,kde-1)
1203    j_st = MAX(jts,jds+1)
1204    j_en = MIN(jte,jde-1)
1205
1206 ! --- distribute fluxes
1207
1208    DO j = j_st,j_en
1209       DO k = k_st,k_en
1210          DO i = i_st,i_en
1211
1212             ! --- set z (in meters above ground)
1213
1214             z_w = z_at_w(i,k,j) - zs(i,j) ! should be zero when k=k_st
1215
1216             ! --- heat flux
1217
1218             fact_g = cp_i * EXP( - alfg_i * z_w )
1219             IF ( z_w < z1can ) THEN
1220                fact_c = cp_i
1221             ELSE
1222                fact_c = cp_i * EXP( - alfc_i * (z_w - z1can) )
1223             END IF
1224             hfx(i,k,j) = fact_g * grnhfx(i,j) + fact_c * canhfx(i,j)
1225
1226 !!            write(msg,2)i,k,j,z_w,grnhfx(i,j),hfx(i,k,j)
1227 !!2           format('hfx:',3i4,6e11.3)
1228 !!            call message(msg)
1229
1230             ! --- vapor flux
1231
1232             fact_g = xlv_i * EXP( - alfg_i * z_w )
1233             IF (z_w < z1can) THEN
1234                fact_c = xlv_i
1235             ELSE
1236                fact_c = xlv_i * EXP( - alfc_i * (z_w - z1can) )
1237             END IF
1238             qfx(i,k,j) = fact_g * grnqfx(i,j) + fact_c * canqfx(i,j)
1239
1240 !!            if(hfx(i,k,j).ne.0. .or. qfx(i,k,j) .ne. 0.)then
1241 !!                write(msg,1)i,k,j,hfx(i,k,j),qfx(i,k,j)
1242 !!1               format('tend:',3i6,2e11.3)
1243 !!                call message(msg)
1244 !            endif
1245
1246          END DO
1247       END DO
1248    END DO
1249
1250 ! --- add flux divergence to tendencies
1251 !
1252 !   multiply by dry air mass (mu) to eliminate the need to
1253 !   call sr. calculate_phy_tend (in dyn_em/module_em.F)
1254
1255    ! print *,'fire_tendency:',i_st,i_en,j_st,j_en
1256    DO j = j_st,j_en
1257       DO k = k_st,k_en-1
1258          DO i = i_st,i_en
1259
1260             rho_i = 1./rho(i,k,j)
1261
1262             rthfrten(i,k,j) = - mu(i,j) * rho_i * (hfx(i,k+1,j)-hfx(i,k,j)) / dz8w(i,k,j)
1263             rqvfrten(i,k,j) = - mu(i,j) * rho_i * (qfx(i,k+1,j)-qfx(i,k,j)) / dz8w(i,k,j)
1264
1265             ! print *,i,j,k,rthfrten(i,k,j)
1266
1267          END DO
1268       END DO
1269    END DO
1270
1271    call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rthfrten,'fire_tendency:rthfrten')
1272    call print_3d_stats(its,ite,kts,kte,jts,jte,ims,ime,kms,kme,jms,jme,rqvfrten,'fire_tendency:rqvfrten')
1273
1274    RETURN
1275
1276 END SUBROUTINE fire_tendency
1277
1278 !
1279 !***
1280 !
1281 subroutine interpolate_atm2fire(id,               & ! for debug output, <= 0 no output
1282     ids,ide, kds,kde, jds,jde,                    & ! atm grid dimensions
1283     ims,ime, kms,kme, jms,jme,                    &
1284     ips,ipe,jps,jpe,                              &
1285     its,ite,jts,jte,                              &
1286     ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
1287     ifms, ifme, jfms, jfme,                       &
1288     ifps, ifpe, jfps, jfpe,                       & ! fire patch bounds
1289     ifts,ifte,jfts,jfte,                          &
1290     ir,jr,                                        & ! atm/fire grid ratio
1291     u_frame, v_frame,                             & ! velocity frame correction
1292     u,v,                                          & ! atm grid arrays in
1293     ph,phb,                                       &
1294     z0,zs,                                        &
1295     uah,vah,                                      &
1296     uf,vf)                                          ! fire grid arrays out
1297
1298 implicit none
1299 ! Jan Mandel, October 2010
1300 !*** purpose: interpolate winds and height
1301
1302 !*** arguments
1303 integer, intent(in)::id
1304 integer, intent(in)::                             &
1305     ids,ide, kds,kde, jds,jde,                    & ! atm domain bounds
1306     ims,ime, kms,kme, jms,jme,                    & ! atm memory bounds
1307     ips,ipe,jps,jpe,                              &
1308     its,ite,jts,jte,                              & ! atm tile bounds
1309     ifds, ifde, jfds, jfde,                       & ! fire domain bounds
1310     ifms, ifme, jfms, jfme,                       & ! fire memory bounds
1311     ifps, ifpe, jfps, jfpe,                       & ! fire patch bounds
1312     ifts,ifte,jfts,jfte,                          & ! fire tile bounds
1313     ir,jr                                         ! atm/fire grid refinement ratio
1314 real, intent(in):: u_frame, v_frame                 ! velocity frame correction
1315 real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::&
1316     u,v,                                          & ! atm wind velocity, staggered
1317     ph,phb                                          ! geopotential
1318 real,intent(in),dimension(ims:ime,jms:jme)::&
1319     z0,                                           & ! roughness height
1320     zs                                              ! terrain height
1321 real,intent(out),dimension(ims:ime,jms:jme)::&
1322     uah,                                           & ! atm wind at fire_wind_height, diagnostics
1323     vah                                              !
1324 real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
1325     uf,vf                                           ! wind velocity fire grid nodes
1326
1327
1328 !*** local
1329 character(len=256)::msg
1330 #define TDIMS its-2,ite+2,jts-2,jte+2
1331 real, dimension(its-2:ite+2,jts-2:jte+2):: ua,va   ! atm winds, interpolated over height, still staggered grid
1332 real, dimension(its-2:ite+2,kds:kde,jts-2:jte+2):: altw,altub,altvb,hgtu,hgtv ! altitudes
1333 integer:: i,j,k,ifts1,ifte1,jfts1,jfte1,ite1,jte1
1334 integer::itst,itet,jtst,jtet,itsu,iteu,jtsu,jteu,itsv,itev,jtsv,jtev
1335 integer::kdmax,its1,jts1,ips1,jps1
1336 integer::itsou,iteou,jtsou,jteou,itsov,iteov,jtsov,jteov
1337 real:: ground,loght,loglast,logz0,logfwh,ht,zr
1338 real::r_nan
1339 integer::i_nan
1340 equivalence (i_nan,r_nan)
1341
1342 !*** executable
1343
1344 ! debug init local arrays
1345 i_nan=2147483647
1346 ua=r_nan
1347 va=r_nan
1348 altw=r_nan
1349 altub=r_nan
1350 hgtu=r_nan
1351 hgtv=r_nan
1352
1353
1354 if(kds.ne.1)then
1355 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1356   write(msg,*)'WARNING: bottom index kds=',kds,' should be 1?'
1357   call message(msg)
1358 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1359 endif
1360
1361 !                            ^ j
1362 !        ------------        |
1363 !        |          |         ----> i
1364 !        u    p     |
1365 !        |          |    nodes in cell (i,j)
1366 !        ------v-----    view from top
1367 !
1368 !             v(ide,jde+1)
1369 !            -------x------
1370 !            |            |
1371 !            |            |
1372 ! u(ide,jde) x      x     x u(ide+1,jde)
1373 !            | p(ide,hde) |
1374 !            |            |   p=ph,phb,z0,...
1375 !            -------x------
1376 !              v(ide,jde)
1377 !
1378 ! staggered values set u(ids:ide+1,jds:jde) v(ids:ide,jds:jde+1)
1379 ! p=ph+phb set at (ids:ide,jds:jde)
1380 ! location of u(i,j) needs p(i-1,j) and p(i,j)
1381 ! location of v(i,j) needs p(i,j-1) and p(i,j)
1382 ! *** NOTE need HALO in ph, phb
1383 ! so we can compute only u(ids+1:ide,jds:jde) v(ids:ide,jds+1,jde)
1384 ! unless we extend p at the boundary
1385 ! but because we care about the fire way in the inside it does not matter
1386 ! if the fire gets close to domain boundary the simulation is over anyway
1387
1388     ite1=snode(ite,ide,1)
1389     jte1=snode(jte,jde,1)
1390     ! do this in any case to check for nans
1391     call print_3d_stats(its,ite1,kds,kde,jts,jte,ims,ime,kms,kme,jms,jme,u,'wind U in')
1392     call print_3d_stats(its,ite,kds,kde,jts,jte1,ims,ime,kms,kme,jms,jme,v,'wind V in')
1393
1394     if(fire_print_msg.gt.0)then
1395 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1396        write(msg,'(a,f7.2,a)')'interpolate_atm2fire: log-interpolation of wind to',fire_wind_height,' m'
1397        call message(msg)
1398 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1399     endif
1400
1401 ! indexing
1402
1403 ! for w
1404     itst=ifval(ids.eq.its,its,its-1)
1405     itet=ifval(ide.eq.ite,ite,ite+1)
1406     jtst=ifval(jds.ge.jts,jts,jts-1)
1407     jtet=ifval(jde.eq.jte,jte,jte+1)
1408
1409 ! for u
1410     itsu=ifval(ids.eq.its,its+1,its)  ! staggered direction
1411     iteu=ifval(ide.eq.ite,ite,ite+1)  ! staggered direction
1412     jtsu=ifval(jds.ge.jts,jts,jts-1)
1413     jteu=ifval(jde.eq.jte,jte,jte+1)
1414
1415 ! for v
1416     jtsv=ifval(jds.eq.jts,jts+1,jts)  ! staggered direction
1417     jtev=ifval(jde.eq.jte,jte,jte+1)  ! staggered direction
1418     itsv=ifval(ids.ge.its,its,its-1)
1419     itev=ifval(ide.eq.ite,ite,ite+1)
1420
1421
1422 if(fire_print_msg.ge.1)then
1423 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1424   write(msg,7001)'atm input  ','tile',its,ite,jts,jte
1425   call message(msg)
1426   write(msg,7001)'altw       ','tile',itst,itet,jtst,jtet
1427   call message(msg)
1428 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1429 endif
1430 7001 format(a,' dimensions ',a4,':',i6,' to ',i6,' by ',i6,' to ',i6)
1431
1432 !**********************************************************
1433 !*                                                        *
1434 !*           find the altitude of the w points            *
1435 !*                                                        *
1436 !**********************************************************
1437 !! in future, replace by z8w & test if the same
1438
1439     kdmax=kde-1   ! max layer to interpolate from, can be less
1440
1441     do j = jtst,jtet
1442       do k=kds,kdmax+1
1443         do i = itst,itet
1444           altw(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g             ! altitude of the bottom w-point
1445         enddo
1446       enddo
1447     enddo
1448
1449 ! values at u points
1450 if(fire_print_msg.ge.1)then
1451 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1452   write(msg,7001)'u interp at','tile',itsu,iteu,jtsu,jteu
1453   call message(msg)
1454 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1455 endif
1456
1457 !**********************************************************
1458 !*                                                        *
1459 !*  interpolate the altitude from w points to u points    *
1460 !*                                                        *
1461 !**********************************************************
1462
1463     do j = jtsu,jteu
1464       do k=kds,kdmax+1
1465         do i = itsu,iteu
1466           altub(i,k,j)= 0.5*(altw(i-1,k,j)+altw(i,k,j))      ! altitude of the bottom point under u-point
1467         enddo
1468       enddo
1469       do k=kds,kdmax
1470         do i = itsu,iteu
1471           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
1472         enddo
1473       enddo
1474     enddo
1475
1476 ! values at v points
1477 if(fire_print_msg.ge.1)then
1478 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1479   write(msg,7001)'v interp at','tile',itsv,itev,jtsv,jtev
1480   call message(msg)
1481 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1482 endif
1483
1484 !**********************************************************
1485 !*                                                        *
1486 !*  interpolate the altitude from w points to v points    *
1487 !*                                                        *
1488 !**********************************************************
1489
1490     do j = jtsv,jtev
1491       do k=kds,kdmax+1
1492         do i = itsv,itev
1493           altvb(i,k,j)= 0.5*(altw(i,k,j-1)+altw(i,k,j))      ! altitude of the bottom point under v-point
1494         enddo
1495       enddo
1496       do k=kds,kdmax
1497         do i = itsv,itev
1498           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
1499         enddo
1500       enddo
1501     enddo
1502
1503 #ifdef DEBUG_OUT
1504         call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,altub,'altub',id)
1505         call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,altvb,'altvb',id)
1506         call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,hgtu,'hgtu',id)
1507         call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,hgtv,'hgtv',id)
1508 #endif
1509
1510     logfwh = log(fire_wind_height)
1511
1512 !**********************************************************
1513 !*                                                        *
1514 !*  interpolate u vertically to fire_wind_height          *
1515 !*                                                        *
1516 !**********************************************************
1517
1518     ! interpolate u, staggered in X
1519
1520     do j = jtsu,jteu            ! compute on domain by 1 smaller, otherwise z0 and ph not available
1521       do i = itsu,iteu        ! compute with halo 2
1522         zr = 0.5*(z0(i,j)+z0(i-1,j)) ! interpolated roughness length under this u point
1523         if(fire_wind_height > zr)then       !
1524           do k=kds,kdmax
1525             ht = hgtu(i,k,j)      ! height of this u point above the ground
1526             if( .not. ht < fire_wind_height) then ! found layer k this point is in
1527               loght = log(ht)
1528               if(k.eq.kds)then               ! first layer, log linear interpolation from 0 at zr
1529                 logz0 = log(zr)
1530                 ua(i,j)= u(i,k,j)*(logfwh-logz0)/(loght-logz0)
1531               else                           ! log linear interpolation
1532                 loglast=log(hgtu(i,k-1,j))
1533                 ua(i,j)= u(i,k-1,j) + (u(i,k,j) - u(i,k-1,j)) * ( logfwh - loglast) / (loght - loglast)
1534               endif
1535               goto 10
1536             endif
1537             if(k.eq.kdmax)then                 ! last layer, still not high enough
1538               ua(i,j)=u(i,k,j)
1539             endif
1540           enddo
1541 10        continue
1542         else  ! roughness higher than the fire wind height
1543           ua(i,j)=0.
1544         endif
1545       enddo
1546     enddo
1547
1548
1549 !**********************************************************
1550 !*                                                        *
1551 !*  interpolate v vertically to fire_wind_height          *
1552 !*                                                        *
1553 !**********************************************************
1554
1555     ! interpolate v, staggered in Y
1556
1557     do j = jtsv,jtev
1558       do i = itsv,itev
1559         zr = 0.5*(z0(i,j-1)+z0(i,j)) ! roughness length under this v point
1560         if(fire_wind_height > zr)then       !
1561           do k=kds,kdmax
1562             ht = hgtv(i,k,j)      ! height of this u point above the ground
1563             if( .not. ht < fire_wind_height) then ! found layer k this point is in
1564               loght = log(ht)
1565               if(k.eq.kds)then               ! first layer, log linear interpolation from 0 at zr
1566                 logz0 = log(zr)
1567                 va(i,j)= v(i,k,j)*(logfwh-logz0)/(loght-logz0)
1568               else                           ! log linear interpolation
1569                 loglast=log(hgtv(i,k-1,j))
1570                 va(i,j)= v(i,k-1,j) + (v(i,k,j) - v(i,k-1,j)) * ( logfwh - loglast) / (loght - loglast)
1571               endif
1572               goto 11
1573             endif
1574             if(k.eq.kdmax)then                 ! last layer, still not high enough
1575               va(i,j)=v(i,k,j)
1576             endif
1577           enddo
1578 11        continue
1579         else  ! roughness higher than the fire wind height
1580           va(i,j)=0.
1581         endif
1582       enddo
1583     enddo
1584
1585 #ifdef DEBUG_OUT
1586 !   store the output for diagnostics
1587     do j = jts,jte1
1588       do i = its,ite1
1589         uah(i,j)=ua(i,j)
1590         vah(i,j)=va(i,j)
1591       enddo
1592     enddo
1593
1594     call write_array_m(its,ite1,jts,jte,ims,ime,jms,jme,uah,'uah_n',id) ! no reflection
1595     call write_array_m(its,ite,jts,jte1,ims,ime,jms,jme,vah,'vah_n',id)
1596 #endif
1597
1598 !**********************************************************
1599 !*                                                        *
1600 !*  interpolate ua,va vertically to the fire mesh         *
1601 !*                                                        *
1602 !**********************************************************
1603
1604
1605     ips1 = ifval(ips.eq.ids,ips+1,ips)
1606     call continue_at_boundary(1,1,0., & ! x direction
1607        TDIMS,                  &! memory dims atm grid tile
1608        ids+1,ide,jds,jde, &     ! domain dims - where u defined
1609        ips1,ipe,jps,jpe, &     ! patch dims
1610        itsu,iteu,jtsu,jteu, & ! tile dims - in nonextended direction one beyond if at patch boundary but not domain
1611        itsou,iteou,jtsou,jteou, & ! out, where set
1612        ua)                           ! array
1613
1614     jps1 = ifval(jps.eq.jds,jps+1,jps)
1615     call continue_at_boundary(1,1,0., & ! y direction
1616        TDIMS,                  & ! memory dims atm grid tile
1617        ids,ide,jds+1,jde, &      ! domain dims - where v wind defined
1618        ips,ipe,jps1,jpe, &        ! patch dims
1619        itsv,itev,jtsv,jtev, & ! tile dims
1620        itsov,iteov,jtsov,jteov, & ! where set
1621        va)                           ! array
1622
1623 !   store the output for diagnostics
1624     do j = jts,jte1
1625       do i = its,ite1
1626         uah(i,j)=ua(i,j)
1627         vah(i,j)=va(i,j)
1628       enddo
1629     enddo
1630
1631 #ifdef DEBUG_OUT
1632         call write_array_m(itsou,iteou,jtsou,jteou,TDIMS,ua,'ua',id)
1633         call write_array_m(itsov,iteov,jtsov,jteov,TDIMS,va,'va',id)
1634 #endif
1635
1636 !$OMP CRITICAL(SFIRE_DRIVER_CRIT)
1637     ! don't have all values valid, don't check
1638      write(msg,12)'atm mesh wind U at',fire_wind_height,' m'
1639      call print_2d_stats(itsou,iteou,jtsou,jteou,TDIMS,ua,msg)
1640      write(msg,12)'atm mesh wind V at',fire_wind_height,' m'
1641      call print_2d_stats(itsov,iteov,jtsov,jteov,TDIMS,va,msg)
1642 12  format(a,f6.2,a)
1643     call print_2d_stats(its,ite1,jts,jte,ims,ime,jms,jme,uah,'UAH')
1644     call print_2d_stats(its,ite,jts,jte1,ims,ime,jms,jme,vah,'VAH')
1645     !call write_array_m(its,ite1,jts,jte,ims,ime,jms,jme,uah,'uah',id)
1646     !call write_array_m(its,ite,jts,jte1,ims,ime,jms,jme,vah,'vah',id)
1647 !$OMP END CRITICAL(SFIRE_DRIVER_CRIT)
1648
1649 !      ---------------
1650 !     | F | F | F | F |   Example of atmospheric and fire grid with
1651 !     |-------|-------|   ir=jr=4.
1652 !     | F | F | F | F |   Winds are given at the midpoints of the sides of the atmosphere grid,
1653 !     ua------z-------|   interpolated to midpoints of the cells of the fine fire grid F.
1654 !     | F | F | F | F |   This is (1,1) cell of atmosphere grid, and [*] is the (1,1) cell of the fire grid.
1655 !     |---------------|   ua(1,1) <--> uf(0.5,2.5)
1656 !     | * | F | F | F |   va(1,1) <--> vf(2.5,0.5)
1657 !      -------va------    za(1,1) <--> zf(2.5,2.5)
1658 !
1659 !   ^ x2
1660 !   |  --------va(1,2)---------
1661 !   | |            |           |   Example of atmospheric and fire grid with
1662 !   | |            |           |   ir=jr=1.
1663 !   | |          za,zf         |   Winds are given at the midpoints of the sides of the atmosphere grid,
1664 !   | ua(1,1)----uf,vf-----ua(2,1) interpolated to midpoints of the cells of the (the same) fire grid
1665 !   | |           (1,1)        |   ua(1,1) <--> uf(0.5,1)
1666 !   | |            |           |   va(1,1) <--> vf(1,0.5)
1667 !   | |            |           |   za(1,1) <--> zf(1,1)
1668 !   |  --------va(1,1)---------
1669 !   |--------------------> x1
1670 !
1671 ! Meshes are aligned by the lower left cell of the domain. Then in the above figure
1672 ! u = node with the ua component of the wind at (ids,jds), midpoint of side
1673 ! v = node with the va component of the wind at (ids,jds), midpoint of side
1674 ! * = fire grid node at (ifds,jfds)
1675 ! z = node with height, midpoint of cell
1676 !
1677 ! ua(ids,jds)=uf(ifds-0.5,jfds+jr*0.5-0.5)         = uf(ifds-0.5,jfds+(jr-1)*0.5)
1678 ! va(ids,jds)=vf(ifds+ir*0.5-0.5,jfds-0.5)         = vf(ifds+(ir-1)*0.5,jfds-0.5)
1679 ! 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)
1680
1681     ! ifts1=max(snode(ifts,ifps,-1),ifds) ! go 1 beyond patch boundary but not at domain boundary
1682     ! ifte1=min(snode(ifte,ifpe,+1),ifde)
1683     ! jfts1=max(snode(jfts,jfps,-1),jfds)
1684     ! jfte1=min(snode(jfte,jfpe,+1),jfde)
1685
1686     call interpolate_2d(  &
1687         TDIMS,                  & ! memory dims atm grid tile
1688         itsou,iteou,jtsou,jteou,& ! where set
1689         ifms,ifme,jfms,jfme,    & ! array dims fire grid
1690         ifts,ifte,jfts,jfte,& ! dimensions on the fire grid to interpolate to
1691         ir,jr,                  & ! refinement ratio
1692         real(ids),real(jds),ifds-0.5,jfds+(jr-1)*0.5, & ! line up by lower left corner of domain
1693         ua,                     & ! in atm grid
1694         uf)                      ! out fire grid
1695
1696     call interpolate_2d(  &
1697         TDIMS,                  & ! memory dims atm grid tile
1698         itsov,iteov,jtsov,jteov,& ! where set
1699         ifms,ifme,jfms,jfme,    & ! array dims fire grid
1700         ifts,ifte,jfts,jfte,& ! dimensions on the fire grid to interpolate to
1701         ir,jr,                  & ! refinement ratio
1702         real(ids),real(jds),ifds+(ir-1)*0.5,jfds-0.5, & ! line up by lower left corner of domain
1703         va,                     & ! in atm grid
1704         vf)                      ! out fire grid
1705
1706 call print_2d_stats_vec(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,vf,'fire wind (m/s)')
1707 ! call print_2d_stats_vec(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,vf,'fire wind extended')
1708 #ifdef DEBUG_OUT
1709         call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,'uf1',id)
1710         call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,vf,'vf1',id)
1711         ! call write_array_m(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,'uf1',id)
1712         ! call write_array_m(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,vf,'vf1',id)
1713 #endif
1714
1715 return
1716
1717 end subroutine interpolate_atm2fire
1718
1719 subroutine apply_windrf(                            &
1720             ifms,ifme,jfms,jfme,                    &
1721             ifts,ifte,jfts,jfte,                    &
1722             nfuel_cat,uf,vf)
1723 integer::                                           &
1724     ifms, ifme, jfms, jfme,                       & ! fire memory bounds
1725     ifts, ifte, jfts, jfte                          ! fire tile bounds
1726 real,intent(in),dimension(ifms:ifme,jfms:jfme)::nfuel_cat
1727 real,intent(inout),dimension(ifms:ifme,jfms:jfme)::uf,vf
1728 !*** local
1729     integer::i,j,k
1730 !*** executable
1731
1732     do j = jfts,jfte
1733       do i = ifts,ifte
1734           k=int( nfuel_cat(i,j) )
1735           if(k.lt.no_fuel_cat)then
1736               uf(i,j)=uf(i,j)*windrf(k)
1737               vf(i,j)=vf(i,j)*windrf(k)
1738           else
1739               uf(i,j)=0.
1740               vf(i,j)=0.
1741           endif
1742       enddo
1743     enddo
1744
1745 end subroutine apply_windrf
1746
1747 !
1748 !***
1749 !
1750
1751 subroutine setup_wind_log_interpolation(          &
1752             ids,ide,  jds,jde,                    & ! atm grid dimensions
1753             ims,ime,  jms,jme,                    &
1754             ips,ipe,jps,jpe,                              &
1755             its,ite,jts,jte,                              &
1756             ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
1757             ifms, ifme, jfms, jfme,                       &
1758             ifts, ifte, jfts, jfte,                       &
1759             ir,jr,                                        & ! atm/fire grid ratio
1760             z0,                                           & ! atm grid arrays in
1761             nfuel_cat,                                    & ! fire array in
1762             fz0,fwh)                                        ! fire arrays out
1763 !*** arguments
1764 integer, intent(in)::                    &
1765     ids,ide, jds,jde,                    & ! atm domain bounds
1766     ims,ime, jms,jme,                    & ! atm memory bounds
1767     ips,ipe,jps,jpe,                     &
1768     its,ite,jts,jte,                     & ! atm tile bounds
1769     ifds, ifde, jfds, jfde,                       & ! fire domain bounds
1770     ifms, ifme, jfms, jfme,                       & ! fire memory bounds
1771     ifts, ifte, jfts, jfte,                       & ! fire tile bounds
1772     ir,jr                                         ! atm/fire grid refinement ratio
1773 real,intent(in),dimension(ims:ime,jms:jme)::z0    ! landuse roughness length
1774 real,intent(in),dimension(ifms:ifme,jfms:jfme)::nfuel_cat   ! fuel category
1775 real,intent(out),dimension(ifms:ifme,jfms:jfme)::&
1776     fz0,                                          & ! roughness height
1777     fwh                                             ! height to read the wind from
1778 !*** local
1779 integer::i,j,ii,jj,k,id=0
1780 character(len=128)::msg
1781 real::r
1782 !*** executable
1783
1784      if(.not.have_fuel_cats)call crash('setup_wind_log_interpolation: fuel categories not yet set')
1785
1786      select case(fire_wind_log_interp)
1787
1788      case(1)
1789        call message('fire_wind_log_interp=1: log interpolation on fire mesh, roughness and wind height from fuel categories')
1790        do j=jfts,jfte
1791           do i=ifts,ifte
1792             k = int(nfuel_cat(i,j))
1793             if(k.ge.no_fuel_cat.and.k.le.no_fuel_cat2)then ! no fuel
1794                 fz0(i,j)=-1.
1795                 fwh(i,j)=-1.
1796             elseif(k < 1 .or. k  > nfuelcats) then
1797 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1798                  write(msg,*)'i,j,nfuel_cat,nfuelcats=',i,j,k,nfuelcats
1799 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1800                  call message(msg)
1801                  call crash('setup_wind_log_interpolation: fuel category out of bounds')
1802             else
1803                 fz0(i,j)=fcz0(k)
1804                 fwh(i,j)=fcwh(k)
1805             endif
1806           enddo
1807        enddo
1808
1809      case(2)
1810         call message('fire_wind_log_interp=2: log interpolation on fire mesh' // &
1811           'piecewise constant roughness from landuse, constant fire_wind_height')
1812         do j=jts,jte
1813            do i=its,ite
1814               do jj=(j-1)*jr+1,(j-1)*jr+jr
1815                  do ii=(i-1)*ir+1,(i-1)*ir+ir
1816                      fz0(ii,jj)=z0(i,j)
1817                  enddo
1818               enddo
1819            enddo
1820        enddo
1821        do j=jfts,jfte
1822           do i=ifts,ifte
1823             k = int(nfuel_cat(i,j))
1824             if(k.lt.no_fuel_cat)then ! no fuel, interpolation does not matter
1825                 fwh(i,j)=fcwh(k)
1826             else
1827                 fz0(i,j)=-1.
1828                 fwh(i,j)=-1.
1829             endif
1830           enddo
1831        enddo
1832
1833      case(3)
1834        call message('fire_wind_log_interp=3: log interpolation on fire mesh,' // &
1835            ' interpolated roughness from landuse, constant fire_wind_height')
1836        call interpolate_z2fire(id,1,                & ! for debug output, <= 0 no output
1837             ids,ide,  jds,jde,                    & ! atm grid dimensions
1838             ims,ime,  jms,jme,                    &
1839             ips,ipe,jps,jpe,                              &
1840             its,ite,jts,jte,                              &
1841             ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
1842             ifms, ifme, jfms, jfme,                       &
1843             ifts,ifte,jfts,jfte,                          &
1844             ir,jr,                                        & ! atm/fire grid ratio
1845             z0,                                       & ! atm grid arrays in
1846             fz0)                                      ! fire grid arrays out
1847        do j=jfts,jfte
1848           do i=ifts,ifte
1849             k = int(nfuel_cat(i,j))
1850             if(k.ne.no_fuel_cat)then ! no fuel, interpolation does not matter
1851                 fwh(i,j)=fcwh(k)
1852             else
1853                 fz0(i,j)=-1.
1854                 fwh(i,j)=-1.
1855             endif
1856           enddo
1857        enddo
1858
1859      case(4)
1860         call message('fire_wind_log_interp=4: log interpolation on atmospheric' // &
1861            ' mesh, roughness from landuse, constant fire_wind_height')
1862         return
1863
1864      case default
1865         !$OMP CRITICAL(SFIRE_ATM_CRIT)
1866         write(msg,*)'setup_wind_log_interpolation: invalid fire_wind_log_interp=',fire_wind_log_interp
1867         !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1868         call crash('msg')
1869
1870      end select
1871
1872      select case(fire_use_windrf)
1873
1874      case(0)
1875        call message('setup_wind_log_interpolation: not using wind reduction factors')
1876
1877      case(1)
1878        call message('setup_wind_log_interpolation: multiplying wind by reduction factors')
1879
1880      case(2)
1881        call message('setup_wind_log_interpolation: resetting wind interpolation height from wind reduction factors')
1882        do j=jfts,jfte
1883           do i=ifts,ifte
1884             k = int(nfuel_cat(i,j))
1885             if(k.ne.no_fuel_cat)then
1886                 fwh(i,j) = fz0(i,j) ** (1.-windrf(k)) * fire_wind_height ** windrf(k) ! GMD paper eq. (26)
1887
1888                 if (.not. fz0(i,j) > 0. .or. .not. fwh(i,j) > fz0(i,j))then
1889 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1890                   write(msg,*)'category ',k,'windrf=',windrf(k),' fire_wind_height=',fire_wind_height
1891                   call message(msg,level=-1)
1892                   write(msg,*)'i=',i,' j=',j,' fz0(i,j)=',fz0(i,j),' fwh(i,j)=',fwh(i,j)
1893                   call message(msg,level=-1)
1894 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1895                   call crash('setup_wind_log_interpolation: must have fwh > fz0 > 0')
1896                 endif
1897
1898             endif
1899           enddo
1900        enddo
1901
1902      case(3)
1903       if(fire_wind_log_interp.eq.2.or.fire_wind_log_interp.eq.3)then
1904        call message('setup_wind_log_interpolation: adjusting wind interpolation height for LANDUSE roughness height')
1905        do j=jfts,jfte
1906           do i=ifts,ifte
1907             k = int(nfuel_cat(i,j))
1908             if(k.lt.no_fuel_cat)then
1909                 r = log(fcwh(k)/fcz0(k))/log(fire_wind_height/fcz0(k))! GMD paper eq. (25)
1910                 fwh(i,j) = fz0(i,j) ** (1.-r) * fire_wind_height ** r ! GMD paper eq. (26)
1911
1912                 if (.not. fz0(i,j) > 0. .or. .not. fwh(i,j) > fz0(i,j))then
1913 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1914                   write(msg,*)'category ',k, 'roughness ',fcz0(k),' midflame height ',fcwh(k),' fire_wind_height=',fire_wind_height
1915                   call message(msg,level=-1)
1916                   write(msg,*)'computed wind reduction factor ',r
1917                   call message(msg,level=-1)
1918                   write(msg,*)'i=',i,' j=',j,' fz0(i,j)=',fz0(i,j),' fwh(i,j)=',fwh(i,j)
1919                   call message(msg,level=-1)
1920 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1921                   call crash('setup_wind_log_interpolation: must have fwh > fz0 > 0')
1922                 endif
1923
1924             endif
1925           enddo
1926        enddo
1927       else
1928        call message('setup_wind_log_interpolation: not using wind reduction factors')
1929       endif
1930
1931      case default
1932         !$OMP CRITICAL(SFIRE_ATM_CRIT)
1933         write(msg,*)'setup_wind_log_interpolation: invalid fire_use_windrf=',fire_use_windrf
1934         !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1935         call crash('msg')
1936
1937      end select
1938
1939 !      consistency check
1940        do j=jfts,jfte
1941           do i=ifts,ifte
1942             k = int(nfuel_cat(i,j))
1943             if(k.lt.no_fuel_cat)then
1944               if (.not. fz0(i,j) > 0. .or. .not. fwh(i,j) > fz0(i,j))then
1945 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1946                  write(msg,*)'i=',i,' j=',j,' fz0(i,j)=',fz0(i,j),' fwh(i,j)=',fwh(i,j)
1947                  call message(msg,level=-1)
1948 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1949                  call crash('setup_wind_log_interpolation: must have fwh > fz0 > 0')
1950               endif
1951             else
1952               if(.not.fwh(i,j)<0.)then
1953 !$OMP CRITICAL(SFIRE_ATM_CRIT)
1954                  write(msg,*)'i=',i,' j=',j,' fz0(i,j)=',fz0(i,j),' fwh(i,j)=',fwh(i,j)
1955                  call message(msg,level=-1)
1956 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
1957                  call crash('setup_wind_log_interpolation: no fuel must be signalled by fwh<0')
1958               endif
1959             endif
1960           enddo
1961        enddo
1962
1963      have_wind_log_interpolation = .true.
1964
1965      call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fz0,'setup_wind_log:fz0')
1966      call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fz0,'setup_wind_log:fwh')
1967
1968 end subroutine setup_wind_log_interpolation
1969
1970 !
1971 !***
1972 !
1973
1974 subroutine interpolate_wind2fire_height(id,       & ! to identify debugging prints and files if needed
1975     ids,ide, kds,kde, jds,jde,                    & ! atm grid dimensions
1976     ims,ime, kms,kme, jms,jme,                    &
1977     ips,ipe,jps,jpe,                              &
1978     its,ite,jts,jte,                              &
1979     ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
1980     ifms, ifme, jfms, jfme,                       &
1981     ifps, ifpe, jfps, jfpe,                       & ! fire patch bounds
1982     ifts,ifte,jfts,jfte,                          &
1983     ir,jr,                                        & ! atm/fire grid ratio
1984     u_frame, v_frame,                             & ! velocity frame correction
1985     u,v,ph,phb,                                   & ! input atmospheric arrays
1986     fz0,fwh,                                      & ! input fire arrays
1987     uf,vf)                                          ! output fire arrays
1988
1989   implicit none
1990 !*** arguments
1991 integer, intent(in):: id,                         & ! debug identification
1992     ids,ide, kds,kde, jds,jde,                    & ! atm domain bounds
1993     ims,ime, kms,kme, jms,jme,                    & ! atm memory bounds
1994     ips,ipe,jps,jpe,                              &
1995     its,ite,jts,jte,                              & ! atm tile bounds
1996     ifds, ifde, jfds, jfde,                       & ! fire domain bounds
1997     ifms, ifme, jfms, jfme,                       & ! fire memory bounds
1998     ifps, ifpe, jfps, jfpe,                       & ! fire patch bounds
1999     ifts,ifte,jfts,jfte,                          & ! fire tile bounds
2000     ir,jr                                         ! atm/fire grid refinement ratio
2001 real, intent(in):: u_frame, v_frame                 ! velocity frame correction
2002 real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::&
2003     u,v,                                          & ! atm wind velocity, staggered
2004     ph,phb                                          ! geopotential
2005 real,intent(in),dimension(ifms:ifme,jfms:jfme)::&
2006     fz0,                                          & ! roughness height
2007     fwh                                             ! height to read the wind from
2008 real,intent(out),dimension(ifms:ifme,jfms:jfme)::&
2009     uf,                                           & ! atm wind at fire_wind_height, diagnostics
2010     vf                                              !
2011
2012 !*** local
2013   integer:: i,j,k,jcb,jcm,icb,icm,kdmax,kmin,kmax
2014   integer::itst,itet,jtst,jtet
2015   integer::iftst,iftet,jftst,jftet
2016   real:: wjcb,wjcm,wicb,wicm,ht,i_g,loght,zr,ht_last,logwh,wh,loght_last,uk,vk,uk1,vk1,z0,logz0
2017   real, dimension (its-1:ite+1,kds:kde,jts-1:jte+1):: z
2018   character(len=128)::msg
2019
2020 !*** executable
2021
2022   if(.not. have_wind_log_interpolation) call crash('interpolate_wind2fire_height: wind_log_interpolation must be set up first')
2023
2024   ! print *,'interpolate_wind2fire_height start, id=',id
2025
2026   kdmax=kde-1                                       ! max layer to use
2027
2028 ! find the altitude of atm cell centers
2029
2030 ! index bounds for cell centers - need to go one beyond if at end of tile but not domain
2031   itst=ifval(ids.eq.its,its,its-1)
2032   itet=ifval(ide.eq.ite,ite,ite+1)
2033   jtst=ifval(jds.ge.jts,jts,jts-1)
2034   jtet=ifval(jde.eq.jte,jte,jte+1)
2035
2036   ! print *,'its, ite, jts, jte    =',its, ite, jts, jte
2037   ! print *,'itst, itet, jtst, jtet=',itst, itet, jtst, jtet
2038
2039 ! get altitudes
2040   i_g = 1./g
2041   do j = jtst,jtet
2042     do k=kds,kdmax+1
2043       do i = itst,itet
2044         z(i,k,j) = (ph(i,k,j)+phb(i,k,j))*i_g      ! altitude of the bottom w-point
2045
2046         ! print *,'i,k,j=',i,k,j,'ph, phb, z=',ph(i,k,j),phb(i,k,j),z(i,k,j)
2047
2048       enddo
2049     enddo
2050     do k=kds,kdmax
2051       do i = itst,itet
2052         z(i,k,j) = (z(i,k,j)+z(i,k+1,j))*0.5 - z(i,kds,j)   ! height of the cell center
2053       enddo
2054     enddo
2055   enddo
2056
2057 ! index bounds for fire mesh cell centers
2058 ! to prevent setting values from uninitialized memory
2059   iftst=ifval(ifds.eq.ifts,ifts+ir/2,ifts)
2060   iftet=ifval(ifde.eq.ifte,ifte-ir/2,ifte)
2061   jftst=ifval(jfds.ge.jfts,jfts+jr/2,jfts)
2062   jftet=ifval(jfde.eq.jfte,jfte-jr/2,jfte)
2063
2064   ! print *,'iftst, iftet, jftst, jftet=',iftst, iftet, jftst, jftet
2065
2066 ! zero out first, to prevent unitialized values on strips along domain boundaries
2067 ! it would be faster but longer code to do just cleanup loop on the strips
2068    do j = jfts,jfte
2069      do i = ifts,ifte
2070        uf(i,j)=0.
2071        vf(i,j)=0.
2072      enddo
2073   enddo
2074
2075 ! vertical and horizontal interpolation
2076
2077   kmin=kde  ! init stats
2078   kmax=kds
2079
2080   loop_j: do j = jftst,jftet
2081     call coarse(j,jr,-2,jcb,wjcb)    ! get interpolation coefficients from the boundary
2082     call coarse(j,jr,ir,jcm,wjcm)    ! get interpolation coefficients from the midpoint
2083     loop_i: do i = iftst,iftet
2084       call coarse(i,ir,-2,icb,wicb)    ! get interpolation coefficients from the boundary
2085       call coarse(i,ir,ir,icm,wicm)    ! get interpolation coefficients from the midpoint
2086       z0 = fz0(i,j)                    ! roughness length
2087       wh = fwh(i,j)                    ! wind height
2088
2089       ! print *,'i=',i,' j=',j,' icb=',icb,' jcb=',jcb,' z0=',z0,' wh=',wh
2090
2091
2092       if( wh > z0 .and. z0 > 0)then
2093
2094       ht_last=z0                     ! initialize starting height of this layer
2095       loop_k: do k=kds,kdmax         ! search for layer k such that ht(k-1)<=wh<ht(k), ht(0)=z0
2096         ! interpolate height from atmospheric cell midpoints
2097         ht=interpolate_h(its-1,ite+1,kds,kde,jts-1,jts+1,icm,k,jcm,wicm,wjcm,z)
2098
2099         ! print *,'i=',i,' j=',j,'k=',k,' ht=',ht
2100
2101         if(.not. ht < wh) exit loop_k ! found layer k this point is in
2102         ht_last = ht
2103       enddo loop_k
2104       if(k .gt. kdmax) then
2105           goto 91  ! run out of vertical levels, this must be wrong
2106       endif
2107       kmin=min(k,kmin)
2108       kmax=max(k,kmax)
2109       ! found layer k, ht_last < wh <= ht
2110       logz0 = log(z0)
2111       logwh= log(wh)
2112       loght_last = log(ht_last)
2113       loght = log(ht)
2114
2115       ! interpolate u at level k from staggered coarse grid: boundary in i, midpoints in j
2116       uk=interpolate_h(ims,ime,kms,kme,jms,jme,icb,k,jcm,wicb,wjcm,u) - u_frame
2117
2118       ! print *,'i=',i,' j=',j,' uk=',uk
2119
2120       ! interpolate v at level k from staggered coarse grid: midpoints in i, boundary in j
2121       vk=interpolate_h(ims,ime,kms,kme,jms,jme,icm,k,jcb,wicm,wjcb,v) - v_frame
2122
2123       ! print *,'i=',i,' j=',j,' vk=',vk
2124
2125       if(k.gt.kds)then           ! interpolate u,v horizontaly at the previous layer,  k-1
2126         ! interpolate u at level k-1 from staggered coarse grid: boundary in i, midpoints in j
2127         uk1=interpolate_h(ims,ime,kms,kme,jms,jme,icb,k-1,jcm,wicb,wjcm,u)
2128
2129         ! print *,'i=',i,' j=',j,' uk1=',uk1
2130
2131         ! interpolate v at level k-1 from staggered coarse grid: midpoints in i, boundary in j
2132         vk1=interpolate_h(ims,ime,kms,kme,jms,jme,icm,k-1,jcb,wicm,wjcb,v)
2133
2134         ! print *,'i=',i,' j=',j,' uk1=',uk1
2135
2136       else ! there is no previous layer, wind at roughness is assumed 0
2137         uk1=0.
2138         vk1=0.
2139       endif
2140
2141       ! log interpolate the wind to wh
2142       uf(i,j)= uk1 + (uk - uk1) * ( logwh - loght_last) / (loght - loght_last)
2143       vf(i,j)= vk1 + (vk - vk1) * ( logwh - loght_last) / (loght - loght_last)
2144
2145       ! print *,'i=',i,' j=',j,' uk=',uk,' vk=',vk,' uk1=',uk1,' vk1=',vk1,' uf=',uf(i,j),' vf=',vf(i,j)
2146
2147     else
2148       ! no fuel, no wind
2149       uf(i,j) = 0.
2150       vf(i,j) = 0.
2151
2152     endif
2153
2154     enddo loop_i
2155   enddo loop_j
2156
2157   ! print *,'interpolate_wind2fire_height complete, id=',id
2158
2159 !$OMP CRITICAL(SFIRE_ATM_CRIT)
2160   write(msg,*)'wind interpolated from layers',kmin,' to ',kmax
2161   call message(msg)
2162 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
2163
2164   return
2165
2166   91 call crash('interpolate_wind2fire_height: fire wind height too large, increase kdmax or atm height')
2167   92 continue
2168 !$OMP CRITICAL(SFIRE_ATM_CRIT)
2169   write(msg,*)'fz0(',i,j,')=',fz0(i,j),'fwh(',i,j,')=',fwh(i,j)
2170   call message(msg)
2171 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
2172   call crash('interpolate_wind2fire_height: must have fire wind height > roughness height > 0')
2173
2174 contains
2175
2176 real function interpolate_h(ims,ime,kms,kme,jms,jme,ic,kc,jc,wic,wjc,a)
2177 !*** purpose: bilinear interpolation from a(ic:ic+1,k,jc:jc+1) with weights wicm wjcm
2178   implicit none
2179 !*** arguments
2180   integer, intent(in)::ims,ime,jms,kms,kme,jme,ic,kc,jc
2181   real, intent(in)::wic,wjc,a(ims:ime,kms:kme,jms:jme)
2182 !*** executable
2183   interpolate_h = &
2184     a(ic,kc,jc)    *wic     *wjc      + &
2185     a(ic,kc,jc+1)  *wic     *(1.-wjc) + &
2186     a(ic+1,kc,jc)  *(1.-wic)*wjc      + &
2187     a(ic+1,kc,jc+1)*(1.-wic)*(1.-wjc)
2188 end function interpolate_h
2189
2190
2191 subroutine coarse(ix,nr,ia,ic,w)
2192 !*** find coarse mesh index and interpolation weight
2193 !*** arguments
2194     implicit none
2195     integer, intent(in)::ix,nr,ia
2196     integer, intent(out)::ic
2197     real, intent(out)::w
2198 !*** description
2199 ! given fine mesh with nr cells for each coarse mesh cell and such that
2200 !   coarse mesh node 1 is aligned with the fine mesh at (na+1)/2
2201 ! for fine mesh node ix find its coarse mesh coordinate c and return
2202 !   ic=floor(c), the index of the nearest coarse mesh node below
2203 !   w =1-(c-ic), the interpolation weight from coarse mesh node ic to fine mesh node ix
2204 !
2205 ! Intended use:
2206 ! fine mesh nodes are always at the middle of their cells
2207 !
2208 ! the alignment when the coarse nodes are on the boundary of coarse cells:
2209 ! |---1---|---2---|.......|--nr---|   fine mesh
2210 ! 1-------------------------------2   coarse mesh
2211 ! ia = -2  because coarse node 1 is aligned with the fine mesh at -1/2 = (-2 + 1)/2
2212 !
2213 ! the alignment when the coarse node is at the midpoint of coarse  cell:
2214 ! |---1---|---2---|---3---|---4---|   fine mesh, here nr=4
2215 ! |---------------1---------------|   coarse mesh
2216 ! ia = nr because coarse node 1 is aligned with the fine mesh at (nr + 1)/2
2217 ! here,  (4 + 1)/2 = 2.5
2218 !
2219
2220 !*** local
2221     real:: c,a
2222
2223 !*** executable
2224
2225   a = (ia + 1)*0.5    ! the location on the fine grid where coarse node 1 is aligned
2226   c = 1 + (ix - a)/nr ! coarse mesh coordinate of ix
2227   ic= floor(c)        ! nearest coarse node to the left
2228   w = (1 + ic) - c    ! interpolation weight, 1-(c-ic)
2229
2230 end subroutine coarse
2231
2232 end subroutine interpolate_wind2fire_height
2233
2234 !#ifdef WRF_CHEM
2235 subroutine fire_emission(                    &
2236              tracer_opt,                     &
2237              ids,ide, kds,kde, jds,jde,      & ! domain dimensions
2238              ims,ime, kms,kme, jms,jme,      &
2239              its,ite, kts,kte, jts,jte,      &
2240              rho,dz8w,                       &
2241              grnhfx,                         & ! input variables from fire model
2242              tracer                          ) ! output emissions array
2243 use module_state_description , only: num_tracer, p_tr17_1
2244 #ifdef WRF_CHEM
2245 use module_state_description , only: p_smoke
2246 #endif
2247     integer, intent(in)::tracer_opt
2248     INTEGER , INTENT(in) :: ids,ide, kds,kde, jds,jde, &
2249                             ims,ime, kms,kme, jms,jme, &
2250                             its,ite, kts,kte, jts,jte
2251
2252     real,intent(in),dimension( ims:ime,jms:jme ) :: grnhfx
2253     real,intent(inout),dimension( ims:ime,kms:kme,jms:jme,num_tracer ) :: tracer
2254     real,intent(in),dimension( ims:ime,kms:kme,jms:jme ) :: rho,dz8w
2255
2256     integer::i,j,k,l
2257     character(len=128)::msg
2258     real::t,s
2259
2260     ! just a dumb placeholder
2261     k=kds      ! dump into surface layer
2262
2263     select case(tracer_opt)
2264     case(0)
2265        return ! no tracers
2266     case(1)
2267 #ifdef WRF_CHEM
2268         l=p_smoke
2269 #else
2270         call crash('fire_emission: tracer_opt=1 requires WRF-Chem')
2271 #endif
2272     case(2)
2273         l=p_tr17_1
2274     case default
2275         call crash('fire_emission: tracer_opt not supported')
2276     end select
2277
2278     if(num_tracer<l)call crash('fire_emission: num_tracer too small')
2279     s=0
2280     do j=jts,jte
2281         do i=its,ite
2282             t = grnhfx(i,j)/(rho(i,k,j)*dz8w(i,k,j))
2283             tracer(i,k,j,l) = tracer(i,k,j,l) + t
2284             s = s + t
2285         enddo
2286     enddo
2287 !$OMP CRITICAL(SFIRE_ATM_CRIT)
2288     if(fire_print_msg >= 1)then
2289         write(msg,'(a,4i6,a,e13.4,a,i4)')'Tile ',its,ite,jts,jte,' added ',s,' total to tracer',l
2290         call message(msg,level=0)
2291     endif
2292 !$OMP END CRITICAL(SFIRE_ATM_CRIT)
2293 end subroutine fire_emission
2294 !#endif
2295
2296
2297 end module module_fr_sfire_atm