standalone/fire_main.F

   1 module module_fire_standalone
   2
   3 use module_fr_sfire_driver, only: set_flags, fire_ignition_convert, &
   4                                   set_fp_from_grid
   5 use module_fr_sfire_util, only: message,crash, &
   6           lines_type, print_2d_stats
   7 use module_fr_sfire_phys, only: fire_params, init_fuel_cats
   8 use module_fr_sfire_model, only: sfire_model
   9 use module_domain, only: domain
  10 use module_configure, only: grid_config_rec_type,read_namelist
  11 use wrf_netcdf, only : grid_info, read_info, &
  12                        create_output_file,write_vars, &
  13                        read_vars, debug_print
  14 implicit none
  15
  16 contains
  17
  18 subroutine sub_main
  19
  20 !*** purpose: standalone driver with compatible files to WRF-Fire
  21
  22 implicit none
  23
  24 !*** local
  25
  26 ! arguments to SFIRE
  27
  28 type(domain)::grid          ! all: state+inputs+outputs, compatible with wrf
  29 TYPE (grid_config_rec_type):: config_flags ! the namelist
  30 integer::  &                ! fire mesh dimensions
  31     ifds,ifde,jfds,jfde, &  ! the physical domain
  32     ifps,ifpe,jfps,jfpe, &  ! patch - assigned to one process. Here the same as domain.
  33     ifts,ifte,jfts,jfte, &  ! memory allocated, needs a strip around the patch
  34     ifms,ifme,jfms,jfme     ! memory allocated, needs a strip around the patch
  35
  36 ! I/O interface
  37 character(len=*),parameter::inputfile='fire_input.nc'
  38 character(len=*),parameter::outputfile='fire_output.nc'
  39 real, pointer, dimension(:,:) ::  uf1, vf1, uf2, vf2, fmc_g1, fmc_g2   ! stored input fields
  40
  41 ! other derived types
  42 type(grid_info)::info                    ! dimensions, grid controls
  43
  44 ! scalars
  45 integer:: nsteps,itimestep,ifun_start,ifun_end,id,ifun,iframe,istep
  46 integer::nhalo=5
  47 double precision:: dt_double,dt_config,duration_s,frame_s  ! may need more accurate time computation to get the number of timesteps right
  48 real:: time_start,dt,t
  49 logical::do_ouput
  50 TYPE(lines_type) :: ignition, hfx
  51 type(fire_params)::fp
  52 logical::restart=.false.,replay=.false.,uniform=.false.
  53 integer::iframe_start,iframe_end
  54 logical::run_fuel_moisture=.false.
  55
  56 !*** executable
  57
  58 call read_namelist(config_flags)           ! read flags from namelist.input
  59 call set_flags(config_flags)               ! copy configuration flags to sfire internal structures
  60
  61 debug_print = config_flags%fire_print_msg.ge.2 ! if we write a lot
  62 debug_print = config_flags%fire_print_msg.ge.1 ! if we write a lot
  63
  64 call read_info(inputfile,info)             ! get dimensions
  65
  66 ! start empty NetCDF file with the dimensions
  67 ! not here, may want to overwrite something
  68 ! call create_output_file(outputfile,info)
  69
  70
  71
  72 ! set dimensions
  73 ifds=1
  74 ifde=info%nfirex
  75 jfds=1
  76 jfde=info%nfirey
  77 ifms=ifds-nhalo
  78 ifme=ifde+nhalo
  79 jfms=ifds-nhalo
  80 jfme=ifde+nhalo
  81 ifps=1
  82 ifpe=ifde
  83 jfps=1
  84 jfpe=jfde
  85 ifts=1
  86 ifte=ifde
  87 jfts=1
  88 jfte=jfde
  89
  90 write(6,2)'fire domain dimensions       ',ifds,ifde,jfds,jfde
  91 write(6,2)'fire memory dimensions       ',ifms,ifme,jfms,jfme
  92 write(6,2)'fire patch  dimensions       ',ifps,ifpe,jfps,jfpe
  93 write(6,2)'fire tile   dimensions       ',ifts,ifte,jfts,jfte
  94 2 format(a,6i6)
  95
  96 ! allocate
  97
  98 ! inputs
  99 call allocate2d(grid%uf,ifms,ifme,jfms,jfme,'uf')              ! fire winds
 100 call allocate2d(grid%vf,ifms,ifme,jfms,jfme,'vf')              ! fire winds
 101 call allocate2d(grid%zsf,ifms,ifme,jfms,jfme,'zsf')             ! terrain height
 102 call allocate2d(grid%dzdxf,ifms,ifme,jfms,jfme,'dzdxf')           ! terrain grad
 103 call allocate2d(grid%dzdyf,ifms,ifme,jfms,jfme,'dzdyf')           ! terrain grad
 104 call allocate2d(grid%fxlong,ifms,ifme,jfms,jfme,'fxlong')          !
 105 call allocate2d(grid%fxlat,ifms,ifme,jfms,jfme,'fxlat')           !
 106 call allocate2d(grid%nfuel_cat,ifms,ifme,jfms,jfme,'nfuel_cat')          !
 107 call allocate2d(grid%fmc_g,ifms,ifme,jfms,jfme,'fmc_g')          !
 108
 109 ! state
 110 call allocate2d(grid%bbb,ifms,ifme,jfms,jfme,'bbb')             ! spread formula coeff
 111 call allocate2d(grid%betafl,ifms,ifme,jfms,jfme,'betafl')          ! spread formula coeff
 112 call allocate2d(grid%phiwc,ifms,ifme,jfms,jfme,'phiwc')           ! spread formula coeff
 113 call allocate2d(grid%phisc,ifms,ifme,jfms,jfme,'phisc')           ! spread formula coeff
 114 call allocate2d(grid%r_0,ifms,ifme,jfms,jfme,'r_0')             ! spread formula coeff
 115 call allocate2d(grid%fgip,ifms,ifme,jfms,jfme,'fgip')            ! spread formula coeff
 116 call allocate2d(grid%ischap,ifms,ifme,jfms,jfme,'ischap')          ! spread formula coeff
 117 call allocate2d(grid%fuel_time,ifms,ifme,jfms,jfme,'fuel_time')        !
 118 call allocate2d(grid%lfn,ifms,ifme,jfms,jfme,'lfn')
 119 call allocate2d(grid%tign_g,ifms,ifme,jfms,jfme,'tign_g')
 120 call allocate2d(grid%fuel_frac,ifms,ifme,jfms,jfme,'fuel_frac')
 121 call allocate2d(grid%fuel_frac_burnt,ifms,ifme,jfms,jfme,'fuel_frac_burnt')
 122 call allocate2d(grid%lfn_out,ifms,ifme,jfms,jfme,'lfn_out')
 123
 124 ! outputs
 125 call allocate2d(grid%fire_area,ifms,ifme,jfms,jfme,'fire_area')
 126 call allocate2d(grid%ros,ifms,ifme,jfms,jfme,'ros')
 127 call allocate2d(grid%flineint,ifms,ifme,jfms,jfme,'flineint')
 128 call allocate2d(grid%flineint2,ifms,ifme,jfms,jfme,'flineint2')
 129 call allocate2d(grid%fgrnhfx,ifms,ifme,jfms,jfme,'fgrnhfx')          !
 130 call allocate2d(grid%fgrnqfx,ifms,ifme,jfms,jfme,'fgrnqfx')          !
 131 call allocate2d(grid%fcanhfx,ifms,ifme,jfms,jfme,'fcanhfx')          !
 132 call allocate2d(grid%fcanqfx,ifms,ifme,jfms,jfme,'fcanqfx')          !
 133 call allocate2d(grid%f_ros,ifms,ifme,jfms,jfme,'f_ros')              !
 134 call allocate2d(grid%f_ros0,ifms,ifme,jfms,jfme,'f_ros0')            !
 135 call allocate2d(grid%f_rosx,ifms,ifme,jfms,jfme,'f_rosx')            !
 136 call allocate2d(grid%f_rosy,ifms,ifme,jfms,jfme,'f_rosy')            !
 137 call allocate2d(grid%f_lineint,ifms,ifme,jfms,jfme,'f_lineint')      !
 138 call allocate2d(grid%f_lineint2,ifms,ifme,jfms,jfme,'f_lineint2')    !
 139 call allocate2d(grid%f_int,ifms,ifme,jfms,jfme,'f_int')              !
 140
 141 ! local
 142 call allocate2d(uf1,ifms,ifme,jfms,jfme,'uf1')              ! fire winds
 143 call allocate2d(vf1,ifms,ifme,jfms,jfme,'vf1')              ! fire winds
 144 call allocate2d(uf2,ifms,ifme,jfms,jfme,'uf2')              ! fire winds
 145 call allocate2d(vf2,ifms,ifme,jfms,jfme,'vf2')              ! fire winds
 146 call allocate2d(fmc_g1,ifms,ifme,jfms,jfme,'fmc_g1')              ! moisture
 147 call allocate2d(fmc_g2,ifms,ifme,jfms,jfme,'fmc_g2')              ! moisture
 148
 149 ! copy pointers to grid fields, to pass to the spread rate calculation
 150 call set_fp_from_grid(grid,fp)
 151 call init_fuel_cats(.true.)
 152
 153 ! time control
 154 ! NOTE: dt in the netcdf input file as returned in info%dt is WRONG !!
 155 dt_config=config_flags%time_step
 156 if(config_flags%time_step_fract_den.ne.0)then
 157   dt_config=dt_config+dble(config_flags%time_step_fract_num)/dble(config_flags%time_step_fract_den)
 158 endif
 159 duration_s = config_flags%run_seconds           &
 160            + 60d0*(config_flags%run_minutes     &
 161            + 60d0*(config_flags%run_hours       &
 162            + 24d0*(config_flags%run_days)))
 163
 164 if(config_flags%history_interval.ne.0)config_flags%history_interval_m=config_flags%history_interval
 165 frame_s = config_flags%history_interval_s           &
 166            + 60d0*(config_flags%history_interval_m     &
 167            + 60d0*(config_flags%history_interval_h       &
 168            + 24d0*(config_flags%history_interval_d)))
 169
 170 write(*,'(a,f10.2,a,f10.6,a,f10.6)')'from namelist.input: history_interval=',frame_s, &
 171      ' time_step=',dt_config,'s time steps in history frame ', frame_s / dt_config
 172
 173 nsteps = nint( frame_s / dt_config ) ! number of time steps for the duration
 174 dt_double = frame_s / nsteps
 175 dt = dt_double
 176
 177 write(*,'(a,f10.2,a,f8.4,a,i6)')'adjusted: ',frame_s,' time_step=',dt,'s time steps in history frame ',nsteps
 178
 179 ! divide up for shared memory parallel execution
 180 !!call set_tiles(1,1,ips,ipe,jps,jpe,grid%num_tiles,grid%i_start,grid%i_end,grid%j_start,grid%j_end)
 181
 182 ! set the scalars in grid type
 183 grid%dt = dt
 184 grid%itimestep=0
 185 grid%xtime=0.
 186 grid%u_frame=0.
 187 grid%v_frame=0.
 188 info%dt = dt     ! dt may be different than it was in the input file
 189
 190 ! start output file
 191 call create_output_file(outputfile,info)
 192
 193
 194 if(info%ntimes.lt.3)then
 195   !call crash('need at least 3 steps')
 196   uniform=.true.
 197   call read_vars(inputfile,info,1,grid)
 198   iframe_start=1
 199   iframe_end=int(duration_s/frame_s)
 200 else
 201   uniform=.false.
 202   call read_vars(inputfile,info,2,grid)
 203   iframe_start=3
 204   write(*,'(a)')'starting from frame 3 because need to interpolate between frames and frame 1 may be invalid'
 205   iframe_end=info%ntimes
 206   uf1=grid%uf
 207   vf1=grid%vf
 208   fmc_g1=grid%fmc_g
 209 endif
 210
 211 print *,'Fire mesh:'
 212 print *,'fxlat  lower bounds:',lbound(grid%fxlat)
 213 print *,'fxlat  upper bounds:',ubound(grid%fxlat)
 214 print *,'fxlat(1,1)=',grid%fxlat(1,1),' fxlat(',ifpe,',',jfpe,')=',grid%fxlat(ifpe,jfpe)
 215 print *,'fxlong lower bounds:',lbound(grid%fxlong)
 216 print *,'fxlong upper bounds:',ubound(grid%fxlong)
 217 print *,'fxlong(1,1)=',grid%fxlong(1,1),' fxlong(',ifpe,',',jfpe,')=',grid%fxlong(ifpe,jfpe)
 218 call print_2d_stats(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,grid%fxlong,'fire:fxlong')
 219 call print_2d_stats(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,grid%fxlat,'fire:fxlat')
 220
 221 ! get ignition data - should have fxlong fxlat now
 222 call fire_ignition_convert (config_flags,ignition,                   &
 223              grid%fxlong, grid%fxlat,                                &
 224              ifds,ifde, jfds,jfde,                                   &
 225              ifms,ifme, jfms,jfme,                                   &
 226              ifps,ifpe, jfps,jfpe )
 227
 228 write(*,'(a,i5)')'number of frames in the file: ntimes=',info%ntimes
 229 write(*,'(2(a,i5))')'using frames from',iframe_start,' to',iframe_end
 230
 231 itimestep = 0
 232 ifun_start=1
 233 do iframe=iframe_start,iframe_end ! interval ending with iframe
 234   if(.not.uniform)then
 235     call read_vars(inputfile,info,iframe,grid)
 236     uf2=grid%uf
 237     vf2=grid%vf
 238     fmc_g2=grid%fmc_g
 239   endif
 240   do istep=1,nsteps
 241     itimestep=info%ntimes * (iframe - 1) + istep
 242     grid%itimestep = itimestep
 243     grid%xtime = itimestep * grid%dt / 60.
 244     id=itimestep
 245     ifun_end=6
 246     ! interpolate time
 247     time_start = dt_double * (nsteps * (iframe - 1) + istep - 1)
 248     ! interpolate wind
 249     if(.not.uniform)then
 250       t = (istep - 1.)/real(nsteps)
 251       write(*,'(a,i4,a,i3,a,i8,a,f10.3,a,f10.3)')'frame',iframe,' istep',istep, &
 252          ' itimestep',itimestep, &
 253          ' start at ',time_start,'s weight ',t
 254       grid%uf = (1. - t)*uf1 + t*uf2
 255       grid%vf = (1. - t)*vf1 + t*vf2
 256       grid%fmc_g = (1. - t)*fmc_g1 + t*fmc_g2
 257     endif
 258
 259     do ifun=ifun_start,ifun_end
 260
 261       if(ifun.eq.4)then
 262         call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fp%fmc_g,'fire:fmc_g')
 263       endif
 264
 265       call sfire_model (                    &
 266         id,                                     & ! unique number for prints and debug
 267         ifun,                                   & ! what to do see below
 268         restart,replay,                         & ! use existing state; prescribe fire arrival time
 269         run_fuel_moisture,                      & ! run the moisture model
 270         config_flags%fire_fuel_read,config_flags%fire_fuel_cat,   & ! legacy initial constant fuel category
 271         ifds,ifde,jfds,jfde,                    & ! fire domain dims - the whole domain
 272         ifms,ifme,jfms,jfme,                    & ! fire memory dims - how declared
 273         ifps,ifpe,jfps,jfpe,                    & ! patch - nodes owned by this process
 274         ifts,ifte,jfts,jfte,                    & ! fire tile dims  - this thread
 275         time_start,dt,                          & ! time and increment
 276         info%fdx,info%fdy,                                & ! fire mesh spacing,
 277         ignition,hfx,                               & ! small array of ignition line descriptions
 278         grid%fxlong,grid%fxlat,                      & ! fire mesh coordinates
 279         grid%fire_hfx,                          & ! given heat flux (experimental)
 280         grid%tign_in,                           & ! ignition time, if given
 281         grid%lfn,grid%lfn_out,grid%tign_g,grid%fuel_frac,grid%fire_area,   & ! state: level function, ign time, fuel left, area burning
 282         grid%fuel_frac_burnt,                   &
 283         grid%fgrnhfx,grid%fgrnqfx,                          & ! output: heat fluxes
 284         grid%ros,grid%flineint,grid%flineint2,                 & ! diagnostic variables
 285         grid%f_ros0,grid%f_rosx,grid%f_rosy,grid%f_ros,             & ! fire risk spread
 286         grid%f_int,grid%f_lineint,grid%f_lineint2,             & ! fire risk intensities
 287         grid%nfuel_cat,                              & ! fuel data per point
 288         grid%fuel_time,grid%fwh,grid%fz0,                      & ! save derived internal data
 289         fp &
 290       )
 291
 292     enddo
 293     ifun_start=3
 294   enddo
 295   call write_vars(outputfile,grid,info,iframe)
 296   if(.not.uniform)then
 297     uf1=uf2
 298     vf1=vf2
 299     fmc_g1=fmc_g2
 300   endif
 301 enddo
 302
 303 end subroutine sub_main
 304
 305
 306 !subroutine model_driver(grid,config_flags)
 307
 308 !
 309 !******************************
 310 !
 311
 312 subroutine set_tiles(itiles,jtiles,ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 313 !*** set tiles for standalone/testing
 314 implicit none
 315 !*** arguments
 316 integer,intent(in)::itiles,jtiles,ips,ipe,jps,jpe
 317 integer,intent(out)::num_tiles
 318 integer,intent(out),dimension(itiles*jtiles)::i_start,i_end,j_start,j_end
 319 !*** local
 320 integer::i,j,istep,jstep,ij
 321 character(len=128)::msg
 322 write(msg,1)'patch',ips,':',ipe,jps,':',jpe
 323 1 format(a,5x,i6,a,2i6,a,i6)
 324 call message(msg,level=-1)
 325 !if(ips.ge.ipe.or.jps.ge.jpe)call crash('bad domain bounds')
 326 num_tiles=itiles*jtiles
 327 istep=(ipe-ips+itiles)/itiles
 328 jstep=(jpe-jps+jtiles)/jtiles
 329 do i=1,itiles
 330     do j=1,jtiles
 331         ij=j+(i-1)*jtiles
 332         i_start(ij)=min(ipe,ips+(i-1)*istep)
 333         i_end(ij)  =min(ipe,ips+(i  )*istep-1)
 334         j_start(ij)=min(jpe,jps+(j-1)*jstep)
 335         j_end(ij)  =min(jpe,jps+(j  )*jstep-1)
 336     enddo
 337 enddo
 338 call check_tiles(ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 339 end subroutine set_tiles
 340
 341
 342 subroutine check_tiles(ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 343 implicit none
 344 !*** purpose: check if tiles fit
 345 !*** arguments
 346 integer,intent(in)::ips,ipe,jps,jpe,num_tiles
 347 integer,intent(in),dimension(num_tiles)::i_start,i_end,j_start,j_end
 348 !*** local
 349 character(len=128)::msg
 350 integer:: ij,ie
 351 !*** executable
 352 if(num_tiles.lt.1)call crash('check_tiles: need at least one tile')
 353 ie=0
 354 do ij=1,num_tiles
 355     if(i_start(ij).lt.ips.or.i_end(ij).gt.ipe &
 356     .or.j_start(ij).lt.jps.or.j_end(ij).gt.jpe)then
 357         write(msg,1)'patch',ips,':',ipe,jps,':',jpe
 358 1       format(a,5x,i6,a,2i6,a,i6)
 359         call message(msg,level=-1)
 360         write(msg,2)'tile',ij,i_start(ij),':',i_end(ij),j_start(ij),':',j_end(ij)
 361 2       format(a,2i6,a,2i6,a,i6)
 362         call message(msg,level=-1)
 363         call crash('bad tile bounds')
 364     endif
 365 enddo
 366 end subroutine check_tiles
 367
 368
 369 subroutine allocate2d(p,ims,ime,jms,jme,s)
 370 !*** allocate a pointer with error checking and initialization
 371 implicit none
 372 !*** arguments
 373 real, pointer, dimension(:,:)::p
 374 integer, intent(in):: ims,ime,jms,jme
 375 character(len=*),intent(in)::s
 376 !*** local
 377 integer::err
 378 !*** executable
 379 if(debug_print)write(6,1) ims,ime,jms,jme,trim(s)
 380 if(associated(p))call crash(trim(s) // ' already allocated')
 381 1 format('allocate2d',2(1x,i6,' :',i6),1x,a)
 382 allocate(p(ims:ime,jms:jme),stat=err)
 383 if(err.ne.0)then
 384    write(6,1)ims,ime,jms,jme,trim(s)
 385    call crash('memory allocation failed')
 386 endif
 387 p=0.
 388 end subroutine allocate2d
 389
 390 subroutine allocate3d(p,ims,ime,jms,jme,kms,kme,s)
 391 !*** allocate a pointer with error checking and initialization
 392 implicit none
 393 !*** arguments
 394 real, pointer, dimension(:,:,:)::p
 395 integer, intent(in):: ims,ime,jms,jme,kms,kme
 396 character(len=*),intent(in)::s
 397 !*** local
 398 integer::err
 399 !*** executable
 400 if(debug_print)write(6,1) ims,ime,jms,jme,kms,kme,trim(s)
 401 1 format('allocate3d',3(1x,i6,' :',i6),1x,a)
 402 if(associated(p))call crash('already allocated')
 403 allocate(p(ims:ime,jms:jme,kms:kme),stat=err)
 404 if(err.ne.0)then
 405    write(6,1)ims,ime,jms,jme,kms,kme,trim(s)
 406    call crash('memory allocation failed')
 407 endif
 408 p=0.
 409 end subroutine allocate3d
 410
 411 end module module_fire_standalone
 412
 413 !
 414 !******************************
 415 !
 416
 417
 418 program fire
 419 use module_fire_standalone, only: sub_main
 420 call  sub_main
 421 end program fire