phys/module_fr_fire_driver.F

   1
   2 ! SFIRE - Spread fire model in WRF-Fire
   3 !
   4 ! This module is the entry point from WRF ARW to the wildland
   5 ! fire module. The call to fire_driver advances the fire module by
   6 ! one timestep. The fire module inputs the wind and outputs
   7 ! temperature and humidity tendencies. The fire module also inputs a
   8 ! number of constant arrays (fuel data, topography). Additional
   9 ! arguments are model state (for data assimilation) and constant arrays
  10 ! the model gives to WRF for safekeeping because it is not allowed
  11 ! to save anything.
  12 !
  13 ! Contributions to this wildland fire module have come from individuals at
  14 ! NCAR, the U.S.D.A. Forest Service, the Australian Bureau of Meteorology,
  15 ! and the University of Colorado at Denver.
  16 !
  17
  18 module module_fr_fire_driver
  19 ! use this module for standalone call, you only need to provide some mock-up wrf modules
  20
  21 use module_fr_fire_model
  22 use module_fr_fire_phys, only : fire_params , init_fuel_cats
  23 use module_fr_fire_util
  24 use module_fr_fire_core, only: ignition_line_type
  25 use module_fr_fire_atm, only: add_fire_tracer_emissions
  26
  27 USE module_domain, only: domain
  28 USE module_configure, only: grid_config_rec_type
  29 use module_model_constants, only: reradius,    & ! 1/earth radiusw
  30                                   g,           & ! gravitation acceleration
  31                                   pi2            ! 2*pi
  32
  33 implicit none
  34
  35
  36 private
  37 public fire_driver_em,use_atm_vars
  38
  39 logical:: use_atm_vars=.true.   !  interpolate wind from atm mesh and average output fluxes back
  40
  41 logical:: fmoist_run, fmoist_interp, fire_run  ! which kind of model to run overall
  42
  43 contains
  44
  45 #define DEBUG_OUT
  46 subroutine fire_driver_em ( grid , config_flags                    &
  47             ,fire_ifun_start,fire_ifun_end,tsteps                   &
  48             ,ids,ide, kds,kde, jds,jde                              &
  49             ,ims,ime, kms,kme, jms,jme                              &
  50             ,ips,ipe, kps,kpe, jps,jpe                              &
  51             ,ifds,ifde, jfds,jfde                                   &
  52             ,ifms,ifme, jfms,jfme                                   &
  53             ,ifps,ifpe, jfps,jfpe                                   &
  54             ,rho,z_at_w,dz8w           &
  55             )
  56
  57 !*** purpose: driver from grid structure
  58
  59 ! Driver layer modules
  60 #ifdef DM_PARALLEL
  61     USE module_dm        , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks
  62     USE module_comm_dm , ONLY : halo_fire_fuel_sub, halo_fire_tign_sub, halo_fire_wind_f_sub, &
  63 halo_fire_wind_a_sub, halo_fire_ph_sub, halo_fire_zsf_sub, halo_fire_longlat_sub, &
  64 halo_fire_phb_sub, halo_fire_z0_sub, halo_fire_lfn_sub, halo_fire_mfg_sub, &
  65 halo_fire_mag_sub
  66 #endif
  67
  68     implicit none
  69 !*** arguments
  70     TYPE(domain) , TARGET :: grid                             ! state
  71     TYPE (grid_config_rec_type) , INTENT(IN)  :: config_flags ! namelist
  72     integer, intent(in)::     fire_ifun_start,fire_ifun_end,tsteps ! driver cycle controls
  73     integer, intent(in):: &
  74              ids,ide, kds,kde, jds,jde,                              &
  75              ims,ime, kms,kme, jms,jme,                              &
  76              ips,ipe, kps,kpe, jps,jpe,                              &
  77              ifds,ifde, jfds,jfde,                                   &
  78              ifms,ifme, jfms,jfme,                                   &
  79              ifps,ifpe, jfps,jfpe
  80     real,dimension(ims:ime, kms:kme, jms:jme),intent(in), optional::rho,z_at_w,dz8w
  81
  82 !*** local
  83     INTEGER:: fire_num_ignitions
  84     integer, parameter::fire_max_ignitions=5
  85     TYPE(ignition_line_type), DIMENSION(fire_max_ignitions):: ignition_line
  86     integer::fire_ifun,ir,jr,fire_ignition_longlat,istep,itimestep
  87     logical::need_lfn_update,restart
  88     real, dimension(ifms:ifme, jfms:jfme)::lfn_out
  89     real:: corner_ll,corner_ul,corner_ur,corner_lr
  90     character(len=128)msg
  91     real:: unit_fxlong ,unit_fxlat
  92     type(fire_params)::fp
  93     real:: time_start, dt
  94     real:: moisture_time
  95     logical:: run_advance_moisture,run_fuel_moisture, moisture_initializing
  96     real::    dt_moisture
  97
  98
  99
 100 !*** executable
 101
 102 ! populate our structures from wrf
 103
 104     ! pointers to be passed to fire rate of spread formulas
 105     fp%vx => grid%uf         ! W-E winds used in fire module
 106     fp%vy => grid%vf         ! S-N winds used in fire module
 107     fp%zsf => grid%zsf       ! terrain height
 108     fp%dzdxf => grid%dzdxf   ! terrain grad
 109     fp%dzdyf => grid%dzdyf   ! terrain grad
 110     fp%bbb => grid%bbb       ! a rate of spread formula coeff
 111     fp%betafl => grid%betafl ! a rate of spread formula variable
 112     fp%phiwc => grid%phiwc   ! a rate of spread formula coeff
 113     fp%r_0 => grid%r_0       ! a rate of spread formula variable
 114     fp%fgip => grid%fgip     ! a rate of spread formula coeff
 115     fp%ischap => grid%ischap ! a rate of spread formula switch
 116     fp%iboros => grid%iboros ! Ib divided by ROS
 117     fp%fmc_g => grid%fmc_g   ! fuel moisture, ground
 118
 119     ! get ignition data
 120     call fire_ignition_convert (config_flags,fire_max_ignitions,fire_ignition_longlat, &
 121         ignition_line,fire_num_ignitions,unit_fxlong,unit_fxlat)
 122
 123     ! copy configuration flags to fire internal structures
 124     call set_flags(config_flags)
 125
 126     ! refinement r
 127     ir=grid%sr_x ! refinement ratio
 128     jr=grid%sr_y
 129     itimestep=grid%itimestep
 130     restart=config_flags%restart
 131
 132     ! time - assume dt does not change
 133     dt = grid%dt
 134     time_start = itimestep * dt
 135
 136     ! **** moisture model
 137
 138     ! decide what to run - moisture, interpolation, or fire model itself
 139     fmoist_run    = config_flags%fmoist_run
 140     fmoist_interp = config_flags%fmoist_interp
 141     if(fire_fmc_read.ne.0.and.fmoist_run)call crash('fmoist_run=T requires fire_fmc_read=0')
 142     fire_run = .not. config_flags%fmoist_only
 143
 144     !decide what to run
 145     moisture_time = time_start
 146     run_advance_moisture = .false. ! default
 147     run_fuel_moisture = .false. ! default
 148     moisture_initializing = fire_ifun_start < 3
 149
 150
 151
 152     if(fmoist_run)then
 153         if(moisture_initializing)then
 154             if(fire_ifun_end>2)call crash('initialization must be run separately')
 155             grid%fmoist_lasttime=moisture_time ! initialize the last time the model has run to start of run
 156             grid%fmoist_nexttime=moisture_time
 157             call message('moisture initialization')
 158             run_advance_moisture = .true.
 159         else ! regular timestep
 160             if(config_flags%fmoist_freq > 0)then  ! regular timestep. go by multiples?
 161                 if(mod(grid%itimestep,config_flags%fmoist_freq) .eq. 0)then
 162                     write(msg,'(a,i10,a,i10)')'moisture model runs because timestep ',grid%itimestep,' is a multiple of ',config_flags%fmoist_freq
 163                     call message(msg)
 164                     run_advance_moisture = .true.
 165                 endif
 166             else
 167                 if(.not. moisture_time  < grid%fmoist_nexttime) then ! no, by time interval
 168                     write(msg,'(a,f12.2,a)')'moisture model runs because time ',grid%fmoist_nexttime,'s has arrived'
 169                     call message(msg)
 170                     run_advance_moisture = .true.
 171                 endif
 172             endif
 173             if(run_advance_moisture)then ! decide on timing
 174                 dt_moisture  = moisture_time - grid%fmoist_lasttime  ! Time since moisture model run the last time. Should be long.
 175                 grid%fmoist_lasttime = moisture_time
 176                 if(config_flags%fmoist_freq > 0)then
 177                     write(msg,'(a,f12.2,a,i10,a)')'moisture time step is ',dt_moisture,'s running every ',config_flags%fmoist_freq,' steps'
 178                     call message(msg)
 179                 else
 180                     grid%fmoist_nexttime = moisture_time + config_flags%fmoist_dt
 181                     write(msg,'(a,f12.2,a,f12.2,a)')'moisture time step is ',dt_moisture,'s next run at ',grid%fmoist_nexttime,'s'
 182                     call message(msg)
 183                 endif
 184                 if(fmoist_interp)then
 185                     call message('moisture interpolation to fuels will run because moisture model does')
 186                     run_fuel_moisture=.true.
 187                 endif
 188             endif
 189         endif
 190     elseif(itimestep.eq.1.and.fmoist_interp)then
 191             call message('initializing, moisture interpolation to fuels will run from input data')
 192             run_fuel_moisture=.true.
 193     endif
 194
 195 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
 196     write(msg,'(a,i1,a,i1,a,i4)') &
 197        'fire_driver_em: ifun from ',fire_ifun_start,' to ',fire_ifun_end,' test steps',tsteps
 198 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
 199     call message(msg)
 200
 201     do istep=0,tsteps ! istep >0 is for testing only, exit after the first call
 202       itimestep = grid%itimestep + istep ! in the first call, do fire_test_steps steps of the fire model
 203
 204       need_lfn_update=.false.
 205       do fire_ifun=fire_ifun_start,fire_ifun_end
 206
 207         ! 1 = run pass 1: interpolate height to zsf=terrain, moisture_model initialize
 208         ! 2 = initialize run pass 2: set fuel data, terrain gradient
 209         ! 3 = initialize timestep: interpolate winds, check for ignition, time step on moisture model
 210         ! 4 = do one timestep
 211         ! 5 = copy timestep output to input
 212         ! 6 = compute output fluxes
 213
 214 #ifdef DM_PARALLEL
 215
 216        if(fire_run)then
 217         if(need_lfn_update)then
 218 !           halo exchange on lfn width 3
 219 #include "HALO_FIRE_LFN.inc"
 220         endif
 221
 222         if(fire_ifun.eq.1)then
 223 !       halo exchange on topography
 224 #include "HALO_FIRE_LONGLAT.inc"
 225 !!            if(fire_topo_from_atm.eq.1)then
 226 !!#include "HALO_FIRE_HT.inc"
 227 !!            endif
 228 ! base geopotential and roughness
 229 #include "HALO_FIRE_PHB.inc"
 230
 231 #include "HALO_FIRE_Z0.inc"
 232
 233         elseif(fire_ifun.eq.2)then
 234 !           halo exchange on zsf width 2
 235 #include "HALO_FIRE_ZSF.inc"
 236
 237         elseif(fire_ifun.eq.3)then
 238 !           halo exchange on atm winds and geopotential, width 1 for interpolation
 239 #include "HALO_FIRE_WIND_A.inc"
 240 #include "HALO_FIRE_PH.inc"
 241
 242         elseif(fire_ifun.eq.4)then
 243 !           halo exchange on fire winds width 2 for a 2-step RK method
 244 #include "HALO_FIRE_WIND_F.inc"
 245
 246             if(run_fuel_moisture)then
 247             ! have interpolated to the fire grid
 248 #include "HALO_FIRE_MFG.inc"
 249             endif
 250
 251         elseif(fire_ifun.eq.6)then
 252 !           computing fuel_left needs ignition time from neighbors
 253 #include "HALO_FIRE_TIGN.inc"
 254
 255         endif
 256        endif
 257 #endif
 258
 259
 260         ! need domain by 1 smaller, in last row.col winds are not set properly
 261         call fire_driver_phys ( &
 262             fire_ifun,need_lfn_update,                  &
 263             ids,ide-1, kds,kde, jds,jde-1,                          &
 264             ims,ime, kms,kme, jms,jme,                          &
 265             ips,min(ipe,ide-1), kps,kpe, jps,min(jpe,jde-1),                          &
 266             ifds,ifde-ir, jfds,jfde-jr,                    &
 267             ifms,ifme, jfms,jfme,                    &
 268             ifps,min(ifpe,ifde-ir), jfps,min(jfpe,jfde-jr),      &
 269             ir,jr,                                      & ! atm/fire grid ratio
 270             grid%num_tiles,                             & ! atm grid tiling
 271             grid%i_start,min(grid%i_end,ide-1),                    &
 272             grid%j_start,min(grid%j_end,jde-1),                    &
 273             time_start,                                 &
 274             itimestep,restart,config_flags%fire_fuel_read,config_flags%fire_fuel_cat, &  ! in scalars
 275             grid%dt,grid%dx,grid%dy,                    &
 276             grid%u_frame,grid%v_frame,                  &
 277             unit_fxlong,unit_fxlat,                           & ! coordinates of grid center
 278             config_flags%fire_ext_grnd,config_flags%fire_ext_crwn,config_flags%fire_crwn_hgt, &
 279             config_flags%fire_wind_height,              &  ! height of wind input to fire spread formula
 280             fire_num_ignitions,                                &
 281             fire_ignition_longlat,      &
 282             ignition_line,              &
 283             grid%u_2,grid%v_2,           &          ! atm arrays in
 284             grid%ph_2,grid%phb,               & ! geopotential
 285             grid%z0,                        & ! roughness height
 286             grid%ht,                        &                         ! terrain height
 287             grid%xlong,grid%xlat,                         & ! coordinates of atm grid centers, for ignition location
 288             grid%lfn, &
 289             grid%lfn_hist, &  ! PAJ: to initialize fire perimeter from obs
 290             config_flags%fire_is_real_perim,  & ! PAJ: use obs fire perimeter.
 291             grid%lfn_0,grid%lfn_1,grid%lfn_2,grid%lfn_s0,grid%lfn_s1,grid%lfn_s2,grid%lfn_s3,grid%flame_length,grid%ros_front, &
 292             grid%tign_g,grid%fuel_frac,          & ! state arrays, fire grid
 293             grid%fire_area,                               & ! redundant, for display, fire grid
 294             grid%burnt_area_dt,                               & ! redundant, for display
 295             lfn_out,                                      & ! work - one timestep
 296             grid%avg_fuel_frac,                           & ! out redundant arrays, atm grid
 297             grid%grnhfx,grid%grnqfx,grid%canhfx,grid%canqfx, & ! out redundant arrays, atm grid
 298             grid%grnhfx_fu,grid%grnqfx_fu, & ! out redundant arrays, atm grid
 299             grid%uah,grid%vah,                          &
 300             grid%fgrnhfx,grid%fgrnqfx,grid%fcanhfx,grid%fcanqfx, & ! out redundant arrays, atm grid
 301             grid%ros,                                   & ! rate of spread
 302             grid%fxlong,grid%fxlat,grid%fz0,                           &
 303             grid%nfuel_cat,                               & ! input, or internal for safekeeping
 304             grid%fuel_time,                      &
 305             config_flags%nfmc,         & ! moisture model variables start
 306             run_advance_moisture,run_fuel_moisture,dt_moisture,     &    ! moisture model control
 307             config_flags%fmep_decay_tlag,                               & ! moisture extended model assim. diffs decay time lag
 308             grid%rainc, grid%rainnc,                       & ! accumulated rain from different sources
 309             grid%t2, grid%q2, grid%psfc,               & ! temperature (K), vapor contents (kg/kg), pressure (Pa) at the surface
 310             grid%rain_old,                   & ! previous value of accumulated rain
 311             grid%t2_old, grid%q2_old, grid%psfc_old,   & ! previous values of the atmospheric state at surface
 312             grid%rh_fire,                    & ! relative humidity, diagnostics
 313             grid%fmc_gc,                      & ! fuel moisture fields updated, by class, assumed set to something reasonable
 314             grid%fmep,                      & ! fuel moisture extended model parameters
 315             grid%fmc_equi,                      & ! fuel moisture fields updated, by class, equilibrium diagnostic
 316             grid%fmc_lag,                      & ! fuel moisture fields updated, by class, tendency diagnostic
 317             fp%fmc_g , &                             !  write-only alias. need to exit before using fp again
 318             fp, &
 319             grid, & ! DME added for halo update
 320             ids,ide,jds,jde,kds,kde, &
 321             ims,ime,jms,jme,kms,kme, &
 322             ips,ipe,jps,jpe,kps,kpe, &
 323             config_flags &
 324         )
 325
 326
 327 #ifdef DM_PARALLEL
 328           if(fire_run)then
 329             if(fire_ifun.eq.2)then
 330 !               halo exchange on all fuel data width 2
 331 #include "HALO_FIRE_FUEL.inc"
 332 !               fire state was initialized
 333                 call message('halo exchange on lfn width 2')
 334 #include "HALO_FIRE_LFN.inc"
 335           endif
 336           if(run_fuel_moisture)then
 337             if(fire_ifun.eq.3)then
 338                      ! prepare for interpolation to the fire grid
 339 #include "HALO_FIRE_MAG.inc"
 340                 endif
 341             endif
 342             endif
 343 #endif
 344
 345 ! DME update tracer emissions
 346      if(fire_ifun.eq.6 .AND. config_flags%tracer_opt.eq.3)then
 347        call add_fire_tracer_emissions(config_flags%tracer_opt,grid%dt,grid%dx,grid%dy, &
 348                  ifms,ifme,jfms,jfme,    &
 349                  ifps,ifpe,jfps,jfpe,    &
 350                  ids,ide,kds,kde,jds,jde,          &
 351                  ims,ime,kms,kme,jms,jme,          &
 352                  ips,ipe,kps,kpe,jps,jpe,          &
 353                  rho,dz8w,                         &
 354                  grid%burnt_area_dt,grid%fgip, &
 355                  grid%tracer,config_flags%fire_tracer_smoke)
 356      endif
 357 ! DME
 358       enddo
 359     enddo
 360     if(tsteps>0)call crash('fire_driver_em: test run of uncoupled fire model completed')
 361
 362 end subroutine fire_driver_em
 363
 364 !
 365 !*******************
 366 !
 367
 368 ! module_fr_fire_driver%%fire_driver
 369 subroutine fire_driver_phys (ifun,need_lfn_update,    &
 370     ids,ide, kds,kde, jds,jde,                    & ! atm grid dimensions
 371     ims,ime, kms,kme, jms,jme,                    &
 372     ips,ipe, kps,kpe, jps,jpe,                    &
 373     ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
 374     ifms, ifme, jfms, jfme,                       &
 375     ifps, ifpe, jfps, jfpe,                       & ! fire patch in - will use smaller
 376     ir,jr,                                        & ! atm/fire grid ratio
 377     num_tiles,i_start,i_end,j_start,j_end,        & ! atm grid tiling
 378     time_start,                                   & ! start of the timestep
 379     itimestep,restart,ifuelread,nfuel_cat0,dt,dx,dy,      & ! in scalars
 380     u_frame,v_frame,                              &
 381     unit_fxlong,unit_fxlat,                       & ! fxlong, fxlat units in m
 382     fire_ext_grnd,fire_ext_crwn,fire_crwn_hgt,    &
 383     fire_wind_height,                             & ! for vertical wind interpolation
 384     num_ignitions,                                &
 385     ignition_longlat,                             &
 386     ignition_line,                                &
 387     u,v,                                       & ! in arrays, atm grid
 388     ph,phb,                                       &
 389     z0,zs,                                        &
 390     xlong,xlat,                                   &
 391     lfn,                                      &
 392     lfn_hist,                                      & ! PAJ to init perimeter from obs
 393     fire_is_real_perim,                            & ! PAJ to init perimeter from obs
 394     lfn_0,lfn_1,lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3,flame_length,ros_front, &
 395     tign,fuel_frac,                           & ! state arrays, fire grid
 396     fire_area,                                    & ! redundant state, fire grid
 397     burnt_area_dt,                                    & ! redundant state, fire grid
 398     lfn_out,                                      & ! out level set function
 399     avg_fuel_frac,                                &
 400     grnhfx,grnqfx,canhfx,canqfx,                  & ! out redundant arrays, atm grid
 401     grnhfx_fu,grnqfx_fu,                 & ! out redundant arrays, atm grid
 402     uah,vah,                                      & ! out atm grid
 403     fgrnhfx,fgrnqfx,fcanhfx,fcanqfx,              & ! out redundant arrays, fire grid
 404     ros,                                          &
 405     fxlong,fxlat,fz0,                                 &
 406     nfuel_cat,                                    & ! in array, data, fire grid, or constant internal
 407     fuel_time,                                & ! save constant internal data, fire grid
 408     nfmc,                                     & ! number of fuel moisture classes
 409     run_advance_moisture,run_fuel_moisture,dt_moisture,& ! moisture model control
 410     fmep_decay_tlag,                              & ! moist. extended model assim. diffs time lag
 411     rainc,rainnc,              & ! accumulated rain from different sources
 412     t2, q2, psfc,               & ! temperature (K), vapor contents (kg/kg), pressure (Pa) at the surface
 413     rain_old,                   & ! previous value of accumulated rain
 414     t2_old, q2_old, psfc_old,   & ! previous values of the atmospheric state at surface
 415     rh_fire,                    & ! relative humidity, diagnostics
 416     fmc_gc,                    &  ! fuel moisture fields updated, by class, assumed set to something reasonable
 417     fmep,                       &  ! fuel moisture extended model parameters
 418     fmc_equi,                    &  ! fuel moisture fields updated, by class equilibrium diagnostic
 419     fmc_lag,                    &  ! fuel moisture fields updated, by class tendency diagnostic
 420     fmc_g, &                          ! fuel moisture, alias of fp%fmc_g
 421     fp,                                           & ! fire params
 422     grid,                                         &
 423     ids_hu,ide_hu,jds_hu,jde_hu,kds_hu,kde_hu, &
 424     ims_hu,ime_hu,jms_hu,jme_hu,kms_hu,kme_hu, &
 425     ips_hu,ipe_hu,jps_hu,jpe_hu,kps_hu,kpe_hu, &
 426     config_flags &
 427     )
 428
 429 USE module_dm, only:wrf_dm_maxval
 430 USE module_domain, only: domain
 431
 432 implicit none
 433
 434     type(domain) , target :: grid                             ! state
 435
 436 !*** arguments
 437
 438 integer, intent(in)::ifun,                        &
 439     ids,ide, kds,kde, jds,jde,                    & ! atm domain bounds
 440     ims,ime, kms,kme, jms,jme,                    & ! atm memory bounds
 441     ips,ipe, kps,kpe, jps,jpe,                    & ! atm patch bounds
 442     ifds, ifde, jfds, jfde,                       & ! fire domain bounds
 443     ifms, ifme, jfms, jfme,                       & ! fire memory bounds
 444     ifps, ifpe, jfps, jfpe,                       & ! fire patch bounds
 445     ir,jr,                                        & ! atm/fire grid refinement ratio
 446     nfmc,                                     & ! number of fuel moisture classes
 447     itimestep,                                    & ! number of this timestep
 448     ifuelread,                                    & ! how to initialize nfuel_cat:
 449                                                        ! -1=not at all, done outside
 450                                                        ! 0=from nfuel_cat0
 451                                                        ! 1=from altitude
 452                                                        ! 2=from file
 453     nfuel_cat0,                                   & ! fuel category to initialize everything to
 454     num_tiles                                       ! number of tiles
 455
 456 integer, intent(in)::ids_hu,ide_hu,jds_hu,jde_hu,kds_hu,kde_hu, &
 457                      ims_hu,ime_hu,jms_hu,jme_hu,kms_hu,kme_hu, &
 458                      ips_hu,ipe_hu,jps_hu,jpe_hu,kps_hu,kpe_hu
 459
 460 type (grid_config_rec_type), intent(in)          :: config_flags
 461
 462 logical, intent(in)::restart
 463
 464
 465 logical, intent(out)::need_lfn_update
 466
 467 integer,dimension(num_tiles),intent(in) :: i_start,i_end,j_start,j_end  ! atm grid tiling
 468
 469 real, intent(in):: &
 470     time_start,                                   & ! start of this time step from the beginning of the simulation
 471     dt,                                           & ! time step
 472     dx,dy,                                        & ! atm grid step
 473     u_frame,v_frame,                              & ! velocity offset
 474     unit_fxlong,unit_fxlat,                       & ! fxlong, fxlat units in m
 475     fire_crwn_hgt,                                & ! lowest height crown fire heat is released (m)
 476     fire_ext_grnd,                                & ! extinction depth of surface fire heat flux (m)
 477     fire_ext_crwn,                                & ! wind height for vertical interploation to fire spread
 478     fire_wind_height
 479
 480
 481 integer, intent(in):: num_ignitions                 ! number of ignitions, can be 0
 482 integer, intent(in):: ignition_longlat       ! if 1 ignition give as long/lat, otherwise as m from lower left corner
 483 TYPE (ignition_line_type), DIMENSION(num_ignitions), intent(out):: ignition_line
 484
 485 real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::u,v, & ! wind velocity (m/s) (staggered atm grid)
 486                               ph, phb                      ! geopotential (w-points atm grid)
 487 real,intent(in),dimension(ims:ime, jms:jme)::   z0, &    ! roughness height
 488                                                 zs       ! terrain height
 489 real,intent(out),dimension(ims:ime,jms:jme)::&
 490     uah,                                           & ! atm wind at fire_wind_height, diagnostics
 491     vah                                              ! atm wind at fire_wind_height, diagnostics
 492
 493 real, dimension(ims:ime, jms:jme), intent(inout)::xlong, xlat ! inout because of extension at bdry
 494
 495 real, intent(inout), dimension(ifms:ifme,jfms:jfme):: &
 496     nfuel_cat                                       ! fuel data; can be also set inside (cell based, fire grid)
 497
 498 real, intent(inout), dimension(ifms:ifme, jfms:jfme)::     &
 499     lfn,tign,fuel_frac,                        &     ! state: level function, ign time, fuel left
 500     lfn_hist,                                  &     ! PAJ: to init obs. fira perimeter
 501     lfn_out                                    ! fire wind velocities
 502
 503 real, intent(inout), dimension(ifms:ifme, jfms:jfme)::     &
 504     lfn_0,lfn_1,lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3,flame_length,ros_front  ! level set function and reinitialization steps
 505
 506 real, intent(out), dimension(ifms:ifme, jfms:jfme)::  &
 507     fire_area                                        ! fraction of each cell burning
 508
 509 real, intent(out), dimension(ifms:ifme, jfms:jfme)::  &
 510     burnt_area_dt
 511 real, intent(out), dimension(ims:ime, jms:jme)::  &
 512    grnhfx_fu,grnqfx_fu
 513
 514 real, intent(out), dimension(ims:ime, jms:jme):: &  ! redundant arrays, for display purposes only (atm grid)
 515     avg_fuel_frac,                               &  ! average fuel fraction
 516     grnhfx,                                      &  ! heat flux from surface fire (W/m^2)
 517     grnqfx,                                      &  ! moisture flux from surface fire (W/m^2)
 518     canhfx,                                      &  ! heat flux from crown fire (W/m^2)
 519     canqfx                                         ! moisture flux from crown fire (W/m^2)
 520
 521 real, intent(out), dimension(ifms:ifme, jfms:jfme):: &  ! redundant arrays, for display only, fire grid
 522     fgrnhfx,                                      &  ! heat flux from surface fire (W/m^2)
 523     fgrnqfx,                                      &  ! moisture flux from surface fire (W/m^2)
 524     fcanhfx,                                      &  ! heat flux from crown fire (W/m^2)
 525     fcanqfx,                                      &   ! moisture flux from crown fire (W/m^2)
 526     ros                                             ! fire rate of spread (m/s)
 527
 528 ! moisture model arguments
 529 logical, intent(in)::run_advance_moisture,run_fuel_moisture
 530 real, intent(in)::dt_moisture
 531 real, intent(in)::fmep_decay_tlag
 532 real, intent(in), dimension(ims:ime,jms:jme):: t2, q2, psfc, rainc, rainnc
 533 real, intent(inout), dimension(ims:ime,jms:jme):: t2_old, q2_old, psfc_old, rain_old
 534 real, intent(out),dimension(ims:ime,jms:jme):: rh_fire
 535 real, intent(inout), dimension(ims:ime,nfmc,jms:jme):: fmc_gc
 536 real, intent(inout), dimension(ims:ime,2,jms:jme):: fmep
 537 real, intent(out), dimension(ims:ime,nfmc,jms:jme):: fmc_equi,fmc_lag
 538 real, intent(inout), dimension(ifms:ifme,jfms:jfme):: fmc_g
 539
 540
 541
 542 !  ***** data (constant in time) *****
 543
 544 real, dimension(ifms:ifme, jfms:jfme), intent(inout)::fxlong,fxlat,fz0 ! fire mesh coordinates
 545 real, intent(out), dimension(ifms:ifme, jfms:jfme)::fuel_time   ! fire params arrays
 546
 547 type(fire_params),intent(inout)::fp
 548
 549 !*** local
 550 real :: dxf,dyf,latm, s
 551 integer :: its,ite,jts,jte,kts,kte, &            ! tile
 552     ij,i,j,k,id,pid,ipe1,jpe1,ite1,jte1, &
 553     ifts,ifte,jfts,jfte                          ! fire tile
 554 character(len=128)::msg
 555 character(len=3)::kk
 556 integer::ignitions_done_tile(num_tiles),ignited_tile(num_ignitions,num_tiles)
 557 integer::ignitions_done,ignited_patch(num_ignitions),idex,jdex
 558
 559     ! PAJ:
 560   logical :: fire_is_real_perim
 561
 562 !*** executable
 563
 564 ! fire mesh step
 565 dxf=dx/ir
 566 dyf=dy/jr
 567
 568
 569
 570 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
 571 write(msg,'(a,2f15.6)')'atmosphere mesh step:',dx,dy
 572 call message(msg)
 573 write(msg,'(a,2f15.6)')'fire mesh step:      ',dxf,dyf
 574 call message(msg)
 575 write(msg,7001)'atm domain      ','ids',ids,ide,jds,jde
 576 call message(msg)
 577 write(msg,7001)'atm memory      ','ims',ims,ime,jms,jme
 578 call message(msg)
 579 write(msg,7001)'atm patch       ','ips',ips,ipe,jps,jpe
 580 call message(msg)
 581 write(msg,7001)'fire domain     ','ifds',ifds,ifde,jfds,jfde
 582 call message(msg)
 583 write(msg,7001)'fire memory     ','ifms',ifms,ifme,jfms,jfme
 584 call message(msg)
 585 write(msg,7001)'fire patch      ','ifps',ifps,ifpe,jfps,jfpe
 586 call message(msg)
 587 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
 588
 589 ! check mesh dimensions
 590 call check_fmesh(ids,ide,ifds,ifde,ir,'id')           ! check if atm and fire grids line up
 591 call check_fmesh(jds,jde,jfds,jfde,jr,'jd')
 592 call check_fmesh(ips,ipe,ifps,ifpe,ir,'ip')
 593 call check_fmesh(jps,jpe,jfps,jfpe,jr,'jp')
 594 call check_mesh_2dim(ips,ipe,jps,jpe,ims,ime,jms,jme)        ! check if atm patch fits in atm array
 595 call check_mesh_2dim(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme) ! check if fire patch fits in fire array
 596 call check_mesh_2dim(ips,ipe,jps,jpe,ids,ide,jds,jde)        ! check if atm patch fits in atm domain
 597 call check_mesh_2dim(ifps,ifpe,jfps,jfpe,ifds,ifde,jfds,jfde) ! check if fire patch fits in fire domain
 598
 599 !$OMP SINGLE
 600 ! init rest of fuel tables with derived constants
 601 if(ifun.eq.1) then
 602    ! properties of fuel categories and moisture classes from namelist.fire
 603    call init_fuel_cats(fmoist_run .or. fmoist_interp) ! common for all threads
 604 endif
 605 !$OMP END SINGLE
 606
 607
 608 pid=0
 609 if(fire_print_file.gt.0)then
 610     if(itimestep.le.fire_print_file.or.mod(itimestep,fire_print_file).eq.0)pid=itimestep ! print 1-fire_print_file then every fire_print_file-th
 611 endif
 612
 613 if(ifun.eq.3)then
 614  call print_chsum(itimestep,ims,ime,kms,kme,jms,jme,ids,ide,kds,kde,jds,jde,ips,ipe,kps,kpe,jps,jpe,1,0,0,u,'u')
 615  call print_chsum(itimestep,ims,ime,kms,kme,jms,jme,ids,ide,kds,kde,jds,jde,ips,ipe,kps,kpe,jps,jpe,0,0,1,v,'v')
 616  call print_chsum(itimestep,ims,ime,kms,kme,jms,jme,ids,ide,kds,kde,jds,jde,ips,ipe,kps,kpe,jps,jpe,0,1,0,ph,'ph')
 617 endif
 618
 619 ! fake atm tile bounds
 620 kts=kps
 621 kte=kpe
 622
 623 ! staggered atm patch bounds
 624 ipe1=ifval(ipe.eq.ide,ipe+1,ipe)
 625 jpe1=ifval(jpe.eq.jde,jpe+1,jpe)
 626
 627 ! set up fire tiles & interpolate to fire grid
 628 !$OMP PARALLEL DO PRIVATE(ij,its,ite,jts,jte,ite1,jte1,ifts,ifte,jfts,jfte,msg,id) &
 629 !$OMP SCHEDULE(STATIC)
 630 do ij=1,num_tiles
 631
 632     id = ifval(pid.ne.0,pid+ij*10000,0) ! for print
 633
 634     ignitions_done_tile(ij)=0
 635
 636     ! set up tile bounds
 637     its = i_start(ij)  ! start atmospheric tile in i
 638     ite = i_end(ij)    ! end atmospheric tile in i
 639     jts = j_start(ij)  ! start atmospheric tile in j
 640     jte = j_end(ij)    ! end atmospheric tile in j
 641     ifts= (its-ids)*ir+ifds       ! start fire tile in i
 642     ifte= (ite-ids+1)*ir+ifds-1   ! end fire tile in i
 643     jfts= (jts-jds)*jr+jfds       ! start fire tile in j
 644     jfte= (jte-jds+1)*jr+jfds-1   ! end fire tile in j
 645
 646 ! staggered atm tile bounds
 647     ite1=ifval(ite.eq.ide,ite+1,ite)
 648     jte1=ifval(jte.eq.jde,jte+1,jte)
 649
 650 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
 651     write(msg,'(a,i3,1x,a,i7,1x,a,i3)')'tile=',ij,' id=',id,' ifun=',ifun
 652     call message(msg)
 653     write(msg,7001)'atm tile   ','its',its,ite,jts,jte
 654     call message(msg)
 655     write(msg,7001)'fire tile  ','ifts',ifts,ifte,jfts,jfte
 656     call message(msg)
 657 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
 658
 659     ! check the tiles
 660     call check_mesh_2dim(its,ite,jts,jte,ips,ipe,jps,jpe)                 ! check if atm tile fits in atm patch
 661     call check_mesh_2dim(ifts,ifte,jfts,jfte,ifps,ifpe,jfps,jfpe)         ! check if fire tile fits in fire patch
 662     call check_mesh_2dim(ifts-2,ifte+2,jfts-2,jfte+2,ifms,ifme,jfms,jfme)! check if fire node tile fits in memory
 663
 664
 665 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
 666     write(msg,'(a,i6,a,2(f15.6,a))')'time step',itimestep,' at',time_start,' duration',dt,'s'
 667     call message(msg)
 668     7001 format(a,' dimensions ',a4,':',i6,' to ',i6,' by ',i6,' to ',i6)
 669     write(msg,'(a,2i9)')'refinement ratio:',ir,jr
 670 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
 671
 672     if(run_advance_moisture)then
 673       if(ifun.eq.3)then
 674
 675       ! one timestep of the moisture model
 676           call message('advance_moisture start')
 677           call advance_moisture(    &
 678               itimestep.eq.1,             & ! initialize?
 679               ims,ime,  jms,jme,          & ! memory dimensions
 680               its,ite,  jts,jte,          & ! tile dimensions
 681               nfmc,                       & ! number of moisture fields
 682               dt_moisture,                & ! moisture model time step
 683               fmep_decay_tlag,            & ! moisture extended model assim. diffs decay tlag
 684               rainc, rainnc,              & ! accumulated rain
 685               t2, q2, psfc,               & ! temperature (K), vapor contents (kg/kg), pressure (Pa) at the surface
 686               rain_old,                   & ! previous value of accumulated rain
 687               t2_old, q2_old, psfc_old,   & ! previous values of the atmospheric state at surface
 688               rh_fire,                    & ! relative humidity, diagnostics
 689               fmc_gc,                     & ! fuel moisture fields updated, by class, assumed set to something reasonable
 690               fmep,                       & ! fuel moisture extended model parameters
 691               fmc_equi,                   &  ! fuel moisture fields updated, by class equilibrium diagnostic
 692               fmc_lag                    &  ! fuel moisture fields updated, by class tendency diagnostic
 693           )
 694           call message('advance_moisture end')
 695       endif
 696     endif
 697
 698     if(fire_run)then
 699
 700      if(ifun.eq.1)then   ! set terrain
 701
 702       if(restart)then
 703
 704           call message('restart - topo initialization skipped')
 705
 706       else
 707 !
 708 !        call print_2d_stats(ips,ipe,jps,jpe,ims,ime,jms,jme,zs,'driver:zs')
 709 !
 710 !        ! interpolate terrain height
 711 !        if(fire_topo_from_atm.eq.1)then
 712 !            call interpolate_z2fire(id,1,                 & ! for debug output, <= 0 no output
 713 !                ids,ide,  jds,jde,                    & ! atm grid dimensions
 714 !                ims,ime,  jms,jme,                    &
 715 !                ips,ipe,jps,jpe,                              &
 716 !                its,ite,jts,jte,                              &
 717 !                ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
 718 !                ifms, ifme, jfms, jfme,                       &
 719 !                ifts,ifte,jfts,jfte,                          &
 720 !                ir,jr,                                        & ! atm/fire grid ratio
 721 !                zs,                                       & ! atm grid arrays in
 722 !                fp%zsf)                                      ! fire grid arrays out
 723 !        else
 724 !!$OMP CRITICAL(FIRE_DRIVER_CRIT)
 725 !           write(msg,'(a,i3,a)')'fire_topo_from_atm=',fire_topo_from_atm,' assuming ZSF set, interpolation skipped'
 726 !!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
 727 !        endif
 728
 729         if(ignition_longlat .eq.0)then
 730             ! set ideal fire mesh coordinates - used for ignition only
 731             ! do not forget to set unit_fxlong, unit_fxlat outside of parallel loop
 732
 733             ! DME: Next call added to fixe a bug when
 734             ! initializing a nested domain (before, fxlong & fxlat where not
 735             ! assigned, they were supposed to be set by a mod in WPS
 736             ! but here we use standard WPS therefore were zero and fire does not ignite)
 737             call set_ideal_coord( dxf,dyf, &
 738                 ifds,ifde,jfds,jfde,  &
 739                 ifms,ifme,jfms,jfme,  &
 740                 ifts,ifte,jfts,jfte,  &
 741                 fxlong,fxlat          )
 742             call interpolate_z2fire(id,             &
 743                 ids,ide,jds,jde,                    &
 744                 ims,ime,jms,jme,                    &
 745                 ips,ipe,jps,jpe,                    &
 746                 its,ite,jts,jte,                    &
 747                 ifds,ifde,jfds,jfde,                &
 748                 ifms,ifme,jfms,jfme,                &
 749                 ifts,ifte,jfts,jfte,                &
 750                 ir,jr,                              &
 751                 z0,                                 &
 752                 fz0,1)
 753             !call set_ideal_coord( dx,dy, &
 754             !    ids,ide,jds,jde,  &
 755             !    ims,ime,jms,jme,  &
 756             !    its,ite,jts,jte,  &
 757             !    xlong,xlat          )
 758         elseif(use_atm_vars)then
 759             ! assume halo xlong xlat
 760             ! interpolate nodal coordinates
 761
 762 #ifdef DEBUG_OUT
 763          call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,xlat,'xlat',id)
 764          call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,xlong,'xlong',id)
 765 #endif
 766         call interpolate_z2fire(id,                 & ! for debug output, <= 0 no output
 767             ids,ide,  jds,jde,                    & ! atm grid dimensions
 768             ims,ime,  jms,jme,                    &
 769             ips,ipe,jps,jpe,                              &
 770             its,ite,jts,jte,                              &
 771             ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
 772             ifms, ifme, jfms, jfme,                       &
 773             ifts,ifte,jfts,jfte,                          &
 774             ir,jr,                                        & ! atm/fire grid ratio
 775             xlat,                                       & ! atm grid arrays in
 776             fxlat,0)                                      ! fire grid arrays out
 777
 778         call interpolate_z2fire(id,                 & ! for debug output, <= 0 no output
 779             ids,ide,  jds,jde,                    & ! atm grid dimensions
 780             ims,ime,  jms,jme,                    &
 781             ips,ipe,jps,jpe,                              &
 782             its,ite,jts,jte,                              &
 783             ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
 784             ifms, ifme, jfms, jfme,                       &
 785             ifts,ifte,jfts,jfte,                          &
 786             ir,jr,                                        & ! atm/fire grid ratio
 787             xlong,                                       & ! atm grid arrays in
 788             fxlong,0)                                      ! fire grid arrays out
 789
 790         call interpolate_z2fire(id,                       &
 791             ids,ide,  jds,jde,                            &
 792             ims,ime,  jms,jme,                            &
 793             ips,ipe,jps,jpe,                              &
 794             its,ite,jts,jte,                              &
 795             ifds, ifde, jfds, jfde,                       &
 796             ifms, ifme, jfms, jfme,                       &
 797             ifts,ifte,jfts,jfte,                          &
 798             ir,jr,                                        &
 799             z0,                                           &
 800             fz0,1)
 801
 802         endif
 803
 804      endif
 805
 806     elseif(ifun.eq.2)then
 807
 808         ! after the loop where zsf created exited and all synced
 809         call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fp%zsf,'driver_phys:zsf')
 810
 811     elseif(ifun.eq.3)then  ! interpolate winds to the fire grid
 812
 813
 814       if(use_atm_vars)then
 815
 816         call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,z0,'z0',id)
 817         call write_array_m3(its,ite1,kts,kde-1,jts,jte,ims,ime,kms,kme,jms,jme,u,'u_2',id)
 818         call write_array_m3(its,ite,kts,kde-1,jts,jte1,ims,ime,kms,kme,jms,jme,v,'v_2',id)
 819         call write_array_m3(its,ite,kts,kde,jts,jte,ims,ime,kms,kme,jms,jme,ph,'ph_2',id)
 820         call write_array_m3(its,ite,kts,kde,jts,jte,ims,ime,kms,kme,jms,jme,phb,'phb',id)
 821         call interpolate_atm2fire(id,                     & ! flag for debug output
 822             fire_wind_height,                             & ! height to interpolate to
 823             ids,ide, kds,kde, jds,jde,                    & ! atm grid dimensions
 824             ims,ime, kms,kme, jms,jme,                    &
 825             ips,ipe, jps,jpe,                             &
 826             its,ite,jts,jte,                              &
 827             ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
 828             ifms, ifme, jfms, jfme,                       &
 829             ifps, ifpe, jfps, jfpe,                       & ! fire patch bounds
 830             ifts, ifte, jfts, jfte,                       &
 831             ir,jr,                                        & ! atm/fire grid ratio
 832             u_frame, v_frame,                             & ! velocity frame correction
 833             u,v,                                          & ! 3D atm grid arrays in
 834             ph,phb,                                       &
 835             z0,zs,                                        & ! 2D atm grid arrays in
 836             uah,vah,                                      & ! 2D atm grid out
 837             fp%vx,fp%vy,fz0)                                    ! fire grid arrays out
 838
 839       endif
 840
 841     ! endif ! ifun
 842     ! endif ! fire_run
 843
 844     elseif(ifun.eq.4)then
 845
 846       ! interpolate and compute weighted average to get the fuel moisture
 847       !! print *,'ifun=4, run_fuel_moisture=',run_fuel_moisture
 848       if(run_fuel_moisture)then
 849         call message('fuel_moisture start')
 850         call fuel_moisture(                &
 851         id,                                     & ! for prints and maybe file names
 852         nfmc,                                &
 853         ids,ide, jds,jde,               & ! atm grid dimensions
 854         ims,ime, jms,jme,           &
 855         ips,ipe,jps,jpe,                &
 856         its,ite,jts,jte,                     &
 857         ifds, ifde, jfds, jfde,         & ! fire grid dimensions
 858         ifms, ifme, jfms, jfme,     &
 859         ifts,ifte,jfts,jfte,                 &
 860         ir,jr,                                   & ! atm/fire grid ratio
 861         nfuel_cat,                         & ! fuel data
 862         fmc_gc,                             & ! moisture contents by class on atmospheric grid
 863         fmc_g                                & ! weighted fuel moisture contents on fire grid
 864         )
 865         call message('fuel_moisture end')
 866       endif ! run_fuel_moisture`
 867
 868      endif ! ifun
 869
 870     endif ! fire_run
 871
 872 !   the following executes in any case
 873
 874     if(fire_run)then
 875       call fire_model (id,ifun,restart,need_lfn_update,  &
 876         run_fuel_moisture,                      & ! if fuel moisture needs to be updated
 877         num_ignitions,                          & ! switches
 878         ifuelread,nfuel_cat0,                   & ! initialize fuel categories
 879         ifds,ifde,jfds,jfde,                    & ! fire domain dims
 880         ifms,ifme,jfms,jfme,                    & ! fire memory dims
 881         ifps,ifpe,jfps,jfpe,                    &
 882         ifts,ifte,jfts,jfte,                    & ! fire patch dims
 883         time_start,dt,                          & ! time and increment
 884         dxf,dyf,                                & ! fire mesh spacing
 885         ignition_line,                          & ! description of ignition lines
 886         ignitions_done_tile(ij),ignited_tile(1,ij),  &
 887         fxlong,fxlat,unit_fxlong,unit_fxlat,      & ! fire mesh coordinates
 888         lfn,                                     & ! state: level function
 889         lfn_hist,                                & ! PAJ: to init obs fire perimeter
 890         fire_is_real_perim,                      & ! PAJ: to init obs fire perimeter
 891         lfn_0,lfn_1,lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3,flame_length,ros_front,  & ! state
 892         lfn_out,tign,fuel_frac,                     & ! state: ign time, fuel left
 893         fire_area,                              & ! output: fraction of cell burning
 894         burnt_area_dt,                              &
 895         fgrnhfx,fgrnqfx,                        & ! output: heat fluxes
 896         ros,                                    & ! output: rate of spread for display
 897         nfuel_cat,                              & ! fuel data per point
 898         fuel_time,                              & ! save derived internal data
 899         fp,                                      & ! fire coefficients
 900         grid,            &  ! DME added for halo update
 901         ids_hu,ide_hu,jds_hu,jde_hu,kds_hu,kde_hu, &
 902         ims_hu,ime_hu,jms_hu,jme_hu,kms_hu,kme_hu, &
 903         ips_hu,ipe_hu,jps_hu,jpe_hu,kps_hu,kpe_hu  &
 904     )
 905
 906     if(ifun.eq.6)then ! heat fluxes into the atmosphere
 907
 908         call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fgrnhfx,'fire_driver:fgrnhfx')
 909         call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fgrnqfx,'fire_driver:fgrnqfx')
 910
 911         ! sum the fluxes over atm cells
 912         if(use_atm_vars)then
 913           call sum_2d_cells(        &
 914             ifms,ifme,jfms,jfme,  &
 915             ifts,ifte,jfts,jfte,  &
 916             fuel_frac,              &
 917             ims, ime, jms, jme,   &
 918             its,ite,jts,jte,      &
 919             avg_fuel_frac)
 920           call sum_2d_cells(        &
 921             ifms,ifme,jfms,jfme,  &
 922             ifts,ifte,jfts,jfte,  &
 923             fgrnhfx,              &
 924             ims, ime, jms, jme,   &
 925             its,ite,jts,jte,      &
 926             grnhfx)
 927 !comment out the next call to get results as before commit 55fd92051196b796891b60cb7ec1c4bdb8800078
 928           call sum_2d_cells(        &
 929             ifms,ifme,jfms,jfme,  &
 930             ifts,ifte,jfts,jfte,  &
 931             fgrnqfx,              &
 932             ims, ime, jms, jme,   &
 933             its,ite,jts,jte,      &
 934             grnqfx)
 935
 936 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
 937           write(msg,'(a,f6.3)')'fire-atmosphere feedback scaling ',fire_atm_feedback
 938 !$OMP end CRITICAL(FIRE_DRIVER_CRIT)
 939           call message(msg)
 940           s = 1./(ir*jr)
 941           do j=jts,jte
 942             do i=its,ite
 943                 ! DME heat fluxes contribution for the case wiythout feedback
 944                 grnhfx_fu(i,j)=grnhfx(i,j)*s
 945                 grnqfx_fu(i,j)=grnqfx(i,j)*s
 946                 ! scale surface fluxes to get the averages
 947                 avg_fuel_frac(i,j)=avg_fuel_frac(i,j)*s
 948                 grnhfx(i,j)=fire_atm_feedback*grnhfx(i,j)*s
 949                 grnqfx(i,j)=fire_atm_feedback*grnqfx(i,j)*s
 950                 ! we do not have canopy fluxes yet...
 951                 canhfx(i,j)=0
 952                 canqfx(i,j)=0
 953             enddo
 954           enddo
 955
 956           call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnhfx,'fire_driver:grnhfx')
 957           call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnqfx,'fire_driver:grnqfx')
 958        endif
 959
 960     endif ! ifun=6
 961   endif ! fire_run
 962 enddo ! tiles
 963 !$OMP END PARALLEL DO
 964
 965 #ifdef DEBUG_OUT
 966 if(ifun.eq.1)then
 967     if(pid.ne.0)then
 968         call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,zs,'zs',pid)
 969         call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fp%zsf,'zsf',pid)
 970     endif
 971 endif
 972 #endif
 973
 974 if (ifun.eq.3)then
 975     ignitions_done=ignitions_done_tile(1) ! all should be the same
 976     do i=1,ignitions_done
 977 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
 978         write(msg,'(2(a,i4,1x))')'fire_driver_phys: checking ignition ',i,' of ',ignitions_done
 979 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
 980         call message(msg)
 981         ignited_patch(i)=0
 982         do ij=1,num_tiles
 983             ignited_patch(i)=ignited_patch(i)+ignited_tile(i,ij)
 984         enddo
 985 #ifdef DM_PARALLEL
 986         call message('fire_driver_phys: checking ignitions, collect counts')
 987         call wrf_dm_maxval(ignited_patch(i),idex,jdex)
 988         call message('fire_driver_phys: collected')
 989 #endif
 990         if(ignited_patch(i).eq.0)then
 991             call crash('fire_driver_phys: Ignition failed, no nodes ignited. Bad coordinates?')
 992         endif
 993     enddo
 994
 995  call print_chsum(itimestep,ims,ime,1,1,jms,jme,ids,ide,1,1,jds,jde,ips,ipe,1,1,jps,jpe,1,0,0,uah,'uah')
 996  call print_chsum(itimestep,ims,ime,1,1,jms,jme,ids,ide,1,1,jds,jde,ips,ipe,1,1,jps,jpe,0,0,1,vah,'vah')
 997  call print_chsum(itimestep,ifms,ifme,1,1,jfms,jfme,ifds,ifde,1,1,jfds,jfde,ifps,ifpe,1,1,jfps,jfpe,0,0,0,fp%vx,'uf')
 998  call print_chsum(itimestep,ifms,ifme,1,1,jfms,jfme,ifds,ifde,1,1,jfds,jfde,ifps,ifpe,1,1,jfps,jfpe,0,0,0,fp%vy,'vf')
 999 #ifdef DEBUG_OUT
1000     if(pid.gt.0)then
1001         call write_array_m(ips,ipe1,jps,jpe,ims,ime,jms,jme,uah,'uah',pid)
1002         call write_array_m(ips,ipe,jps,jpe1,ims,ime,jms,jme,vah,'vah',pid)
1003         call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnhfx,'grnhfx',pid)
1004         call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnqfx,'grnqfx',pid)
1005         call write_array_m3(ips,ipe1,kds,kde+1,jps,jpe,ims,ime,kms,kme,jms,jme,u,'u',pid)
1006         call write_array_m3(ips,ipe,kds,kde+1,jps,jpe1,ims,ime,kms,kme,jms,jme,v,'v',pid)
1007         call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fp%vx,'uf',pid)
1008         call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fp%vy,'vf',pid)
1009     endif
1010 #endif
1011 endif
1012
1013 if(ifun.eq.5)then
1014 #ifdef DEBUG_OUT
1015     if(pid.gt.0)then
1016         call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,lfn,'lfn',pid)
1017         call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,tign,'tign',pid)
1018     endif
1019 #endif
1020 endif
1021
1022 if(ifun.eq.6)then
1023     call print_chsum(itimestep,ifms,ifme,1,1,jfms,jfme,ifds,ifde,1,1,jfds,jfde,ifps,ifpe,1,1,jfps,jfpe,0,0,0,fgrnhfx,'fgrnhfx')
1024     call print_chsum(itimestep,ifms,ifme,1,1,jfms,jfme,ifds,ifde,1,1,jfds,jfde,ifps,ifpe,1,1,jfps,jfpe,0,0,0,fgrnqfx,'fgrnqfx')
1025     call print_chsum(itimestep,ims,ime,1,1,jms,jme,ids,ide,1,1,jds,jde,ips,ipe,1,1,jps,jpe,0,0,0,grnhfx,'grnhfx')
1026     call print_chsum(itimestep,ims,ime,1,1,jms,jme,ids,ide,1,1,jds,jde,ips,ipe,1,1,jps,jpe,0,0,0,grnqfx,'grnqfx')
1027 #ifdef DEBUG_OUT
1028     if(pid.gt.0)then
1029         call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnhfx,'grnhfx',pid)
1030         call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnqfx,'grnqfx',pid)
1031         call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fuel_frac,'fuel_frac',pid)
1032         call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fgrnhfx,'fgrnhfx',pid)
1033         call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fgrnqfx,'fgrnqfx',pid)
1034     endif
1035 #endif
1036 endif
1037
1038 end subroutine fire_driver_phys
1039 !
1040 !*******************
1041 !
1042
1043 subroutine fire_ignition_convert (config_flags,fire_max_ignitions,fire_ignition_longlat, &
1044     ignition_line,fire_num_ignitions,unit_fxlong,unit_fxlat)
1045     USE module_configure, only : grid_config_rec_type
1046     implicit none
1047 ! create ignition arrays from scalar flags
1048 !*** arguments
1049     TYPE (grid_config_rec_type) , INTENT(IN)          :: config_flags
1050     integer, intent(in)::fire_max_ignitions
1051     TYPE (ignition_line_type), DIMENSION(fire_max_ignitions), intent(out):: ignition_line ! any values from input discarded
1052     integer, intent(out)::fire_num_ignitions,fire_ignition_longlat
1053     real, intent(out)::unit_fxlong,unit_fxlat
1054 !*** local
1055     integer::i
1056     logical:: real,ideal
1057     real::lat_ctr,lon_ctr
1058 !*** executable
1059     ! this is only until I figure out how to input arrays through the namelist...
1060     if(fire_max_ignitions.lt.5)call crash('fire_max_ignitions too small')
1061     ! figure out which kind of coordinates from the first given
1062     ideal=config_flags%fire_ignition_start_x1 .ne.0. .or. config_flags%fire_ignition_start_y1 .ne. 0.
1063     real=config_flags%fire_ignition_start_lon1 .ne. 0. .or. config_flags%fire_ignition_start_lat1 .ne. 0.
1064     if(ideal)call message('Using ideal ignition coordinates, m from the lower left domain corner')
1065     if(real)call message('Using real ignition coordinates, longitude and latitude')
1066     if(ideal.and.real)call crash('Only one of the ideal or real coordinates may be given')
1067
1068     fire_ignition_longlat=0  ! default, if no ignition
1069     if(ideal)then
1070        ! use values from _x and _y variables
1071        fire_ignition_longlat=0
1072        ignition_line(1)%start_x=config_flags%fire_ignition_start_x1
1073        ignition_line(1)%start_y=config_flags%fire_ignition_start_y1
1074        ignition_line(1)%end_x=config_flags%fire_ignition_end_x1
1075        ignition_line(1)%end_y=config_flags%fire_ignition_end_y1
1076        ignition_line(2)%start_x=config_flags%fire_ignition_start_x2
1077        ignition_line(2)%start_y=config_flags%fire_ignition_start_y2
1078        ignition_line(2)%end_x=config_flags%fire_ignition_end_x2
1079        ignition_line(2)%end_y=config_flags%fire_ignition_end_y2
1080        ignition_line(3)%start_x=config_flags%fire_ignition_start_x3
1081        ignition_line(3)%start_y=config_flags%fire_ignition_start_y3
1082        ignition_line(3)%end_x=config_flags%fire_ignition_end_x3
1083        ignition_line(3)%end_y=config_flags%fire_ignition_end_y3
1084        ignition_line(4)%start_x=config_flags%fire_ignition_start_x4
1085        ignition_line(4)%start_y=config_flags%fire_ignition_start_y4
1086        ignition_line(4)%end_x=config_flags%fire_ignition_end_x4
1087        ignition_line(4)%end_y=config_flags%fire_ignition_end_y4
1088        ignition_line(5)%start_x=config_flags%fire_ignition_start_x5
1089        ignition_line(5)%start_y=config_flags%fire_ignition_start_y5
1090        ignition_line(5)%end_x=config_flags%fire_ignition_end_x5
1091        ignition_line(5)%end_y=config_flags%fire_ignition_end_y5
1092     endif
1093     if(real)then
1094         ! use values from _long and _lat
1095        fire_ignition_longlat=1
1096        ignition_line(1)%start_x=config_flags%fire_ignition_start_lon1
1097        ignition_line(1)%start_y=config_flags%fire_ignition_start_lat1
1098        ignition_line(1)%end_x=config_flags%fire_ignition_end_lon1
1099        ignition_line(1)%end_y=config_flags%fire_ignition_end_lat1
1100        ignition_line(2)%start_x=config_flags%fire_ignition_start_lon2
1101        ignition_line(2)%start_y=config_flags%fire_ignition_start_lat2
1102        ignition_line(2)%end_x=config_flags%fire_ignition_end_lon2
1103        ignition_line(2)%end_y=config_flags%fire_ignition_end_lat2
1104        ignition_line(3)%start_x=config_flags%fire_ignition_start_lon3
1105        ignition_line(3)%start_y=config_flags%fire_ignition_start_lat3
1106        ignition_line(3)%end_x=config_flags%fire_ignition_end_lon3
1107        ignition_line(3)%end_y=config_flags%fire_ignition_end_lat3
1108        ignition_line(4)%start_x=config_flags%fire_ignition_start_lon4
1109        ignition_line(4)%start_y=config_flags%fire_ignition_start_lat4
1110        ignition_line(4)%end_x=config_flags%fire_ignition_end_lon4
1111        ignition_line(4)%end_y=config_flags%fire_ignition_end_lat4
1112        ignition_line(5)%start_x=config_flags%fire_ignition_start_lon5
1113        ignition_line(5)%start_y=config_flags%fire_ignition_start_lat5
1114        ignition_line(5)%end_x=config_flags%fire_ignition_end_lon5
1115        ignition_line(5)%end_y=config_flags%fire_ignition_end_lat5
1116     endif
1117     ! common to both cases
1118        ignition_line(1)%ros=config_flags%fire_ignition_ros1
1119        ignition_line(1)%radius=config_flags%fire_ignition_radius1
1120        ignition_line(1)%start_time=config_flags%fire_ignition_start_time1
1121        ignition_line(1)%end_time=config_flags%fire_ignition_end_time1
1122        ignition_line(2)%ros=config_flags%fire_ignition_ros2
1123        ignition_line(2)%radius=config_flags%fire_ignition_radius2
1124        ignition_line(2)%start_time=config_flags%fire_ignition_start_time2
1125        ignition_line(2)%end_time=config_flags%fire_ignition_end_time2
1126        ignition_line(3)%ros=config_flags%fire_ignition_ros3
1127        ignition_line(3)%radius=config_flags%fire_ignition_radius3
1128        ignition_line(3)%start_time=config_flags%fire_ignition_start_time3
1129        ignition_line(3)%end_time=config_flags%fire_ignition_end_time3
1130        ignition_line(4)%ros=config_flags%fire_ignition_ros4
1131        ignition_line(4)%radius=config_flags%fire_ignition_radius4
1132        ignition_line(4)%start_time=config_flags%fire_ignition_start_time4
1133        ignition_line(4)%end_time=config_flags%fire_ignition_end_time4
1134        ignition_line(5)%ros=config_flags%fire_ignition_ros5
1135        ignition_line(5)%radius=config_flags%fire_ignition_radius5
1136        ignition_line(5)%start_time=config_flags%fire_ignition_start_time5
1137        ignition_line(5)%end_time=config_flags%fire_ignition_end_time5
1138
1139     !
1140         fire_num_ignitions=0
1141         do i=1,min(5,config_flags%fire_num_ignitions)
1142             ! count the ignitions
1143             if(ignition_line(i)%radius.gt.0.)fire_num_ignitions=i
1144             ! expand ignition data given as zero
1145             if(ignition_line(i)%end_x.eq.0.)ignition_line(i)%end_x=ignition_line(i)%start_x
1146             if(ignition_line(i)%end_y.eq.0.)ignition_line(i)%end_y=ignition_line(i)%start_y
1147             if(ignition_line(i)%end_time.eq.0.)ignition_line(i)%end_time=ignition_line(i)%start_time
1148         enddo
1149
1150     if(fire_ignition_longlat .eq. 0)then
1151        ! ideal
1152        !  ignition is in m
1153        unit_fxlong=1.
1154        unit_fxlat=1.
1155        ! will set fire mesh coordinates to uniform mesh below
1156     else
1157        ! real
1158        lat_ctr=config_flags%cen_lat
1159        lon_ctr=config_flags%cen_lon
1160        ! 1 degree in m (approximate OK)
1161        unit_fxlat=pi2/(360.*reradius)  ! earth circumference in m / 360 degrees
1162        unit_fxlong=cos(lat_ctr*pi2/360.)*unit_fxlat  ! latitude
1163     endif
1164
1165 end subroutine fire_ignition_convert
1166
1167 !
1168 !*****************************
1169 !
1170
1171 subroutine interpolate_atm2fire(id,               & ! for debug output, <= 0 no output
1172     fire_wind_height,                             & ! interpolation height
1173     ids,ide, kds,kde, jds,jde,                    & ! atm grid dimensions
1174     ims,ime, kms,kme, jms,jme,                    &
1175     ips,ipe,jps,jpe,                              &
1176     its,ite,jts,jte,                              &
1177     ifds, ifde, jfds, jfde,                       & ! fire grid dimensions
1178     ifms, ifme, jfms, jfme,                       &
1179     ifps, ifpe, jfps, jfpe,                       & ! fire patch bounds
1180     ifts,ifte,jfts,jfte,                          &
1181     ir,jr,                                        & ! atm/fire grid ratio
1182     u_frame, v_frame,                             & ! velocity frame correction
1183     u,v,                                          & ! atm grid arrays in
1184     ph,phb,                                       &
1185     z0,zs,                                        &
1186     uah,vah,                                      &
1187     uf,vf,z0f)                                          ! fire grid arrays out
1188
1189 implicit none
1190 !*** purpose: interpolate winds and height
1191
1192 !*** arguments
1193 integer, intent(in)::id
1194 real, intent(in):: fire_wind_height                 ! height above the terrain for vertical interpolation
1195 integer, intent(in)::                             &
1196     ids,ide, kds,kde, jds,jde,                    & ! atm domain bounds
1197     ims,ime, kms,kme, jms,jme,                    & ! atm memory bounds
1198     ips,ipe,jps,jpe,                              &
1199     its,ite,jts,jte,                              & ! atm tile bounds
1200     ifds, ifde, jfds, jfde,                       & ! fire domain bounds
1201     ifms, ifme, jfms, jfme,                       & ! fire memory bounds
1202     ifps, ifpe, jfps, jfpe,                       & ! fire patch bounds
1203     ifts,ifte,jfts,jfte,                          & ! fire tile bounds
1204     ir,jr                                         ! atm/fire grid refinement ratio
1205 real, intent(in):: u_frame, v_frame                 ! velocity frame correction
1206 real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::&
1207     u,v,                                          & ! atm wind velocity, staggered
1208     ph,phb                                          ! geopotential
1209 real,intent(in),dimension(ims:ime,jms:jme)::&
1210     z0,                                           & ! roughness height
1211     zs                                              ! terrain height
1212 real,intent(out),dimension(ims:ime,jms:jme)::&
1213     uah,                                           & ! atm wind at fire_wind_height, diagnostics
1214     vah                                              !
1215 real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
1216     uf,vf                                           ! wind velocity fire grid nodes
1217 real,intent(in),dimension(ifms:ifme,jfms:jfme)::z0f ! roughness length in fire grid
1218
1219
1220 !*** local
1221 character(len=256)::msg
1222 #define TDIMS its-2,ite+2,jts-2,jte+2
1223 real, dimension(its-2:ite+2,jts-2:jte+2):: ua,va   ! atm winds, interpolated over height, still staggered grid
1224 real, dimension(its-2:ite+2,kds:kde,jts-2:jte+2):: altw,altub,altvb,hgtu,hgtv ! altitudes
1225 integer:: i,j,k,ifts1,ifte1,jfts1,jfte1,ite1,jte1
1226 integer::itst,itet,jtst,jtet,itsu,iteu,jtsu,jteu,itsv,itev,jtsv,jtev
1227 integer::kdmax,its1,jts1,ips1,jps1
1228 integer::itsou,iteou,jtsou,jteou,itsov,iteov,jtsov,jteov
1229 real:: ground,loght,loglast,logz0,logfwh,ht,zr
1230 real::r_nan
1231 integer::i_nan
1232 equivalence (i_nan,r_nan)
1233 real::fire_wind_height_local,z0fc
1234 real::ust_d,wsf,wsf1,uf_temp,vf_temp
1235 real,parameter::vk_kappa=0.4
1236
1237 !*** executable
1238
1239 ! debug init local arrays
1240 i_nan=2147483647
1241 ua=r_nan
1242 va=r_nan
1243 altw=r_nan
1244 altub=r_nan
1245 hgtu=r_nan
1246 hgtv=r_nan
1247
1248
1249 if(kds.ne.1)then
1250 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
1251   write(msg,*)'WARNING: bottom index kds=',kds,' should be 1?'
1252   call message(msg)
1253 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
1254 endif
1255
1256 !                            ^ j
1257 !        ------------        |
1258 !        |          |         ----> i
1259 !        u    p     |
1260 !        |          |    nodes in cell (i,j)
1261 !        ------v-----    view from top
1262 !
1263 !             v(ide,jde+1)
1264 !            -------x------
1265 !            |            |
1266 !            |            |
1267 ! u(ide,jde) x      x     x u(ide+1,jde)
1268 !            | p(ide,hde) |
1269 !            |            |   p=ph,phb,z0,...
1270 !            -------x------
1271 !              v(ide,jde)
1272 !
1273 ! staggered values set u(ids:ide+1,jds:jde) v(ids:ide,jds:jde+1)
1274 ! p=ph+phb set at (ids:ide,jds:jde)
1275 ! location of u(i,j) needs p(i-1,j) and p(i,j)
1276 ! location of v(i,j) needs p(i,j-1) and p(i,j)
1277 ! *** NOTE need HALO in ph, phb
1278 ! so we can compute only u(ids+1:ide,jds:jde) v(ids:ide,jds+1,jde)
1279 ! unless we extend p at the boundary
1280 ! but because we care about the fire way in the inside it does not matter
1281 ! if the fire gets close to domain boundary the simulation is over anyway
1282
1283     ite1=snode(ite,ide,1)
1284     jte1=snode(jte,jde,1)
1285     ! do this in any case to check for nans
1286     call print_3d_stats(its,ite1,kds,kde,jts,jte,ims,ime,kms,kme,jms,jme,u,'wind U in')
1287     call print_3d_stats(its,ite,kds,kde,jts,jte1,ims,ime,kms,kme,jms,jme,v,'wind V in')
1288
1289     if(fire_print_msg.gt.0)then
1290 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
1291        write(msg,'(a,f7.2,a)')'interpolate_atm2fire: log-interpolation of wind to',fire_wind_height,' m'
1292        call message(msg)
1293 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
1294     endif
1295
1296     ! indexing
1297
1298     ! file for u
1299     itst=ifval(ids.eq.its,its,its-1)
1300     itet=ifval(ide.eq.ite,ite,ite+1)
1301     jtst=ifval(jds.ge.jts,jts,jts-1)
1302     jtet=ifval(jde.eq.jte,jte,jte+1)
1303
1304 if(fire_print_msg.ge.1)then
1305 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
1306   write(msg,7001)'atm input  ','tile',its,ite,jts,jte
1307   call message(msg)
1308   write(msg,7001)'altw       ','tile',itst,itet,jtst,jtet
1309   call message(msg)
1310 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
1311 endif
1312 7001 format(a,' dimensions ',a4,':',i6,' to ',i6,' by ',i6,' to ',i6)
1313
1314     kdmax=kde-1   ! max layer to interpolate from, can be less
1315
1316     do j = jtst,jtet
1317       do k=kds,kdmax+1
1318         do i = itst,itet
1319           altw(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g             ! altitude of the bottom w-point
1320         enddo
1321       enddo
1322     enddo
1323
1324 ! values at u points
1325     itsu=ifval(ids.eq.its,its+1,its)  ! staggered direction
1326     iteu=ifval(ide.eq.ite,ite,ite+1)  ! staggered direction
1327     jtsu=ifval(jds.ge.jts,jts,jts-1)
1328     jteu=ifval(jde.eq.jte,jte,jte+1)
1329
1330 if(fire_print_msg.ge.1)then
1331 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
1332   write(msg,7001)'u interp at','tile',itsu,iteu,jtsu,jteu
1333   call message(msg)
1334 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
1335 endif
1336
1337     do j = jtsu,jteu
1338       do k=kds,kdmax+1
1339         do i = itsu,iteu
1340           altub(i,k,j)= 0.5*(altw(i-1,k,j)+altw(i,k,j))      ! altitude of the bottom point under u-point
1341         enddo
1342       enddo
1343       do k=kds,kdmax
1344         do i = itsu,iteu
1345           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
1346         enddo
1347       enddo
1348     enddo
1349
1350 ! values at v points
1351     jtsv=ifval(jds.eq.jts,jts+1,jts)  ! staggered direction
1352     jtev=ifval(jde.eq.jte,jte,jte+1)  ! staggered direction
1353     itsv=ifval(ids.ge.its,its,its-1)
1354     itev=ifval(ide.eq.ite,ite,ite+1)
1355
1356 if(fire_print_msg.ge.1)then
1357 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
1358   write(msg,7001)'v interp at','tile',itsv,itev,jtsv,jtev
1359   call message(msg)
1360 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
1361 endif
1362     do j = jtsv,jtev
1363       do k=kds,kdmax+1
1364         do i = itsv,itev
1365           altvb(i,k,j)= 0.5*(altw(i,k,j-1)+altw(i,k,j))      ! altitude of the bottom point under v-point
1366         enddo
1367       enddo
1368       do k=kds,kdmax
1369         do i = itsv,itev
1370           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
1371         enddo
1372       enddo
1373     enddo
1374
1375 #ifdef DEBUG_OUT
1376         call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,altub,'altub',id)
1377         call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,altvb,'altvb',id)
1378         call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,hgtu,'hgtu',id)
1379         call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,hgtv,'hgtv',id)
1380 #endif
1381 ! DME
1382 if (fire_lsm_zcoupling) then
1383   logfwh = log(fire_lsm_zcoupling_ref)
1384   fire_wind_height_local = fire_lsm_zcoupling_ref
1385 else
1386   logfwh = log(fire_wind_height)
1387   fire_wind_height_local = fire_wind_height
1388 endif
1389
1390     ! interpolate u, staggered in X
1391
1392     do j = jtsu,jteu            ! compute on domain by 1 smaller, otherwise z0 and ph not available
1393 !      do i = itsu,iteu        ! compute with halo 2
1394       do i = itsu,iteu
1395         zr = 0.5*(z0(i,j)+z0(i-1,j)) ! interpolated roughness length under this u point
1396         if(fire_wind_height_local > zr)then       !
1397           do k=kds,kdmax
1398             ht = hgtu(i,k,j)      ! height of this u point above the ground
1399             if( .not. ht < fire_wind_height_local) then ! found layer k this point is in
1400               loght = log(ht)
1401               if(k.eq.kds)then               ! first layer, log linear interpolation from 0 at zr
1402                 logz0 = log(zr)
1403                 ua(i,j)= u(i,k,j)*(logfwh-logz0)/(loght-logz0)
1404               else                           ! log linear interpolation
1405                 loglast=log(hgtu(i,k-1,j))
1406                 ua(i,j)= u(i,k-1,j) + (u(i,k,j) - u(i,k-1,j)) * ( logfwh - loglast) / (loght - loglast)
1407               endif
1408               goto 10
1409             endif
1410             if(k.eq.kdmax)then                 ! last layer, still not high enough
1411               ua(i,j)=u(i,k,j)
1412             endif
1413           enddo
1414 10        continue
1415         else  ! roughness higher than the fire wind height
1416           ua(i,j)=0.
1417         endif
1418       enddo
1419     enddo
1420
1421     do j = jtsu,jteu
1422       ua(itsu-1,j)=ua(itsu,j)
1423     enddo
1424
1425
1426     ! interpolate v, staggered in Y
1427
1428     do j = jtsv,jtev
1429       do i = itsv,itev
1430         zr = 0.5*(z0(i,j-1)+z0(i,j)) ! roughness length under this v point
1431         if(fire_wind_height_local > zr)then       !
1432           do k=kds,kdmax
1433             ht = hgtv(i,k,j)      ! height of this u point above the ground
1434             if( .not. ht < fire_wind_height_local) then ! found layer k this point is in
1435               loght = log(ht)
1436               if(k.eq.kds)then               ! first layer, log linear interpolation from 0 at zr
1437                 logz0 = log(zr)
1438                 va(i,j)= v(i,k,j)*(logfwh-logz0)/(loght-logz0)
1439               else                           ! log linear interpolation
1440                 loglast=log(hgtv(i,k-1,j))
1441                 va(i,j)= v(i,k-1,j) + (v(i,k,j) - v(i,k-1,j)) * ( logfwh - loglast) / (loght - loglast)
1442               endif
1443               goto 11
1444             endif
1445             if(k.eq.kdmax)then                 ! last layer, still not high enough
1446               va(i,j)=v(i,k,j)
1447             endif
1448           enddo
1449 11        continue
1450         else  ! roughness higher than the fire wind height
1451           va(i,j)=0.
1452         endif
1453       enddo
1454     enddo
1455
1456     do i = itsv,itev
1457       va(i,jtsv-1)=va(i,jtsv)
1458     enddo
1459
1460 #ifdef DEBUG_OUT
1461 !   store the output for diagnostics
1462     do j = jts,jte1
1463       do i = its,ite1
1464         uah(i,j)=ua(i,j)
1465         vah(i,j)=va(i,j)
1466       enddo
1467     enddo
1468
1469     call write_array_m(its,ite1,jts,jte,ims,ime,jms,jme,uah,'uah_n',id) ! no reflection
1470     call write_array_m(its,ite,jts,jte1,ims,ime,jms,jme,vah,'vah_n',id)
1471 #endif
1472
1473     ips1 = ifval(ips.eq.ids,ips+1,ips)
1474
1475     call continue_at_boundary(1,1,0., & ! x direction
1476        TDIMS,                  &! memory dims atm grid tile
1477        ids+1,ide,jds,jde, &     ! domain dims - where u defined
1478        ips1,ipe,jps,jpe, &     ! patch dims
1479        itsu,iteu,jtsu,jteu, & ! tile dims - in nonextended direction one beyond if at patch boundary but not domain
1480        itsou,iteou,jtsou,jteou, & ! out, where set
1481        ua)                           ! array
1482
1483     jps1 = ifval(jps.eq.jds,jps+1,jps)
1484
1485     call continue_at_boundary(1,1,0., & ! y direction
1486        TDIMS,                  & ! memory dims atm grid tile
1487        ids,ide,jds+1,jde, &      ! domain dims - where v wind defined
1488        ips,ipe,jps1,jpe, &        ! patch dims
1489        itsv,itev,jtsv,jtev, & ! tile dims
1490        itsov,iteov,jtsov,jteov, & ! where set
1491        va)                           ! array
1492
1493 !   store the output for diagnostics
1494     do j = jts,jte1
1495       do i = its,ite1
1496         uah(i,j)=ua(i,j)
1497         vah(i,j)=va(i,j)
1498       enddo
1499     enddo
1500
1501 #ifdef DEBUG_OUT
1502         call write_array_m(itsou,iteou,jtsou,jteou,TDIMS,ua,'ua',id)
1503         call write_array_m(itsov,iteov,jtsov,jteov,TDIMS,va,'va',id)
1504 #endif
1505
1506 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
1507     ! don't have all values valid, don't check
1508      write(msg,12)'atm mesh wind U at',fire_wind_height,' m'
1509      call print_2d_stats(itsou,iteou,jtsou,jteou,TDIMS,ua,msg)
1510      write(msg,12)'atm mesh wind V at',fire_wind_height,' m'
1511      call print_2d_stats(itsov,iteov,jtsov,jteov,TDIMS,va,msg)
1512 12  format(a,f6.2,a)
1513     call print_2d_stats(its,ite1,jts,jte,ims,ime,jms,jme,uah,'UAH')
1514     call print_2d_stats(its,ite,jts,jte1,ims,ime,jms,jme,vah,'VAH')
1515     !call write_array_m(its,ite1,jts,jte,ims,ime,jms,jme,uah,'uah',id)
1516     !call write_array_m(its,ite,jts,jte1,ims,ime,jms,jme,vah,'vah',id)
1517 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
1518
1519 !      ---------------
1520 !     | F | F | F | F |   Example of atmospheric and fire grid with
1521 !     |-------|-------|   ir=jr=4.
1522 !     | F | F | F | F |   Winds are given at the midpoints of the sides of the atmosphere grid,
1523 !     ua------z-------|   interpolated to midpoints of the cells of the fine fire grid F.
1524 !     | F | F | F | F |   This is (1,1) cell of atmosphere grid, and [*] is the (1,1) cell of the fire grid.
1525 !     |---------------|   ua(1,1) <--> uf(0.5,2.5)
1526 !     | * | F | F | F |   va(1,1) <--> vf(2.5,0.5)
1527 !      -------va------    za(1,1) <--> zf(2.5,2.5)
1528 !
1529 !   ^ x2
1530 !   |  --------va(1,2)---------
1531 !   | |            |           |   Example of atmospheric and fire grid with
1532 !   | |            |           |   ir=jr=1.
1533 !   | |          za,zf         |   Winds are given at the midpoints of the sides of the atmosphere grid,
1534 !   | ua(1,1)----uf,vf-----ua(2,1) interpolated to midpoints of the cells of the (the same) fire grid
1535 !   | |           (1,1)        |   ua(1,1) <--> uf(0.5,1)
1536 !   | |            |           |   va(1,1) <--> vf(1,0.5)
1537 !   | |            |           |   za(1,1) <--> zf(1,1)
1538 !   |  --------va(1,1)---------
1539 !   |--------------------> x1
1540 !
1541 ! Meshes are aligned by the lower left cell of the domain. Then in the above figure
1542 ! u = node with the ua component of the wind at (ids,jds), midpoint of side
1543 ! v = node with the va component of the wind at (ids,jds), midpoint of side
1544 ! * = fire grid node at (ifds,jfds)
1545 ! z = node with height, midpoint of cell
1546 !
1547 ! ua(ids,jds)=uf(ifds-0.5,jfds+jr*0.5-0.5)         = uf(ifds-0.5,jfds+(jr-1)*0.5)
1548 ! va(ids,jds)=vf(ifds+ir*0.5-0.5,jfds-0.5)         = vf(ifds+(ir-1)*0.5,jfds-0.5)
1549 ! 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)
1550
1551     ! ifts1=max(snode(ifts,ifps,-1),ifds) ! go 1 beyond patch boundary but not at domain boundary
1552     ! ifte1=min(snode(ifte,ifpe,+1),ifde)
1553     ! jfts1=max(snode(jfts,jfps,-1),jfds)
1554     ! jfte1=min(snode(jfte,jfpe,+1),jfde)
1555
1556     call interpolate_2d(  &
1557         TDIMS,                  & ! memory dims atm grid tile
1558         itsou,iteou,jtsou,jteou,& ! where set
1559         ifms,ifme,jfms,jfme,    & ! array dims fire grid
1560         ifts,ifte,jfts,jfte,& ! dimensions on the fire grid to interpolate to
1561         ir,jr,                  & ! refinement ratio
1562         real(ids),real(jds),ifds-0.5,jfds+(jr-1)*0.5, & ! line up by lower left corner of domain
1563         ua,                     & ! in atm grid
1564         uf)                      ! out fire grid
1565
1566     call interpolate_2d(  &
1567         TDIMS,                  & ! memory dims atm grid tile
1568         itsov,iteov,jtsov,jteov,& ! where set
1569         ifms,ifme,jfms,jfme,    & ! array dims fire grid
1570         ifts,ifte,jfts,jfte,& ! dimensions on the fire grid to interpolate to
1571         ir,jr,                  & ! refinement ratio
1572         real(ids),real(jds),ifds+(ir-1)*0.5,jfds-0.5, & ! line up by lower left corner of domain
1573         va,                     & ! in atm grid
1574         vf)                      ! out fire grid
1575
1576 ! DME here code to extrapolate mid-flame height velocity -> fire_lsm_zcoupling = .true.
1577 if (fire_lsm_zcoupling) then
1578   do j = jfts,jfte
1579     do i = ifts,ifte
1580       uf_temp=uf(i,j)
1581       vf_temp=vf(i,j)
1582       wsf=max(sqrt(uf_temp**2.+vf_temp**2.),0.1)
1583       z0fc=z0f(i,j)
1584       ust_d=wsf*vk_kappa/log(fire_lsm_zcoupling_ref/z0fc)
1585       wsf1=(ust_d/vk_kappa)*log((fire_wind_height+z0fc)/z0fc)
1586       uf(i,j)=wsf1*uf_temp/wsf
1587       vf(i,j)=wsf1*vf_temp/wsf
1588     enddo
1589   enddo
1590 endif
1591
1592
1593 call print_2d_stats_vec(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,vf,'fire wind (m/s)')
1594 ! call print_2d_stats_vec(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,vf,'fire wind extended')
1595 #ifdef DEBUG_OUT
1596         call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,'uf1',id)
1597         call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,vf,'vf1',id)
1598         ! call write_array_m(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,'uf1',id)
1599         ! call write_array_m(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,vf,'vf1',id)
1600 #endif
1601 return
1602
1603 end subroutine interpolate_atm2fire
1604
1605 !
1606 !***
1607 !
1608
1609 subroutine check_fmesh(ids,ide,ifds,ifde,ir,s)
1610 !*** purpose: check if fire and atm meshes line up
1611 implicit none
1612 !*** arguments
1613 integer, intent(in)::ids,ide,ifds,ifde,ir
1614 character(len=*),intent(in)::s
1615 !*** local
1616 character(len=128)msg
1617 !*** executable
1618 if ((ide-ids+1)*ir.ne.(ifde-ifds+1))then
1619 !$OMP CRITICAL(FIRE_DRIVER_CRIT)
1620     write(msg,1)s,ids,ide,ifds,ifde,ir
1621 1   format('module_fr_fire_driver: incompatible bounds ',a,' atm ',i5,':',i5,' fire ',i5,':',i5,' ratio ',i3)
1622 !$OMP END CRITICAL(FIRE_DRIVER_CRIT)
1623     call crash(msg)
1624 endif
1625 end subroutine check_fmesh
1626
1627 !
1628 !*****************************
1629 !
1630 subroutine set_flags(config_flags)
1631 implicit none
1632 TYPE (grid_config_rec_type) , INTENT(IN)          :: config_flags
1633 ! copy flags from wrf to module_fr_fire_util
1634 ! for instructions how to add a flag see the top of module_fr_fire_util.F
1635
1636
1637 fire_print_msg          = config_flags%fire_print_msg
1638 fire_print_file         = config_flags%fire_print_file
1639 fuel_left_method        = config_flags%fire_fuel_left_method
1640 fuel_left_irl           = config_flags%fire_fuel_left_irl
1641 fuel_left_jrl           = config_flags%fire_fuel_left_jrl
1642 fire_const_time         = config_flags%fire_const_time
1643 fire_const_grnhfx       = config_flags%fire_const_grnhfx
1644 fire_const_grnqfx       = config_flags%fire_const_grnqfx
1645 fire_atm_feedback       = config_flags%fire_atm_feedback
1646 boundary_guard          = config_flags%fire_boundary_guard
1647 fire_grows_only         = config_flags%fire_grows_only
1648 fire_upwinding          = config_flags%fire_upwinding
1649 fire_upwind_split       = config_flags%fire_upwind_split
1650 fire_viscosity          = config_flags%fire_viscosity
1651 fire_lfn_ext_up         = config_flags%fire_lfn_ext_up
1652 fire_test_steps         = config_flags%fire_test_steps
1653 fire_advection          = config_flags%fire_advection
1654 fire_lsm_reinit         = config_flags%fire_lsm_reinit
1655 fire_lsm_reinit_iter    = config_flags%fire_lsm_reinit_iter
1656 fire_upwinding_reinit   = config_flags%fire_upwinding_reinit
1657 fire_lsm_band_ngp       = config_flags%fire_lsm_band_ngp
1658 fire_lsm_zcoupling      = config_flags%fire_lsm_zcoupling
1659 fire_lsm_zcoupling_ref  = config_flags%fire_lsm_zcoupling_ref
1660 fire_viscosity_bg       = config_flags%fire_viscosity_bg
1661 fire_viscosity_band     = config_flags%fire_viscosity_band
1662 fire_viscosity_ngp      = config_flags%fire_viscosity_ngp
1663 fire_slope_factor       = config_flags%fire_slope_factor
1664 fire_fmc_read           = config_flags%fire_fmc_read
1665
1666
1667
1668
1669 end subroutine set_flags
1670
1671 end module module_fr_fire_driver