!
! This module is the entry point from WRF ARW to the wildland
! fire module. The call to fire_driver advances the fire module by
! one timestep. The fire module inputs the wind and outputs
! temperature and humidity tendencies. The fire module also inputs a
! number of constant arrays (fuel data, topography). Additional
! arguments are model state (for data assimilation) and constant arrays
! the model gives to WRF for safekeeping because it is not allowed
! to save anything.
!
! Contributions to this wildland fire module have come from individuals at
! NCAR, the U.S.D.A. Forest Service, the Australian Bureau of Meteorology,
! and the University of Colorado at Denver.
!
module module_fr_fire_driver 1
! use this module for standalone call, you only need to provide some mock-up wrf modules
use module_fr_fire_model
use module_fr_fire_phys, only : fire_params , init_fuel_cats
use module_fr_fire_util
use module_fr_fire_core, only: ignition_line_type
USE module_domain, only: domain
USE module_configure, only: grid_config_rec_type
use module_model_constants, only: reradius, & ! 1/earth radiusw
g, & ! gravitation acceleration
pi2 ! 2*pi
implicit none
private
public fire_driver_em,use_atm_vars
logical:: use_atm_vars=.true. ! interpolate wind from atm mesh and average output fluxes back
contains
#define DEBUG_OUT
subroutine fire_driver_em ( grid , config_flags & 2,7
,fire_ifun_start,fire_ifun_end,tsteps &
,ids,ide, kds,kde, jds,jde &
,ims,ime, kms,kme, jms,jme &
,ips,ipe, kps,kpe, jps,jpe &
,ifds,ifde, jfds,jfde &
,ifms,ifme, jfms,jfme &
,ifps,ifpe, jfps,jfpe )
!*** purpose: driver from grid structure
! Driver layer modules
#ifdef DM_PARALLEL
USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks
USE module_comm_dm , ONLY : halo_fire_fuel_sub, halo_fire_tign_sub, halo_fire_wind_f_sub, &
halo_fire_wind_a_sub, halo_fire_ph_sub, halo_fire_zsf_sub, halo_fire_longlat_sub, &
halo_fire_phb_sub, halo_fire_z0_sub, halo_fire_lfn_sub
#endif
implicit none
!*** arguments
TYPE(domain) , TARGET :: grid ! state
TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags ! namelist
integer, intent(in):: fire_ifun_start,fire_ifun_end,tsteps ! driver cycle controls
integer, intent(in):: &
ids,ide, kds,kde, jds,jde, &
ims,ime, kms,kme, jms,jme, &
ips,ipe, kps,kpe, jps,jpe, &
ifds,ifde, jfds,jfde, &
ifms,ifme, jfms,jfme, &
ifps,ifpe, jfps,jfpe
!*** local
INTEGER:: fire_num_ignitions
integer, parameter::fire_max_ignitions=5
TYPE(ignition_line_type), DIMENSION(fire_max_ignitions):: ignition_line
integer::fire_ifun,ir,jr,fire_ignition_longlat,istep,itimestep
logical::need_lfn_update,restart
real, dimension(ifms:ifme, jfms:jfme)::lfn_out
real:: corner_ll,corner_ul,corner_ur,corner_lr
character(len=128)msg
real:: unit_fxlong ,unit_fxlat
type(fire_params)::fp
!*** executable
! populate our structures from wrf
! pointers to be passed to fire rate of spread formulas
fp%vx => grid%uf ! W-E winds used in fire module
fp%vy => grid%vf ! S-N winds used in fire module
fp%zsf => grid%zsf ! terrain height
fp%dzdxf => grid%dzdxf ! terrain grad
fp%dzdyf => grid%dzdyf ! terrain grad
fp%bbb => grid%bbb ! a rate of spread formula coeff
fp%betafl => grid%betafl ! a rate of spread formula variable
fp%phiwc => grid%phiwc ! a rate of spread formula coeff
fp%r_0 => grid%r_0 ! a rate of spread formula variable
fp%fgip => grid%fgip ! a rate of spread formula coeff
fp%ischap => grid%ischap ! a rate of spread formula switch
! get ignition data
call fire_ignition_convert (config_flags,fire_max_ignitions,fire_ignition_longlat, &
ignition_line,fire_num_ignitions,unit_fxlong,unit_fxlat)
! copy configuration flags to fire internal structures
call set_flags(config_flags)
! refinement r
ir=grid%sr_x ! refinement ratio
jr=grid%sr_y
itimestep=grid%itimestep
restart=config_flags%restart
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,'(a,i1,a,i1,a,i4)') &
'fire_driver_em: ifun from ',fire_ifun_start,' to ',fire_ifun_end,' test steps',tsteps
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
call message(msg)
do istep=0,tsteps ! istep >0 is for testing only, exit after the first call
itimestep = grid%itimestep + istep ! in the first call, do fire_test_steps steps of the fire model
need_lfn_update=.false.
do fire_ifun=fire_ifun_start,fire_ifun_end
! 1 = initialize run pass 1: interpolate height to zsf=terrain
! 2 = initialize run pass 2: set fuel data, terrain gradient
! 3 = initialize timestep: interpolate winds, check for ignition
! 4 = do one timestep
! 5 = copy timestep output to input
! 6 = compute output fluxes
#ifdef DM_PARALLEL
if(need_lfn_update)then
! halo exchange on lfn width 2
#include "HALO_FIRE_LFN.inc"
endif
if(fire_ifun.eq.1)then
! halo exchange on topography
#include "HALO_FIRE_LONGLAT.inc"
!! if(fire_topo_from_atm.eq.1)then
!!#include "HALO_FIRE_HT.inc"
!! endif
! base geopotential and roughness
#include "HALO_FIRE_PHB.inc"
#include "HALO_FIRE_Z0.inc"
elseif(fire_ifun.eq.2)then
! halo exchange on zsf width 2
#include "HALO_FIRE_ZSF.inc"
elseif(fire_ifun.eq.3)then
! halo exchange on atm winds and geopotential, width 1 for interpolation
#include "HALO_FIRE_WIND_A.inc"
#include "HALO_FIRE_PH.inc"
elseif(fire_ifun.eq.4)then
! halo exchange on fire winds width 2 for a 2-step RK method
#include "HALO_FIRE_WIND_F.inc"
elseif(fire_ifun.eq.6)then
! computing fuel_left needs ignition time from neighbors
#include "HALO_FIRE_TIGN.inc"
endif
#endif
! need domain by 1 smaller, in last row.col winds are not set properly
call fire_driver_phys ( &
fire_ifun,need_lfn_update, &
ids,ide-1, kds,kde, jds,jde-1, &
ims,ime, kms,kme, jms,jme, &
ips,min(ipe,ide-1), kps,kpe, jps,min(jpe,jde-1), &
ifds,ifde-ir, jfds,jfde-jr, &
ifms,ifme, jfms,jfme, &
ifps,min(ifpe,ifde-ir), jfps,min(jfpe,jfde-jr), &
ir,jr, & ! atm/fire grid ratio
grid%num_tiles, & ! atm grid tiling
grid%i_start,min(grid%i_end,ide-1), &
grid%j_start,min(grid%j_end,jde-1), &
itimestep,restart,config_flags%fire_fuel_read,config_flags%fire_fuel_cat, & ! in scalars
grid%dt,grid%dx,grid%dy, &
grid%u_frame,grid%v_frame, &
unit_fxlong,unit_fxlat, & ! coordinates of grid center
config_flags%fire_ext_grnd,config_flags%fire_ext_crwn,config_flags%fire_crwn_hgt, &
config_flags%fire_wind_height, & ! height of wind input to fire spread formula
fire_num_ignitions, &
fire_ignition_longlat, &
ignition_line, &
grid%u_2,grid%v_2, & ! atm arrays in
grid%ph_2,grid%phb, & ! geopotential
grid%z0, & ! roughness height
grid%ht, & ! terrain height
grid%xlong,grid%xlat, & ! coordinates of atm grid centers, for ignition location
grid%lfn,grid%tign_g,grid%fuel_frac, & ! state arrays, fire grid
grid%fire_area, & ! redundant, for display, fire grid
lfn_out, & ! work - one timestep
grid%avg_fuel_frac, & ! out redundant arrays, atm grid
grid%grnhfx,grid%grnqfx,grid%canhfx,grid%canqfx, & ! out redundant arrays, atm grid
grid%uah,grid%vah, &
grid%fgrnhfx,grid%fgrnqfx,grid%fcanhfx,grid%fcanqfx, & ! out redundant arrays, atm grid
grid%ros, & ! rate of spread
grid%fxlong,grid%fxlat, &
grid%nfuel_cat, & ! input, or internal for safekeeping
grid%fuel_time, &
fp &
)
#ifdef DM_PARALLEL
if(fire_ifun.eq.2)then
! halo exchange on all fuel data width 2
#include "HALO_FIRE_FUEL.inc"
endif
#endif
enddo
enddo
if(tsteps>0)call crash('fire_driver_em: test run of uncoupled fire model completed')
end subroutine fire_driver_em
!
!*******************
!
! module_fr_fire_driver%%fire_driver
subroutine fire_driver_phys (ifun,need_lfn_update, & 1,84
ids,ide, kds,kde, jds,jde, & ! atm grid dimensions
ims,ime, kms,kme, jms,jme, &
ips,ipe, kps,kpe, jps,jpe, &
ifds, ifde, jfds, jfde, & ! fire grid dimensions
ifms, ifme, jfms, jfme, &
ifps, ifpe, jfps, jfpe, & ! fire patch in - will use smaller
ir,jr, & ! atm/fire grid ratio
num_tiles,i_start,i_end,j_start,j_end, & ! atm grid tiling
itimestep,restart,ifuelread,nfuel_cat0,dt,dx,dy, & ! in scalars
u_frame,v_frame, &
unit_fxlong,unit_fxlat, & ! fxlong, fxlat units in m
fire_ext_grnd,fire_ext_crwn,fire_crwn_hgt, &
fire_wind_height, & ! for vertical wind interpolation
num_ignitions, &
ignition_longlat, &
ignition_line, &
u,v, & ! in arrays, atm grid
ph,phb, &
z0,zs, &
xlong,xlat, &
lfn,tign,fuel_frac, & ! state arrays, fire grid
fire_area, & ! redundant state, fire grid
lfn_out, & ! out level set function
avg_fuel_frac, &
grnhfx,grnqfx,canhfx,canqfx, & ! out redundant arrays, atm grid
uah,vah, & ! out atm grid
fgrnhfx,fgrnqfx,fcanhfx,fcanqfx, & ! out redundant arrays, fire grid
ros, &
fxlong,fxlat, & !
nfuel_cat, & ! in array, data, fire grid, or constant internal
fuel_time, & ! save constant internal data, fire grid
fp & ! fire params
)
USE module_dm, only:wrf_dm_maxval
implicit none
!*** arguments
integer, intent(in)::ifun, &
ids,ide, kds,kde, jds,jde, & ! atm domain bounds
ims,ime, kms,kme, jms,jme, & ! atm memory bounds
ips,ipe, kps,kpe, jps,jpe, & ! atm patch bounds
ifds, ifde, jfds, jfde, & ! fire domain bounds
ifms, ifme, jfms, jfme, & ! fire memory bounds
ifps, ifpe, jfps, jfpe, & ! fire patch bounds
ir,jr, & ! atm/fire grid refinement ratio
itimestep, & ! number of this timestep
ifuelread, & ! how to initialize nfuel_cat:
! -1=not at all, done outside
! 0=from nfuel_cat0
! 1=from altitude
! 2=from file
nfuel_cat0, & ! fuel category to initialize everything to
num_tiles ! number of tiles
logical, intent(in)::restart
logical, intent(out)::need_lfn_update
integer,dimension(num_tiles),intent(in) :: i_start,i_end,j_start,j_end ! atm grid tiling
real, intent(in):: &
dt, & ! time step
dx,dy, & ! atm grid step
u_frame,v_frame, & ! velocity offset
unit_fxlong,unit_fxlat, & ! fxlong, fxlat units in m
fire_crwn_hgt, & ! lowest height crown fire heat is released (m)
fire_ext_grnd, & ! extinction depth of surface fire heat flux (m)
fire_ext_crwn, & ! wind height for vertical interploation to fire spread
fire_wind_height
integer, intent(in):: num_ignitions ! number of ignitions, can be 0
integer, intent(in):: ignition_longlat ! if 1 ignition give as long/lat, otherwise as m from lower left corner
TYPE (ignition_line_type), DIMENSION(num_ignitions), intent(out):: ignition_line
real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::u,v, & ! wind velocity (m/s) (staggered atm grid)
ph, phb ! geopotential (w-points atm grid)
real,intent(in),dimension(ims:ime, jms:jme):: z0, & ! roughness height
zs ! terrain height
real,intent(out),dimension(ims:ime,jms:jme)::&
uah, & ! atm wind at fire_wind_height, diagnostics
vah ! atm wind at fire_wind_height, diagnostics
real, dimension(ims:ime, jms:jme), intent(inout)::xlong, xlat ! inout because of extension at bdry
real, intent(inout), dimension(ifms:ifme,jfms:jfme):: &
nfuel_cat ! fuel data; can be also set inside (cell based, fire grid)
real, intent(inout), dimension(ifms:ifme, jfms:jfme):: &
lfn,tign,fuel_frac, & ! state: level function, ign time, fuel left
lfn_out ! fire wind velocities
real, intent(out), dimension(ifms:ifme, jfms:jfme):: &
fire_area ! fraction of each cell burning
real, intent(out), dimension(ims:ime, jms:jme):: & ! redundant arrays, for display purposes only (atm grid)
avg_fuel_frac, & ! average fuel fraction
grnhfx, & ! heat flux from surface fire (W/m^2)
grnqfx, & ! moisture flux from surface fire (W/m^2)
canhfx, & ! heat flux from crown fire (W/m^2)
canqfx ! moisture flux from crown fire (W/m^2)
real, intent(out), dimension(ifms:ifme, jfms:jfme):: & ! redundant arrays, for display only, fire grid
fgrnhfx, & ! heat flux from surface fire (W/m^2)
fgrnqfx, & ! moisture flux from surface fire (W/m^2)
fcanhfx, & ! heat flux from crown fire (W/m^2)
fcanqfx, & ! moisture flux from crown fire (W/m^2)
ros ! fire rate of spread (m/s)
! ***** data (constant in time) *****
real, dimension(ifms:ifme, jfms:jfme), intent(inout)::fxlong,fxlat ! fire mesh coordinates
real, intent(out), dimension(ifms:ifme, jfms:jfme)::fuel_time ! fire params arrays
type(fire_params),intent(inout)::fp
!*** local
real :: dxf,dyf,time_start,latm, s
integer :: its,ite,jts,jte,kts,kte, & ! tile
ij,i,j,k,id,pid,ipe1,jpe1,ite1,jte1, &
ifts,ifte,jfts,jfte ! fire tile
character(len=128)::msg
character(len=3)::kk
integer::ignitions_done_tile(num_tiles),ignited_tile(num_ignitions,num_tiles)
integer::ignitions_done,ignited_patch(num_ignitions),idex,jdex
!*** executable
! time - assume dt does not change
time_start = itimestep * dt
! fire mesh step
dxf=dx/ir
dyf=dy/jr
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,'(a,2f15.6)')'atmosphere mesh step:',dx,dy
call message(msg)
write(msg,'(a,2f15.6)')'fire mesh step: ',dxf,dyf
call message(msg)
write(msg,7001)'atm domain ','ids',ids,ide,jds,jde
call message(msg)
write(msg,7001)'atm memory ','ims',ims,ime,jms,jme
call message(msg)
write(msg,7001)'atm patch ','ips',ips,ipe,jps,jpe
call message(msg)
write(msg,7001)'fire domain ','ifds',ifds,ifde,jfds,jfde
call message(msg)
write(msg,7001)'fire memory ','ifms',ifms,ifme,jfms,jfme
call message(msg)
write(msg,7001)'fire patch ','ifps',ifps,ifpe,jfps,jfpe
call message(msg)
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
! check mesh dimensions
call check_fmesh(ids,ide,ifds,ifde,ir,'id') ! check if atm and fire grids line up
call check_fmesh(jds,jde,jfds,jfde,jr,'jd')
call check_fmesh(ips,ipe,ifps,ifpe,ir,'ip')
call check_fmesh(jps,jpe,jfps,jfpe,jr,'jp')
call check_mesh_2dim(ips,ipe,jps,jpe,ims,ime,jms,jme) ! check if atm patch fits in atm array
call check_mesh_2dim(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme) ! check if fire patch fits in fire array
call check_mesh_2dim(ips,ipe,jps,jpe,ids,ide,jds,jde) ! check if atm patch fits in atm domain
call check_mesh_2dim(ifps,ifpe,jfps,jfpe,ifds,ifde,jfds,jfde) ! check if fire patch fits in fire domain
!$OMP SINGLE
! init rest of fuel tables with derived constants
if(ifun.eq.1) then
call init_fuel_cats ! common for all threads
endif
!$OMP END SINGLE
pid=0
if(fire_print_file.gt.0)then
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
endif
if(ifun.eq.3)then
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')
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')
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')
endif
! fake atm tile bounds
kts=kps
kte=kpe
! staggered atm patch bounds
ipe1=ifval(ipe.eq.ide,ipe+1,ipe)
jpe1=ifval(jpe.eq.jde,jpe+1,jpe)
! set up fire tiles & interpolate to fire grid
!$OMP PARALLEL DO PRIVATE(ij,its,ite,jts,jte,ite1,jte1,ifts,ifte,jfts,jfte,msg,id) &
!$OMP SCHEDULE(STATIC)
do ij=1,num_tiles
id = ifval(pid.ne.0,pid+ij*10000,0) ! for print
ignitions_done_tile(ij)=0
! set up tile bounds
its = i_start(ij) ! start atmospheric tile in i
ite = i_end(ij) ! end atmospheric tile in i
jts = j_start(ij) ! start atmospheric tile in j
jte = j_end(ij) ! end atmospheric tile in j
ifts= (its-ids)*ir+ifds ! start fire tile in i
ifte= (ite-ids+1)*ir+ifds-1 ! end fire tile in i
jfts= (jts-jds)*jr+jfds ! start fire tile in j
jfte= (jte-jds+1)*jr+jfds-1 ! end fire tile in j
! staggered atm tile bounds
ite1=ifval(ite.eq.ide,ite+1,ite)
jte1=ifval(jte.eq.jde,jte+1,jte)
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,'(a,i3,1x,a,i7,1x,a,i3)')'tile=',ij,' id=',id,' ifun=',ifun
call message(msg)
write(msg,7001)'atm tile ','its',its,ite,jts,jte
call message(msg)
write(msg,7001)'fire tile ','ifts',ifts,ifte,jfts,jfte
call message(msg)
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
! check the tiles
call check_mesh_2dim(its,ite,jts,jte,ips,ipe,jps,jpe) ! check if atm tile fits in atm patch
call check_mesh_2dim(ifts,ifte,jfts,jfte,ifps,ifpe,jfps,jfpe) ! check if fire tile fits in fire patch
call check_mesh_2dim(ifts-2,ifte+2,jfts-2,jfte+2,ifms,ifme,jfms,jfme)! check if fire node tile fits in memory
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,'(a,i6,a,2(f15.6,a))')'time step',itimestep,' at',time_start,' duration',dt,'s'
call message(msg)
7001 format(a,' dimensions ',a4,':',i6,' to ',i6,' by ',i6,' to ',i6)
write(msg,'(a,2i9)')'refinement ratio:',ir,jr
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
if(ifun.eq.1)then ! set terrain
if(restart)then
call message('restart - topo initialization skipped')
else
!
! call print_2d_stats(ips,ipe,jps,jpe,ims,ime,jms,jme,zs,'driver:zs')
!
! ! interpolate terrain height
! if(fire_topo_from_atm.eq.1)then
! call interpolate_z2fire(id, & ! for debug output, <= 0 no output
! ids,ide, jds,jde, & ! atm grid dimensions
! ims,ime, jms,jme, &
! ips,ipe,jps,jpe, &
! its,ite,jts,jte, &
! ifds, ifde, jfds, jfde, & ! fire grid dimensions
! ifms, ifme, jfms, jfme, &
! ifts,ifte,jfts,jfte, &
! ir,jr, & ! atm/fire grid ratio
! zs, & ! atm grid arrays in
! fp%zsf) ! fire grid arrays out
! else
!!$OMP CRITICAL(FIRE_DRIVER_CRIT)
! write(msg,'(a,i3,a)')'fire_topo_from_atm=',fire_topo_from_atm,' assuming ZSF set, interpolation skipped'
!!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
! endif
if(ignition_longlat .eq.0)then
! set ideal fire mesh coordinates - used for ignition only
! do not forget to set unit_fxlong, unit_fxlat outside of parallel loop
!call set_ideal_coord( dxf,dyf, &
! ifds,ifde,jfds,jfde, &
! ifms,ifme,jfms,jfme, &
! ifts,ifte,jfts,jfte, &
! fxlong,fxlat )
!call set_ideal_coord( dx,dy, &
! ids,ide,jds,jde, &
! ims,ime,jms,jme, &
! its,ite,jts,jte, &
! xlong,xlat )
elseif(use_atm_vars)then
! assume halo xlong xlat
! interpolate nodal coordinates
#ifdef DEBUG_OUT
call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,xlat,'xlat',id)
call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,xlong,'xlong',id)
#endif
call interpolate_z2fire(id, & ! for debug output, <= 0 no output
ids,ide, jds,jde, & ! atm grid dimensions
ims,ime, jms,jme, &
ips,ipe,jps,jpe, &
its,ite,jts,jte, &
ifds, ifde, jfds, jfde, & ! fire grid dimensions
ifms, ifme, jfms, jfme, &
ifts,ifte,jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
xlat, & ! atm grid arrays in
fxlat) ! fire grid arrays out
call interpolate_z2fire(id, & ! for debug output, <= 0 no output
ids,ide, jds,jde, & ! atm grid dimensions
ims,ime, jms,jme, &
ips,ipe,jps,jpe, &
its,ite,jts,jte, &
ifds, ifde, jfds, jfde, & ! fire grid dimensions
ifms, ifme, jfms, jfme, &
ifts,ifte,jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
xlong, & ! atm grid arrays in
fxlong) ! fire grid arrays out
endif
endif
elseif(ifun.eq.2)then
! after the loop where zsf created exited and all synced
call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fp%zsf,'driver_phys:zsf')
elseif(ifun.eq.3)then ! interpolate winds to the fire grid
if(use_atm_vars)then
call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,z0,'z0',id)
call write_array_m3(its,ite1,kts,kde-1,jts,jte,ims,ime,kms,kme,jms,jme,u,'u_2',id)
call write_array_m3(its,ite,kts,kde-1,jts,jte1,ims,ime,kms,kme,jms,jme,v,'v_2',id)
call write_array_m3(its,ite,kts,kde,jts,jte,ims,ime,kms,kme,jms,jme,ph,'ph_2',id)
call write_array_m3(its,ite,kts,kde,jts,jte,ims,ime,kms,kme,jms,jme,phb,'phb',id)
call interpolate_atm2fire(id, & ! flag for debug output
fire_wind_height, & ! height to interpolate to
ids,ide, kds,kde, jds,jde, & ! atm grid dimensions
ims,ime, kms,kme, jms,jme, &
ips,ipe, jps,jpe, &
its,ite,jts,jte, &
ifds, ifde, jfds, jfde, & ! fire grid dimensions
ifms, ifme, jfms, jfme, &
ifps, ifpe, jfps, jfpe, & ! fire patch bounds
ifts, ifte, jfts, jfte, &
ir,jr, & ! atm/fire grid ratio
u_frame, v_frame, & ! velocity frame correction
u,v, & ! 3D atm grid arrays in
ph,phb, &
z0,zs, & ! 2D atm grid arrays in
uah,vah, & ! 2D atm grid out
fp%vx,fp%vy) ! fire grid arrays out
endif
endif
! the following executes in any case
call fire_model (id,ifun,restart,need_lfn_update, &
num_ignitions, & ! switches
ifuelread,nfuel_cat0, & ! initialize fuel categories
ifds,ifde,jfds,jfde, & ! fire domain dims
ifms,ifme,jfms,jfme, & ! fire memory dims
ifps,ifpe,jfps,jfpe, &
ifts,ifte,jfts,jfte, & ! fire patch dims
time_start,dt, & ! time and increment
dxf,dyf, & ! fire mesh spacing
ignition_line, & ! description of ignition lines
ignitions_done_tile(ij),ignited_tile(1,ij), &
fxlong,fxlat,unit_fxlong,unit_fxlat, & ! fire mesh coordinates
lfn,lfn_out,tign,fuel_frac, & ! state: level function, ign time, fuel left
fire_area, & ! output: fraction of cell burning
fgrnhfx,fgrnqfx, & ! output: heat fluxes
ros, & ! output: rate of spread for display
nfuel_cat, & ! fuel data per point
fuel_time, & ! save derived internal data
fp & ! fire coefficients
)
if(ifun.eq.6)then ! heat fluxes into the atmosphere
call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fgrnhfx,'fire_driver:fgrnhfx')
call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fgrnqfx,'fire_driver:fgrnqfx')
! sum the fluxes over atm cells
if(use_atm_vars)then
call sum_2d_cells( &
ifms,ifme,jfms,jfme, &
ifts,ifte,jfts,jfte, &
fuel_frac, &
ims, ime, jms, jme, &
its,ite,jts,jte, &
avg_fuel_frac)
call sum_2d_cells( &
ifms,ifme,jfms,jfme, &
ifts,ifte,jfts,jfte, &
fgrnhfx, &
ims, ime, jms, jme, &
its,ite,jts,jte, &
grnhfx)
!comment out the next call to get results as before commit 55fd92051196b796891b60cb7ec1c4bdb8800078
call sum_2d_cells( &
ifms,ifme,jfms,jfme, &
ifts,ifte,jfts,jfte, &
fgrnqfx, &
ims, ime, jms, jme, &
its,ite,jts,jte, &
grnqfx)
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,'(a,f6.3)')'fire-atmosphere feedback scaling ',fire_atm_feedback
!$OMP end CRITICAL(FIRE_DRIVER_CRIT)
call message(msg)
s = 1./(ir*jr)
do j=jts,jte
do i=its,ite
! scale surface fluxes to get the averages
avg_fuel_frac(i,j)=avg_fuel_frac(i,j)*s
grnhfx(i,j)=fire_atm_feedback*grnhfx(i,j)*s
grnqfx(i,j)=fire_atm_feedback*grnqfx(i,j)*s
! we do not have canopy fluxes yet...
canhfx(i,j)=0
canqfx(i,j)=0
enddo
enddo
call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnhfx,'fire_driver:grnhfx')
call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme,grnqfx,'fire_driver:grnqfx')
endif
endif ! ifun=6
enddo ! tiles
!$OMP END PARALLEL DO
#ifdef DEBUG_OUT
if(ifun.eq.1)then
if(pid.ne.0)then
call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,zs,'zs',pid)
call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fp%zsf,'zsf',pid)
endif
endif
#endif
if (ifun.eq.3)then
ignitions_done=ignitions_done_tile(1) ! all should be the same
do i=1,ignitions_done
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,'(2(a,i4,1x))')'fire_driver_phys: checking ignition ',i,' of ',ignitions_done
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
call message(msg)
ignited_patch(i)=0
do ij=1,num_tiles
ignited_patch(i)=ignited_patch(i)+ignited_tile(i,ij)
enddo
#ifdef DM_PARALLEL
call message('fire_driver_phys: checking ignitions, collect counts')
call wrf_dm_maxval(ignited_patch(i),idex,jdex)
call message('fire_driver_phys: collected')
#endif
if(ignited_patch(i).eq.0)then
call crash('fire_driver_phys: Ignition failed, no nodes ignited. Bad coordinates?')
endif
enddo
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')
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')
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')
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')
#ifdef DEBUG_OUT
if(pid.gt.0)then
call write_array_m(ips,ipe1,jps,jpe,ims,ime,jms,jme,uah,'uah',pid)
call write_array_m(ips,ipe,jps,jpe1,ims,ime,jms,jme,vah,'vah',pid)
call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnhfx,'grnhfx',pid)
call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnqfx,'grnqfx',pid)
call write_array_m3(ips,ipe1,kds,kde+1,jps,jpe,ims,ime,kms,kme,jms,jme,u,'u',pid)
call write_array_m3(ips,ipe,kds,kde+1,jps,jpe1,ims,ime,kms,kme,jms,jme,v,'v',pid)
call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fp%vx,'uf',pid)
call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fp%vy,'vf',pid)
endif
#endif
endif
if(ifun.eq.5)then
#ifdef DEBUG_OUT
if(pid.gt.0)then
call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,lfn,'lfn',pid)
call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,tign,'tign',pid)
endif
#endif
endif
if(ifun.eq.6)then
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')
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')
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')
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')
#ifdef DEBUG_OUT
if(pid.gt.0)then
call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnhfx,'grnhfx',pid)
call write_array_m(ips,ipe,jps,jpe,ims,ime,jms,jme,grnqfx,'grnqfx',pid)
call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fuel_frac,'fuel_frac',pid)
call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fgrnhfx,'fgrnhfx',pid)
call write_array_m(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,fgrnqfx,'fgrnqfx',pid)
endif
#endif
endif
end subroutine fire_driver_phys
!
!*******************
!
subroutine fire_ignition_convert (config_flags,fire_max_ignitions,fire_ignition_longlat, & 1,5
ignition_line,fire_num_ignitions,unit_fxlong,unit_fxlat)
USE module_configure, only : grid_config_rec_type
implicit none
! create ignition arrays from scalar flags
!*** arguments
TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags
integer, intent(in)::fire_max_ignitions
TYPE (ignition_line_type), DIMENSION(fire_max_ignitions), intent(out):: ignition_line ! any values from input discarded
integer, intent(out)::fire_num_ignitions,fire_ignition_longlat
real, intent(out)::unit_fxlong,unit_fxlat
!*** local
integer::i
logical:: real,ideal
real::lat_ctr,lon_ctr
!*** executable
! this is only until I figure out how to input arrays through the namelist...
if(fire_max_ignitions.lt.5)call crash('fire_max_ignitions too small')
! figure out which kind of coordinates from the first given
ideal=config_flags%fire_ignition_start_x1 .ne.0. .or. config_flags%fire_ignition_start_y1 .ne. 0.
real=config_flags%fire_ignition_start_lon1 .ne. 0. .or. config_flags%fire_ignition_start_lat1 .ne. 0.
if(ideal)call message('Using ideal ignition coordinates, m from the lower left domain corner')
if(real)call message('Using real ignition coordinates, longitude and latitude')
if(ideal.and.real)call crash('Only one of the ideal or real coordinates may be given')
fire_ignition_longlat=0 ! default, if no ignition
if(ideal)then
! use values from _x and _y variables
fire_ignition_longlat=0
ignition_line(1)%start_x=config_flags%fire_ignition_start_x1
ignition_line(1)%start_y=config_flags%fire_ignition_start_y1
ignition_line(1)%end_x=config_flags%fire_ignition_end_x1
ignition_line(1)%end_y=config_flags%fire_ignition_end_y1
ignition_line(2)%start_x=config_flags%fire_ignition_start_x2
ignition_line(2)%start_y=config_flags%fire_ignition_start_y2
ignition_line(2)%end_x=config_flags%fire_ignition_end_x2
ignition_line(2)%end_y=config_flags%fire_ignition_end_y2
ignition_line(3)%start_x=config_flags%fire_ignition_start_x3
ignition_line(3)%start_y=config_flags%fire_ignition_start_y3
ignition_line(3)%end_x=config_flags%fire_ignition_end_x3
ignition_line(3)%end_y=config_flags%fire_ignition_end_y3
ignition_line(4)%start_x=config_flags%fire_ignition_start_x4
ignition_line(4)%start_y=config_flags%fire_ignition_start_y4
ignition_line(4)%end_x=config_flags%fire_ignition_end_x4
ignition_line(4)%end_y=config_flags%fire_ignition_end_y4
ignition_line(5)%start_x=config_flags%fire_ignition_start_x5
ignition_line(5)%start_y=config_flags%fire_ignition_start_y5
ignition_line(5)%end_x=config_flags%fire_ignition_end_x5
ignition_line(5)%end_y=config_flags%fire_ignition_end_y5
endif
if(real)then
! use values from _long and _lat
fire_ignition_longlat=1
ignition_line(1)%start_x=config_flags%fire_ignition_start_lon1
ignition_line(1)%start_y=config_flags%fire_ignition_start_lat1
ignition_line(1)%end_x=config_flags%fire_ignition_end_lon1
ignition_line(1)%end_y=config_flags%fire_ignition_end_lat1
ignition_line(2)%start_x=config_flags%fire_ignition_start_lon2
ignition_line(2)%start_y=config_flags%fire_ignition_start_lat2
ignition_line(2)%end_x=config_flags%fire_ignition_end_lon2
ignition_line(2)%end_y=config_flags%fire_ignition_end_lat2
ignition_line(3)%start_x=config_flags%fire_ignition_start_lon3
ignition_line(3)%start_y=config_flags%fire_ignition_start_lat3
ignition_line(3)%end_x=config_flags%fire_ignition_end_lon3
ignition_line(3)%end_y=config_flags%fire_ignition_end_lat3
ignition_line(4)%start_x=config_flags%fire_ignition_start_lon4
ignition_line(4)%start_y=config_flags%fire_ignition_start_lat4
ignition_line(4)%end_x=config_flags%fire_ignition_end_lon4
ignition_line(4)%end_y=config_flags%fire_ignition_end_lat4
ignition_line(5)%start_x=config_flags%fire_ignition_start_lon5
ignition_line(5)%start_y=config_flags%fire_ignition_start_lat5
ignition_line(5)%end_x=config_flags%fire_ignition_end_lon5
ignition_line(5)%end_y=config_flags%fire_ignition_end_lat5
endif
! common to both cases
ignition_line(1)%ros=config_flags%fire_ignition_ros1
ignition_line(1)%radius=config_flags%fire_ignition_radius1
ignition_line(1)%start_time=config_flags%fire_ignition_start_time1
ignition_line(1)%end_time=config_flags%fire_ignition_end_time1
ignition_line(2)%ros=config_flags%fire_ignition_ros2
ignition_line(2)%radius=config_flags%fire_ignition_radius2
ignition_line(2)%start_time=config_flags%fire_ignition_start_time2
ignition_line(2)%end_time=config_flags%fire_ignition_end_time2
ignition_line(3)%ros=config_flags%fire_ignition_ros3
ignition_line(3)%radius=config_flags%fire_ignition_radius3
ignition_line(3)%start_time=config_flags%fire_ignition_start_time3
ignition_line(3)%end_time=config_flags%fire_ignition_end_time3
ignition_line(4)%ros=config_flags%fire_ignition_ros4
ignition_line(4)%radius=config_flags%fire_ignition_radius4
ignition_line(4)%start_time=config_flags%fire_ignition_start_time4
ignition_line(4)%end_time=config_flags%fire_ignition_end_time4
ignition_line(5)%ros=config_flags%fire_ignition_ros5
ignition_line(5)%radius=config_flags%fire_ignition_radius5
ignition_line(5)%start_time=config_flags%fire_ignition_start_time5
ignition_line(5)%end_time=config_flags%fire_ignition_end_time5
!
fire_num_ignitions=0
do i=1,min(5,config_flags%fire_num_ignitions)
! count the ignitions
if(ignition_line(i)%radius.gt.0.)fire_num_ignitions=i
! expand ignition data given as zero
if(ignition_line(i)%end_x.eq.0.)ignition_line(i)%end_x=ignition_line(i)%start_x
if(ignition_line(i)%end_y.eq.0.)ignition_line(i)%end_y=ignition_line(i)%start_y
if(ignition_line(i)%end_time.eq.0.)ignition_line(i)%end_time=ignition_line(i)%start_time
enddo
if(fire_ignition_longlat .eq. 0)then
! ideal
! ignition is in m
unit_fxlong=1.
unit_fxlat=1.
! will set fire mesh coordinates to uniform mesh below
else
! real
lat_ctr=config_flags%cen_lat
lon_ctr=config_flags%cen_lon
! 1 degree in m (approximate OK)
unit_fxlat=pi2/(360.*reradius) ! earth circumference in m / 360 degrees
unit_fxlong=cos(lat_ctr*pi2/360.)*unit_fxlat ! latitude
endif
end subroutine fire_ignition_convert
!
!*****************************
!
subroutine interpolate_z2fire(id, & ! for debug output, <= 0 no output 2,6
ids,ide, jds,jde, & ! atm grid dimensions
ims,ime, jms,jme, &
ips,ipe,jps,jpe, &
its,ite,jts,jte, &
ifds, ifde, jfds, jfde, & ! fire grid dimensions
ifms, ifme, jfms, jfme, &
ifts,ifte,jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
zs, & ! atm grid arrays in
zsf) ! fire grid arrays out
implicit none
!*** purpose: interpolate height
!*** arguments
integer, intent(in)::id, &
ids,ide, jds,jde, & ! atm domain bounds
ims,ime,jms,jme, & ! atm memory bounds
ips,ipe,jps,jpe, &
its,ite,jts,jte, & ! atm tile bounds
ifds, ifde, jfds, jfde, & ! fire domain bounds
ifms, ifme, jfms, jfme, & ! fire memory bounds
ifts,ifte,jfts,jfte, & ! fire tile bounds
ir,jr ! atm/fire grid refinement ratio
real, intent(in), dimension(ims:ime, jms:jme):: zs ! terrain height at atm cell centers & ! terrain height
real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
zsf ! terrain height fire grid nodes
!*** local
real, dimension(its-2:ite+2,jts-2:jte+2):: za ! terrain height
integer:: i,j,jts1,jte1,its1,ite1,jfts1,jfte1,ifts1,ifte1,itso,jtso,iteo,jteo
! terrain height
jts1=max(jts-1,jds) ! lower loop limit by one less when at end of domain
its1=max(its-1,ids) ! ASSUMES THE HALO IS THERE if patch != domain
jte1=min(jte+1,jde)
ite1=min(ite+1,ide)
do j = jts1,jte1
do i = its1,ite1
! copy to local array
za(i,j)=zs(i,j)
enddo
enddo
call continue_at_boundary(1,1,0., & ! do x direction or y direction
its-2,ite+2,jts-2,jte+2, & ! memory dims
ids,ide,jds,jde, & ! domain dims - winds defined up to +1
ips,ipe,jps,jpe, & ! patch dims - winds defined up to +1
its1,ite1,jts1,jte1, & ! tile dims
itso,jtso,iteo,jteo, &
za) ! array
! interpolate to tile plus strip along domain boundary if at boundary
jfts1=snode(jfts,jfds,-1) ! lower loop limit by one less when at end of domain
ifts1=snode(ifts,ifds,-1)
jfte1=snode(jfte,jfde,+1)
ifte1=snode(ifte,ifde,+1)
call interpolate_2d( &
its-2,ite+2,jts-2,jte+2, & ! memory dims atm grid tile
its1-1,ite1+1,jts1-1,jte1+1, & ! where atm grid values set
ifms,ifme,jfms,jfme, & ! array dims fire grid
ifts1,ifte1,jfts1,jfte1, & ! dimensions fire grid tile
ir,jr, & ! refinement ratio
real(ids),real(jds),ifds+(ir-1)*0.5,jfds+(jr-1)*0.5, & ! line up by lower left corner of domain
za, & ! in atm grid
zsf) ! out fire grid
end subroutine interpolate_z2fire
!
!*****************************
!
subroutine interpolate_atm2fire(id, & ! for debug output, <= 0 no output 1,43
fire_wind_height, & ! interpolation height
ids,ide, kds,kde, jds,jde, & ! atm grid dimensions
ims,ime, kms,kme, jms,jme, &
ips,ipe,jps,jpe, &
its,ite,jts,jte, &
ifds, ifde, jfds, jfde, & ! fire grid dimensions
ifms, ifme, jfms, jfme, &
ifps, ifpe, jfps, jfpe, & ! fire patch bounds
ifts,ifte,jfts,jfte, &
ir,jr, & ! atm/fire grid ratio
u_frame, v_frame, & ! velocity frame correction
u,v, & ! atm grid arrays in
ph,phb, &
z0,zs, &
uah,vah, &
uf,vf) ! fire grid arrays out
implicit none
!*** purpose: interpolate winds and height
!*** arguments
integer, intent(in)::id
real, intent(in):: fire_wind_height ! height above the terrain for vertical interpolation
integer, intent(in):: &
ids,ide, kds,kde, jds,jde, & ! atm domain bounds
ims,ime, kms,kme, jms,jme, & ! atm memory bounds
ips,ipe,jps,jpe, &
its,ite,jts,jte, & ! atm tile bounds
ifds, ifde, jfds, jfde, & ! fire domain bounds
ifms, ifme, jfms, jfme, & ! fire memory bounds
ifps, ifpe, jfps, jfpe, & ! fire patch bounds
ifts,ifte,jfts,jfte, & ! fire tile bounds
ir,jr ! atm/fire grid refinement ratio
real, intent(in):: u_frame, v_frame ! velocity frame correction
real,intent(in),dimension(ims:ime,kms:kme,jms:jme)::&
u,v, & ! atm wind velocity, staggered
ph,phb ! geopotential
real,intent(in),dimension(ims:ime,jms:jme)::&
z0, & ! roughness height
zs ! terrain height
real,intent(out),dimension(ims:ime,jms:jme)::&
uah, & ! atm wind at fire_wind_height, diagnostics
vah !
real,intent(out), dimension(ifms:ifme,jfms:jfme)::&
uf,vf ! wind velocity fire grid nodes
!*** local
character(len=256)::msg
#define TDIMS its-2,ite+2,jts-2,jte+2
real, dimension(its-2:ite+2,jts-2:jte+2):: ua,va ! atm winds, interpolated over height, still staggered grid
real, dimension(its-2:ite+2,kds:kde,jts-2:jte+2):: altw,altub,altvb,hgtu,hgtv ! altitudes
integer:: i,j,k,ifts1,ifte1,jfts1,jfte1,ite1,jte1
integer::itst,itet,jtst,jtet,itsu,iteu,jtsu,jteu,itsv,itev,jtsv,jtev
integer::kdmax,its1,jts1,ips1,jps1
integer::itsou,iteou,jtsou,jteou,itsov,iteov,jtsov,jteov
real:: ground,loght,loglast,logz0,logfwh,ht,zr
real::r_nan
integer::i_nan
equivalence (i_nan,r_nan)
!*** executable
! debug init local arrays
i_nan=2147483647
ua=r_nan
va=r_nan
altw=r_nan
altub=r_nan
hgtu=r_nan
hgtv=r_nan
if(kds.ne.1)then
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,*)'WARNING: bottom index kds=',kds,' should be 1?'
call message(msg)
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
endif
! ^ j
! ------------ |
! | | ----> i
! u p |
! | | nodes in cell (i,j)
! ------v----- view from top
!
! v(ide,jde+1)
! -------x------
! | |
! | |
! u(ide,jde) x x x u(ide+1,jde)
! | p(ide,hde) |
! | | p=ph,phb,z0,...
! -------x------
! v(ide,jde)
!
! staggered values set u(ids:ide+1,jds:jde) v(ids:ide,jds:jde+1)
! p=ph+phb set at (ids:ide,jds:jde)
! location of u(i,j) needs p(i-1,j) and p(i,j)
! location of v(i,j) needs p(i,j-1) and p(i,j)
! *** NOTE need HALO in ph, phb
! so we can compute only u(ids+1:ide,jds:jde) v(ids:ide,jds+1,jde)
! unless we extend p at the boundary
! but because we care about the fire way in the inside it does not matter
! if the fire gets close to domain boundary the simulation is over anyway
ite1=snode(ite,ide,1)
jte1=snode(jte,jde,1)
! do this in any case to check for nans
call print_3d_stats(its,ite1,kds,kde,jts,jte,ims,ime,kms,kme,jms,jme,u,'wind U in')
call print_3d_stats(its,ite,kds,kde,jts,jte1,ims,ime,kms,kme,jms,jme,v,'wind V in')
if(fire_print_msg.gt.0)then
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,'(a,f7.2,a)')'interpolate_atm2fire: log-interpolation of wind to',fire_wind_height,' m'
call message(msg)
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
endif
! indexing
! file for u
itst=ifval(ids.eq.its,its,its-1)
itet=ifval(ide.eq.ite,ite,ite+1)
jtst=ifval(jds.ge.jts,jts,jts-1)
jtet=ifval(jde.eq.jte,jte,jte+1)
if(fire_print_msg.ge.1)then
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,7001)'atm input ','tile',its,ite,jts,jte
call message(msg)
write(msg,7001)'altw ','tile',itst,itet,jtst,jtet
call message(msg)
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
endif
7001 format(a,' dimensions ',a4,':',i6,' to ',i6,' by ',i6,' to ',i6)
kdmax=kde-1 ! max layer to interpolate from, can be less
do j = jtst,jtet
do k=kds,kdmax+1
do i = itst,itet
altw(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g ! altitude of the bottom w-point
enddo
enddo
enddo
! values at u points
itsu=ifval(ids.eq.its,its+1,its) ! staggered direction
iteu=ifval(ide.eq.ite,ite,ite+1) ! staggered direction
jtsu=ifval(jds.ge.jts,jts,jts-1)
jteu=ifval(jde.eq.jte,jte,jte+1)
if(fire_print_msg.ge.1)then
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,7001)'u interp at','tile',itsu,iteu,jtsu,jteu
call message(msg)
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
endif
do j = jtsu,jteu
do k=kds,kdmax+1
do i = itsu,iteu
altub(i,k,j)= 0.5*(altw(i-1,k,j)+altw(i,k,j)) ! altitude of the bottom point under u-point
enddo
enddo
do k=kds,kdmax
do i = itsu,iteu
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
enddo
enddo
enddo
! values at v points
jtsv=ifval(jds.eq.jts,jts+1,jts) ! staggered direction
jtev=ifval(jde.eq.jte,jte,jte+1) ! staggered direction
itsv=ifval(ids.ge.its,its,its-1)
itev=ifval(ide.eq.ite,ite,ite+1)
if(fire_print_msg.ge.1)then
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,7001)'v interp at','tile',itsv,itev,jtsv,jtev
call message(msg)
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
endif
do j = jtsv,jtev
do k=kds,kdmax+1
do i = itsv,itev
altvb(i,k,j)= 0.5*(altw(i,k,j-1)+altw(i,k,j)) ! altitude of the bottom point under v-point
enddo
enddo
do k=kds,kdmax
do i = itsv,itev
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
enddo
enddo
enddo
#ifdef DEBUG_OUT
call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,altub,'altub',id)
call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,altvb,'altvb',id)
call write_array_m3(itsu,iteu,kds,kdmax,jtsu,jteu,its-2,ite+2,kds,kde,jts-2,jte+2,hgtu,'hgtu',id)
call write_array_m3(itsv,itev,kds,kdmax,jtsv,jtev,its-2,ite+2,kds,kde,jts-2,jte+2,hgtv,'hgtv',id)
#endif
logfwh = log(fire_wind_height)
! interpolate u, staggered in X
do j = jtsu,jteu ! compute on domain by 1 smaller, otherwise z0 and ph not available
do i = itsu,iteu ! compute with halo 2
zr = 0.5*(z0(i,j)+z0(i-1,j)) ! interpolated roughness length under this u point
if(fire_wind_height > zr)then !
do k=kds,kdmax
ht = hgtu(i,k,j) ! height of this u point above the ground
if( .not. ht < fire_wind_height) then ! found layer k this point is in
loght = log(ht)
if(k.eq.kds)then ! first layer, log linear interpolation from 0 at zr
logz0 = log(zr)
ua(i,j)= u(i,k,j)*(logfwh-logz0)/(loght-logz0)
else ! log linear interpolation
loglast=log(hgtu(i,k-1,j))
ua(i,j)= u(i,k-1,j) + (u(i,k,j) - u(i,k-1,j)) * ( logfwh - loglast) / (loght - loglast)
endif
goto 10
endif
if(k.eq.kdmax)then ! last layer, still not high enough
ua(i,j)=u(i,k,j)
endif
enddo
10 continue
else ! roughness higher than the fire wind height
ua(i,j)=0.
endif
enddo
enddo
! interpolate v, staggered in Y
do j = jtsv,jtev
do i = itsv,itev
zr = 0.5*(z0(i,j-1)+z0(i,j)) ! roughness length under this v point
if(fire_wind_height > zr)then !
do k=kds,kdmax
ht = hgtv(i,k,j) ! height of this u point above the ground
if( .not. ht < fire_wind_height) then ! found layer k this point is in
loght = log(ht)
if(k.eq.kds)then ! first layer, log linear interpolation from 0 at zr
logz0 = log(zr)
va(i,j)= v(i,k,j)*(logfwh-logz0)/(loght-logz0)
else ! log linear interpolation
loglast=log(hgtv(i,k-1,j))
va(i,j)= v(i,k-1,j) + (v(i,k,j) - v(i,k-1,j)) * ( logfwh - loglast) / (loght - loglast)
endif
goto 11
endif
if(k.eq.kdmax)then ! last layer, still not high enough
va(i,j)=v(i,k,j)
endif
enddo
11 continue
else ! roughness higher than the fire wind height
va(i,j)=0.
endif
enddo
enddo
#ifdef DEBUG_OUT
! store the output for diagnostics
do j = jts,jte1
do i = its,ite1
uah(i,j)=ua(i,j)
vah(i,j)=va(i,j)
enddo
enddo
call write_array_m(its,ite1,jts,jte,ims,ime,jms,jme,uah,'uah_n',id) ! no reflection
call write_array_m(its,ite,jts,jte1,ims,ime,jms,jme,vah,'vah_n',id)
#endif
ips1 = ifval(ips.eq.ids,ips+1,ips)
call continue_at_boundary(1,1,0., & ! x direction
TDIMS, &! memory dims atm grid tile
ids+1,ide,jds,jde, & ! domain dims - where u defined
ips1,ipe,jps,jpe, & ! patch dims
itsu,iteu,jtsu,jteu, & ! tile dims - in nonextended direction one beyond if at patch boundary but not domain
itsou,iteou,jtsou,jteou, & ! out, where set
ua) ! array
jps1 = ifval(jps.eq.jds,jps+1,jps)
call continue_at_boundary(1,1,0., & ! y direction
TDIMS, & ! memory dims atm grid tile
ids,ide,jds+1,jde, & ! domain dims - where v wind defined
ips,ipe,jps1,jpe, & ! patch dims
itsv,itev,jtsv,jtev, & ! tile dims
itsov,iteov,jtsov,jteov, & ! where set
va) ! array
! store the output for diagnostics
do j = jts,jte1
do i = its,ite1
uah(i,j)=ua(i,j)
vah(i,j)=va(i,j)
enddo
enddo
#ifdef DEBUG_OUT
call write_array_m(itsou,iteou,jtsou,jteou,TDIMS,ua,'ua',id)
call write_array_m(itsov,iteov,jtsov,jteov,TDIMS,va,'va',id)
#endif
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
! don't have all values valid, don't check
write(msg,12)'atm mesh wind U at',fire_wind_height,' m'
call print_2d_stats(itsou,iteou,jtsou,jteou,TDIMS,ua,msg)
write(msg,12)'atm mesh wind V at',fire_wind_height,' m'
call print_2d_stats(itsov,iteov,jtsov,jteov,TDIMS,va,msg)
12 format(a,f6.2,a)
call print_2d_stats(its,ite1,jts,jte,ims,ime,jms,jme,uah,'UAH')
call print_2d_stats(its,ite,jts,jte1,ims,ime,jms,jme,vah,'VAH')
!call write_array_m(its,ite1,jts,jte,ims,ime,jms,jme,uah,'uah',id)
!call write_array_m(its,ite,jts,jte1,ims,ime,jms,jme,vah,'vah',id)
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
! ---------------
! | F | F | F | F | Example of atmospheric and fire grid with
! |-------|-------| ir=jr=4.
! | F | F | F | F | Winds are given at the midpoints of the sides of the atmosphere grid,
! ua------z-------| interpolated to midpoints of the cells of the fine fire grid F.
! | F | F | F | F | This is (1,1) cell of atmosphere grid, and [*] is the (1,1) cell of the fire grid.
! |---------------| ua(1,1) <--> uf(0.5,2.5)
! | * | F | F | F | va(1,1) <--> vf(2.5,0.5)
! -------va------ za(1,1) <--> zf(2.5,2.5)
!
! ^ x2
! | --------va(1,2)---------
! | | | | Example of atmospheric and fire grid with
! | | | | ir=jr=1.
! | | za,zf | Winds are given at the midpoints of the sides of the atmosphere grid,
! | ua(1,1)----uf,vf-----ua(2,1) interpolated to midpoints of the cells of the (the same) fire grid
! | | (1,1) | ua(1,1) <--> uf(0.5,1)
! | | | | va(1,1) <--> vf(1,0.5)
! | | | | za(1,1) <--> zf(1,1)
! | --------va(1,1)---------
! |--------------------> x1
!
! Meshes are aligned by the lower left cell of the domain. Then in the above figure
! u = node with the ua component of the wind at (ids,jds), midpoint of side
! v = node with the va component of the wind at (ids,jds), midpoint of side
! * = fire grid node at (ifds,jfds)
! z = node with height, midpoint of cell
!
! ua(ids,jds)=uf(ifds-0.5,jfds+jr*0.5-0.5) = uf(ifds-0.5,jfds+(jr-1)*0.5)
! va(ids,jds)=vf(ifds+ir*0.5-0.5,jfds-0.5) = vf(ifds+(ir-1)*0.5,jfds-0.5)
! 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)
! ifts1=max(snode(ifts,ifps,-1),ifds) ! go 1 beyond patch boundary but not at domain boundary
! ifte1=min(snode(ifte,ifpe,+1),ifde)
! jfts1=max(snode(jfts,jfps,-1),jfds)
! jfte1=min(snode(jfte,jfpe,+1),jfde)
call interpolate_2d( &
TDIMS, & ! memory dims atm grid tile
itsou,iteou,jtsou,jteou,& ! where set
ifms,ifme,jfms,jfme, & ! array dims fire grid
ifts,ifte,jfts,jfte,& ! dimensions on the fire grid to interpolate to
ir,jr, & ! refinement ratio
real(ids),real(jds),ifds-0.5,jfds+(jr-1)*0.5, & ! line up by lower left corner of domain
ua, & ! in atm grid
uf) ! out fire grid
call interpolate_2d( &
TDIMS, & ! memory dims atm grid tile
itsov,iteov,jtsov,jteov,& ! where set
ifms,ifme,jfms,jfme, & ! array dims fire grid
ifts,ifte,jfts,jfte,& ! dimensions on the fire grid to interpolate to
ir,jr, & ! refinement ratio
real(ids),real(jds),ifds+(ir-1)*0.5,jfds-0.5, & ! line up by lower left corner of domain
va, & ! in atm grid
vf) ! out fire grid
call print_2d_stats_vec(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,vf,'fire wind (m/s)')
! call print_2d_stats_vec(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,vf,'fire wind extended')
#ifdef DEBUG_OUT
call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,uf,'uf1',id)
call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,vf,'vf1',id)
! call write_array_m(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,uf,'uf1',id)
! call write_array_m(ifts1,ifte1,jfts1,jfte1,ifms,ifme,jfms,jfme,vf,'vf1',id)
#endif
return
end subroutine interpolate_atm2fire
!
!***
!
subroutine check_fmesh(ids,ide,ifds,ifde,ir,s) 4,1
!*** purpose: check if fire and atm meshes line up
implicit none
!*** arguments
integer, intent(in)::ids,ide,ifds,ifde,ir
character(len=*),intent(in)::s
!*** local
character(len=128)msg
!*** executable
if ((ide-ids+1)*ir.ne.(ifde-ifds+1))then
!$OMP CRITICAL(FIRE_DRIVER_CRIT)
write(msg,1)s,ids,ide,ifds,ifde,ir
1 format('module_fr_fire_driver: incompatible bounds ',a,' atm ',i5,':',i5,' fire ',i5,':',i5,' ratio ',i3)
!$OMP END CRITICAL(FIRE_DRIVER_CRIT)
call crash(msg)
endif
end subroutine check_fmesh
!
!*****************************
!
subroutine set_flags(config_flags) 1
implicit none
TYPE (grid_config_rec_type) , INTENT(IN) :: config_flags
! copy flags from wrf to module_fr_fire_util
! for instructions how to add a flag see the top of module_fr_fire_util.F
fire_print_msg = config_flags%fire_print_msg
fire_print_file = config_flags%fire_print_file
fuel_left_method = config_flags%fire_fuel_left_method
fuel_left_irl = config_flags%fire_fuel_left_irl
fuel_left_jrl = config_flags%fire_fuel_left_jrl
fire_const_time = config_flags%fire_const_time
fire_const_grnhfx = config_flags%fire_const_grnhfx
fire_const_grnqfx = config_flags%fire_const_grnqfx
fire_atm_feedback = config_flags%fire_atm_feedback
boundary_guard = config_flags%fire_boundary_guard
fire_back_weight = config_flags%fire_back_weight
fire_grows_only = config_flags%fire_grows_only
fire_upwinding = config_flags%fire_upwinding
fire_upwind_split = config_flags%fire_upwind_split
fire_viscosity = config_flags%fire_viscosity
fire_lfn_ext_up = config_flags%fire_lfn_ext_up
fire_test_steps = config_flags%fire_test_steps
!fire_topo_from_atm = config_flags%fire_topo_from_atm
fire_advection = config_flags%fire_advection
end subroutine set_flags
end module module_fr_fire_driver