standalone/fire.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,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.,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
  63 call read_info(inputfile,info)             ! get dimensions
  64
  65 ! start empty NetCDF file with the dimensions
  66 call create_output_file(outputfile,info)
  67
  68 ! get ignition data
  69 call fire_ignition_convert (config_flags,ignition,                   &
  70              grid%fxlong, grid%fxlat,                                &
  71              ifds,ifde, jfds,jfde,                                   &
  72              ifms,ifme, jfms,jfme,                                   &
  73              ifps,ifpe, jfps,jfpe )
  74
  75
  76 ! set dimensions
  77 ifds=1
  78 ifde=info%nfirex
  79 jfds=1
  80 jfde=info%nfirey
  81 ifms=ifds-nhalo
  82 ifme=ifde+nhalo
  83 jfms=ifds-nhalo
  84 jfme=ifde+nhalo
  85 ifps=1
  86 ifpe=ifde
  87 jfps=1
  88 jfpe=jfde
  89 ifts=1
  90 ifte=ifde
  91 jfts=1
  92 jfte=jfde
  93
  94 write(6,2)'fire domain dimensions       ',ifds,ifde,jfds,jfde
  95 write(6,2)'fire memory dimensions       ',ifms,ifme,jfms,jfme
  96 write(6,2)'fire patch  dimensions       ',ifps,ifpe,jfps,jfpe
  97 write(6,2)'fire tile   dimensions       ',ifts,ifte,jfts,jfte
  98 2 format(a,6i6)
  99
 100 ! allocate
 101
 102 ! inputs
 103 call allocate2d(grid%uf,ifms,ifme,jfms,jfme,'uf')              ! fire winds
 104 call allocate2d(grid%vf,ifms,ifme,jfms,jfme,'vf')              ! fire winds
 105 call allocate2d(grid%zsf,ifms,ifme,jfms,jfme,'zsf')             ! terrain height
 106 call allocate2d(grid%dzdxf,ifms,ifme,jfms,jfme,'dzdxf')           ! terrain grad
 107 call allocate2d(grid%dzdyf,ifms,ifme,jfms,jfme,'dzdyf')           ! terrain grad
 108 call allocate2d(grid%fxlong,ifms,ifme,jfms,jfme,'fxlong')          !
 109 call allocate2d(grid%fxlat,ifms,ifme,jfms,jfme,'fxlat')           !
 110 call allocate2d(grid%nfuel_cat,ifms,ifme,jfms,jfme,'nfuel_cat')          !
 111 call allocate2d(grid%fmc_g,ifms,ifme,jfms,jfme,'fmc_g')          !
 112
 113 ! state
 114 call allocate2d(grid%bbb,ifms,ifme,jfms,jfme,'bbb')             ! spread formula coeff
 115 call allocate2d(grid%betafl,ifms,ifme,jfms,jfme,'betafl')          ! spread formula coeff
 116 call allocate2d(grid%phiwc,ifms,ifme,jfms,jfme,'phiwc')           ! spread formula coeff
 117 call allocate2d(grid%phisc,ifms,ifme,jfms,jfme,'phisc')           ! spread formula coeff
 118 call allocate2d(grid%r_0,ifms,ifme,jfms,jfme,'r_0')             ! spread formula coeff
 119 call allocate2d(grid%fgip,ifms,ifme,jfms,jfme,'fgip')            ! spread formula coeff
 120 call allocate2d(grid%ischap,ifms,ifme,jfms,jfme,'ischap')          ! spread formula coeff
 121 call allocate2d(grid%fuel_time,ifms,ifme,jfms,jfme,'fuel_time')        !
 122 call allocate2d(grid%lfn,ifms,ifme,jfms,jfme,'lfn')
 123 call allocate2d(grid%tign_g,ifms,ifme,jfms,jfme,'tign_g')
 124 call allocate2d(grid%fuel_frac,ifms,ifme,jfms,jfme,'fuel_frac')
 125 call allocate2d(grid%fuel_frac_burnt,ifms,ifme,jfms,jfme,'fuel_frac_burnt')
 126 call allocate2d(grid%lfn_out,ifms,ifme,jfms,jfme,'lfn_out')
 127
 128 ! outputs
 129 call allocate2d(grid%fire_area,ifms,ifme,jfms,jfme,'fire_area')
 130 call allocate2d(grid%ros,ifms,ifme,jfms,jfme,'ros')
 131 call allocate2d(grid%flineint,ifms,ifme,jfms,jfme,'flineint')
 132 call allocate2d(grid%flineint2,ifms,ifme,jfms,jfme,'flineint2')
 133 call allocate2d(grid%fgrnhfx,ifms,ifme,jfms,jfme,'fgrnhfx')          !
 134 call allocate2d(grid%fgrnqfx,ifms,ifme,jfms,jfme,'fgrnqfx')          !
 135 call allocate2d(grid%fcanhfx,ifms,ifme,jfms,jfme,'fcanhfx')          !
 136 call allocate2d(grid%fcanqfx,ifms,ifme,jfms,jfme,'fcanqfx')          !
 137 call allocate2d(grid%f_ros,ifms,ifme,jfms,jfme,'f_ros')              !
 138 call allocate2d(grid%f_ros0,ifms,ifme,jfms,jfme,'f_ros0')            !
 139 call allocate2d(grid%f_rosx,ifms,ifme,jfms,jfme,'f_rosx')            !
 140 call allocate2d(grid%f_rosy,ifms,ifme,jfms,jfme,'f_rosy')            !
 141 call allocate2d(grid%f_lineint,ifms,ifme,jfms,jfme,'f_lineint')      !
 142 call allocate2d(grid%f_lineint2,ifms,ifme,jfms,jfme,'f_lineint2')    !
 143 call allocate2d(grid%f_int,ifms,ifme,jfms,jfme,'f_int')              !
 144
 145 ! local
 146 call allocate2d(uf1,ifms,ifme,jfms,jfme,'uf1')              ! fire winds
 147 call allocate2d(vf1,ifms,ifme,jfms,jfme,'vf1')              ! fire winds
 148 call allocate2d(uf2,ifms,ifme,jfms,jfme,'uf2')              ! fire winds
 149 call allocate2d(vf2,ifms,ifme,jfms,jfme,'vf2')              ! fire winds
 150 call allocate2d(fmc_g1,ifms,ifme,jfms,jfme,'fmc_g1')              ! moisture
 151 call allocate2d(fmc_g2,ifms,ifme,jfms,jfme,'fmc_g2')              ! moisture
 152
 153 ! copy pointers to grid fields, to pass to the spread rate calculation
 154 call set_fp_from_grid(grid,fp)
 155 call init_fuel_cats(.true.)
 156
 157 ! time control
 158 ! NOTE: dt in the netcdf input file as returned in info%dt is WRONG !!
 159 dt_double=config_flags%time_step
 160 if(config_flags%time_step_fract_den.ne.0)then
 161   dt_double=dt_double+dble(config_flags%time_step_fract_num)/dble(config_flags%time_step_fract_den)
 162 endif
 163 duration_s = config_flags%run_seconds           &
 164            + 60d0*(config_flags%run_minutes     &
 165            + 60d0*(config_flags%run_hours       &
 166            + 24d0*(config_flags%run_days)))
 167
 168 if(config_flags%history_interval.ne.0)config_flags%history_interval_m=config_flags%history_interval
 169 frame_s = config_flags%history_interval_s           &
 170            + 60d0*(config_flags%history_interval_m     &
 171            + 60d0*(config_flags%history_interval_h       &
 172            + 24d0*(config_flags%history_interval_d)))
 173
 174 nsteps = nint( frame_s / dt_double ) ! number of time steps for the duration
 175 dt_double = frame_s / nsteps
 176 dt = dt_double
 177
 178 write(*,'(a,f10.3,a,i6,a,f10.3,a)')'frame ',frame_s,'s ',nsteps,' time steps at ',dt_double,'s'
 179
 180 ! divide up for shared memory parallel execution
 181 !!call set_tiles(1,1,ips,ipe,jps,jpe,grid%num_tiles,grid%i_start,grid%i_end,grid%j_start,grid%j_end)
 182
 183 ! set the scalars in grid type
 184 grid%dt = dt
 185 grid%itimestep=0
 186 grid%u_frame=0.
 187 grid%v_frame=0.
 188
 189 ! start output file
 190 !! call create_output_file(outputfile,info)
 191
 192 if(info%ntimes.lt.3)then
 193   !write(*,'(a,i5)')'ntimes=',info%ntimes
 194   !call crash('need at least 3 steps')
 195   uniform=.true.
 196   call read_vars(inputfile,info,1,grid)
 197   iframe_start=1
 198   iframe_end=int(duration_s/frame_s)
 199 else
 200   uniform=.false.
 201   call read_vars(inputfile,info,2,grid)
 202   iframe_start=3
 203   iframe_end=info%ntimes
 204   uf1=grid%uf
 205   vf1=grid%vf
 206   fmc_g1=grid%fmc_g
 207 endif
 208 itimestep = 0
 209 ifun_start=1
 210 do iframe=iframe_start,iframe_end ! interval ending with iframe
 211   if(.not.uniform)then
 212     call read_vars(inputfile,info,iframe,grid)
 213     uf2=grid%uf
 214     vf2=grid%vf
 215     fmc_g2=grid%fmc_g
 216   endif
 217   do istep=1,nsteps
 218     itimestep=info%ntimes * (iframe - 1) + istep
 219     grid%itimestep = itimestep
 220     id=itimestep
 221     ifun_end=6
 222     ! interpolate time
 223     time_start = dt_double * (nsteps * (iframe - 1) + istep - 1)
 224     ! interpolate wind
 225     if(.not.uniform)then
 226       t = (istep - 1.)/real(nsteps)
 227       write(*,'(a,i4,a,i3,a,i8,a,f10.3,a,f10.3)')'frame',iframe,' step',istep,' id',id, &
 228          ' start at ',time_start,'s t=',t
 229       grid%uf = (1. - t)*uf1 + t*uf2
 230       grid%vf = (1. - t)*vf1 + t*vf2
 231       grid%fmc_g = (1. - t)*fmc_g1 + t*fmc_g2
 232     endif
 233
 234     do ifun=ifun_start,ifun_end
 235
 236       if(ifun.eq.4)then
 237         call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fp%fmc_g,'fire:fmc_g')
 238       endif
 239
 240       call sfire_model (                    &
 241         id,                                     & ! unique number for prints and debug
 242         ifun,                                   & ! what to do see below
 243         restart,                                & ! use existing state
 244         run_fuel_moisture,                      & ! run the moisture model
 245         config_flags%fire_fuel_read,config_flags%fire_fuel_cat,   & ! legacy initial constant fuel category
 246         ifds,ifde,jfds,jfde,                    & ! fire domain dims - the whole domain
 247         ifms,ifme,jfms,jfme,                    & ! fire memory dims - how declared
 248         ifps,ifpe,jfps,jfpe,                    & ! patch - nodes owned by this process
 249         ifts,ifte,jfts,jfte,                    & ! fire tile dims  - this thread
 250         time_start,dt,                          & ! time and increment
 251         info%fdx,info%fdy,                                & ! fire mesh spacing,
 252         ignition,hfx,                               & ! small array of ignition line descriptions
 253         grid%fxlong,grid%fxlat,                      & ! fire mesh coordinates
 254         grid%fire_hfx,                          & ! given heat flux (experimental)
 255         grid%lfn,grid%lfn_out,grid%tign_g,grid%fuel_frac,grid%fire_area,   & ! state: level function, ign time, fuel left, area burning
 256         grid%fuel_frac_burnt,                   &
 257         grid%fgrnhfx,grid%fgrnqfx,                          & ! output: heat fluxes
 258         grid%ros,grid%flineint,grid%flineint2,                 & ! diagnostic variables
 259         grid%f_ros0,grid%f_rosx,grid%f_rosy,grid%f_ros,             & ! fire risk spread
 260         grid%f_int,grid%f_lineint,grid%f_lineint2,             & ! fire risk intensities
 261         grid%nfuel_cat,                              & ! fuel data per point
 262         grid%fuel_time,grid%fwh,grid%fz0,                      & ! save derived internal data
 263         fp &
 264       )
 265
 266     enddo
 267     ifun_start=3
 268   enddo
 269   call write_vars(outputfile,grid,info,iframe)
 270   if(.not.uniform)then
 271     uf1=uf2
 272     vf1=vf2
 273     fmc_g1=fmc_g2
 274   endif
 275 enddo
 276
 277 end subroutine sub_main
 278
 279
 280 !subroutine model_driver(grid,config_flags)
 281
 282 !
 283 !******************************
 284 !
 285
 286 subroutine set_tiles(itiles,jtiles,ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 287 !*** set tiles for standalone/testing
 288 implicit none
 289 !*** arguments
 290 integer,intent(in)::itiles,jtiles,ips,ipe,jps,jpe
 291 integer,intent(out)::num_tiles
 292 integer,intent(out),dimension(itiles*jtiles)::i_start,i_end,j_start,j_end
 293 !*** local
 294 integer::i,j,istep,jstep,ij
 295 character(len=128)::msg
 296 write(msg,1)'patch',ips,':',ipe,jps,':',jpe
 297 1 format(a,5x,i6,a,2i6,a,i6)
 298 call message(msg,level=-1)
 299 !if(ips.ge.ipe.or.jps.ge.jpe)call crash('bad domain bounds')
 300 !num_tiles=itiles*jtiles
 301 !istep=(ipe-ips+itiles)/itiles
 302 !jstep=(jpe-jps+jtiles)/jtiles
 303 !do i=1,itiles
 304 !    do j=1,jtiles
 305 !        ij=j+(i-1)*jtiles
 306 !        i_start(ij)=min(ipe,ips+(i-1)*istep)
 307 !        i_end(ij)  =min(ipe,ips+(i  )*istep-1)
 308 !        j_start(ij)=min(jpe,jps+(j-1)*jstep)
 309 !        j_end(ij)  =min(jpe,jps+(j  )*jstep-1)
 310 !    enddo
 311 !enddo
 312 !call check_tiles(ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 313 end subroutine set_tiles
 314
 315
 316 subroutine check_tiles(ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 317 implicit none
 318 !*** purpose: check if tiles fit
 319 !*** arguments
 320 integer,intent(in)::ips,ipe,jps,jpe,num_tiles
 321 integer,intent(in),dimension(num_tiles)::i_start,i_end,j_start,j_end
 322 !*** local
 323 character(len=128)::msg
 324 integer:: ij,ie
 325 !*** executable
 326 if(num_tiles.lt.1)call crash('check_tiles: need at least one tile')
 327 ie=0
 328 do ij=1,num_tiles
 329     if(i_start(ij).lt.ips.or.i_end(ij).gt.ipe &
 330     .or.j_start(ij).lt.jps.or.j_end(ij).gt.jpe)then
 331         write(msg,1)'patch',ips,':',ipe,jps,':',jpe
 332 1       format(a,5x,i6,a,2i6,a,i6)
 333         call message(msg,level=-1)
 334         write(msg,2)'tile',ij,i_start(ij),':',i_end(ij),j_start(ij),':',j_end(ij)
 335 2       format(a,2i6,a,2i6,a,i6)
 336         call message(msg,level=-1)
 337         call crash('bad tile bounds')
 338     endif
 339 enddo
 340 end subroutine check_tiles
 341
 342
 343 subroutine allocate2d(p,ims,ime,jms,jme,s)
 344 !*** allocate a pointer with error checking and initialization
 345 implicit none
 346 !*** arguments
 347 real, pointer, dimension(:,:)::p
 348 integer, intent(in):: ims,ime,jms,jme
 349 character(len=*),intent(in)::s
 350 !*** local
 351 integer::err
 352 !*** executable
 353 if(debug_print)write(6,1) ims,ime,jms,jme,trim(s)
 354 if(associated(p))call crash('already allocated')
 355 1 format('allocate2d',2(1x,i6,' :',i6),1x,a)
 356 allocate(p(ims:ime,jms:jme),stat=err)
 357 if(err.ne.0)then
 358    write(6,1)ims,ime,jms,jme,trim(s)
 359    call crash('memory allocation failed')
 360 endif
 361 p=0.
 362 end subroutine allocate2d
 363
 364 subroutine allocate3d(p,ims,ime,jms,jme,kms,kme,s)
 365 !*** allocate a pointer with error checking and initialization
 366 implicit none
 367 !*** arguments
 368 real, pointer, dimension(:,:,:)::p
 369 integer, intent(in):: ims,ime,jms,jme,kms,kme
 370 character(len=*),intent(in)::s
 371 !*** local
 372 integer::err
 373 !*** executable
 374 if(debug_print)write(6,1) ims,ime,jms,jme,kms,kme,trim(s)
 375 1 format('allocate3d',3(1x,i6,' :',i6),1x,a)
 376 if(associated(p))call crash('already allocated')
 377 allocate(p(ims:ime,jms:jme,kms:kme),stat=err)
 378 if(err.ne.0)then
 379    write(6,1)ims,ime,jms,jme,kms,kme,trim(s)
 380    call crash('memory allocation failed')
 381 endif
 382 p=0.
 383 end subroutine allocate3d
 384
 385 end module module_fire_standalone
 386
 387 !
 388 !******************************
 389 !
 390
 391
 392 program fire
 393 use module_fire_standalone, only: sub_main
 394 call  sub_main
 395 end program fire