wrfv2_fire/phys/module_fr_sfire_model.F

   1 !
   2 #define DEBUG_OUT
   3
   4 module module_fr_sfire_model
   5
   6 use module_fr_sfire_core
   7 use module_fr_sfire_util
   8 use module_fr_sfire_phys
   9
  10 contains
  11
  12 subroutine sfire_model (                    &
  13     id,                                     & ! unique number for prints and debug
  14     ifun,                                   & ! what to do see below
  15     restart,                                &
  16     need_lfn_update,                          & ! if lfn needs to be synced between tiles
  17     num_ignitions,                          & ! number of ignitions before advancing
  18     ifuelread,nfuel_cat0,                   & ! initialize fuel categories
  19     ifds,ifde,jfds,jfde,                    & ! fire domain dims - the whole domain
  20     ifms,ifme,jfms,jfme,                    & ! fire memory dims - how declared
  21     ifps,ifpe,jfps,jfpe,                    & ! patch - nodes owned by this process
  22     ifts,ifte,jfts,jfte,                    & ! fire tile dims  - this thread
  23     time_start,dt,                          & ! time and increment
  24     fdx,fdy,                                & ! fire mesh spacing,
  25     ignition_line,                          & ! small array of ignition line descriptions
  26     ignitions_done,ignited_tile,            &
  27     coord_xf,coord_yf,unit_xf,unit_yf,      & ! fire mesh coordinates
  28     lfn,lfn_out,tign,fuel_frac,fire_area,   & ! state: level function, ign time, fuel left, area burning
  29     grnhfx,grnqfx,                          & ! output: heat fluxes
  30     ros,                                    & ! output: rate of spread
  31     nfuel_cat,                              & ! fuel data per point
  32     fuel_time,                              & ! save derived internal data
  33     fp &
  34 )
  35
  36 ! This subroutine implements the fire spread model.
  37 ! All quantities are on the fire grid. It inputs
  38 ! winds given on the nodes of the fire grid
  39 ! and outputs the heat fluxes on the cells of the fire grid.
  40 ! This subroutine has no knowledge of any atmospheric model.
  41 ! This code was written to conform with the WRF parallelism model, however it
  42 ! does not depend on it. It can be called with domain equal to tile.
  43 ! Wind and height must be given on 1 more node beyond the domain bounds.
  44 ! The subroutine changes only array entries of the arguments in the tile.
  45 ! Upon exit with ifun=2 (time step), lfn_out is to be copied into lfn by the caller.
  46 ! When this subroutine is used on separate tiles that make a domain the value, the
  47 ! it uses lfn on a strip of width 2 from neighboring tiles.
  48 !
  49 ! All computation is done on one tile.
  50 !
  51 ! This subroutine is intended to be called in a loop like
  52 !
  53 !
  54 ! do ifun=1,6 (if initizalize run, otherwise 3,6)
  55 !   start parallel loop over tiles
  56 !       if ifun=1, set z and fuel data
  57 !       if ifun=3, set the wind arrays
  58 !       call sfire_model(....)
  59 !   end parallel loop over tiles
  60 !
  61 !   if need_lfn_update, halo exchange on lfn width 2
  62 !
  63 !   if ifun=0
  64 !       halo exchange on z width 2
  65 !       halo exchange on fuel data width 1
  66 !   endif
  67 !
  68 !   if ifun=3, halo exchange on winds width 2
  69 !
  70 ! enddo
  71
  72 implicit none
  73
  74 !*** arguments
  75
  76 ! control switches
  77 integer, intent(in) :: id
  78 integer, intent(in) :: ifun                 ! 1 = initialize run pass 1
  79                                             ! 2 = initialize run pass 2
  80                                             ! 3 = initialize timestep
  81                                             ! 4 = do one timestep
  82                                             ! 5 = copy timestep output to input
  83                                             ! 6 = compute output fluxes
  84 logical, intent(in):: restart               ! if true, use existing state
  85 logical, intent(out)::need_lfn_update       ! if true, halo update on lfn afterwards
  86 ! scalar data
  87 integer, intent(in) :: num_ignitions        ! number of ignition lines
  88 integer, intent(in) :: ifuelread,nfuel_cat0 ! for set_fire_params
  89 integer, intent(in) :: ifds,ifde,jfds,jfde,&  ! fire domain bounds
  90         ifps,ifpe,jfps,jfpe                ! patch - nodes owned by this process
  91 integer, intent(in) :: ifts,ifte,jfts,jfte  ! fire tile bounds
  92 integer, intent(in) :: ifms,ifme,jfms,jfme  ! fire memory array bounds
  93 REAL,INTENT(in) :: time_start,dt            ! starting time, time step
  94 REAL,INTENT(in) :: fdx,fdy                  ! spacing of the fire mesh
  95 ! array data
  96 type(ignition_line_type), dimension (num_ignitions), intent(in):: ignition_line ! descriptions of ignition lines
  97 integer, intent(out):: ignited_tile(num_ignitions),ignitions_done
  98 real, dimension(ifms:ifme, jfms:jfme), intent(in):: &
  99     coord_xf,coord_yf                       !  node coordinates
 100 real, intent(in):: unit_xf,unit_yf          !  coordinate units in m
 101
 102 ! state
 103 REAL, INTENT(inout), dimension(ifms:ifme,jfms:jfme):: &
 104     lfn   , &                               ! level function: fire is where lfn<0 (node)
 105     tign  , &                               ! absolute time of ignition (node)
 106     fuel_frac                               ! fuel fraction (node), currently redundant
 107
 108 REAL, INTENT(out), dimension(ifms:ifme,jfms:jfme):: &
 109     fire_area                               ! fraction of each cell burning
 110
 111 ! output
 112 REAL, INTENT(out), dimension(ifms:ifme,jfms:jfme):: &
 113     lfn_out, &                              !
 114     grnhfx,grnqfx, &                        ! heat fluxes J/m^2/s  (cell)
 115     ros                                     ! output: rate of spread
 116
 117 ! constant arrays - set at initialization
 118 real, intent(inout), dimension(ifms:ifme, jfms:jfme)::nfuel_cat ! cell based, data, constant
 119 real,intent(inout),dimension(ifms:ifme,jfms:jfme):: fuel_time
 120 type(fire_params),intent(inout)::fp
 121
 122 !*** local
 123
 124 integer :: xifms,xifme,xjfms,xjfme  ! memory bounds for pass-through arguments to normal spread
 125 real, dimension(ifts:ifte,jfts:jfte)::fuel_frac_burnt,fuel_frac_end
 126 integer::ignited,ig,i,j,itso,iteo,jtso,jteo
 127 real::tbound,err,erri,errj,maxgrad,grad,tfa,thf,mhf,tqf,mqf,aw,mw
 128 character(len=128)::msg
 129 logical:: freeze_fire
 130 integer:: stat_lev=1
 131
 132 !*** executable
 133
 134 call check_mesh_2dim(ifts-1,ifte+1,jfts-1,jfte+1,ifms,ifme,jfms,jfme)
 135
 136 xifms=ifms  ! dimensions for the include file
 137 xifme=ifme
 138 xjfms=jfms
 139 xjfme=jfme
 140
 141
 142 ! init flags
 143 need_lfn_update=.false.
 144 ignitions_done=0
 145 freeze_fire = fire_const_time > 0. .and. time_start < fire_const_time
 146
 147 if(ifun.eq.1)then       ! do nothing, init pass 1 is outside only
 148 elseif(ifun.eq.2)then
 149         ! initialize all arrays that the model will not change later
 150
 151         ! assuming halo on zsf done
 152         ! extrapolate on 1 row of cells beyond the domain boundary
 153         ! including on the halo regions
 154
 155         call continue_at_boundary(1,1,0., & ! do x direction or y direction
 156             ifms,ifme,jfms,jfme,           &                ! memory dims
 157             ifds,ifde,jfds,jfde, &                     ! domain dims
 158             ifps,ifpe,jfps,jfpe, &            ! patch dims - winds defined up to +1
 159             ifts,ifte,jfts,jfte, &                ! tile dims
 160             itso,iteo,jtso,jteo, &              ! where set now
 161             fp%zsf)                               ! array
 162
 163 !       compute the gradients once for all
 164         err=0.
 165         maxgrad=0.
 166         do j=jfts,jfte
 167             do i=ifts,ifte
 168                 erri = fp%dzdxf(i,j) - (fp%zsf(i+1,j)-fp%zsf(i-1,j))/(2.*fdx)
 169                 errj = fp%dzdyf(i,j) - (fp%zsf(i,j+1)-fp%zsf(i,j-1))/(2.*fdy)
 170                 err=max(err,abs(erri),abs(errj))
 171                 grad=sqrt(fp%dzdxf(i,j)**2+fp%dzdyf(i,j)**2)
 172                 maxgrad=max(maxgrad,grad)
 173             enddo
 174         enddo
 175 !$OMP CRITICAL(SFIRE_MODEL_CRIT)
 176         write(msg,*)'max gradient ',maxgrad,' max error against zsf',err
 177 !$OMP END CRITICAL(SFIRE_MODEL_CRIT)
 178         call message(msg)
 179
 180         if(.not.restart)call set_nfuel_cat( &
 181             ifms,ifme,jfms,jfme, &
 182             ifts,ifte,jfts,jfte, &
 183             ifuelread,nfuel_cat0,&
 184             fp%zsf,nfuel_cat)            ! better not use the extrapolated zsf!!
 185
 186         ! uses nfuel_cat to set the other fuel data arrays
 187         ! needs zsf on halo width 1 to compute the terrain gradient
 188         if(.not.restart)call set_fire_params(   &
 189             ifds,ifde,jfds,jfde, &
 190             ifms,ifme,jfms,jfme, &
 191             ifts,ifte,jfts,jfte, &
 192             fdx,fdy,nfuel_cat0,  &
 193             nfuel_cat,fuel_time, &
 194             fp  &
 195 )
 196
 197         ! initialize model state to no fire
 198         if(.not.restart)then
 199             call init_no_fire  ( &
 200             ifds,ifde,jfds,jfde, &
 201             ifms,ifme,jfms,jfme, &
 202             ifts,ifte,jfts,jfte, &
 203             fdx,fdy,time_start,  &
 204             fuel_frac,fire_area,lfn,tign)
 205
 206             need_lfn_update=.true. ! because we have set lfn
 207
 208         endif
 209
 210 elseif(ifun.eq.3)then   ! ignition if so specified
 211
 212
 213 elseif (ifun.eq.4) then  ! do the timestep
 214
 215     if(fire_print_msg.ge.stat_lev)then
 216       aw=fun_real(RNRM_SUM,  &
 217         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 218         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 219         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 220         0,0,0,       &                            ! staggering
 221         fp%vx,fp%vy)/((ifde-ifds+1)*(jfde-jfds+1))
 222       mw=fun_real(RNRM_MAX,  &
 223         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 224         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 225         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 226         0,0,0,       &                            ! staggering
 227         fp%vx,fp%vy)
 228 !$OMP MASTER
 229       write(msg,91)time_start,'Average wind        ',aw,'m/s'
 230       call message(msg,stat_lev)
 231       write(msg,91)time_start,'Maximum wind        ',mw,'m/s'
 232       call message(msg,stat_lev)
 233 !$OMP END MASTER
 234     endif
 235
 236 !   compute fuel fraction at start
 237 !    call fuel_left( &
 238 !        ifms,ifme,jfms,jfme, &
 239 !        ifts,ifte,jfts,jfte, &
 240 !        ifms,ifme,jfms,jfme, &
 241 !        lfn,tign,fuel_time,time_start,fuel_frac,fire_area) ! fuel frac is shared
 242
 243     call print_2d_stats(ifts,ifte,jfts,jfte, &
 244                    ifms,ifme,jfms,jfme, &
 245                    fuel_frac,'model: fuel_frac start')
 246
 247     ! advance the model from time_start to time_start+dt
 248     ! return the fuel fraction burnt this call in each fire cell
 249     ! will call module_fr_sfire_speed::normal_spread for propagation speed
 250     ! We cannot simply compute the spread rate here because that will change with the
 251     ! angle of the wind and the direction of propagation, thus it is done in subroutine
 252     ! normal_spread at each fire time step. Instead, we pass arguments that
 253     ! the speed function may use as fp.
 254
 255 !   propagate level set function in time
 256 !   set lfn_out tign
 257 !   lfn does not change, tign has no halos
 258
 259     if(.not. freeze_fire)then
 260
 261     call prop_ls(id,     &
 262         ifds,ifde,jfds,jfde,                      & ! fire domain dims - the whole domain
 263         ifms,ifme,jfms,jfme,                      &
 264         ifps,ifpe,jfps,jfpe, &                ! patch - nodes owned by this process
 265         ifts,ifte,jfts,jfte,                      &
 266         time_start,dt,fdx,fdy,tbound,  &
 267         lfn,lfn_out,tign,ros, fp &
 268     )
 269     else
 270         call message('sfire_model: EXPERIMENTAL: skipping fireline propagation')
 271
 272     endif
 273
 274 elseif (ifun.eq.5) then ! copy the result of timestep back to input
 275     ! this cannot be done in the time step itself because of race condition
 276     ! some thread may still be using lfn as input in their tile halo
 277
 278     if(.not. freeze_fire)then
 279
 280     do j=jfts,jfte
 281         do i=ifts,ifte
 282             lfn(i,j)=lfn_out(i,j)
 283             ! if want to try timestep again treat tign the same way here
 284             ! even if tign does not need a halo
 285         enddo
 286     enddo
 287
 288     endif
 289
 290     ! check for ignitions
 291     do ig = 1,num_ignitions
 292
 293 !  for now, check for ignition every time step...
 294 !        if(ignition_line(ig)%end_time>=time_start.and.ignition_line(ig)%start_time<time_start+dt)then
 295             call ignite_fire(                             &
 296                 ifds,ifde,jfds,jfde,                      & ! fire domain dims - the whole domain
 297                 ifms,ifme,jfms,jfme,                      &
 298                 ifts,ifte,jfts,jfte,                      &
 299                 ignition_line(ig),                        &
 300                 time_start,time_start+dt,                 &
 301                 coord_xf,coord_yf,unit_xf,unit_yf,        &
 302                 lfn,tign,ignited)
 303
 304             ignitions_done=ignitions_done+1
 305             ignited_tile(ignitions_done)=ignited
 306
 307 !            need_lfn_update=.true. ! if ignition, lfn changed
 308 #ifdef DEBUG_OUT
 309             call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn,'lfn_ig',id)
 310             call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,coord_xf,'coord_xf_ig',id)
 311             call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,coord_yf,'coord_yf_ig',id)
 312 #endif
 313 !        endif
 314
 315     enddo
 316
 317     call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme, &
 318                    lfn,'sfire_model: lfn out')
 319
 320
 321     need_lfn_update=.true. ! duh
 322
 323 elseif (ifun.eq.6) then ! timestep postprocessing
 324
 325   if(.not. freeze_fire)then
 326
 327     ! compute the heat fluxes from the fuel burned
 328     ! needs lfn and tign from neighbors so halo must be updated before
 329     call fuel_left(&
 330         ifms,ifme,jfms,jfme, &
 331         ifts,ifte,jfts,jfte, &
 332         ifts,ifte,jfts,jfte, &
 333         lfn,tign,fuel_time,time_start+dt,fuel_frac_end,fire_area) !fuel_frac_end is private and tile based
 334
 335     call print_2d_stats(ifts,ifte,jfts,jfte, &
 336                    ifts,ifte,jfts,jfte, &
 337                    fuel_frac_end,'model: fuel_frac end')
 338
 339     do j=jfts,jfte
 340         do i=ifts,ifte
 341             fuel_frac_burnt(i,j)=fuel_frac(i,j)-fuel_frac_end(i,j) ! fuel lost this timestep
 342             fuel_frac(i,j)=fuel_frac_end(i,j) ! copy new value to state array
 343         enddo
 344     enddo
 345
 346     call print_2d_stats(ifts,ifte,jfts,jfte, &
 347                    ifts,ifte,jfts,jfte, &
 348                    fuel_frac_burnt,'model: fuel_frac burned')
 349
 350     call heat_fluxes(dt,                          &
 351         ifms,ifme,jfms,jfme,                      &
 352         ifts,ifte,jfts,jfte,                      &
 353         ifts,ifte,jfts,jfte,                      &  ! fuel_frac_burned is tile dimensioned
 354         fp%fgip,                                     &
 355         fuel_frac_burnt,                          & !
 356         grnhfx,grnqfx)                              !out
 357
 358     if(fire_print_msg.ge.stat_lev)then
 359       tfa=fun_real(REAL_SUM,  &
 360         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 361         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 362         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 363         0,0,0,       &                            ! staggering
 364         fire_area,fire_area) * fdx * fdy
 365       thf=fun_real(REAL_SUM,  &
 366         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 367         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 368         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 369         0,0,0,       &                            ! staggering
 370         grnhfx,grnhfx) * fdx * fdy
 371       mhf=fun_real(REAL_MAX,  &
 372         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 373         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 374         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 375         0,0,0,       &                            ! staggering
 376         grnhfx,grnhfx)
 377       tqf=fun_real(REAL_SUM,  &
 378         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 379         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 380         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 381         0,0,0,       &                            ! staggering
 382         grnqfx,grnqfx) * fdx * fdy
 383       mqf=fun_real(REAL_MAX,  &
 384         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 385         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 386         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 387         0,0,0,       &                            ! staggering
 388         grnqfx,grnqfx)
 389 !$OMP MASTER
 390       write(msg,91)time_start,'Fire area           ',tfa,'m^2'
 391       call message(msg,stat_lev)
 392       write(msg,91)time_start,'Heat output         ',thf,'W'
 393       call message(msg,stat_lev)
 394       write(msg,91)time_start,'Max heat flux       ',mhf,'W/m^2'
 395       call message(msg,stat_lev)
 396       write(msg,91)time_start,'Latent heat output  ',tqf,'W'
 397       call message(msg,stat_lev)
 398       write(msg,91)time_start,'Max latent heat flux',mqf,'W/m^2'
 399       call message(msg,stat_lev)
 400 !$OMP END MASTER
 401 91  format('Time ',f11.3,' s ',a,e12.3,1x,a)
 402     endif
 403
 404
 405   else
 406      call message('sfire_model: EXPERIMENTAL: skipping fuel burnt computation')
 407
 408      if (fire_const_grnhfx >= 0. .and. fire_const_grnqfx >= 0.) then
 409
 410 !$OMP CRITICAL(SFIRE_MODEL_CRIT)
 411         write(msg,'(a,2e12.3,a)')'sfire_model: EXPERIMENTAL output constant heat flux', &
 412            fire_const_grnhfx, fire_const_grnqfx, ' W/s'
 413 !$OMP END CRITICAL(SFIRE_MODEL_CRIT)
 414         call message(msg)
 415
 416         do j=jfts,jfte
 417             do i=ifts,ifte
 418                 grnhfx(i,j)=fire_const_grnhfx
 419                 grnqfx(i,j)=fire_const_grnqfx
 420             enddo
 421         enddo
 422
 423       endif
 424
 425    endif
 426
 427     call print_2d_stats(ifts,ifte,jfts,jfte, &
 428                    ifms,ifme,jfms,jfme, &
 429                    grnhfx,'model: heat flux(J/m^2/s)')
 430
 431 else
 432 !$OMP CRITICAL(SFIRE_MODEL_CRIT)
 433     write(msg,*)'sfire_model: bad ifun=',ifun
 434 !$OMP END CRITICAL(SFIRE_MODEL_CRIT)
 435     call crash(msg)
 436 endif
 437
 438 end subroutine sfire_model
 439
 440 !
 441 !*****************
 442 !
 443
 444 end module module_fr_sfire_model