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 implicit none
  11
  12 contains
  13
  14 subroutine sfire_model (                    &
  15     id,                                     & ! unique number for prints and debug
  16     ifun,                                   & ! what to do see below
  17     restart,replay,                         & ! use existing state, prescribe spread
  18     run_fuel_moisture,                      & ! if need update fuel moisture in pass 4
  19     ifuelread,nfuel_cat0,                   & ! initialize fuel categories
  20     ifds,ifde,jfds,jfde,                    & ! fire domain dims - the whole domain
  21     ifms,ifme,jfms,jfme,                    & ! fire memory dims - how declared
  22     ifps,ifpe,jfps,jfpe,                    & ! patch - nodes owned by this process
  23     ifts,ifte,jfts,jfte,                    & ! fire tile dims  - this thread
  24     time_start,dt,                          & ! time and increment
  25     fdx,fdy,                                & ! fire mesh spacing,
  26     ignition,hfx,                           & ! small array of ignition line descriptions
  27     coord_xf,coord_yf,                      & ! fire mesh coordinates
  28     fire_hfx,                               & ! input: given heat flux, or set inside
  29     tign_in,                                & ! ignition time, if given
  30     lfn,lfn_out,tign,fuel_frac,fire_area,   & ! state: level function, ign time, fuel left, area burning
  31     fuel_frac_burnt,                        & ! output: fuel fraction burnt in this timestep
  32     fgrnhfx,fgrnqfx,                          & ! output: heat fluxes
  33     ros,flineint,flineint2,                 & ! diagnostic variables
  34     f_ros0,f_rosx,f_rosy,f_ros,             & ! fire risk spread
  35     f_int,f_lineint,f_lineint2,             & ! fire risk intensities
  36     nfuel_cat,                              & ! fuel data per point
  37     fuel_time,fwh,fz0,                      & ! save derived internal data
  38     fp &
  39 )
  40
  41 ! This subroutine implements the fire spread model.
  42 ! All quantities are on the fire grid. It inputs
  43 ! winds given on the nodes of the fire grid
  44 ! and outputs the heat fluxes on the cells of the fire grid.
  45 ! This subroutine has no knowledge of any atmospheric model.
  46 ! This code was written to conform with the WRF parallelism model, however it
  47 ! does not depend on it. It can be called with domain equal to tile.
  48 ! Wind and height must be given on 1 more node beyond the domain bounds.
  49 ! The subroutine changes only array entries of the arguments in the tile.
  50 ! Upon exit with ifun=2 (time step), lfn_out is to be copied into lfn by the caller.
  51 ! When this subroutine is used on separate tiles that make a domain the value, the
  52 ! it uses lfn on a strip of width 2 from neighboring tiles.
  53 !
  54 ! All computation is done on one tile.
  55 !
  56 ! This subroutine is intended to be called in a loop like
  57 !
  58 !
  59 ! do ifun=1,6 (if initizalize run, otherwise 3,6)
  60 !   start parallel loop over tiles
  61 !       if ifun=1, set z and fuel data
  62 !       if ifun=3, set the wind arrays
  63 !       call sfire_model(....)
  64 !   end parallel loop over tiles
  65 !
  66 !
  67 !   if ifun=0
  68 !       halo exchange on z width 2
  69 !       halo exchange on fuel data width 1
  70 !   endif
  71 !
  72 !   if ifun=3, halo exchange on winds width 2
  73 !
  74 ! enddo
  75
  76 implicit none
  77
  78 !*** arguments
  79
  80 ! control switches
  81 integer, intent(in) :: id
  82 integer, intent(in) :: ifun                 ! 1 = initialize run pass 1
  83                                             ! 2 = initialize run pass 2
  84                                             ! 3 = initialize timestep
  85                                             ! 4 = do one timestep
  86                                             ! 5 = copy timestep output to input
  87                                             ! 6 = compute output fluxes
  88 logical, intent(in):: restart               ! if true, use existing state
  89 logical, intent(in):: replay                ! if true, use tign_g for level set
  90 logical, intent(in)::run_fuel_moisture      !
  91 ! scalar data
  92 integer, intent(in) :: ifuelread,nfuel_cat0 ! for set_fire_params
  93 integer, intent(in) :: ifds,ifde,jfds,jfde,&  ! fire domain bounds
  94         ifps,ifpe,jfps,jfpe                ! patch - nodes owned by this process
  95 integer, intent(in) :: ifts,ifte,jfts,jfte  ! fire tile bounds
  96 integer, intent(in) :: ifms,ifme,jfms,jfme  ! fire memory array bounds
  97 REAL,INTENT(in) :: time_start,dt            ! starting time, time step
  98 REAL,INTENT(in) :: fdx,fdy                  ! spacing of the fire mesh
  99 ! array data
 100 type(lines_type), intent(in):: ignition,hfx  ! descriptions of ignition lines and hfx lines
 101 real, dimension(ifms:ifme, jfms:jfme), intent(in):: &
 102     coord_xf,coord_yf                      ! node coordinates
 103 real, dimension(ifms:ifme, jfms:jfme), intent(inout):: &
 104     fire_hfx                                ! given heat flux
 105
 106 real, dimension(ifms:ifme, jfms:jfme), intent(inout):: tign_in ! given ignition times
 107 ! state
 108 REAL, INTENT(inout), dimension(ifms:ifme,jfms:jfme):: &
 109     lfn   , &                               ! level function: fire is where lfn<0 (node)
 110     tign  , &                               ! absolute time of ignition (node)
 111     fuel_frac                               ! fuel fraction (node), currently redundant
 112
 113 REAL, INTENT(out), dimension(ifms:ifme,jfms:jfme):: &
 114     fire_area, &                            ! fraction of each cell burning
 115     fuel_frac_burnt                         ! fuel fraction burned in this timestep
 116
 117 ! output
 118 REAL, INTENT(out), dimension(ifms:ifme,jfms:jfme):: &
 119     lfn_out, &                              !
 120     fgrnhfx,fgrnqfx, &                        ! heat fluxes J/m^2/s  (cell)
 121     ros,flineint,flineint2,                 & ! diagnostic variables
 122     f_ros0,f_rosx,f_rosy,f_ros,             & ! fire risk spread
 123     f_int,f_lineint,f_lineint2                ! fire risk intensities
 124
 125
 126 ! constant arrays - set at initialization
 127 real, intent(inout), dimension(ifms:ifme, jfms:jfme)::nfuel_cat ! cell based, data, constant
 128 real,intent(inout),dimension(ifms:ifme,jfms:jfme):: fuel_time,fwh,fz0
 129 type(fire_params),intent(inout)::fp
 130
 131 !*** local
 132
 133 integer :: xifms,xifme,xjfms,xjfme  ! memory bounds for pass-through arguments to normal spread
 134 real, dimension(ifts:ifte,jfts:jfte)::fuel_frac_end
 135 integer::ignited,ig,i,j,itso,iteo,jtso,jteo
 136 real::tbound,err,erri,errj,maxgrad,grad,tfa,thf,mhf,tqf,mqf,aw,mw,t
 137 character(len=128)::msg
 138 logical:: freeze_fire
 139 real::fireline_mask=0.
 140
 141 !*** executable
 142
 143 call check_mesh_2dim(ifts-1,ifte+1,jfts-1,jfte+1,ifms,ifme,jfms,jfme)
 144
 145 xifms=ifms  ! dimensions for the include file
 146 xifme=ifme
 147 xjfms=jfms
 148 xjfme=jfme
 149
 150
 151 ! init flags
 152 freeze_fire = fire_hfx_given .ne.  0
 153
 154
 155 if(ifun.eq.1)then       ! do nothing, init pass 1 is outside only
 156 ! !$OMP SINGLE
 157 !! done in driver now
 158 !   call init_fuel_cats  ! initialize fuel subsystem
 159 ! !$OMP END SINGLE
 160         if(replay) then
 161         ! when starting replay, recompute lfn and fuel_frac
 162           call message('replay, setting the level set function',level=1)
 163           do j=jfts,jfte
 164             do i=ifts,ifte
 165               lfn(i,j) = tign(i,j) - time_start ! <0 if burning at the end of the time step
 166             enddo
 167           enddo
 168           call check_lfn_tign('replay, lfn from tign',time_start,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn,tign)
 169         endif
 170 elseif(ifun.eq.2)then
 171         ! initialize all arrays that the model will not change later
 172
 173         ! assuming halo on zsf done
 174         ! extrapolate on 1 row of cells beyond the domain boundary
 175         ! including on the halo regions
 176
 177         call continue_at_boundary(1,1,0., & ! do x direction or y direction
 178             ifms,ifme,jfms,jfme,           &                ! memory dims
 179             ifds,ifde,jfds,jfde, &                     ! domain dims
 180             ifps,ifpe,jfps,jfpe, &            ! patch dims - winds defined up to +1
 181             ifts,ifte,jfts,jfte, &                ! tile dims
 182             itso,iteo,jtso,jteo, &              ! where set now
 183             fp%zsf)                               ! array
 184
 185 !       compute the gradients once for all
 186         err=0.
 187         maxgrad=0.
 188         do j=jfts,jfte
 189             do i=ifts,ifte
 190                 erri = fp%dzdxf(i,j) - (fp%zsf(i+1,j)-fp%zsf(i-1,j))/(2.*fdx)
 191                 errj = fp%dzdyf(i,j) - (fp%zsf(i,j+1)-fp%zsf(i,j-1))/(2.*fdy)
 192                 err=max(err,abs(erri),abs(errj))
 193                 grad=sqrt(fp%dzdxf(i,j)**2+fp%dzdyf(i,j)**2)
 194                 maxgrad=max(maxgrad,grad)
 195             enddo
 196         enddo
 197 !$OMP CRITICAL(SFIRE_MODEL_CRIT)
 198         write(msg,*)'max gradient ',maxgrad,' max error against zsf',err
 199 !$OMP END CRITICAL(SFIRE_MODEL_CRIT)
 200         call message(msg)
 201
 202         call set_nfuel_cat(      &       ! also on restart
 203             ifms,ifme,jfms,jfme, &
 204             ifts,ifte,jfts,jfte, &
 205             ifuelread,nfuel_cat0,&
 206             fp%zsf,nfuel_cat)            ! better not use the extrapolated zsf!!
 207
 208         ! uses nfuel_cat to set the other fuel data arrays
 209         ! needs zsf on halo width 1 to compute the terrain gradient
 210         call set_fire_params(    &       ! also on restart
 211             ifds,ifde,jfds,jfde, &
 212             ifms,ifme,jfms,jfme, &
 213             ifts,ifte,jfts,jfte, &
 214             fdx,fdy,nfuel_cat0,  &
 215             nfuel_cat,fuel_time, &
 216             fp)
 217
 218         ! initialize model state to no fire / set tign_in
 219         if(.not.restart)then
 220             call init_no_fire  ( &
 221             ifds,ifde,jfds,jfde, &
 222             ifms,ifme,jfms,jfme, &
 223             ifts,ifte,jfts,jfte, &
 224             fdx,fdy,time_start,dt,  &
 225             fuel_frac,fire_area,lfn,tign_in,tign)
 226
 227         endif
 228
 229         call check_lfn_tign('after initialization:',time_start,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn,tign)
 230
 231         if(replay) then
 232           call message('replay, recomputing fuel fraction')
 233           call fuel_left( &
 234             ifds,ifde,jfds,jfde, &
 235             ifms,ifme,jfms,jfme, &
 236             ifts,ifte,jfts,jfte, &
 237             ifms,ifme,jfms,jfme, &
 238             lfn,tign,fuel_time,time_start,fuel_frac,fire_area) !fuel_frac is global
 239           call check_lfn_tign('recomputed fuel fraction',time_start,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn,tign)
 240         endif
 241
 242 elseif(ifun.eq.3)then   ! ignition if so specified
 243
 244
 245 elseif (ifun.eq.4) then  ! do the timestep
 246
 247     ! #35 call check_lfn_tign('time step start',time_start,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn, tign)
 248     if(run_fuel_moisture)then
 249         ! uses nfuel_cat to set the other fuel data arrays
 250         ! needs zsf on halo width 1 to compute the terrain gradient
 251         call set_fire_params(    &       ! also on restart
 252             ifds,ifde,jfds,jfde, &
 253             ifms,ifme,jfms,jfme, &
 254             ifts,ifte,jfts,jfte, &
 255             fdx,fdy,nfuel_cat0,  &
 256             nfuel_cat,fuel_time, &
 257             fp)
 258     endif
 259
 260     if(fire_print_msg.ge.stat_lev)then
 261       aw=fun_real(RNRM_SUM,  &
 262         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 263         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 264         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 265         0,0,0,       &                            ! staggering
 266         fp%vx,fp%vy)/((ifde-ifds+1)*(jfde-jfds+1))
 267       mw=fun_real(RNRM_MAX,  &
 268         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 269         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 270         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 271         0,0,0,       &                            ! staggering
 272         fp%vx,fp%vy)
 273 !$OMP MASTER
 274       write(msg,91)time_start,'Average surface wind',aw,'m/s'
 275       call message(msg,stat_lev)
 276       write(msg,91)time_start,'Maximum surface wind',mw,'m/s'
 277       call message(msg,stat_lev)
 278 !$OMP END MASTER
 279     endif
 280
 281 !   compute fuel fraction at start
 282 !    call fuel_left( &
 283 !        ifms,ifme,jfms,jfme, &
 284 !        ifts,ifte,jfts,jfte, &
 285 !        ifms,ifme,jfms,jfme, &
 286 !        lfn,tign,fuel_time,time_start,fuel_frac,fire_area) ! fuel frac is shared
 287
 288     call print_2d_stats(ifts,ifte,jfts,jfte, &
 289                    ifms,ifme,jfms,jfme, &
 290                    fuel_frac,'model: fuel_frac start')
 291
 292     ! advance the model from time_start to time_start+dt
 293     ! return the fuel fraction burnt this call in each fire cell
 294     ! will call module_fr_sfire_speed::normal_spread for propagation speed
 295     ! We cannot simply compute the spread rate here because that will change with the
 296     ! angle of the wind and the direction of propagation, thus it is done in subroutine
 297     ! normal_spread at each fire time step. Instead, we pass arguments that
 298     ! the speed function may use as fp.
 299
 300 !   propagate level set function in time
 301 !   set lfn_out tign
 302 !   lfn does not change, tign has no halos
 303
 304     if(.not. freeze_fire)then
 305       if(.not.replay) then
 306         call prop_ls(id,1,     &
 307         ifds,ifde,jfds,jfde,                      & ! fire domain dims - the whole domain
 308         ifms,ifme,jfms,jfme,                      &
 309         ifps,ifpe,jfps,jfpe, &                ! patch - nodes owned by this process
 310         ifts,ifte,jfts,jfte,                      &
 311         time_start,dt,fdx,fdy,tbound,  &
 312         lfn,lfn_out,tign,ros, fp &
 313         )
 314       else
 315         do j=jfts,jfte
 316           do i=ifts,ifte
 317             lfn_out(i,j) = tign(i,j) - (time_start + dt) ! <0 if burning at the end of the time step
 318           enddo
 319         enddo
 320       endif
 321     else
 322         call message('sfire_model: EXPERIMENTAL: skipping fireline propagation')
 323
 324     endif
 325
 326     call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme, ros,'model: ros')
 327
 328 elseif (ifun.eq.5) then ! copy the result of timestep back to input
 329     ! this cannot be done in the time step itself because of race condition
 330     ! some thread may still be using lfn as input in their tile halo
 331
 332     ! fix lfn_out if needed and compute tign
 333     if(.not. freeze_fire)then
 334       if(.not.replay) then
 335         call prop_ls(id,2,     &
 336         ifds,ifde,jfds,jfde,                      & ! fire domain dims - the whole domain
 337         ifms,ifme,jfms,jfme,                      &
 338         ifps,ifpe,jfps,jfpe, &                ! patch - nodes owned by this process
 339         ifts,ifte,jfts,jfte,                      &
 340         time_start,dt,fdx,fdy,tbound,  &
 341         lfn,lfn_out,tign,ros, fp &
 342         )
 343       endif
 344     endif
 345
 346     call check_lfn_tign('time step end',time_start+dt,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn_out, tign)
 347
 348     if(.not. freeze_fire)then
 349     do j=jfts,jfte
 350         do i=ifts,ifte
 351             lfn(i,j)=lfn_out(i,j)
 352             ! if want to try timestep again treat tign the same way here
 353             ! even if tign does not need a halo
 354         enddo
 355     enddo
 356
 357     call check_lfn_tign('before ignition',time_start+dt,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn, tign)
 358
 359     ! check for ignitions
 360     if( fire_tign_in_time > 0.) then
 361          if (ignition%num_lines >  0) then
 362               call crash('ignition from lines and from tign_in are not compatible')
 363          endif
 364          call ignite_from_tign_in(                                   &
 365                         ifds,ifde,jfds,jfde,                      & ! fire domain dims - the whole domain
 366                         ifms,ifme,jfms,jfme,                      &
 367                         ifts,ifte,jfts,jfte,                      &
 368                         time_start,time_start+dt,                 &
 369                         tign_in,                                  &
 370                         lfn,tign,ignited)
 371          call check_lfn_tign('after ignite_from_tign_in',time_start+dt,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn, tign)
 372     endif
 373
 374
 375     do ig = 1,ignition%num_lines
 376
 377 !  for now, check for ignition every time step...
 378 !        if(ignition%line(ig)%end_time>=time_start.and.ignition%line(ig)%start_time<time_start+dt)then
 379             call ignite_fire(                             &
 380                 ifds,ifde,jfds,jfde,                      & ! fire domain dims - the whole domain
 381                 ifms,ifme,jfms,jfme,                      &
 382                 ifts,ifte,jfts,jfte,                      &
 383                 ignition%line(ig),                        &
 384                 time_start,time_start+dt,                 &
 385                 coord_xf,coord_yf,ignition%unit_fxlong,ignition%unit_fxlat,        &
 386                 lfn,tign,ignited)
 387
 388 #ifdef DEBUG_OUT
 389             call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn,'lfn_ig',id)
 390             call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,coord_xf,'coord_xf_ig',id)
 391             call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,coord_yf,'coord_yf_ig',id)
 392 #endif
 393 !        endif
 394
 395     enddo
 396     else
 397         call message('sfire_model: EXPERIMENTAL: skipping ignition')
 398     endif
 399
 400     call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme, &
 401                    lfn,'sfire_model: lfn out')
 402
 403     call check_lfn_tign('after ignition',time_start+dt,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn, tign)
 404
 405 elseif (ifun.eq.6) then ! timestep postprocessing
 406     ! have halo on lfn now
 407
 408     ! diagnostics
 409     call fire_intensity(fp,                   &  ! fuel properties
 410         ifms,ifme,jfms,jfme,                      &  ! memory dims
 411         ifts,ifte,jfts,jfte,                      &  ! tile dims
 412         ifms,ifme,jfms,jfme,                      &  ! ros dims
 413         ros,nfuel_cat,                            & !in
 414         flineint,flineint2)                       ! fireline intensities out
 415
 416     if(fireline_mask < 0.) then
 417       do j=jfts,jfte
 418         do i=ifts,ifte
 419           ! if the sign of lfn is the same as all of its neighbors or we are at domain boundary, we are not near the fireline
 420           if( (lfn(i-1,j-1)>0. .and. lfn(i-1,j)>0. .and. lfn(i,j-1)>0. .and. lfn(i,j)>0.  .and. &
 421             lfn(i+1,j+1)>0. .and. lfn(i+1,j)>0. .and. lfn(i,j+1)>0. ) .or.                   &
 422            (lfn(i-1,j-1)<0. .and. lfn(i-1,j)<0. .and. lfn(i,j-1)<0. .and. lfn(i,j)<0.  .and. &
 423             lfn(i+1,j+1)<0. .and. lfn(i+1,j)<0. .and. lfn(i,j+1)<0. ) .or.                   &
 424             i.eq.ifds .or. i .eq. ifde .or. j.eq.jfds .or. j.eq.jfde) then
 425                ros(i,j)=fireline_mask
 426                flineint(i,j)=fireline_mask
 427                flineint2(i,j)=fireline_mask
 428            endif
 429         enddo
 430       enddo
 431     endif
 432
 433   call fire_risk(fp,                              &
 434         ifms,ifme,jfms,jfme,                      &  ! memory dims
 435         ifts,ifte,jfts,jfte,                      &  ! tile dims
 436         nfuel_cat,                                &  !
 437         f_ros0,f_rosx,f_rosy,f_ros,               &  ! fire spread
 438         f_int,f_lineint,f_lineint2)                  ! fire intensities for danger rating
 439
 440
 441   select case(fire_hfx_given)
 442
 443     case(0)   ! normal
 444
 445     ! compute the heat fluxes from the fuel burned
 446     ! needs lfn and tign from neighbors so halo must be updated before
 447 !$OMP CRITICAL(SFIRE_MODEL_CRIT)
 448     write(msg,*)'time_start=',time_start,' dt=',dt,' before fuel_left'
 449 !$OMP END CRITICAL(SFIRE_MODEL_CRIT)
 450     call message(msg)
 451     call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn,'model: lfn')
 452     call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,tign,'model: tign')
 453     call check_lfn_tign('before fuel_left',time_start+dt,ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn, tign)
 454
 455     call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fuel_time,'model: fuel_time')
 456
 457     if(fire_update_fuel_frac.eq.1)then
 458         ! run model the normal way
 459         call fuel_left(&
 460             ifds,ifde,jfds,jfde, &
 461             ifms,ifme,jfms,jfme, &
 462             ifts,ifte,jfts,jfte, &
 463             ifts,ifte,jfts,jfte, &
 464             lfn,tign,fuel_time,time_start+dt,fuel_frac_end,fire_area) !fuel_frac_end is private and tile based
 465
 466         call print_2d_stats(ifts,ifte,jfts,jfte,ifts,ifte,jfts,jfte,fuel_frac_end,'model: fuel_frac end')
 467
 468         do j=jfts,jfte
 469             do i=ifts,ifte
 470                 t = min(fuel_frac(i,j),fuel_frac_end(i,j)) ! do not allow fuel fraction to increase, in case of approximation error
 471                 fuel_frac_burnt(i,j)=fuel_frac(i,j)-t ! fuel lost this timestep
 472                 fuel_frac(i,j)=t ! copy new value to state array
 473             enddo
 474         enddo
 475     elseif(fire_update_fuel_frac.eq.2)then
 476         do j=jfts,jfte
 477             do i=ifts,ifte
 478                 if(lfn(i,j)<0.)then
 479                     fuel_frac_burnt(i,j) = dt / fuel_time(i,j)
 480                 else
 481                     fuel_frac_burnt(i,j) = 0.
 482                 endif
 483             enddo
 484         enddo
 485     else
 486         call crash('fire_update_fuel_frac value not supported')
 487     endif
 488
 489     call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,fuel_frac_burnt,'model: fuel_frac burned')
 490
 491     call heat_fluxes(dt,fp,                       &
 492         ifms,ifme,jfms,jfme,                      &
 493         ifts,ifte,jfts,jfte,                      &
 494         ifms,ifme,jfms,jfme,                      &  ! fuel_frac_burned has standard memory dimensions
 495         fp%fgip,                                     &
 496         fuel_frac_burnt,                          & !
 497         fgrnhfx,fgrnqfx)                              !out
 498
 499     case (1, 2)
 500 !$OMP CRITICAL(SFIRE_MODEL_CRIT)
 501        write(msg,*)"model: expecting fire_hfx to be set in WRF, from wrfinput or wrfrst files"
 502        call message(msg)
 503 !$OMP END CRITICAL(SFIRE_MODEL_CRIT)
 504
 505        do j=jfts,jfte
 506            do i=ifts,ifte
 507               fgrnhfx(i,j) = (1. - fire_hfx_latent_part)*fire_hfx(i,j)
 508               fgrnqfx(i,j) =       fire_hfx_latent_part *fire_hfx(i,j)
 509            enddo
 510        enddo
 511
 512     case (3)
 513
 514        call message('artificial heat flux from parameters given in namelist.input')
 515
 516        call param_hfx( time_start, &
 517                             ifms,ifme,jfms,jfme, &
 518                             ifts,ifte,jfts,jfte, &
 519                             coord_xf,coord_yf,   &
 520                             hfx,                 &
 521                             fire_area,fgrnhfx,fgrnqfx)
 522
 523    case default
 524         call crash('bad fire_hfx_given')
 525    end select
 526
 527    ! this should run in any case
 528
 529     if(fire_print_msg.ge.stat_lev)then
 530       tfa=fun_real(REAL_SUM,  &
 531         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 532         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 533         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 534         0,0,0,       &                            ! staggering
 535         fire_area,fire_area) * fdx * fdy
 536       thf=fun_real(REAL_SUM,  &
 537         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 538         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 539         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 540         0,0,0,       &                            ! staggering
 541         fgrnhfx,fgrnhfx) * fdx * fdy
 542       mhf=fun_real(REAL_MAX,  &
 543         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 544         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 545         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 546         0,0,0,       &                            ! staggering
 547         fgrnhfx,fgrnhfx)
 548       tqf=fun_real(REAL_SUM,  &
 549         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 550         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 551         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 552         0,0,0,       &                            ! staggering
 553         fgrnqfx,fgrnqfx) * fdx * fdy
 554       mqf=fun_real(REAL_MAX,  &
 555         ifms,ifme,1,1,jfms,jfme, &                ! memory dims
 556         ifds,ifde,1,1,jfds,jfde, &                ! domain dims
 557         ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
 558         0,0,0,       &                            ! staggering
 559         fgrnqfx,fgrnqfx)
 560 !$OMP MASTER
 561       write(msg,91)time_start,'Fire area           ',tfa,'m^2'
 562       call message(msg,stat_lev)
 563       write(msg,91)time_start,'Heat output         ',thf,'W'
 564       call message(msg,stat_lev)
 565       write(msg,91)time_start,'Max heat flux       ',mhf,'W/m^2'
 566       call message(msg,stat_lev)
 567       write(msg,91)time_start,'Latent heat output  ',tqf,'W'
 568       call message(msg,stat_lev)
 569       write(msg,91)time_start,'Max latent heat flux',mqf,'W/m^2'
 570       call message(msg,stat_lev)
 571 !$OMP END MASTER
 572 91  format('Time ',f11.3,' s ',a,e12.3,1x,a)
 573     endif
 574
 575
 576    call print_2d_stats(ifts,ifte,jfts,jfte, &
 577                    ifms,ifme,jfms,jfme, &
 578                    fgrnhfx,'model: heat flux(J/m^2/s)')
 579
 580 else
 581 !$OMP CRITICAL(SFIRE_MODEL_CRIT)
 582     write(msg,*)'sfire_model: bad ifun=',ifun
 583 !$OMP END CRITICAL(SFIRE_MODEL_CRIT)
 584     call crash(msg)
 585 endif
 586
 587 end subroutine sfire_model
 588
 589        subroutine param_hfx( time_now,&
 590                             ifms,ifme,jfms,jfme, &
 591                             ifts,ifte,jfts,jfte, &
 592                             coord_xf,coord_yf,   &
 593                             hfx,                 &
 594                             fire_area,fgrnhfx,fgrnqfx)
 595 !***   generate artifical heat flux from a parametric description
 596 !***   arguments
 597        real, intent(in)::time_now
 598        integer, intent(in):: &
 599                             ifms,ifme,jfms,jfme, &
 600                             ifts,ifte,jfts,jfte
 601        type(lines_type), intent(in)::hfx
 602        real, dimension(ifms:ifme,jfms:jfme), intent(in)::coord_xf,coord_yf ! nodal coordinates
 603        real, dimension(ifms:ifme,jfms:jfme), intent(out)::fire_area,fgrnhfx,fgrnqfx ! the imposed heat flux
 604        character(len=128):: msg
 605 !***   local
 606        integer::i,j,k,nfa,ncells
 607        real:: d2,ax,ay,hfrac,fa,thfx,t,r,radius
 608        real, parameter:: sigma_mult=3.   ! 3 gives drop to 1% in trans. time from gaussian kernel
 609        real:: maxspeed=100  ! max speed how the circle can move
 610
 611        do j=jfts,jfte  ! zero out the outputs
 612            do i=ifts,ifte
 613               fire_area(i,j)=0
 614               fgrnhfx(i,j) = 0.
 615               fgrnqfx(i,j) = 0.
 616            enddo
 617        enddo
 618
 619        do k=1,hfx%num_lines
 620           if(hfx%line(k)%radius > 0.)then
 621              ! processing heatflux line i
 622              ! find the time multiplier
 623              t = max(hfx%line(k)%start_time - time_now, time_now - hfx%line(k)%end_time)
 624              if(t > 0.)then  ! postitive distance from the time interval
 625                  r = t / hfx%line(k)%trans_time  ! position in the transition - 1 is at transition distance
 626                  hfrac = exp(-(sigma_mult * r)**2/2.)  ! gaussian kernel
 627              else
 628                  hfrac = 1.0
 629              endif
 630
 631
 632              ! find the coordinates of the center of the heat flux circle now
 633              ax = safe_prop(time_now, &
 634                               hfx%line(k)%start_time,&
 635                               hfx%line(k)%end_time,&
 636                               hfx%line(k)%start_x,&
 637                               hfx%line(k)%end_x, &
 638                               hfx%unit_fxlong*maxspeed)
 639              ay = safe_prop(time_now,&
 640                               hfx%line(k)%start_time,&
 641                               hfx%line(k)%end_time,&
 642                               hfx%line(k)%start_y,&
 643                               hfx%line(k)%end_y, &
 644                               hfx%unit_fxlat*maxspeed)
 645
 646              radius=hfx%line(k)%radius
 647
 648 !$OMP CRITICAL(SFIRE_MODEL_CRIT)
 649              write(msg,*)'hfx line ',i,' at ',time_now,'s ',hfrac,' of max ', hfx%line(k)%hfx_value
 650              call message(msg)
 651              write(msg,*)'center ',ax,ay,' radius ',radius
 652              call message(msg)
 653 !$OMP END CRITICAL(SFIRE_MODEL_CRIT)
 654
 655              nfa=0
 656              ncells=0
 657              do j=jfts,jfte
 658                  do i=ifts,ifte
 659                      ! distance squared from the center
 660                      d2 = (hfx%unit_fxlong*(ax - coord_xf(i,j)))**2 +  (hfx%unit_fxlat*(ay - coord_yf(i,j)))**2
 661                      if(d2 < radius**2)then
 662                          fa=1.
 663                      else
 664                          fa=0.
 665                      endif
 666                      ! set heat fluxes
 667                      thfx= hfx%line(k)%hfx_value * hfrac * fa  ! total heat flux at this point
 668                      fgrnhfx(i,j)= fgrnhfx(i,j) + (1.-fire_hfx_latent_part) * thfx
 669                      fgrnqfx(i,j)= fgrnqfx(i,j) + fire_hfx_latent_part * thfx
 670                      ! set fire area
 671                      fire_area(i,j) = max(fire_area(i,j),fa)
 672                      ! set stats
 673                      nfa=nfa+fa;
 674                      ncells=ncells+1
 675                  enddo
 676              enddo
 677
 678 !$OMP CRITICAL(SFIRE_MODEL_CRIT)
 679              write(msg,*)'Number of cells in fire area in this tile ',nfa,' ',(100.*nfa)/ncells,' %'
 680              call message(msg)
 681 !$OMP END CRITICAL(SFIRE_MODEL_CRIT)
 682
 683          endif
 684       enddo
 685  end subroutine param_hfx
 686
 687
 688
 689 real function safe_prop(t,t1,t2,x1,x2,ms)
 690 !  return x between x1 and x2 in the same proportion as is t between t1 and t2, safe in the case when t1=t2 and x1=x2
 691 !  safe_prop = x1 + (t-t1)*(x2-x1)/(t2-t1)
 692 !  future: abs((x2-x1)/(t2-t1)) capped at ms but still return x1 when t=t1 and x2 when t=t2
 693 real, intent(in)::t,t1,t2,x1,x2,ms
 694 real:: p,x
 695        if(t2 < t1)call crash('safe_prop: must have t2>t1')
 696        if(.not.ms>0.)call crash('safe_prop: must have ms>0')
 697        if(t1 .eq. t2)then
 698            if(x1.eq.x2)then
 699                x=x1
 700            else
 701                call crash('safe_prop: infinite speed')
 702            endif
 703        else
 704            p = (t - t1)/(t2 - t1)   ! 0 at t1, 1 at t2
 705            x = x1*(1.-p) + x2*p
 706        endif
 707        safe_prop=x
 708 end function safe_prop
 709 !
 710 !*****************
 711 !
 712
 713 end module module_fr_sfire_model