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.,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_double=config_flags%time_step
 156 if(config_flags%time_step_fract_den.ne.0)then
 157   dt_double=dt_double+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 nsteps = nint( frame_s / dt_double ) ! number of time steps for the duration
 171 dt_double = frame_s / nsteps
 172 dt = dt_double
 173
 174 write(*,'(a,f10.3,a,i6,a,f10.3,a)')'fire.exe: frame ',frame_s,'s ',nsteps,' time steps at ',dt_double,'s'
 175
 176 ! divide up for shared memory parallel execution
 177 !!call set_tiles(1,1,ips,ipe,jps,jpe,grid%num_tiles,grid%i_start,grid%i_end,grid%j_start,grid%j_end)
 178
 179 ! set the scalars in grid type
 180 grid%dt = dt
 181 grid%itimestep=0
 182 grid%xtime=0.
 183 grid%u_frame=0.
 184 grid%v_frame=0.
 185 info%dt = dt     ! dt may be different than it was in the input file
 186
 187 ! start output file
 188 call create_output_file(outputfile,info)
 189
 190 if(info%ntimes.lt.3)then
 191   !write(*,'(a,i5)')'ntimes=',info%ntimes
 192   !call crash('need at least 3 steps')
 193   uniform=.true.
 194   call read_vars(inputfile,info,1,grid)
 195   iframe_start=1
 196   iframe_end=int(duration_s/frame_s)
 197 else
 198   uniform=.false.
 199   call read_vars(inputfile,info,2,grid)
 200   iframe_start=3
 201   iframe_end=info%ntimes
 202   uf1=grid%uf
 203   vf1=grid%vf
 204   fmc_g1=grid%fmc_g
 205 endif
 206
 207 print *,'fxlat  lower bounds:',lbound(grid%fxlat)
 208 print *,'fxlat  upper bounds:',ubound(grid%fxlat)
 209 print *,'fxlat(1,1)=',grid%fxlat(1,1),' fxlat(',ifpe,',',jfpe,')=',grid%fxlat(ifpe,jfpe)
 210 print *,'fxlong lower bounds:',lbound(grid%fxlong)
 211 print *,'fxlong upper bounds:',ubound(grid%fxlong)
 212 print *,'fxlong(1,1)=',grid%fxlong(1,1),' fxlong(',ifpe,',',jfpe,')=',grid%fxlong(ifpe,jfpe)
 213 call print_2d_stats(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,grid%fxlong,'fire:fxlong')
 214 call print_2d_stats(ifps,ifpe,jfps,jfpe,ifms,ifme,jfms,jfme,grid%fxlat,'fire:fxlat')
 215
 216 ! get ignition data - should have fxlong fxlat now
 217 call fire_ignition_convert (config_flags,ignition,                   &
 218              grid%fxlong, grid%fxlat,                                &
 219              ifds,ifde, jfds,jfde,                                   &
 220              ifms,ifme, jfms,jfme,                                   &
 221              ifps,ifpe, jfps,jfpe )
 222
 223 itimestep = 0
 224 ifun_start=1
 225 do iframe=iframe_start,iframe_end ! interval ending with iframe
 226   if(.not.uniform)then
 227     call read_vars(inputfile,info,iframe,grid)
 228     uf2=grid%uf
 229     vf2=grid%vf
 230     fmc_g2=grid%fmc_g
 231   endif
 232   do istep=1,nsteps
 233     itimestep=info%ntimes * (iframe - 1) + istep
 234     grid%itimestep = itimestep
 235     grid%xtime = itimestep * grid%dt / 60.
 236     id=itimestep
 237     ifun_end=6
 238     ! interpolate time
 239     time_start = dt_double * (nsteps * (iframe - 1) + istep - 1)
 240     ! interpolate wind
 241     if(.not.uniform)then
 242       t = (istep - 1.)/real(nsteps)
 243       write(*,'(a,i4,a,i3,a,i8,a,f10.3,a,f10.3)')'frame',iframe,' step',istep,' id',id, &
 244          ' start at ',time_start,'s t=',t
 245       grid%uf = (1. - t)*uf1 + t*uf2
 246       grid%vf = (1. - t)*vf1 + t*vf2
 247       grid%fmc_g = (1. - t)*fmc_g1 + t*fmc_g2
 248     endif
 249
 250     do ifun=ifun_start,ifun_end
 251
 252       if(ifun.eq.4)then
 253         call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fp%fmc_g,'fire:fmc_g')
 254       endif
 255
 256       call sfire_model (                    &
 257         id,                                     & ! unique number for prints and debug
 258         ifun,                                   & ! what to do see below
 259         restart,replay,                         & ! use existing state; prescribe fire arrival time
 260         run_fuel_moisture,                      & ! run the moisture model
 261         config_flags%fire_fuel_read,config_flags%fire_fuel_cat,   & ! legacy initial constant fuel category
 262         ifds,ifde,jfds,jfde,                    & ! fire domain dims - the whole domain
 263         ifms,ifme,jfms,jfme,                    & ! fire memory dims - how declared
 264         ifps,ifpe,jfps,jfpe,                    & ! patch - nodes owned by this process
 265         ifts,ifte,jfts,jfte,                    & ! fire tile dims  - this thread
 266         time_start,dt,                          & ! time and increment
 267         info%fdx,info%fdy,                                & ! fire mesh spacing,
 268         ignition,hfx,                               & ! small array of ignition line descriptions
 269         grid%fxlong,grid%fxlat,                      & ! fire mesh coordinates
 270         grid%fire_hfx,                          & ! given heat flux (experimental)
 271         grid%tign_in,                           & ! ignition time, if given
 272         grid%lfn,grid%lfn_out,grid%tign_g,grid%fuel_frac,grid%fire_area,   & ! state: level function, ign time, fuel left, area burning
 273         grid%fuel_frac_burnt,                   &
 274         grid%fgrnhfx,grid%fgrnqfx,                          & ! output: heat fluxes
 275         grid%ros,grid%flineint,grid%flineint2,                 & ! diagnostic variables
 276         grid%f_ros0,grid%f_rosx,grid%f_rosy,grid%f_ros,             & ! fire risk spread
 277         grid%f_int,grid%f_lineint,grid%f_lineint2,             & ! fire risk intensities
 278         grid%nfuel_cat,                              & ! fuel data per point
 279         grid%fuel_time,grid%fwh,grid%fz0,                      & ! save derived internal data
 280         fp &
 281       )
 282
 283     enddo
 284     ifun_start=3
 285   enddo
 286   call write_vars(outputfile,grid,info,iframe)
 287   if(.not.uniform)then
 288     uf1=uf2
 289     vf1=vf2
 290     fmc_g1=fmc_g2
 291   endif
 292 enddo
 293
 294 end subroutine sub_main
 295
 296
 297 !subroutine model_driver(grid,config_flags)
 298
 299 !
 300 !******************************
 301 !
 302
 303 subroutine set_tiles(itiles,jtiles,ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 304 !*** set tiles for standalone/testing
 305 implicit none
 306 !*** arguments
 307 integer,intent(in)::itiles,jtiles,ips,ipe,jps,jpe
 308 integer,intent(out)::num_tiles
 309 integer,intent(out),dimension(itiles*jtiles)::i_start,i_end,j_start,j_end
 310 !*** local
 311 integer::i,j,istep,jstep,ij
 312 character(len=128)::msg
 313 write(msg,1)'patch',ips,':',ipe,jps,':',jpe
 314 1 format(a,5x,i6,a,2i6,a,i6)
 315 call message(msg,level=-1)
 316 !if(ips.ge.ipe.or.jps.ge.jpe)call crash('bad domain bounds')
 317 !num_tiles=itiles*jtiles
 318 !istep=(ipe-ips+itiles)/itiles
 319 !jstep=(jpe-jps+jtiles)/jtiles
 320 !do i=1,itiles
 321 !    do j=1,jtiles
 322 !        ij=j+(i-1)*jtiles
 323 !        i_start(ij)=min(ipe,ips+(i-1)*istep)
 324 !        i_end(ij)  =min(ipe,ips+(i  )*istep-1)
 325 !        j_start(ij)=min(jpe,jps+(j-1)*jstep)
 326 !        j_end(ij)  =min(jpe,jps+(j  )*jstep-1)
 327 !    enddo
 328 !enddo
 329 !call check_tiles(ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 330 end subroutine set_tiles
 331
 332
 333 subroutine check_tiles(ips,ipe,jps,jpe,num_tiles,i_start,i_end,j_start,j_end)
 334 implicit none
 335 !*** purpose: check if tiles fit
 336 !*** arguments
 337 integer,intent(in)::ips,ipe,jps,jpe,num_tiles
 338 integer,intent(in),dimension(num_tiles)::i_start,i_end,j_start,j_end
 339 !*** local
 340 character(len=128)::msg
 341 integer:: ij,ie
 342 !*** executable
 343 if(num_tiles.lt.1)call crash('check_tiles: need at least one tile')
 344 ie=0
 345 do ij=1,num_tiles
 346     if(i_start(ij).lt.ips.or.i_end(ij).gt.ipe &
 347     .or.j_start(ij).lt.jps.or.j_end(ij).gt.jpe)then
 348         write(msg,1)'patch',ips,':',ipe,jps,':',jpe
 349 1       format(a,5x,i6,a,2i6,a,i6)
 350         call message(msg,level=-1)
 351         write(msg,2)'tile',ij,i_start(ij),':',i_end(ij),j_start(ij),':',j_end(ij)
 352 2       format(a,2i6,a,2i6,a,i6)
 353         call message(msg,level=-1)
 354         call crash('bad tile bounds')
 355     endif
 356 enddo
 357 end subroutine check_tiles
 358
 359
 360 subroutine allocate2d(p,ims,ime,jms,jme,s)
 361 !*** allocate a pointer with error checking and initialization
 362 implicit none
 363 !*** arguments
 364 real, pointer, dimension(:,:)::p
 365 integer, intent(in):: ims,ime,jms,jme
 366 character(len=*),intent(in)::s
 367 !*** local
 368 integer::err
 369 !*** executable
 370 if(debug_print)write(6,1) ims,ime,jms,jme,trim(s)
 371 if(associated(p))call crash(trim(s) // ' already allocated')
 372 1 format('allocate2d',2(1x,i6,' :',i6),1x,a)
 373 allocate(p(ims:ime,jms:jme),stat=err)
 374 if(err.ne.0)then
 375    write(6,1)ims,ime,jms,jme,trim(s)
 376    call crash('memory allocation failed')
 377 endif
 378 p=0.
 379 end subroutine allocate2d
 380
 381 subroutine allocate3d(p,ims,ime,jms,jme,kms,kme,s)
 382 !*** allocate a pointer with error checking and initialization
 383 implicit none
 384 !*** arguments
 385 real, pointer, dimension(:,:,:)::p
 386 integer, intent(in):: ims,ime,jms,jme,kms,kme
 387 character(len=*),intent(in)::s
 388 !*** local
 389 integer::err
 390 !*** executable
 391 if(debug_print)write(6,1) ims,ime,jms,jme,kms,kme,trim(s)
 392 1 format('allocate3d',3(1x,i6,' :',i6),1x,a)
 393 if(associated(p))call crash('already allocated')
 394 allocate(p(ims:ime,jms:jme,kms:kme),stat=err)
 395 if(err.ne.0)then
 396    write(6,1)ims,ime,jms,jme,kms,kme,trim(s)
 397    call crash('memory allocation failed')
 398 endif
 399 p=0.
 400 end subroutine allocate3d
 401
 402 end module module_fire_standalone
 403
 404 !
 405 !******************************
 406 !
 407
 408
 409 program fire
 410 use module_fire_standalone, only: sub_main
 411 call  sub_main
 412 end program fire