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