wrfv2_fire/phys/module_fr_sfire_core.F

   1 !
   2 #define DEBUG_OUT
   3 #define DEBUG_PRINT
   4 !#define FUEL_LEFT
   5 !#define DEBUG_OUT_FUEL_LEFT
   6
   7 module module_fr_sfire_core
   8
   9 use module_fr_sfire_phys, only: fire_params , fire_ros
  10 use module_fr_sfire_util
  11
  12 ! The mathematical core of the fire spread model. No physical constants here.
  13 !
  14 ! subroutine sfire_core: only this routine should be called from the outside.
  15 ! subroutine fuel_left:  compute remaining fuel from time of ignition.
  16 ! subroutine prop_ls: propagation of curve in normal direction.
  17
  18 ! describe one ignition line
  19 type ignition_line_type
  20   REAL  ros, &          ! subscale rate of spread during the ignition process
  21         stop_time, &    ! when the ignition process stops from ignition start (s)
  22         wind_red,  &    ! wind reduction factor at the ignition line
  23         wrdist,   &     ! distance from the ignition line when the wind reduction stops
  24         wrupwind, &     ! use distance interpolated between 0. = nearest 1. = upwind
  25         start_x, &      ! x coordinate of the ignition line start point (m, or long/lat)
  26         start_y, &      ! y coordinate of the ignition line start point
  27         end_x, &        ! x coordinate of the ignition line end point
  28         end_y, &        ! y coordinate of the ignition line end point
  29         start_time, &   ! ignition time for the start point from simulation start (s)
  30         end_time, &     ! ignition time for the end poin from simulation start (s)
  31         radius          ! all within this radius ignites immediately
  32 end type ignition_line_type
  33
  34 contains
  35
  36 !
  37 !****************************************
  38 !
  39
  40 subroutine init_no_fire(&
  41     ifds,ifde,jfds,jfde, &
  42     ifms,ifme,jfms,jfme, &
  43     ifts,ifte,jfts,jfte, &
  44     fdx,fdy,time_now,    & ! scalars in
  45     fuel_frac,fire_area,lfn,tign)    ! arrays out
  46 implicit none
  47
  48 !*** purpose: initialize model to no fire
  49
  50 !*** arguments
  51 integer, intent(in):: ifds,ifde,jfds,jfde   ! fire domain bounds
  52 integer, intent(in):: ifts,ifte,jfts,jfte   ! fire tile bounds
  53 integer, intent(in):: ifms,ifme,jfms,jfme   ! array bounds
  54 real, intent(in) :: fdx,fdy,time_now        ! mesh spacing, time
  55 real, intent(out), dimension (ifms:ifme,jfms:jfme) :: &
  56                    fuel_frac,fire_area,lfn,tign       ! model state
  57
  58 !*** calls
  59 intrinsic epsilon
  60
  61 !*** local
  62 integer:: i,j
  63 real lfn_init,time_init
  64
  65 lfn_init = 2*max((ifde-ifds+1)*fdx,(jfde-jfds+1)*fdy)      ! more than domain diameter
  66 time_init=time_now + max(time_now,1.0)*epsilon(time_now) ! a bit in future
  67
  68 do j=jfts,jfte
  69     do i=ifts,ifte
  70         fuel_frac(i,j)=1.          ! fuel at start is 1 by definition
  71         fire_area(i,j)=0.          ! nothing burning
  72         tign(i,j) = time_init      ! ignition in future
  73         lfn(i,j) = lfn_init        ! no fire
  74     enddo
  75 enddo
  76 call message('init_no_fire: state set to no fire')
  77
  78 end subroutine init_no_fire
  79
  80 !
  81 !******************
  82 !
  83
  84
  85 subroutine ignite_fire( ifds,ifde,jfds,jfde,                    & ! fire domain dims - the whole domain
  86                         ifms,ifme,jfms,jfme,                      &
  87                         ifts,ifte,jfts,jfte,                      &
  88                         ignition_line,                            &
  89                         start_ts,end_ts,                    &
  90                         coord_xf,coord_yf,                &
  91                         unit_xf,unit_yf,                  &
  92                         lfn,tign,ignited)
  93 implicit none
  94
  95 !*** purpose: ignite a circular/line fire
  96
  97 !*** description
  98 ! ignite fire in the region within radius r from the line (sx,sy) to (ex,ey).
  99 ! the coordinates of nodes are given as the arrays coord_xf and coord_yf
 100 ! r is given in m
 101 ! one unit of coord_xf is unit_xf m
 102 ! one unit of coord_yf is unit_yf m
 103 ! so a node (i,j) will be ignited iff for some (x,y) on the line
 104 ! || ( (coord_xf(i,j) - x)*unit_xf , (coord_yf(i,j) - y)*unit_yf ) || <= r
 105
 106
 107 !*** arguments
 108 integer, intent(in):: ifds,ifde,jfds,jfde   ! fire domain bounds
 109 integer, intent(in):: ifts,ifte,jfts,jfte   ! fire tile bounds
 110 integer, intent(in):: ifms,ifme,jfms,jfme   ! array bounds
 111 type(ignition_line_type), intent(in):: ignition_line    ! description of the ignition line
 112 real, intent(in):: start_ts,end_ts          ! the time step start and end
 113 real, dimension(ifms:ifme, jfms:jfme), intent(in):: &
 114     coord_xf,coord_yf                       !  node coordinates
 115 real, intent(in):: unit_xf,unit_yf          !  coordinate units in m
 116 real, intent(inout), dimension (ifms:ifme,jfms:jfme) :: &
 117                    lfn, tign                ! level function, ignition time (state)
 118 integer, intent(out):: ignited              ! number of nodes newly ignited
 119
 120 !*** local
 121 integer:: i,j
 122 real::lfn_new,time_ign,ax,ay,rels,rele,d
 123 real:: sx,sy                    ! start of ignition line, from lower left corner
 124 real:: ex,ey                    ! end of ignition line, or zero
 125 real:: st,et                    ! start and end of time of the ignition line
 126 character(len=128):: msg
 127 real::cx2,cy2,dmax,axmin,axmax,aymin,aymax,dmin
 128 real:: start_x,start_y          ! start of ignition line, from lower left corner
 129 real:: end_x,end_y              ! end of ignition line, or zero
 130 real:: radius                   ! all within the radius of the line will ignite
 131 real:: start_time,end_time      ! the ignition time for the start and the end of the line
 132 real:: ros,tos                  ! ignition rate and time of spread
 133
 134 !*** executable
 135
 136 ! copy ignition line fields to local variables
 137 start_x    = ignition_line%start_x ! x coordinate of the ignition line start point (m, or long/lat)
 138 start_y    = ignition_line%start_y ! y coordinate of the ignition line start point
 139 end_x      = ignition_line%end_x   ! x coordinate of the ignition line end point
 140 end_y      = ignition_line%end_y   ! y coordinate of the ignition line end point
 141 start_time = ignition_line%start_time ! ignition time for the start point from simulation start (s)
 142 end_time   = ignition_line%end_time! ignition time for the end poin from simulation start (s)
 143 radius     = ignition_line%radius  ! all within this radius ignites immediately
 144 ros        = ignition_line%ros     ! rate of spread
 145
 146
 147 tos        = radius/ros            ! time of spread to the given radius
 148 st         = start_time            ! the start time of ignition considered in this time step
 149 et         = min(end_ts,end_time)  ! the end time of the ignition segment in this time step
 150
 151 ! this should be called whenever (start_ts, end_ts) \subset (start_time, end_time + tos)
 152 if(start_ts>et+tos .or. end_ts<st)return   ! too late or too early, nothing to do
 153
 154 if(start_time < end_time)then  ! we really want to test start_time .ne. end_time, but avoiding test of floats on equality
 155         ! segment of nonzero length
 156         !rels =  (st - start_time) / (end_time - start_time)  ! relative position of st in the segment (start,end)
 157         !sx = start_x + rels * (end_x - start_x)
 158         !sy = start_y + rels * (end_y - start_y)
 159         rels = 0.
 160         sx = start_x
 161         sy = start_y
 162         rele =  (et - start_time) / (end_time - start_time)    ! relative position of et in the segment (start,end)
 163         ex = start_x + rele * (end_x - start_x)
 164         ey = start_y + rele * (end_y - start_y)
 165 else
 166         ! just a point
 167         rels = 0.
 168         rele = 1.
 169         sx = start_x
 170         sy = start_y
 171         ex = end_x
 172         ey = end_y
 173 endif
 174
 175
 176 cx2=unit_xf*unit_xf
 177 cy2=unit_yf*unit_yf
 178
 179 axmin=coord_xf(ifts,jfts)
 180 aymin=coord_yf(ifts,jfts)
 181 axmax=coord_xf(ifte,jfte)
 182 aymax=coord_yf(ifte,jfte)
 183 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 184 write(msg,'(a,2f11.6,a,2f11.6)')'IGN from ',sx,sy,' to ',ex,ey
 185 call message(msg)
 186 write(msg,'(a,2f10.2,a,2f10.2,a)')'IGN timestep [',start_ts,end_ts,'] in [',start_time,end_time,']'
 187 call message(msg)
 188 write(msg,'(a,2g13.6,a,2g13.6)')'IGN tile coord from  ',axmin,aymin,' to ',axmax,aymax
 189 call message(msg)
 190 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 191 ignited=0
 192 dmax=0
 193 dmin=huge(dmax)
 194 11      format('IGN ',6(a,g17.7,1x))
 195 12      format('IGN ',4(a,2g13.7,1x))
 196 do j=jfts,jfte
 197     do i=ifts,ifte
 198         ax=coord_xf(i,j)
 199         ay=coord_yf(i,j)
 200
 201         ! get d= distance from the nearest point on the ignition segment
 202         ! and time_ign = the ignition time there
 203         call nearest(d,time_ign,ax,ay,sx,sy,st,ex,ey,et,cx2,cy2)
 204         dmax=max(d,dmax)
 205         dmin=min(d,dmin)
 206
 207         lfn_new=d - min( radius, ros*(end_ts - time_ign) )  ! lft at end_ts
 208         if(fire_print_msg.ge.3)then
 209 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 210             write(msg,*)'IGN1 i,j=',i,j,' lfn(i,j)=',lfn(i,j),' tign(i,j)=',tign(i,j)
 211             call message(msg)
 212             write(msg,*)'IGN2 i,j=',i,j,' lfn_new= ',lfn_new, ' time_ign= ',time_ign,' d=',d
 213             call message(msg)
 214 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 215         endif
 216         if(.not.lfn_new>0.) then
 217             ignited=ignited+1   ! count
 218         endif
 219         if(lfn(i,j)>0. .and. .not. lfn_new > 0.) then
 220             tign(i,j)=time_ign + d/ros  ! newly ignited now
 221             if(fire_print_msg.ge.3)then
 222 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 223                 write(msg,'(a,2i6,a,2g13.6,a,f10.2,a,2f10.2,a)')'IGN ignited cell ',i,j,' at',ax,ay, &
 224                     ' time',tign(i,j),' in [',start_ts,end_ts,']'
 225                 call message(msg)
 226                 write(msg,'(a,g10.3,a,f10.2,a,2f10.2,a)')'IGN distance',d,' from ignition line at',time_ign,' in [',st,et,']'
 227                 call message(msg)
 228 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 229             endif
 230             if(tign(i,j) < start_ts .or. tign(i,j) > end_ts )then
 231 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 232                 write(msg,'(a,2i6,a,f11.6,a,2f11.6,a)')'WARNING ',i,j, &
 233                 ' fixing ignition time ',tign(i,j),' outside of the time step [',start_ts,end_ts,']'
 234                 call message (msg)
 235 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 236                 tign(i,j) = min(max(tign(i,j),start_ts),end_ts)
 237             endif
 238         endif
 239         lfn(i,j)=min(lfn(i,j),lfn_new)  ! update the level set function
 240         if(fire_print_msg.ge.3)then
 241 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 242             write(msg,*)'IGN3 i,j=',i,j,' lfn(i,j)=',lfn(i,j),' tign(i,j)=',tign(i,j)
 243             call message(msg)
 244 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 245         endif
 246     enddo
 247 enddo
 248 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 249 write(msg,'(a,2g13.2,a,g10.2,a,g10.2)')'IGN units ',unit_xf,unit_yf,' m max dist ',dmax,' min',dmin
 250 call message(msg)
 251 write(msg,'(a,f6.1,a,f8.1,a,i10)')'IGN radius ',radius,' time of spread',tos,' ignited nodes',ignited
 252 call message(msg)
 253 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 254
 255 return
 256 99 continue
 257
 258 end subroutine ignite_fire
 259
 260 ! called from the inside of a loop, inline if possible
 261 !DEC$ ATTRIBUTES FORCEINLINE
 262 SUBROUTINE nearest(d,t,ax,ay,sx,sy,st,ex,ey,et,cx2,cy2)
 263         implicit none
 264 !***    arguments
 265         real, intent(out):: d,t
 266         real, intent(in):: ax,ay,sx,sy,st,ex,ey,et,cx2,cy2
 267         ! input:
 268         ! ax, ay       coordinates of point a
 269         ! sx,sy,ex,ey  coordinates of endpoints of segment [x,y]
 270         ! st,et        values at the endpoints x,y
 271         ! cx2,cy2      x and y unit squared for computing distance
 272
 273         ! output
 274         ! d            the distance of a and the nearest point z on the segment [x,y]
 275         ! t            linear interpolation from the values st,et to the point z
 276         !
 277         ! method: compute d as the distance (ax,ay) from the midpoint (mx,my)
 278         ! minus a correction (because of rounding errors): |a-c|^2 = |a-m|^2 - |m-c|^2
 279         ! when |m-c| >= |s-e|/2 the nearest point is one of the endpoints
 280         ! the computation work also for the case when s=e exactly or approximately
 281         !
 282         !
 283         !           a
 284         !          /| \
 285         !     s---m-c--e
 286         !
 287         ! |m-c| = |a-m| cos (a-m,e-s)
 288         !       = |a-m| (a-m).(e-s))/(|a-m|*|e-s|)
 289 !***    local
 290         real:: mx,my,dam2,dames,am_es,cos2,dmc2,mcrel,mid_t,dif_t,des2,cx,cy
 291         character(len=128):: msg
 292 !***    executable
 293
 294 11      format('IGN ',6(a,g17.7,1x))
 295 12      format('IGN ',4(a,2g13.7,1x))
 296
 297         ! midpoint m = (mx,my)
 298         mx = (sx + ex)*0.5
 299         my = (sy + ey)*0.5
 300         dam2=(ax-mx)*(ax-mx)*cx2+(ay-my)*(ay-my)*cy2      ! |a-m|^2
 301         des2 = (ex-sx)*(ex-sx)*cx2+(ey-sy)*(ey-sy)*cy2          ! des2 = |e-s|^2
 302         dames = dam2*des2
 303         am_es=(ax-mx)*(ex-sx)*cx2+(ay-my)*(ey-sy)*cy2       ! am_es = (a-m).(e-s)
 304         if(dames>0)then
 305             cos2 = (am_es*am_es)/dames                  ! cos2 = cos^2 (a-m,e-s)
 306         else ! point a already is the midpoint
 307             cos2 = 0.
 308         endif
 309         dmc2 = dam2*cos2                                ! dmc2 = |m-c|^2
 310         if(4.*dmc2 < des2)then                          ! if |m-c|<=|e-s|/2
 311             ! d = sqrt(max(dam2 - dmc2,0.))               ! d=|a-m|^2 - |m-c|^2, guard rounding
 312             mcrel = sign(sqrt(4.*dmc2/des2),am_es)      ! relative distance of c from m
 313         elseif(am_es>0)then                             ! if cos > 0, closest is e
 314             mcrel = 1.0
 315         else                                            ! closest is s
 316             mcrel = -1.0
 317         endif
 318         cx = (ex + sx)*0.5 + mcrel*(ex - sx)*0.5     ! interpolate to c by going from m
 319         cy = (ey + sy)*0.5 + mcrel*(ey - sy)*0.5     ! interpolate to c by going from m
 320         d=sqrt((ax-cx)*(ax-cx)*cx2+(ay-cy)*(ay-cy)*cy2) ! |a-c|^2
 321         t = (et + st)*0.5 + mcrel*(et - st)*0.5     ! interpolate to c by going from m
 322         if(fire_print_msg.ge.3)then
 323 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 324             write(msg,12)'find nearest to [',ax,ay,'] from [',sx,sy,'] [',ex,ey,']' ! DEB
 325             call message(msg)
 326             write(msg,12)'end times',st,et,' scale squared',cx2,cy2 ! DEB
 327             call message(msg)
 328             write(msg,11)'nearest at mcrel=',mcrel,'from the midpoint, t=',t ! DEB
 329             call message(msg)
 330             write(msg,12)'nearest is [',cx,cy,'] d=',d ! DEB
 331             call message(msg)
 332             write(msg,11)'dam2=',dam2,'des2=',des2,'dames=',dames
 333             call message(msg)
 334             write(msg,11)'am_es=',am_es,'cos2=',cos2,'dmc2=',dmc2 ! DEB
 335             call message(msg)
 336 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 337         endif
 338 END SUBROUTINE nearest
 339
 340
 341 !
 342 !**********************
 343 !
 344
 345 subroutine fuel_left(&
 346     ims,ime,jms,jme, &
 347     its,ite,jts,jte, &
 348     ifs,ife,jfs,jfe, &
 349     lfn, tign, fuel_time, time_now, fuel_frac, fire_area)
 350 implicit none
 351
 352 !*** purpose: determine fraction of fuel remaining
 353 !*** NOTE: because variables are cell centered, need halo/sync width 1 before
 354
 355 !*** arguments
 356
 357 integer, intent(in) :: its,ite,jts,jte,ims,ime,jms,jme,ifs,ife,jfs,jfe
 358 real, intent(in), dimension(ims:ime,jms:jme)::lfn,tign,fuel_time
 359 real, intent(in):: time_now
 360 real, intent(out), dimension(ifs:ife,jfs:jfe)::fuel_frac
 361 real, intent(out), dimension(ims:ime,jms:jme):: fire_area
 362
 363 ! ims,ime,jms,jme   in   memory dimensions
 364 ! its,ite,jts,jte   in   tile dimensions (cells where fuel_frac computed)
 365 ! ifs,ife,jfs,jfe   in   fuel_frac memory dimensions
 366 ! lfn               in   level function, at nodes at midpoints of cells
 367 ! tign              in   ignition time, at nodes at nodes at midpoints of cells
 368 ! fuel_time         in   time constant of fuel, per cell
 369 ! time_now          in   time now
 370 ! fuel_frac         out  fraction of fuel remaining, per cell
 371 ! fire_area         out  fraction of cell area on fire
 372
 373 !*** local
 374
 375 integer::i,j,ir,jr,icl,jcl,isubcl,jsubcl,i2,j2,ii,jj
 376 real::fmax,frat,helpsum1,helpsum2,fuel_left_ff,fire_area_ff,rx,ry,tignf(2,2)
 377 ! help variables instead of arrays fuel_left_f and fire_area_f
 378 real::lffij,lffi1j,lffij1,lffi1j1,tifij,tifi1j,tifij1,tifi1j1,tx,ty,txx,tyy
 379 ! variables for calculation instead of lff(i,j) and tif(i,j)is lffij,tifij etc..#define IFCELLS (ite-its+1)*fuel_left_irl
 380 #define JFCELLS (jte-jts+1)*fuel_left_jrl
 381 character(len=128)::msg
 382 integer::m,omp_get_thread_num
 383
 384
 385 call check_mesh_2dim(its-1,ite+1,jts-1,jte+1,ims,ime,jms,jme)
 386 call check_mesh_2dim(its,ite,jts,jte,ifs,ife,jfs,jfe)
 387
 388 ! refinement
 389 ir=fuel_left_irl
 390 jr=fuel_left_jrl
 391
 392 if ((ir.ne.2).or.(jr.ne.2)) then
 393    call crash('fuel_left: ir.ne.2 or jr.ne.2 ')
 394 endif
 395
 396 rx=1./ir
 397 ry=1./jr
 398
 399 ! interpolate level set function to finer grid
 400 #ifdef DEBUG_OUT_FUEL_LEFT
 401     call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,lff,'lff',0)
 402     call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,tif,'tif',0)
 403 #endif
 404
 405 !
 406 ! example for ir=2:
 407 !
 408 !                     |      coarse cell      |
 409 !      its-1                     its                                   ite             ite+1
 410 ! -------X------------|-----.-----X-----.-----|--........----|----------X----------|---------X
 411 !           fine node 1           2           3                   2*(ite-its+1)
 412 !                fine cell  1           2                             cell 2*(ite-its+1)
 413
 414
 415
 416 !  Loop over cells in Tile
 417 !  Changes made by Volodymyr Kondratenko 09/24/2009
 418 do icl=its,ite
 419   do jcl=jts,jte
 420     helpsum1=0
 421     helpsum2=0
 422 !   Loop over subcells in cell #(icl,jcl)
 423     do isubcl=1,ir
 424       do jsubcl=1,jr
 425         i=(icl-its)*ir+isubcl
 426         j=(jcl-jts)*jr+jsubcl
 427
 428 !       Direct calculation tif and lff, avoiding arrays, just for case ir=jr=2
 429         if ((isubcl.eq.1).and.(jsubcl.eq.1)) then
 430            i2=icl-1
 431            j2=jcl-1
 432            ty=0.5
 433            tx=0.5
 434            tyy=1.0
 435            txx=1.0
 436         else if ((isubcl.eq.2).and.(jsubcl.eq.1)) then
 437            i2=icl
 438            j2=jcl-1
 439            ty=0.5
 440            tx=0
 441            tyy=1.0
 442            txx=0.5
 443         else if ((isubcl.eq.1).and.(jsubcl.eq.2)) then
 444            i2=icl-1
 445            j2=jcl
 446            tx=0.5
 447            ty=0
 448            txx=1.0
 449            tyy=0.5
 450         else if ((isubcl.eq.2).and.(jsubcl.eq.2)) then
 451            i2=icl
 452            j2=jcl
 453            tx=0
 454            ty=0
 455            txx=0.5
 456            tyy=0.5
 457         else
 458            call crash('fuel_left: isubcl,jsubcl should be only 1 or 2')
 459         endif
 460
 461 ! calculation lff and tif in 4 endpoints of subcell
 462         lffij=                             &
 463                   (1-tx)*(1-ty)*lfn(i2,j2)      &
 464              +    (1-tx)*ty  *lfn(i2,j2+1)      &
 465              +     tx*(1-ty)*lfn(i2+1,j2)       &
 466              +       tx*ty  *lfn(i2+1,j2+1)
 467         lffi1j=                            &
 468                     (1-txx)*(1-ty)*lfn(i2,j2)   &
 469              +      (1-txx)*ty  *lfn(i2,j2+1)   &
 470              +      (txx)*(1-ty)*lfn(i2+1,j2)   &
 471              +      (txx)*ty  *lfn(i2+1,j2+1)
 472         lffij1=                            &
 473                     (1-tx)*(1-tyy)*lfn(i2,j2)   &
 474              +      (1-tx)*(tyy)  *lfn(i2,j2+1) &
 475              +      tx*(1-tyy)*lfn(i2+1,j2)     &
 476              +      tx*(tyy)  *lfn(i2+1,j2+1)
 477         lffi1j1 =                               &
 478                       (1-txx)*(1-tyy)*lfn(i2,j2)     &
 479              +      (1-txx)*(tyy)  *lfn(i2,j2+1)   &
 480              +      (txx)*(1-tyy)*lfn(i2+1,j2)     &
 481              +      (txx)*(tyy)  *lfn(i2+1,j2+1)
 482
 483         ! get ready to fix up tign values to be interpolated
 484         do ii=1,2
 485           do jj=1,2
 486             tignf(ii,jj)=tign(i2+ii-1,j2+jj-1)
 487           enddo
 488         enddo
 489         tifij=                                 &
 490                    (1-tx)*(1-ty)*tignf(1,1)        &
 491              +     (1-tx)*ty*tignf(1,1+1)          &
 492              +     tx*(1-ty)*tignf(1+1,1)          &
 493              +     tx*ty*tignf(1+1,1+1)
 494         tifi1j=                               &
 495                    (1-txx)*(1-ty)*tignf(1,1)      &
 496              +     (1-txx)*ty*tignf(1,1+1)        &
 497              +     (txx)*(1-ty)*tignf(1+1,1)      &
 498              +     (txx)*(ty)*tignf(1+1,1+1)
 499         tifij1=                               &
 500                    (1-tx)*(1-tyy)*tignf(1,1)      &
 501              +     (1-tx)*(tyy)*tignf(1,1+1)      &
 502              +      tx*(1-tyy)*tignf(1+1,1)       &
 503              +      tx*(tyy)*tignf(1+1,1+1)
 504         tifi1j1=                               &
 505                    (1-txx)*(1-tyy)*tignf(1,1)     &
 506              +     (1-txx)*(tyy)*tignf(1,1+1)     &
 507              +     (txx)*(1-tyy)*tignf(1+1,1)     &
 508              +     (txx)*(tyy)*tignf(1+1,1+1)
 509
 510
 511         if(fuel_left_method.eq.1)then
 512           call fuel_left_cell_1( fuel_left_ff, fire_area_ff, &
 513              lffij,lffij1,lffi1j,lffi1j1,&
 514              tifij,tifij1,tifi1j,tifi1j1,&
 515              time_now, fuel_time(icl,jcl))
 516         elseif(fuel_left_method.eq.2)then
 517           fire_area_ff=0  ! initialize to something until computed in fuel_left_f(i,j)
 518           fuel_left_ff=fuel_left_cell_2( &
 519              lffij,lffij1,lffi1j,lffi1j1,&
 520              tifij,tifij1,tifi1j,tifi1j1,&
 521              time_now, fuel_time(icl,jcl))
 522         else
 523           call crash('fuel_left: unknown fuel_left_method')
 524         endif
 525
 526         ! consistency check
 527         if(fire_area_ff.lt.-1e-6 .or.  &
 528           (fire_area_ff.eq.0. .and. fuel_left_ff.lt.1.-1e-6))then
 529 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 530            write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
 531               ' of refined mesh fuel burnt',1-fuel_left_ff,' fire area',fire_area_ff
 532 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 533            call crash(msg)
 534         endif
 535
 536         helpsum1=helpsum1+fuel_left_ff
 537         helpsum2=helpsum2+fire_area_ff
 538       enddo
 539     enddo
 540     fuel_frac(icl,jcl)=helpsum1
 541     fire_area(icl,jcl)=helpsum2
 542   enddo
 543 enddo
 544
 545
 546
 547
 548 #ifdef DEBUG_OUT_FUEL_LEFT
 549     call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fire_area,'fire_area',0)
 550     call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fuel_frac,'fuel_frac',0)
 551     call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fuel_left_f,'fuel_left_f',0)
 552     call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fire_area_f,'fire_area_f',0)
 553 #endif
 554
 555 ! finish the averaging
 556 do j=jts,jte
 557     do i=its,ite
 558         fuel_frac(i,j) = fuel_frac(i,j) /(ir*jr) ! multiply by weight for averaging over coarse cell
 559         fire_area(i,j) = fire_area(i,j) /(ir*jr) !
 560     enddo
 561 enddo
 562
 563 ! consistency check after sum
 564 fmax=0
 565 do j=jts,jte
 566     do i=its,ite
 567        if(fire_area(i,j).eq.0.)then
 568            if(fuel_frac(i,j).lt.1.-1e-6)then
 569 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 570                write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
 571                    ' fuel burnt',1-fuel_frac(i,j),' but fire area',fire_area(i,j)
 572 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 573                call crash(msg)
 574            endif
 575        else
 576            frat=(1-fuel_frac(i,j))/fire_area(i,j)
 577            fmax=max(fmax,frat)
 578        endif
 579     enddo
 580 enddo
 581 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 582 write(msg,'(a,4i6,a,f10.7)')'fuel_left: tile',its,ite,jts,jte,' max fuel burnt/area',fmax
 583 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 584 call message(msg)
 585 return
 586
 587
 588 end subroutine fuel_left
 589
 590 !
 591 !************************
 592 !
 593
 594 subroutine fuel_left_cell_1( fuel_frac_left, fire_frac_area, &
 595     lfn00,lfn01,lfn10,lfn11, &
 596     tign00,tign01,tign10,tign11,&
 597     time_now, fuel_time_cell)
 598 !*** purpose: compute the fuel fraction left in one cell
 599 implicit none
 600 !*** arguments
 601 real, intent(out):: fuel_frac_left, fire_frac_area !
 602 real, intent(in)::lfn00,lfn01,lfn10,lfn11    ! level set function at 4 corners of the cell
 603 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the  4 corners of the cell
 604 real, intent(in)::time_now                   ! the time now
 605 real, intent(in)::fuel_time_cell            ! time to burns off to 1/e
 606 !*** Description
 607 ! The area burning is given by the condition L <= 0, where the function P is
 608 ! interpolated from lfn(i,j)
 609 !
 610 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
 611 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
 612 ! when lfn(i,j)=0.
 613 !
 614 ! The function computes an approxmation  of the integral
 615 !
 616 !
 617 !                                  /\
 618 !                                  |
 619 ! fuel_frac_left  =      1   -     | 1 -  exp(-T(x,y)/fuel_time_cell)) dxdy
 620 !                                  |
 621 !                                 \/
 622 !                                0<x<1
 623 !                                0<y<1
 624 !                             L(x,y)<=0
 625 !
 626 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
 627 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
 628 ! so dx=1 and dy=1 assumed.
 629 !
 630 ! Example:
 631 !
 632 !        lfn<0         lfn>0
 633 !      (0,1) -----O--(1,1)            O = points on the fireline, T=tnow
 634 !            |      \ |               A = the burning area for computing
 635 !            |       \|                        fuel_frac(i,j)
 636 !            |   A    O
 637 !            |        |
 638 !            |        |
 639 !       (0,0)---------(1,0)
 640 !       lfn<0          lfn<0
 641 !
 642 ! Approximations allowed:
 643 ! The fireline can be approximated by straight line(s).
 644 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
 645 !
 646 ! Requirements:
 647 ! 1. The output should be a continuous function of the arrays lfn and
 648 !  tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
 649 ! 2. The output should be invariant to the symmetries of the input in each cell.
 650 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
 651 ! 4. The result should be at least 1st order accurate in the sense that it is
 652 !    exact if the time from ignition is a linear function.
 653 !
 654 ! If time from ignition is approximated by polynomial in the burnt
 655 ! region of the cell, this is integral of polynomial times exponential
 656 ! over a polygon, which can be computed exactly.
 657 !
 658 ! Requirement 4 is particularly important when there is a significant decrease
 659 ! of the fuel fraction behind the fireline on the mesh scale, because the
 660 ! rate of fuel decrease right behind the fireline is much larger
 661 ! (exponential...). This will happen when
 662 !
 663 ! change of time from ignition within one mesh cell / fuel_time_cell is not << 1
 664 !
 665 ! This is the same as
 666 !
 667 !               mesh cell size
 668 !  X =    -------------------------      is not << 1
 669 !       fireline speed * fuel_time_cell
 670 !
 671 !
 672 ! When X is large then the fuel burnt in one timestep in the cell is
 673 ! approximately proportional to length of  fireline in that cell.
 674 !
 675 ! When X is small then the fuel burnt in one timestep in the cell is
 676 ! approximately proportional to the area of the burning region.
 677 !
 678
 679 !*** calls
 680 intrinsic tiny
 681
 682 !*** local
 683 real::ps,aps,area,ta,out
 684 real::t00,t01,t10,t11
 685 real,parameter::safe=tiny(aps)
 686 character(len=128)::msg
 687
 688 ! the following algorithm is a very crude approximation
 689
 690 ! minus time since ignition, 0 if no ignition yet
 691 ! it is possible to have 0 in fire region when ignitin time falls in
 692 ! inside the time step because lfn is updated at the beginning of the time step
 693
 694 t00=tign00-time_now
 695 if(lfn00>0. .or. t00>0.)t00=0.
 696 t01=tign01-time_now
 697 if(lfn01>0. .or. t01>0.)t01=0.
 698 t10=tign10-time_now
 699 if(lfn10>0. .or. t10>0.)t10=0.
 700 t11=tign11-time_now
 701 if(lfn11>0. .or. t11>0.)t11=0.
 702
 703 ! approximate burning area, between 0 and 1
 704 ps = lfn00+lfn01+lfn10+lfn11
 705 aps = abs(lfn00)+abs(lfn01)+abs(lfn10)+abs(lfn11)
 706 aps=max(aps,safe)
 707 area =(-ps/aps+1.)/2.
 708 area = max(area,0.) ! make sure area is between 0 and 1
 709 area = min(area,1.)
 710
 711 ! average negative time since ignition
 712 ta=0.25*(t00+t01+t10+t11)
 713
 714 ! exp decay in the burning area
 715 out=1.
 716 !if(area>0.)out=1. - area*(1. - exp(ta/fuel_time_cell))
 717 if(area>0)out=area*exp(ta/fuel_time_cell) + (1. - area)
 718
 719 if(out>1.)then
 720 !$OMP CRITICAL(SFIRE_CORE_CRIT)
 721     write(msg,*)'out=',out,'>1 area=',area,' ta=',ta
 722     call message(msg)
 723     write(msg,*)'tign=',tign00,tign01,tign10,tign11,' time_now=',time_now
 724 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
 725     call message(msg)
 726     !call message('WARNING: fuel_left_cell_1: fuel fraction > 1')
 727     call crash('fuel_left_cell_1: fuel fraction > 1')
 728 endif
 729
 730 !out = max(out,0.) ! make sure out is between 0 and 1
 731 !out = min(out,1.)
 732
 733 fuel_frac_left = out
 734 fire_frac_area = area
 735
 736 end subroutine fuel_left_cell_1
 737
 738 !
 739 !****************************************
 740 !
 741
 742 real function fuel_left_cell_2(  &
 743     lfn00,lfn01,lfn10,lfn11, &
 744     tign00,tign01,tign10,tign11,&
 745     time_now, fuel_time_cell)
 746 !*** purpose: compute the fuel fraction left in one cell
 747 implicit none
 748 !*** arguments
 749 real, intent(in)::lfn00,lfn01,lfn10,lfn11    ! level set function at 4 corners of the cell
 750 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the  4 corners of the cell
 751 real, intent(in)::time_now                   ! the time now
 752 real, intent(in)::fuel_time_cell            ! time to burns off to 1/e
 753 !*** Description
 754 ! The area burning is given by the condition L <= 0, where the function P is
 755 ! interpolated from lfn(i,j)
 756 !
 757 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
 758 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
 759 ! when lfn(i,j)=0.
 760 !
 761 ! The function computes an approxmation  of the integral
 762 !
 763 !
 764 !                                  /\
 765 !                                  |
 766 ! fuel_frac_left  =      1   -     | 1 -  exp(-T(x,y)/fuel_time_cell)) dxdy
 767 !                                  |
 768 !                                 \/
 769 !                                0<x<1
 770 !                                0<y<1
 771 !                             L(x,y)<=0
 772 !
 773 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
 774 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
 775 ! so dx=1 and dy=1 assumed.
 776 !
 777 ! Example:
 778 !
 779 !        lfn<0         lfn>0
 780 !      (0,1) -----O--(1,1)            O = points on the fireline, T=tnow
 781 !            |      \ |               A = the burning area for computing
 782 !            |       \|                        fuel_frac(i,j)
 783 !            |   A    O
 784 !            |        |
 785 !            |        |
 786 !       (0,0)---------(1,0)
 787 !       lfn<0          lfn<0
 788 !
 789 ! Approximations allowed:
 790 ! The fireline can be approximated by straight line(s).
 791 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
 792 !
 793 ! Requirements:
 794 ! 1. The output should be a continuous function of the arrays lfn and
 795 !  tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
 796 ! 2. The output should be invariant to the symmetries of the input in each cell.
 797 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
 798 ! 4. The result should be at least 1st order accurate in the sense that it is
 799 !    exact if the time from ignition is a linear function.
 800 !
 801 ! If time from ignition is approximated by polynomial in the burnt
 802 ! region of the cell, this is integral of polynomial times exponential
 803 ! over a polygon, which can be computed exactly.
 804 !
 805 ! Requirement 4 is particularly important when there is a significant decrease
 806 ! of the fuel fraction behind the fireline on the mesh scale, because the
 807 ! rate of fuel decrease right behind the fireline is much larger
 808 ! (exponential...). This will happen when
 809 !
 810 ! change of time from ignition within one mesh cell * fuel speed is not << 1
 811 !
 812 ! This is the same as
 813 !
 814 !         mesh cell size*fuel_speed
 815 !         -------------------------      is not << 1
 816 !             fireline speed
 817 !
 818 !
 819 ! When X is large then the fuel burnt in one timestep in the cell is
 820 ! approximately proportional to length of  fireline in that cell.
 821 !
 822 ! When X is small then the fuel burnt in one timestep in the cell is
 823 ! approximately proportional to the area of the burning region.
 824 !
 825
 826 #ifndef FUEL_LEFT
 827 call crash('fuel_left_cell_2: not implemented, please use fire_fuel_left_method=1')
 828 fuel_left_cell_2=0.  ! to avoid compiler warning about value not set
 829 end function fuel_left_cell_2
 830 #else
 831
 832 !*** calls
 833 intrinsic tiny
 834
 835 !*** local
 836 real::ps,aps,area,ta,out
 837 real::t00,t01,t10,t11
 838 real,parameter::safe=tiny(aps)
 839 character(len=128)::msg
 840
 841 !*** local
 842 integer::i,j,k
 843
 844 ! least squares
 845 integer::mmax,nb,nmax,pmax,nin,nout
 846 parameter(mmax=3,nb=64,nmax=8,pmax=8)
 847 integer lda, ldb, lwork, info
 848 parameter (lda=nmax, ldb=nmax, lwork=mmax+nmax+nb*(nmax+pmax))
 849 integer n,m,p
 850 real,dimension(lda,mmax):: mA
 851 real,dimension(nmax):: vecD
 852 real,dimension(lwork):: WORK
 853 real,dimension(pmax):: vecY
 854 real,dimension(mmax):: vecX
 855 real,dimension(ldb,pmax)::mB
 856 real,dimension(mmax)::u
 857
 858 real::tweight,tdist
 859 integer::kk,ll,ss
 860 real::rnorm
 861 real,dimension(8,2)::xylist,xytlist
 862 real,dimension(8)::tlist,llist,xt
 863 real,dimension(5)::xx,yy
 864 real,dimension(5)::lfn,tign
 865
 866 integer:: npoint
 867 real::tt,x0,y0,xts,xte,yts,yte,xt1,xt2
 868 real::lfn0,lfn1,dist,nr,c,s,errQ,ae,ce,ceae,a0,a1,a2,d,cet
 869 real::s1,s2,s3
 870 real::upper,lower,ah,ch,aa,cc,aupp,cupp,alow,clow
 871 real,dimension(2,2)::mQ
 872 real,dimension(2)::ut
 873
 874 !calls
 875 intrinsic epsilon
 876
 877 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
 878
 879
 880 ! external functions
 881 real::snrm2
 882 double precision :: dnrm2
 883 external dnrm2
 884 external snrm2
 885 ! external subroutines
 886 external sggglm
 887 external dggglm
 888 !executable statements
 889
 890 ! a very crude approximation - replace by a better code
 891 ! call check_mesh_2dim(ids,ide+1,jds,jde+1,ims,ime,jms,jme)
 892 t00=tign00-time_now
 893 if(lfn00>=0. .or. t00>0.)t00=0.
 894 t01=tign01-time_now
 895 if(lfn01>=0. .or. t01>0.)t01=0.
 896 t10=tign10-time_now
 897 if(lfn10>=0. .or. t10>0.)t10=0.
 898 t11=tign11-time_now
 899 if(lfn11>=0. .or. t11>0.)t11=0.
 900
 901 !if (t00<0 .or. t01 <0 .or. t10<0 .or. t11<0) then
 902 !   print *,'tign00=',tign00,'tign10=',tign10,&
 903 !          'tign01=',tign01,'tign11=',tign11
 904 !end if
 905
 906
 907
 908 !*** case0 Do nothing
 909 if ( lfn00>=0 .and. lfn10>=0 .and. lfn01>=0 .and. lfn11>=0 ) then
 910     out = 1.0 !  fuel_left, no burning
 911 !*** case4 all four coners are burning
 912 else if (lfn00<=0 .and. lfn10<=0 .and. lfn01<=0 .and. lfn11<=0) then
 913     ! least squares, A matrix for points
 914     mA(1,1)=0.0
 915     mA(2,1)=1.0
 916     mA(3,1)=0.0
 917     mA(4,1)=1.0
 918     mA(1,2)=0.0
 919     mA(2,2)=0.0
 920     mA(3,2)=1.0
 921     mA(4,2)=1.0
 922     mA(1,3)=1.0
 923     mA(2,3)=1.0
 924     mA(3,3)=1.0
 925     mA(4,3)=1.0
 926     ! D vector, time from ignition
 927     vecD(1)=time_now-tign00
 928     vecD(2)=time_now-tign10
 929     vecD(3)=time_now-tign01
 930     vecD(4)=time_now-tign11
 931     ! B matrix, weights
 932     do kk=1,4
 933     do ll=1,4
 934         mB(kk,ll)=0.0
 935     end do
 936         mB(kk,kk)=2.0
 937     end do
 938     ! set the m,n,p
 939     n=4 ! rows of matrix A and B
 940     m=3 ! columns of matrix A
 941     p=4 ! columns of matrix B
 942     ! call least squqres in LAPACK
 943     call SGGGLM(N,M,P,mA,LDA,mB,LDB,vecD,vecX,vecY, &
 944                 WORK,LWORK,INFO)
 945     rnorm=snrm2(p,vecY,1)
 946     ! integrate
 947     u(1)=-vecX(1)/fuel_time_cell
 948     u(2)=-vecX(2)/fuel_time_cell
 949     u(3)=-vecX(3)/fuel_time_cell
 950     !fuel_burn(i,j)=1-exp(u(3))*intexp(u(1)*dx)*intexp(u(2)*dy)
 951     s1=u(1)
 952     s2=u(2)
 953     out=1-exp(u(3))*intexp(s1)*intexp(s2)
 954     !print *,'intexp
 955     if ( out<0 .or. out>1.0 ) then
 956         print *,'case4, out should be between 0 and 1'
 957     end if
 958 !*** case 1,2,3
 959 else
 960     ! set xx, yy for the coner points
 961     ! move these values out of i and j loop to speed up
 962     xx(1) = -0.5
 963     xx(2) =  0.5
 964     xx(3) =  0.5
 965     xx(4) = -0.5
 966     xx(5) = -0.5
 967     yy(1) = -0.5
 968     yy(2) = -0.5
 969     yy(3) =  0.5
 970     yy(4) =  0.5
 971     yy(5) = -0.5
 972     lfn(1)=lfn00
 973     lfn(2)=lfn10
 974     lfn(3)=lfn11
 975     lfn(4)=lfn01
 976     lfn(5)=lfn00
 977     tign(1)=t00
 978     tign(2)=t10
 979     tign(3)=t11
 980     tign(4)=t01
 981     tign(5)=t00
 982     npoint = 0 ! number of points in polygon
 983     !print *,'time_now=',time_now
 984     !print *,'lfn00=',lfn00,'lfn10=',lfn10,&
 985     !        'lfn01=',lfn01,'lfn11=',lfn11
 986     !print *,'t00=',t00,'t10=',t10,&
 987     !        't01=',t01,'t11=',t11
 988
 989     do k=1,4
 990         lfn0=lfn(k  )
 991         lfn1=lfn(k+1)
 992         if ( lfn0 <= 0.0 ) then
 993             npoint = npoint + 1
 994             xylist(npoint,1)=xx(k)
 995             xylist(npoint,2)=yy(k)
 996             tlist(npoint)=-tign(k)
 997             llist(npoint)=lfn0
 998         end if
 999         if ( lfn0*lfn1 < 0 ) then
1000             npoint = npoint + 1
1001             tt=lfn0/(lfn0-lfn1)
1002             x0=xx(k)+( xx(k+1)-xx(k) )*tt
1003             y0=yy(k)+( yy(k+1)-yy(k) )*tt
1004             xylist(npoint,1)=x0
1005             xylist(npoint,2)=y0
1006             tlist(npoint)=0 ! on fireline
1007             llist(npoint)=0
1008         end if
1009     end do
1010
1011     ! make the list circular
1012     tlist(npoint+1)=tlist(1)
1013     llist(npoint+1)=llist(1)
1014     xylist(npoint+1,1)=xylist(1,1)
1015     xylist(npoint+1,2)=xylist(1,2)
1016     !* least squares, A matrix for points
1017     do kk=1,npoint
1018         mA(kk,1)=xylist(kk,1)
1019         mA(kk,2)=xylist(kk,2)
1020         mA(kk,3)=1.0
1021         vecD(kk)=tlist(kk) ! D vector,time from ignition
1022     end do
1023     ! B matrix, weights
1024     do kk=1,ldb
1025     do ll=1,pmax
1026         mB(kk,ll)=0.0 ! clear
1027     end do
1028     end do
1029
1030     do kk=1,npoint
1031         mb(kk,kk) = 10 ! large enough
1032         do ll=1,npoint
1033             if ( kk .ne. ll ) then
1034                 dist = sqrt( (xylist(kk,1)-xylist(ll,1))**2+ &
1035                              (xylist(kk,2)-xylist(ll,2))**2 )
1036                 mB(kk,kk)=min( mB(kk,kk) , dist )
1037             end if
1038         end do !ll
1039         mB(kk,kk)=mB(kk,kk)+1.
1040     end do ! kk
1041     ! set the m,n,p
1042     n=npoint ! rows of matrix A and B
1043     m=3 ! columns of matrix A
1044     p=npoint ! columns of matrix B
1045     !* call least squqres in LAPACK
1046     call SGGGLM(N,M,P,mA,LDA,mB,LDB,vecD,vecX,vecY, &
1047                         WORK,LWORK,INFO)
1048     !print *,'after LS in case3'
1049     !print *,'vecX from LS',vecX
1050     !print *,'tign inputed',tign00,tign10,tign11,tign01
1051     rnorm=snrm2(p,vecY,1)
1052     u(1)=vecX(1)
1053     u(2)=vecX(2)
1054     u(3)=vecX(3)
1055     ! rotate to gradient on x only
1056     nr = sqrt(u(1)**2+u(2)**2)
1057     if(.not.nr.gt.eps)then
1058         out=1.
1059         goto 900
1060     endif
1061     c = u(1)/nr
1062     s = u(2)/nr
1063     mQ(1,1)=c
1064     mQ(1,2)=s
1065     mQ(2,1)=-s
1066     mQ(2,2)=c
1067     ! mat vec multiplication
1068     call matvec(mQ,2,2,u,3,ut,2,2,2)
1069     errQ = ut(2) ! should be zero
1070     ae = -ut(1)/fuel_time_cell
1071     ce = -u(3)/fuel_time_cell
1072     cet=ce!keep ce
1073     call matmatp(xylist,8,2,mQ,2,2,xytlist,8,2,npoint+1,2,2)
1074     call sortxt( xytlist, 8,2, xt,8,npoint )
1075     out=0.0
1076     aupp=0.0
1077     cupp=0.0
1078     alow=0.0
1079     clow=0.0
1080     do k=1,npoint-1
1081         xt1=xt(k)
1082         xt2=xt(k+1)
1083         upper=0
1084         lower=0
1085         ah=0
1086         ch=0
1087         if ( xt2-xt1 > eps*100 ) then
1088
1089             do ss=1,npoint
1090                 xts=xytlist(ss,1)
1091                 yts=xytlist(ss,2)
1092                 xte=xytlist(ss+1,1)
1093                 yte=xytlist(ss+1,2)
1094
1095                 if ( (xts>xt1 .and. xte>xt1) .or. &
1096                      (xts<xt2 .and. xte<xt2) ) then
1097                     aa = 0 ! do nothing
1098                     cc = 0
1099                 else
1100                     aa = (yts-yte)/(xts-xte)
1101                     cc = (xts*yte-xte*yts)/(xts-xte)
1102                     if (xte<xts) then
1103                         aupp = aa
1104                         cupp = cc
1105                         ah=ah+aa
1106                         ch=ch+cc
1107                         upper=upper+1
1108                     else
1109                         alow = aa
1110                         clow = cc
1111                         lower=lower+1
1112                     end if
1113                 end if!(xts>xt1 .and. xte>xt1)
1114             end do ! ss
1115             ce=cet !use stored ce
1116             if (ae*xt1+ce > 0 ) then
1117               ce=ce-(ae*xt1+ce)!shift small amounts exp(-**)
1118             end if
1119             if (ae*xt2+ce > 0) then
1120             ce=ce-(ae*xt2+ce)
1121             end if
1122
1123             ah = aupp-alow
1124             ch = cupp-clow
1125             ! integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2
1126             ! numerically sound for ae->0, ae -> infty
1127             ! this can be important for different model scales
1128             ! esp. if someone runs the model in single precision!!
1129             ! s1=int((ah*x+ch),x,xt1,xt2)
1130             s1 = (xt2-xt1)*((1./2.)*ah*(xt2+xt1)+ch)
1131             ! s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)
1132             ceae=ce/ae;
1133             s2 = -ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1))
1134             ! s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)
1135             ! s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)
1136             ! expand in Taylor series around ae=0
1137             ! collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)
1138             ! =(1/8*xt1^4+1/3*xt1^3*ceae+1/4*xt1^2*ceae^2-1/8*xt2^4-1/3*xt2^3*ceae-1/4*xt2^2*ceae^2)*ae^2
1139             !     + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae
1140             !     + 1/2*xt1^2-1/2*xt2^2
1141             !
1142             ! coefficient at ae^2 in the expansion, after some algebra
1143             a2=(xt1-xt2)*((1./4.)*(xt1+xt2)*ceae**2+(1./3.)* &
1144                (xt1**2+xt1*xt2+xt2**2)*ceae+(1./8.)* &
1145                (xt1**3+xt1*(xt2**2)+xt1**2*xt2+xt2**3))
1146             d=(ae**4)*a2
1147
1148             if (abs(d)>eps) then
1149             ! since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae
1150             ! for ae, ce -> 0 rounding error approx eps/ae^2
1151                 s3=( exp(ae*(xt1+ceae))*(ae*xt1-1)-&
1152                      exp(ae*(xt2+ceae))*(ae*xt2-1) )/(ae**2)
1153
1154             !we do not worry about rounding as xt1 -> xt2, then s3 -> 0
1155             else
1156                 ! coefficient at ae^1 in the expansion
1157                 a1=(xt1-xt2)*((1./2.)*ceae*(xt1+xt2)+(1./3.)*&
1158                    (xt1**2+xt1*xt2+xt2**2))
1159                 ! coefficient at ae^0 in the expansion for ae->0
1160                 a0=(1./2.)*(xt1-xt2)*(xt1+xt2)
1161                 s3=a0+a1*ae+a2*ae**2; ! approximate the integral
1162                             end if
1163
1164             s3=ah*s3
1165             out=out+s1+s2+s3
1166             out=1-out !fuel_left
1167             if(out<0 .or. out>1) then
1168                 print *,':fuel_fraction should be between 0 and 1'
1169                 !print *, 'eps= ', eps
1170                 !print *, 'xt1= ', xt1, 'xt2= ', xt2
1171                 !print *,'ae= ',ae,'ce= ',ce,'ah= ',ah,'ch= ',ch
1172                 !print *,'a0= ', a0,'a1= ', a1,'a2= ', a2
1173                 print *,'s1= ', s1,'s2= ', s2,'s3= ', s3
1174                 print *,':fuel_fraction =',out
1175             end if!print
1176
1177         end if
1178     end do ! k
1179
1180 end if ! if case0, elseif case4 ,else case123
1181
1182 900 continue
1183 if(out>1. .or. out<0.)call crash('fuel_left_cell_2: fuel fraction out of bounds [0,1]')
1184 fuel_left_cell_2 = out
1185 end function fuel_left_cell_2
1186
1187 !
1188 !****************************************
1189 !
1190 real function intexp(ab)
1191 implicit none
1192 real::ab
1193 !calls
1194 intrinsic epsilon
1195
1196 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
1197
1198 !eps = 2.2204*(10.0**(-8))!from matlab
1199 if ( eps < abs(ab)**3/6. ) then
1200     intexp=(exp(ab)-1)/ab
1201   else
1202     intexp=1+ab/2.
1203 end if
1204 end function
1205 !
1206 !****************************************
1207 !
1208 subroutine sortxt(xytlist,nrow,ncolumn,xt,nxt,nvec)
1209 implicit none
1210 integer::nrow,ncolumn,nxt,nvec
1211 real,dimension(nrow,ncolumn)::xytlist
1212 real,dimension(nxt)::xt
1213
1214 integer::i,j
1215 real::temp
1216
1217 do i=1,nvec
1218   xt(i)=xytlist(i,1)
1219 end do
1220
1221 do i=1,nvec-1
1222   do j=i+1,nvec
1223     if ( xt(i) > xt(j) ) then
1224          temp = xt(i)
1225          xt(i)=xt(j)
1226          xt(j)=temp
1227     end if
1228   end do
1229 end do
1230
1231 end subroutine !sortxt
1232 !
1233 !****************************************
1234 !
1235 subroutine matvec(A,m,n,V,nv,out,nout,nrow,ncolumn)
1236 implicit none
1237 integer::m,n,nv,nout,nrow,ncolumn
1238 real,dimension(m,n)::A   ! allocated m by n
1239 real,dimension(nv)::V    ! allocated nv
1240 real,dimension(nout)::out! allocated nout
1241
1242 integer::i,j
1243
1244 do i=1,nrow
1245   out(i)=0.0
1246   do j=1,ncolumn
1247     out(i)=out(i)+A(i,j)*V(j)
1248   end do
1249 end do
1250 end subroutine
1251 !
1252 !****************************************
1253 !
1254 subroutine matmatp(A,mA,nA,B,mB,nB,C,mC,nC,nrow,ncolumn,nP)
1255 implicit none
1256 integer::mA,nA,mB,nB,mC,nC,nrow,ncolumn,nP
1257 real,dimension(mA,nA)::A   ! allocated m by n
1258 real,dimension(mB,nB)::B   ! allocated m by n
1259 real,dimension(mC,nC)::C   ! allocated m by n
1260 integer::i,j,k
1261 do i=1,nrow
1262   do j=1,ncolumn
1263     C(i,j)=0.0
1264   do k=1,nP
1265     C(i,j)=C(i,j)+A(i,k)*B(j,k) ! B'
1266   end do
1267 end do
1268 end do
1269 end subroutine
1270
1271 !
1272 !****************************************
1273 !
1274 #endif
1275
1276 subroutine prop_ls( id, &                                ! for debug
1277                 ids,ide,jds,jde, &                       ! domain dims
1278                 ims,ime,jms,jme, &                       ! memory dims
1279                 ips,ipe,jps,jpe, &                ! patch - nodes owned by this process
1280                 its,ite,jts,jte, &                       ! tile dims
1281                 ts,dt,dx,dy,     &                       ! scalars in
1282                 tbound,          &                       ! scalars out
1283                 lfn_in,lfn_out,tign,ros,  &              ! arrays inout
1284                 fp               &
1285                    )
1286 implicit none
1287
1288 !*** purpose: advance level function in time
1289
1290 !*** description
1291 !
1292 ! Propagation of closed curve by a level function method. The level function
1293 ! lfn is defined by its values at the nodes of a rectangular grid.
1294 ! The area where lfn < 0 is inside the curve. The curve is
1295 ! described implicitly by lfn=0. Points where the curve intersects gridlines
1296 ! can be found by linear interpolation from nodes.
1297 !
1298 ! The level function is advanced from time ts to time ts + dt.
1299 !
1300 ! The level function should be initialized to (an approximation of) the signed
1301 ! distance from the curve. If the initial curve is a circle, the initial level
1302 ! function is simply the distance from the center minus the radius.
1303 !
1304 ! The curve moves outside with speed given by function speed_func.
1305 !
1306 ! Method: Godunov/ENO method for the normal motion. The timestep is checked for
1307 ! CFL condition. For a straight segment in a constant field and locally linear
1308 ! level function, the method reduces to the exact normal motion. The advantage of
1309 ! the level set method is that it treats automatically special cases such as
1310 ! the curve approaching itself and merging components of the area inside the curve.
1311 !
1312 ! Based on S. Osher and R. Fedkiw, Level set methods and dynamic implicit surfaces,
1313 ! Springer, 2003, Sec. 6.4, as implemented in toolboxLS for Matlab by
1314 ! I. Mitchell, A toolbox of Level Set Methods (Version 1.1), TR-2007-11,
1315 ! Dept. Computer Science, University of British Columbia, 2007
1316 ! http://www.cs.ubc.ca/\~mitchell/Toolbo\LS
1317 !
1318
1319 !*** arguments
1320
1321 ! id                in    unique identification for prints and dumps
1322 ! ids,ide,jds,jde   in    domain dimensions
1323 ! ims,ime,jms,jme   in    memory dimensions
1324 ! its,ite,jts,jte   in    tile dimensions
1325 ! ts                in    start time
1326 ! dt                in    time step
1327 ! dx,dy             in    grid spacing
1328 ! tbound            out   bound on stable time step from CFL condition, if tbound>=dt then OK
1329 ! lfn_in,lfn_out    inout,out the level set function at nodes
1330 ! tign              inout the ignition time at nodes
1331
1332 ! The dimensions are cell-based, the nodal value is associated with the south-west corner.
1333 ! The whole computation is on domain indices ids:ide+1,jds:jde+1.
1334 !
1335 ! The region where new lfn and tign are computed is the tile its:ite,jts:jte
1336 ! except when the tile is at domain upper boundary, an extra band of points is added:
1337 ! if ite=ide then region goes up to ite+1, if jte=jde then region goes up to jte+1.
1338
1339 ! The time step requires values from 2 rows of nodes beyond the region except when at the
1340 ! domain boundary one-sided derivatives are used. This is implemented by extending the input
1341 ! beyond the domain boundary so sufficient memory bounds must be allocated.
1342 ! The update on all tiles can be done in parallel. To avoid the race condition (different regions
1343 ! of the same array updated by different threads), the in and out versions of the
1344 ! arrays lft and tign are distinct. If the time step dt is larger
1345 ! that the returned tbound, the routine should be called again with timestep td<=tbound, and then
1346 ! having distinct inputs and outputs comes handy.
1347
1348 !*** calls
1349 !
1350 ! tend_ls
1351 !
1352
1353 integer,intent(in)::id,ims,ime,jms,jme,ids,ide,jds,jde,its,ite,jts,jte,ips,ipe,jps,jpe
1354 real,dimension(ims:ime,jms:jme),intent(inout)::lfn_in,tign
1355 real,dimension(ims:ime,jms:jme),intent(out)::lfn_out,ros
1356 real,intent(in)::dx,dy,ts,dt
1357 real,intent(out)::tbound
1358 type(fire_params),intent(in)::fp
1359
1360 !*** local
1361 ! arrays
1362 #define IMTS its-1
1363 #define IMTE ite+1
1364 #define JMTS jts-1
1365 #define JMTE jte+1
1366 real,dimension(IMTS:IMTE,JMTS:JMTE):: tend, lfn1 ! region-sized with halo
1367 ! scalars
1368 real::grad2,rr,tbound2,a,a1 ! a=0 euler, a=0.5 heun
1369
1370 real::gradx,grady,aspeed,err,aerr,time_now
1371 integer::ihs,ihe,jhs,jhe
1372 integer::ihs2,ihe2,jhs2,jhe2
1373 integer::i,j,its1,ite1,jts1,jte1,k,kk,id1
1374 character(len=128)::msg
1375 integer::nfirenodes,nfireline
1376 real::sum_err,min_err,max_err,sum_aerr,min_aerr,max_aerr
1377
1378 ! constants
1379 integer,parameter :: mstep=1000, printl=1
1380 real, parameter:: zero=0.,one=1.,eps=epsilon(zero),tol=100*eps, &
1381     safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
1382
1383 ! f90 intrinsic function
1384
1385 intrinsic max,min,sqrt,nint,epsilon,tiny,huge
1386
1387 !*** executable
1388
1389 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1390 write(msg,'(a8,i5,a6,i5,3(a1,i5))')'prop_ls:',id,' tile ',its,':',ite,',',jts,':',jte
1391 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1392 call message(msg)
1393
1394     a=fire_back_weight ! from module_fr_sfire_util
1395     a1=1. - a
1396
1397     ! tend = F(lfn)
1398
1399     ihs2=max(its-2,ids)   ! need lfn two beyond the tile but not outside the domain
1400     ihe2=min(ite+2,ide)
1401     jhs2=max(jts-2,jds)
1402     jhe2=min(jte+2,jde)
1403
1404     ihs=max(its-1,ids)   ! compute tend one beyond the tile but not outside the domain
1405     ihe=min(ite+1,ide)
1406     jhs=max(jts-1,jds)
1407     jhe=min(jte+1,jde)
1408
1409 #ifdef DEBUG_OUT
1410     call write_array_m(ihs,ihe,jhs,jhe,ims,ime,jms,jme,lfn_in,'lfn_in',id)
1411 #endif
1412
1413     ! check array dimensions
1414     call check_mesh_2dim(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme)
1415     call print_2d_stats(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme, &
1416                    lfn_in,'prop_ls: lfn in')
1417
1418     ! NOTE: tend_ls will extrapolate to one node strip at domain boundaries
1419     ! so that it can compute gradient at domain boundaries.
1420     ! To avoid copying, lfn_in is declared inout.
1421     ! At tile boundaries that are not domain boundaries values of lfn_in two nodes
1422     ! outside of the tile are needed.
1423     id1 = id  ! for debug prints
1424     if(id1.ne.0)id1=id1+1000
1425     call  tend_ls( id1, &
1426     ims,ime,jms,jme, &                       ! memory dims for lfn_in
1427     IMTS,IMTE,JMTS,JMTE, &                   ! memory dims for tend
1428     ids,ide,jds,jde, &                       ! domain dims - where lfn exists
1429     ips,ipe,jps,jpe, &                       ! patch - nodes owned by this process
1430     ihs,ihe,jhs,jhe, &                       ! where tend computed
1431     ims,ime,jms,jme, &                       ! memory dims for ros
1432     its,ite,jts,jte, &                       ! tile dims - where to set ros
1433     ts,dt,dx,dy,      &                      ! scalars in
1434     lfn_in, &                                ! arrays in
1435     tbound, &                                ! scalars out
1436     tend, ros, &                              ! arrays out
1437     fp         &                             ! params
1438 )
1439
1440 #ifdef DEBUG_OUT
1441     call write_array_m(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE,tend,'tend1',id)
1442 #endif
1443
1444     ! Euler method, the half-step, same region as ted
1445     do j=jhs,jhe
1446         do i=ihs,ihe
1447             lfn1(i,j) = lfn_in(i,j) + dt*tend(i,j)
1448         enddo
1449     enddo
1450
1451     call print_2d_stats(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE, &
1452                    lfn1,'prop_ls: lfn1')
1453     ! tend = F(lfn1) on the tile (not beyond)
1454
1455     if(id1.ne.0)id1=id1+1000
1456     call  tend_ls( id1,&
1457     IMTS,IMTE,JMTS,JMTE, &                   ! memory dims for lfn
1458     IMTS,IMTE,JMTS,JMTE, &                   ! memory dims for tend
1459     ids,ide,jds,jde,     &                   ! domain dims - where lfn exists
1460     ips,ipe,jps,jpe, &                       ! patch - nodes owned by this process
1461     its,ite,jts,jte, &                       ! tile dims - where is tend computed
1462     ims,ime,jms,jme, &                       ! memory dims for ros
1463     its,ite,jts,jte, &                       ! tile dims - where is ros set
1464     ts+dt,dt,dx,dy,      &                   ! scalars in
1465     lfn1, &                                  ! arrays in
1466     tbound2, &                               ! scalars out
1467     tend,ros, &                               ! arrays out
1468     fp &
1469 )
1470
1471 #ifdef DEBUG_OUT
1472     call write_array_m(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'tend2',id)
1473 #endif
1474
1475     call print_2d_stats(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'prop_ls: tend2')
1476
1477     tbound=min(tbound,tbound2)
1478
1479 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1480     write(msg,'(a,f10.2,4(a,f7.2))')'prop_ls: time',ts,' dt=',dt,' bound',min(tbound,999.99), &
1481         ' dx=',dx,' dy=',dy
1482 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1483     call message(msg)
1484     if(dt>tbound)then
1485 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1486         write(msg,'(2(a,f10.2))')'prop_ls: WARNING: time step ',dt, &
1487         ' > bound =',tbound
1488 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1489         call message(msg)
1490     endif
1491
1492     ! combine lfn1 and lfn_in + dt*tend -> lfn_out
1493
1494     do j=jts,jte
1495         do i=its,ite
1496             lfn_out(i,j) = a1*lfn1(i,j) + a*(lfn_in(i,j) + dt*tend(i,j))
1497             lfn_out(i,j) = min(lfn_out(i,j),lfn_in(i,j)) ! fire area can only increase
1498         enddo
1499     enddo
1500
1501     ! compute ignition time by interpolation
1502     ! the node was not burning at start but it is burning at end
1503     ! interpolate from the level functions at start and at end
1504     ! lfn_in   is the level set function value at time ts
1505     ! lfn_out  is the level set function value at time ts+dt
1506     ! 0        is the level set function value at time tign(i,j)
1507     ! thus assuming the level function is approximately linear =>
1508     ! tign(i,j)= ts + ((ts + td) - ts) * lfn_in / (lfn_in - lfn_out)
1509     !        = ts + dt * lfn_in / (lfn_in - lfn_out)
1510
1511     time_now=ts+dt
1512     time_now = time_now + abs(time_now)*epsilon(time_now)*2.
1513     do j=jts,jte
1514         do i=its,ite
1515             ! interpolate the cross-over time
1516             if (.not. lfn_out(i,j)>0 .and. lfn_in(i,j)>0)then
1517                 tign(i,j) = ts + dt * lfn_in(i,j) / (lfn_in(i,j) - lfn_out(i,j))
1518             endif
1519             ! set the ignition time outside of burning region
1520             if(lfn_out(i,j)>0.)tign(i,j)=time_now
1521         enddo
1522     enddo
1523
1524     ! check local speed error and stats
1525     ! may not work correctly in parallel
1526     ! init stats
1527     nfirenodes=0
1528     nfireline=0
1529     sum_err=0.
1530     min_err=rmax
1531     max_err=rmin
1532     sum_aerr=0.
1533     min_aerr=rmax
1534     max_aerr=rmin
1535     its1=its+1
1536     jts1=jts+1
1537     ite1=ite-1
1538     jte1=jte-1
1539     ! loop over right inside of the domain
1540     ! cannot use values outside of the domain, would have to reflect and that
1541     ! would change lfn_out
1542     ! cannot use values outside of tile, not synchronized yet
1543     ! so in parallel mode the statistics is not accurate
1544     do j=jts1,jte1
1545         do i=its1,ite1
1546             if(lfn_out(i,j)>0.0)then   ! a point out of burning region
1547                 if(lfn_out(i+1,j)<=0.or.lfn_out(i,j+1)<=0.or. & ! neighbor in burning region
1548                    lfn_out(i-1,j)<=0.or.lfn_out(i,j-1)<=0)then ! point next to fireline
1549                    gradx=(lfn_out(i+1,j)-lfn_out(i-1,j))/(2.0*dx) ! central differences
1550                    grady=(lfn_out(i,j+1)-lfn_out(i,j-1))/(2.0*dy)
1551                    grad2=sqrt(gradx*gradx+grady*grady)
1552                    aspeed = (lfn_in(i,j)-lfn_out(i,j))/(dt*max(grad2,rmin))
1553                     rr = speed_func(gradx,grady,dx,dy,i,j,fp)
1554                    err=aspeed-rr
1555                    sum_err=sum_err+err
1556                    min_err=min(min_err,err)
1557                    max_err=max(max_err,err)
1558                    aerr=abs(err)
1559                    sum_aerr=sum_aerr+aerr
1560                    min_aerr=min(min_aerr,aerr)
1561                    max_aerr=max(max_aerr,aerr)
1562                    nfireline=nfireline+1
1563                 endif
1564             else
1565                 nfirenodes=nfirenodes+1
1566             endif
1567         enddo
1568     enddo
1569 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1570     write(msg,'(2(a,i6,f8.4))')'prop_ls: nodes burning',nfirenodes, &
1571         (100.*nfirenodes)/((ite1-its1+1)*(jte1-jts1+1)),'% next to fireline',nfireline
1572 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1573     call message(msg)
1574     if(nfireline>0)then
1575         call print_stat_line('speed error',its1,ite1,jts1,jte1,min_err,max_err,sum_err/nfireline)
1576         call print_stat_line('abs(speed error)',its1,ite1,jts1,jte1,min_aerr,max_aerr,sum_aerr/nfireline)
1577     endif
1578
1579     ! check if the fire did not get to the domain boundary
1580     do k=-1,1,2
1581         ! do kk=1,(boundary_guard*(ide-ids+1))/100  ! in %
1582         do kk=1,boundary_guard   ! measured in cells
1583             i=ids+k*kk
1584             if(i.ge.its.and.i.le.ite)then
1585                 do j=jts,jte
1586                     if(lfn_out(i,j)<=0.)goto 9
1587                 enddo
1588             endif
1589     enddo
1590         ! do kk=1,(boundary_guard*(jde-jds+1))/100
1591         do kk=1,boundary_guard    ! measured in cells
1592             j=jds+k*kk
1593             if(j.ge.jts.and.j.le.jte)then
1594                 do i=its,ite
1595                     if(lfn_out(i,j)<=0.)goto 9
1596                 enddo
1597             endif
1598         enddo
1599     enddo
1600     goto 10
1601 9   continue
1602 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1603     write(msg,'(a,i2,a,2i8)')'prop_ls: fire',boundary_guard, &
1604         ' cells from domain boundary at node ',i,j
1605 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1606     call message(msg)
1607     call crash('prop_ls: increase the fire region')
1608 10  continue
1609
1610     call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme, &
1611                    lfn_out,'prop_ls: lfn out')
1612
1613 end subroutine prop_ls
1614
1615 !
1616 !*****************************
1617 !
1618
1619 subroutine tend_ls( id, &
1620     lims,lime,ljms,ljme, &                   ! memory dims for lfn
1621     tims,time,tjms,tjme, &                   ! memory dims for tend
1622     ids,ide,jds,jde, &                       ! domain - nodes where lfn defined
1623     ips,ipe,jps,jpe, &                       ! patch - nodes owned by this process
1624     ints,inte,jnts,jnte, &                   ! region - nodes where tend computed
1625     ims,ime,jms,jme, &                       ! memory dims for ros
1626     its,ite,jts,jte, &                       ! tile dims - where is ros set
1627     t,dt,dx,dy,      &                       ! scalars in
1628     lfn, &                                   ! arrays in
1629     tbound, &                                ! scalars out
1630     tend, ros,  &                              ! arrays out
1631     fp &
1632 )
1633
1634 implicit none
1635 ! purpose
1636 ! compute the right hand side of the level set equation
1637
1638 !*** arguments
1639 integer,intent(in)::id,lims,lime,ljms,ljme,tims,time,tjms,tjme
1640 integer,intent(in)::ims,ime,jms,jme,its,ite,jts,jte
1641 integer, intent(in)::ids,ide,jds,jde,ints,inte,jnts,jnte,ips,ipe,jps,jpe
1642 real,intent(in)::t                                     ! time
1643 real,intent(in)::dt,dx,dy                                 ! mesh step
1644 real,dimension(lims:lime,ljms:ljme),intent(inout)::lfn ! level set function
1645 real,intent(out)::tbound                               ! max allowed time step
1646 real,dimension(tims:time,tjms:tjme),intent(out)::tend  ! tendency (rhs of the level set pde)
1647 real,dimension(ims:ime,jms:jme),intent(out)::ros  ! rate of spread
1648 type(fire_params),intent(in)::fp
1649
1650 !*** local
1651 real:: te,diffLx,diffLy,diffRx,diffRy, &
1652    diffCx,diffCy,diff2x,diff2y,grad,rr, &
1653    ros_base,ros_wind,ros_slope,ros_back,advx,advy,scale,nvx,nvy, &
1654    speed,tanphi
1655 integer::i,j,itso,iteo,jtso,jteo
1656 character(len=128)msg
1657
1658 ! constants
1659 real, parameter:: eps=epsilon(0.0)
1660 !intrinsic epsilon
1661 real, parameter:: zero=0.,one=1.,tol=100*eps, &
1662     safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
1663
1664
1665 ! f90 intrinsic function
1666
1667 intrinsic max,min,sqrt,nint,tiny,huge
1668
1669
1670 #ifdef DEBUG_OUT
1671 real,dimension(tims:time,tjms:tjme)::rra,grada,speeda,tanphia
1672 #endif
1673
1674 !*** executable
1675
1676     ! check array dimensions
1677     call check_mesh_2dim(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme)
1678     call check_mesh_2dim(ints,inte,jnts,jnte,tims,time,tjms,tjme)
1679
1680     call continue_at_boundary(1,1,fire_lfn_ext_up, &   !extend by extrapolation but never down
1681     lims,lime,ljms,ljme, &                ! memory dims
1682     ids,ide,jds,jde, &                    ! domain - nodes where lfn defined
1683     ips,ipe,jps,jpe, &                    ! patch - nodes owned by this process
1684     ints,inte,jnts,jnte, &                ! tile - nodes update by this thread
1685     itso,iteo,jtso,jteo, &                ! where set now
1686     lfn)                                  ! array
1687
1688     call print_2d_stats(itso,iteo,jtso,jteo,lims,lime,ljms,ljme, &
1689                    lfn,'tend_ls: lfn cont')
1690
1691 #ifdef DEBUG_OUT
1692     call write_array_m(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme,lfn,'tend_lfn_in',id)
1693 #endif
1694
1695     tbound=0
1696     do j=jnts,jnte
1697         do i=ints,inte
1698             ! one sided differences
1699             diffRx = (lfn(i+1,j)-lfn(i,j))/dx
1700             diffLx = (lfn(i,j)-lfn(i-1,j))/dx
1701             diffRy = (lfn(i,j+1)-lfn(i,j))/dy
1702             diffLy = (lfn(i,j)-lfn(i,j-1))/dy
1703             diffCx = diffLx+diffRx   !  TWICE CENTRAL DIFFERENCE
1704             diffCy = diffLy+diffRy
1705
1706             !upwinding - select right or left derivative
1707             select case(fire_upwinding)
1708             case(0)  ! none
1709                 grad=sqrt(diffCx**2 + diffCy**2)
1710             case(1) ! standard
1711                 diff2x=select_upwind(diffLx,diffRx)
1712                 diff2y=select_upwind(diffLy,diffRy)
1713                 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1714             case(2) ! godunov per osher/fedkiw
1715                 diff2x=select_godunov(diffLx,diffRx)
1716                 diff2y=select_godunov(diffLy,diffRy)
1717                 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1718             case(3) ! eno
1719                 diff2x=select_eno(diffLx,diffRx)
1720                 diff2y=select_eno(diffLy,diffRy)
1721                 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1722             case(4) ! Sethian - twice stronger pushdown of bumps
1723                 grad=sqrt(max(diffLx,0.)**2+min(diffRx,0.)**2   &
1724                         + max(diffLy,0.)**2+min(diffRy,0.)**2)
1725             case default
1726                 grad=0.
1727             end select
1728
1729             ! normal direction, from central differences
1730             scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps)
1731             nvx=diffCx/scale
1732             nvy=diffCy/scale
1733
1734             ! wind speed in direction of spread
1735             speed =  fp%vx(i,j)*nvx + fp%vy(i,j)*nvy
1736
1737
1738             ! get rate of spread from wind speed and slope
1739
1740             call fire_ros(ros_base,ros_wind,ros_slope, &
1741             nvx,nvy,i,j,fp)
1742
1743             rr=ros_base + ros_wind + ros_slope
1744             if(fire_grows_only.gt.0)rr=max(rr,0.)
1745
1746             ! set ros for output
1747             if(i.ge.its.and.i.le.ite.and.j.ge.jts.and.j.le.jte)ros(i,j)=rr
1748
1749             if(fire_upwind_split.eq.0)then
1750
1751                 ! get rate of spread
1752                 te = -rr*grad   ! normal term
1753
1754             else
1755
1756                 ! normal direction backing rate only
1757                 te = - ros_base*grad
1758
1759                 ! advection in wind direction
1760                 if (abs(speed)> eps) then
1761                     advx=fp%vx(i,j)*ros_wind/speed
1762                     advy=fp%vy(i,j)*ros_wind/speed
1763                 else
1764                     advx=0
1765                     advy=0
1766                 endif
1767
1768                 tanphi =  fp%dzdxf(i,j)*nvx + fp%dzdyf(i,j)*nvy
1769                 ! advection from slope direction
1770                 if(abs(tanphi)>eps) then
1771                     advx=advx+fp%dzdxf(i,j)*ros_slope/tanphi
1772                     advy=advy+fp%dzdyf(i,j)*ros_slope/tanphi
1773                 endif
1774
1775                 if(fire_upwind_split.eq.1)then
1776
1777                     ! one-sided upwinding
1778                     te = te - max(advx,0.)*diffLx - min(advx,0.)*diffRy &
1779                             - max(advy,0.)*diffLy - min(advy,0.)*diffRy
1780
1781
1782                 elseif(fire_upwind_split.eq.2)then
1783
1784                     ! Lax-Friedrichs
1785                     call crash('prop_ls: bad fire_upwind_split, Lax-Friedrichs not done yet')
1786
1787                 else
1788
1789                     call crash('prop_ls: bad fire_upwind_split')
1790
1791                 endif
1792             endif
1793
1794             ! cfl condition
1795             if (grad > 0.) then
1796                  tbound = max(tbound,rr*(abs(diff2x)/dx+abs(diff2y)/dy)/grad)
1797             endif
1798
1799             ! add numerical viscosity
1800             te=te + fire_viscosity*abs(rr)*((diffRx-diffLx)+(diffRy-diffLy))
1801
1802             tend(i,j)=te
1803 #ifdef DEBUG_OUT
1804             rra(i,j)=rr
1805             grada(i,j)=grad
1806             speeda(i,j)=speed
1807             tanphia(i,j)=tanphi
1808 #endif
1809             !write(msg,*)i,j,grad,dzdx,dzdy
1810             !call message(msg)
1811
1812             !if(abs(te)>1e4)then
1813             !    write(msg,'(a,2i5,e12.3)')'tend_ls: tend out of bounds at ',i,j,te
1814             !    call crash(msg)
1815             !endif
1816         enddo
1817     enddo
1818
1819 #ifdef DEBUG_OUT
1820     call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,rra,'rr',id)
1821     call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,grada,'grad',id)
1822     call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,speeda,'speed',id)
1823     call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,tanphia,'tanphi',id)
1824     call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,tend,'tend',id)
1825 #endif
1826
1827     call print_2d_stats(ints,inte,jnts,jnte,tims,time,tjms,tjme, &
1828                    tend,'tend_ls: tend out')
1829
1830     ! the final CFL bound
1831     tbound = 1/(tbound+tol)
1832
1833 end subroutine tend_ls
1834
1835 !
1836 !**************************
1837 !
1838
1839 real function select_upwind(diffLx,diffRx)
1840 implicit none
1841 real, intent(in):: diffLx, diffRx
1842 real diff2x
1843
1844 ! upwind differences, L or R if bith same sign, otherwise zero
1845
1846 diff2x=0
1847 if (diffLx>0.and.diffRx>0.)diff2x=diffLx
1848 if (diffLx<0.and.diffRx<0.)diff2x=diffRx
1849
1850 select_upwind=diff2x
1851 end function select_upwind
1852
1853
1854 !
1855 !**************************
1856 !
1857
1858
1859 real function select_godunov(diffLx,diffRx)
1860 implicit none
1861 real, intent(in):: diffLx, diffRx
1862 real diff2x,diffCx
1863
1864 ! Godunov scheme: upwind differences, L or R or none
1865 ! always test on > or < never = , much faster because of IEEE
1866 ! central diff >= 0 => take left diff if >0, ortherwise 0
1867 ! central diff <= 0 => take right diff if <0, ortherwise 0
1868
1869 diff2x=0
1870 diffCx=diffRx+diffLx
1871 if (diffLx>0.and..not.diffCx<0)diff2x=diffLx
1872 if (diffRx<0.and.     diffCx<0)diff2x=diffRx
1873
1874 select_godunov=diff2x
1875 end function select_godunov
1876
1877 !
1878 !**************************
1879 !
1880
1881 real function select_eno(diffLx,diffRx)
1882 implicit none
1883 real, intent(in):: diffLx, diffRx
1884 real diff2x
1885
1886 ! 1st order ENO scheme
1887
1888 if    (.not.diffLx>0 .and. .not.diffRx>0)then
1889     diff2x=diffRx
1890 elseif(.not.diffLx<0 .and. .not.diffRx<0)then
1891     diff2x=diffLx
1892 elseif(.not.diffLx<0 .and. .not.diffRx>0)then
1893     if(.not. abs(diffRx) < abs(diffLx))then
1894         diff2x=diffRx
1895     else
1896         diff2x=diffLx
1897     endif
1898 else
1899     diff2x=0.
1900 endif
1901
1902 select_eno=diff2x
1903 end function select_eno
1904
1905 !
1906 !**************************
1907 !
1908
1909 real function speed_func(diffCx,diffCy,dx,dy,i,j,fp)
1910 !*** purpose
1911 !    the level set method speed function
1912 implicit none
1913 !*** arguments
1914 real, intent(in)::diffCx,diffCy  ! x and y coordinates of the direction of propagation
1915 real, intent(in)::dx,dy  ! x and y coordinates of the direction of propagation
1916 integer, intent(in)::i,j         ! indices of the node to compute the speed at
1917 type(fire_params),intent(in)::fp
1918 !*** local
1919 real::scale,nvx,nvy,r
1920 real::ros_base , ros_wind , ros_slope
1921 real, parameter:: eps=epsilon(0.0)
1922 !*** executable
1923             ! normal direction, from central differences
1924             scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps)
1925             nvx=diffCx/scale
1926             nvy=diffCy/scale
1927
1928             ! get rate of spread from wind speed and slope
1929
1930             call fire_ros(ros_base,ros_wind,ros_slope, &
1931             nvx,nvy,i,j,fp)
1932
1933             r=ros_base + ros_wind + ros_slope
1934             if(fire_grows_only.gt.0)r=max(r,0.)
1935             speed_func=r
1936
1937 end function speed_func
1938
1939 end module module_fr_sfire_core