2 !*** Jan Mandel August-October 2007 email: jmandel@ucar.edu or Jan.Mandel@gmail.com
4 ! With contributions by Minjeong Kim.
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.
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.
24 !****************************************
27 subroutine init_no_fire(&
28 ifds,ifde,jfds,jfde, &
29 ifms,ifme,jfms,jfme, &
30 ifts,ifte,jfts,jfte, &
31 fdx,fdy,time_now, & ! scalars in
32 fuel_frac,fire_area,lfn,tign) ! arrays out
35 !*** purpose: initialize model to no fire
38 integer, intent(in):: ifds,ifde,jfds,jfde ! fire domain bounds
39 integer, intent(in):: ifts,ifte,jfts,jfte ! fire tile bounds
40 integer, intent(in):: ifms,ifme,jfms,jfme ! array bounds
41 real, intent(in) :: fdx,fdy,time_now ! mesh spacing, time
42 real, intent(out), dimension (ifms:ifme,jfms:jfme) :: &
43 fuel_frac,fire_area,lfn,tign ! model state
50 real lfn_init,time_init
52 lfn_init = 2*max((ifde-ifds+1)*fdx,(jfde-jfds+1)*fdy) ! more than domain diameter
53 time_init=time_now + max(time_now,1.0)*epsilon(time_now) ! a bit in future
57 fuel_frac(i,j)=1. ! fuel at start is 1 by definition
58 fire_area(i,j)=0. ! nothing burning
59 tign(i,j) = time_init ! ignition in future
60 lfn(i,j) = lfn_init ! no fire
63 call message('init_model_no_fire: state set to no fire')
65 end subroutine init_no_fire
72 subroutine ignite_fire( ifds,ifde,jfds,jfde, & ! fire domain dims - the whole domain
73 ifms,ifme,jfms,jfme, &
74 ifts,ifte,jfts,jfte, &
75 start_x,start_y,end_x,end_y, &
83 !*** purpose: ignite a circular/line fire
86 ! ignite fire in the region within radius r from the line (sx,sy) to (ex,ey).
87 ! the coordinates of nodes are given as the arrays coord_xf and coord_yf
89 ! one unit of coord_xf is unit_xf m
90 ! one unit of coord_yf is unit_yf m
91 ! so a node (i,j) will be ignited iff for some (x,y) on the line
92 ! || ( (coord_xf(i,j) - x)*unit_xf , (coord_yf(i,j) - y)*unit_yf ) || <= r
96 integer, intent(in):: ifds,ifde,jfds,jfde ! fire domain bounds
97 integer, intent(in):: ifts,ifte,jfts,jfte ! fire tile bounds
98 integer, intent(in):: ifms,ifme,jfms,jfme ! array bounds
99 real, intent(in):: start_x,start_y ! start of ignition line, from lower left corner
100 real, intent(in):: end_x,end_y ! end of ignition line, or zero
101 real, intent(in):: r ! all within the radius of the line will ignite
102 real, intent(in):: start_t,end_t ! the ignition time for the start and the end of the line
103 real, intent(in):: start_ts,end_ts ! the time step start and end
104 real, dimension(ifms:ifme, jfms:jfme), intent(in):: &
105 coord_xf,coord_yf ! node coordinates
106 real, intent(in):: unit_xf,unit_yf ! coordinate units in m
107 real, intent(inout), dimension (ifms:ifme,jfms:jfme) :: &
108 lfn, tign ! level function, ignition time (state)
109 integer, intent(out):: ignited ! number of nodes newly ignited
113 real::mx,my,ax,ay,dam2,d,dames,des2,am_es,cos2,lfn_new,dmc2,time_ign,rels,rele,mid_t,cos_ame,dif_th
114 real:: sx,sy ! start of ignition line, from lower left corner
115 real:: ex,ey ! end of ignition line, or zero
116 real:: st,et ! start and end of time of the ignition line
117 character(len=128):: msg
120 st=max(start_ts,start_t) ! the start time of ignition in this time step
121 et=min(end_ts,end_t) ! the end time of ignition in this time step
123 if(st>et)return ! no ignition in this time step, nothing to do
126 ! find the points where the ignition line starts and ends this time step
127 ! compute in a way that is stable when end_t - start_t is small
128 if(start_t < end_t)then
129 rels = (st - start_t) / (end_t - start_t)
130 sx = start_x + rels * (end_x - start_x)
131 sy = start_y + rels * (end_y - start_y)
132 rele = (et - end_t) / (end_t - start_t)
133 ex = end_x + rele * (end_x - start_x)
134 ey = end_y + rele * (end_y - start_y)
142 mid_t = (end_t + start_t)*0.5 ! ignition time in the middle
143 dif_th = (end_t - start_t)*0.5
150 ! midpoint m = (mx,my)
155 ! coordinates of the node (i,j), the lower left corner of the domain is (0 0)
156 ! ax = fdx*(i - ifds + 0.5)
157 ! ay = fdy*(j - jfds + 0.5)
160 ! the following computation of distance is also for the case
161 ! when s=e exactly or approximately
162 dam2=(ax-mx)*(ax-mx)*cx2+(ay-my)*(ay-my)*cy2 ! |a-m|^2
163 ! compute distance of a=(ax,ay) and the nearest point on the segment
164 ! [(sx,sy), (ex,ey)] as the distance of (ax,ay) from the midpoint (mx,my)
165 ! minus a correction (because of rounding errors): |a-c|^2 = |a-m|^2 - |m-c|^2
166 ! when |m-c| >= |s-e|/2 the nearest point is one of the endpoints
172 ! |m-c| = |a-m| cos (a-m,e-s)
173 ! = |a-m| (a-m).(e-s))/(|a-m|*|e-s|)
174 des2 = (ex-sx)*(ex-sx)*cx2+(ey-sy)*(ey-sy)*cy2 ! des2 = |e-s|^2
176 am_es=(ax-mx)*(ex-sx)*cx2+(ay-my)*(ey-sy)*cy2 ! am_es = (a-m).(e-s)
178 cos2 = (am_es*am_es)/dames ! cos2 = cos^2 (a-m,e-s)
179 else ! point a already is the midpoint
182 dmc2 = dam2*cos2 ! dmc2 = |m-c|^2
183 if(4.*dmc2 <= des2)then ! if |m-c|<=|e-s|/2
184 d = sqrt(max(dam2 - dmc2,0.)) ! d=|a-m|^2 - |m-c|^2, guard rounding
185 elseif(am_es>0)then ! if cos > 0, closest is e
186 d = sqrt((ax-ex)*(ax-ex)*cx2+(ay-ey)*(ay-ey)*cy2) ! |a-e|
188 d = sqrt((ax-sx)*(ax-sx)*cx2+(ay-sy)*(ay-sy)*cy2) ! |a-s|
193 ignited=ignited+1 ! count
195 if(lfn(i,j)>0 .and. lfn_new<=0) then
196 cos_ame = sqrt(cos2) ! relative distance of c from m
197 time_ign = mid_t + sign(cos_ame,am_es)*dif_th ! ignition time at c by going from m
198 tign(i,j)=time_ign ! newly ignited now
200 lfn(i,j)=min(lfn(i,j),lfn_new) ! update the level set function
203 ! write(msg,'(2i4,10f12.6)')i,j,d,sx,ex,sy,ey
205 ! write(msg,'(10f12.6)')ax,ay,dmc2,des2,cos2
207 ! write(msg,'(10f16.6)')unit_xf,unit_yf,cx2,cy2,dam2
209 ! write(msg,'(10f16.6)')des2,dames,am_es,mx,my
213 !$OMP CRITICAL(SFIRE_CORE_CRIT)
214 write(msg,'(a,2f11.6,a,2f11.6)')'ignite_fire: from',sx,sy,' to ',ex,ey
216 write(msg,'(a,2f11.2,a,f8.1,a)')'units ',unit_xf,unit_yf,' m max dist ',dmax,' m'
218 write(msg,'(a,f4.1,a,f8.1,a,i10)')' radius ',r,' time',time_ign,' ignited nodes',ignited
220 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
221 end subroutine ignite_fire
224 !**********************
227 subroutine fuel_left(&
231 lfn, tign, fuel_time, time_now, fuel_frac, fire_area)
234 !*** purpose: determine fraction of fuel remaining
235 !*** NOTE: because variables are cell centered, need halo/sync width 1 before
237 !*** Jan Mandel August 2007 email: jmandel@ucar.edu or Jan.Mandel@gmail.com
241 integer, intent(in) :: its,ite,jts,jte,ims,ime,jms,jme,ifs,ife,jfs,jfe
242 real, intent(in), dimension(ims:ime,jms:jme)::lfn,tign,fuel_time
243 real, intent(in):: time_now
244 real, intent(out), dimension(ifs:ife,jfs:jfe)::fuel_frac
245 real, intent(out), dimension(ims:ime,jms:jme):: fire_area
247 ! ims,ime,jms,jme in memory dimensions
248 ! its,ite,jts,jte in tile dimensions (cells where fuel_frac computed)
249 ! ifs,ife,jfs,jfe in fuel_frac memory dimensions
250 ! lfn in level function, at nodes at midpoints of cells
251 ! tign in ignition time, at nodes at nodes at midpoints of cells
252 ! fuel_time in time constant of fuel, per cell
253 ! time_now in time now
254 ! fuel_frac out fraction of fuel remaining, per cell
255 ! fire_area out fraction of cell area on fire
259 integer::i,j,ir,jr,icl,jcl,isubcl,jsubcl,i2,j2,ii,jj
260 real::fmax,frat,helpsum1,helpsum2,fuel_left_ff,fire_area_ff,rx,ry,tignf(2,2)
261 ! help variables instead of arrays fuel_left_f and fire_area_f
262 real::lffij,lffi1j,lffij1,lffi1j1,tifij,tifi1j,tifij1,tifi1j1,tx,ty,txx,tyy
263 ! variables for calculation instead of lff(i,j) and tif(i,j)is lffij,tifij etc..#define IFCELLS (ite-its+1)*fuel_left_irl
264 #define JFCELLS (jte-jts+1)*fuel_left_jrl
265 character(len=128)::msg
266 integer::m,omp_get_thread_num
269 call check_mesh_2dim(its-1,ite+1,jts-1,jte+1,ims,ime,jms,jme)
270 call check_mesh_2dim(its,ite,jts,jte,ifs,ife,jfs,jfe)
276 if ((ir.ne.2).or.(jr.ne.2)) then
277 call crash('fuel_left: ir.ne.2 or jr.ne.2 ')
283 ! interpolate level set function to finer grid
284 #ifdef DEBUG_OUT_FUEL_LEFT
285 call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,lff,'lff',0)
286 call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,tif,'tif',0)
293 ! its-1 its ite ite+1
294 ! -------X------------|-----.-----X-----.-----|--........----|----------X----------|---------X
295 ! fine node 1 2 3 2*(ite-its+1)
296 ! fine cell 1 2 cell 2*(ite-its+1)
300 ! Loop over cells in Tile
301 ! Changes made by Volodymyr Kondratenko 09/24/2009
306 ! Loop over subcells in cell #(icl,jcl)
309 i=(icl-its)*ir+isubcl
310 j=(jcl-jts)*jr+jsubcl
312 ! Direct calculation tif and lff, avoiding arrays, just for case ir=jr=2
313 if ((isubcl.eq.1).and.(jsubcl.eq.1)) then
320 else if ((isubcl.eq.2).and.(jsubcl.eq.1)) then
327 else if ((isubcl.eq.1).and.(jsubcl.eq.2)) then
334 else if ((isubcl.eq.2).and.(jsubcl.eq.2)) then
342 call crash('fuel_left: isubcl,jsubcl should be only 1 or 2')
345 ! calculation lff and tif in 4 endpoints of subcell
347 (1-tx)*(1-ty)*lfn(i2,j2) &
348 + (1-tx)*ty *lfn(i2,j2+1) &
349 + tx*(1-ty)*lfn(i2+1,j2) &
350 + tx*ty *lfn(i2+1,j2+1)
352 (1-txx)*(1-ty)*lfn(i2,j2) &
353 + (1-txx)*ty *lfn(i2,j2+1) &
354 + (txx)*(1-ty)*lfn(i2+1,j2) &
355 + (txx)*ty *lfn(i2+1,j2+1)
357 (1-tx)*(1-tyy)*lfn(i2,j2) &
358 + (1-tx)*(tyy) *lfn(i2,j2+1) &
359 + tx*(1-tyy)*lfn(i2+1,j2) &
360 + tx*(tyy) *lfn(i2+1,j2+1)
362 (1-txx)*(1-tyy)*lfn(i2,j2) &
363 + (1-txx)*(tyy) *lfn(i2,j2+1) &
364 + (txx)*(1-tyy)*lfn(i2+1,j2) &
365 + (txx)*(tyy) *lfn(i2+1,j2+1)
367 ! get ready to fix up tign values to be interpolated
370 tignf(ii,jj)=tign(i2+ii-1,j2+jj-1)
374 (1-tx)*(1-ty)*tignf(1,1) &
375 + (1-tx)*ty*tignf(1,1+1) &
376 + tx*(1-ty)*tignf(1+1,1) &
377 + tx*ty*tignf(1+1,1+1)
379 (1-txx)*(1-ty)*tignf(1,1) &
380 + (1-txx)*ty*tignf(1,1+1) &
381 + (txx)*(1-ty)*tignf(1+1,1) &
382 + (txx)*(ty)*tignf(1+1,1+1)
384 (1-tx)*(1-tyy)*tignf(1,1) &
385 + (1-tx)*(tyy)*tignf(1,1+1) &
386 + tx*(1-tyy)*tignf(1+1,1) &
387 + tx*(tyy)*tignf(1+1,1+1)
389 (1-txx)*(1-tyy)*tignf(1,1) &
390 + (1-txx)*(tyy)*tignf(1,1+1) &
391 + (txx)*(1-tyy)*tignf(1+1,1) &
392 + (txx)*(tyy)*tignf(1+1,1+1)
395 if(fuel_left_method.eq.1)then
396 call fuel_left_cell_1( fuel_left_ff, fire_area_ff, &
397 lffij,lffij1,lffi1j,lffi1j1,&
398 tifij,tifij1,tifi1j,tifi1j1,&
399 time_now, fuel_time(icl,jcl))
400 elseif(fuel_left_method.eq.2)then
401 fire_area_ff=0 ! initialize to something until computed in fuel_left_f(i,j)
402 fuel_left_ff=fuel_left_cell_2( &
403 lffij,lffij1,lffi1j,lffi1j1,&
404 tifij,tifij1,tifi1j,tifi1j1,&
405 time_now, fuel_time(icl,jcl))
407 call crash('fuel_left: unknown fuel_left_method')
411 if(fire_area_ff.lt.-1e-6 .or. &
412 (fire_area_ff.eq.0. .and. fuel_left_ff.lt.1.-1e-6))then
413 !$OMP CRITICAL(SFIRE_CORE_CRIT)
414 write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
415 ' of refined mesh fuel burnt',1-fuel_left_ff,' fire area',fire_area_ff
416 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
420 helpsum1=helpsum1+fuel_left_ff
421 helpsum2=helpsum2+fire_area_ff
424 fuel_frac(icl,jcl)=helpsum1
425 fire_area(icl,jcl)=helpsum2
432 #ifdef DEBUG_OUT_FUEL_LEFT
433 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fire_area,'fire_area',0)
434 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fuel_frac,'fuel_frac',0)
435 call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fuel_left_f,'fuel_left_f',0)
436 call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fire_area_f,'fire_area_f',0)
439 ! finish the averaging
442 fuel_frac(i,j) = fuel_frac(i,j) /(ir*jr) ! multiply by weight for averaging over coarse cell
443 fire_area(i,j) = fire_area(i,j) /(ir*jr) !
447 ! consistency check after sum
451 if(fire_area(i,j).eq.0.)then
452 if(fuel_frac(i,j).lt.1.-1e-6)then
453 !$OMP CRITICAL(SFIRE_CORE_CRIT)
454 write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
455 ' fuel burnt',1-fuel_frac(i,j),' but fire area',fire_area(i,j)
456 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
460 frat=(1-fuel_frac(i,j))/fire_area(i,j)
465 !$OMP CRITICAL(SFIRE_CORE_CRIT)
466 write(msg,'(a,4i6,a,f10.7)')'fuel_left: tile',its,ite,jts,jte,' max fuel burnt/area',fmax
467 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
472 end subroutine fuel_left
475 !************************
478 subroutine fuel_left_cell_1( fuel_frac_left, fire_frac_area, &
479 lfn00,lfn01,lfn10,lfn11, &
480 tign00,tign01,tign10,tign11,&
481 time_now, fuel_time_cell)
482 !*** purpose: compute the fuel fraction left in one cell
485 real, intent(out):: fuel_frac_left, fire_frac_area !
486 real, intent(in)::lfn00,lfn01,lfn10,lfn11 ! level set function at 4 corners of the cell
487 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the 4 corners of the cell
488 real, intent(in)::time_now ! the time now
489 real, intent(in)::fuel_time_cell ! time to burns off to 1/e
491 ! The area burning is given by the condition L <= 0, where the function P is
492 ! interpolated from lfn(i,j)
494 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
495 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
498 ! The function computes an approxmation of the integral
503 ! fuel_frac_left = 1 - | 1 - exp(-T(x,y)/fuel_time_cell)) dxdy
510 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
511 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
512 ! so dx=1 and dy=1 assumed.
517 ! (0,1) -----O--(1,1) O = points on the fireline, T=tnow
518 ! | \ | A = the burning area for computing
519 ! | \| fuel_frac(i,j)
523 ! (0,0)---------(1,0)
526 ! Approximations allowed:
527 ! The fireline can be approximated by straight line(s).
528 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
531 ! 1. The output should be a continuous function of the arrays lfn and
532 ! tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
533 ! 2. The output should be invariant to the symmetries of the input in each cell.
534 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
535 ! 4. The result should be at least 1st order accurate in the sense that it is
536 ! exact if the time from ignition is a linear function.
538 ! If time from ignition is approximated by polynomial in the burnt
539 ! region of the cell, this is integral of polynomial times exponential
540 ! over a polygon, which can be computed exactly.
542 ! Requirement 4 is particularly important when there is a significant decrease
543 ! of the fuel fraction behind the fireline on the mesh scale, because the
544 ! rate of fuel decrease right behind the fireline is much larger
545 ! (exponential...). This will happen when
547 ! change of time from ignition within one mesh cell / fuel_time_cell is not << 1
549 ! This is the same as
552 ! X = ------------------------- is not << 1
553 ! fireline speed * fuel_time_cell
556 ! When X is large then the fuel burnt in one timestep in the cell is
557 ! approximately proportional to length of fireline in that cell.
559 ! When X is small then the fuel burnt in one timestep in the cell is
560 ! approximately proportional to the area of the burning region.
567 real::ps,aps,area,ta,out
568 real::t00,t01,t10,t11
569 real,parameter::safe=tiny(aps)
570 character(len=128)::msg
572 ! the following algorithm is a very crude approximation
574 ! minus time since ignition, 0 if no ignition yet
575 ! it is possible to have 0 in fire region when ignitin time falls in
576 ! inside the time step because lfn is updated at the beginning of the time step
579 if(lfn00>0. .or. t00>0.)t00=0.
581 if(lfn01>0. .or. t01>0.)t01=0.
583 if(lfn10>0. .or. t10>0.)t10=0.
585 if(lfn11>0. .or. t11>0.)t11=0.
587 ! approximate burning area, between 0 and 1
588 ps = lfn00+lfn01+lfn10+lfn11
589 aps = abs(lfn00)+abs(lfn01)+abs(lfn10)+abs(lfn11)
591 area =(-ps/aps+1.)/2.
592 area = max(area,0.) ! make sure area is between 0 and 1
595 ! average negative time since ignition
596 ta=0.25*(t00+t01+t10+t11)
598 ! exp decay in the burning area
600 !if(area>0.)out=1. - area*(1. - exp(ta/fuel_time_cell))
601 if(area>0)out=area*exp(ta/fuel_time_cell) + (1. - area)
604 !$OMP CRITICAL(SFIRE_CORE_CRIT)
605 write(msg,*)'out=',out,'>1 area=',area,' ta=',ta
607 write(msg,*)'tign=',tign00,tign01,tign10,tign11,' time_now=',time_now
608 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
610 !call message('WARNING: fuel_left_cell_1: fuel fraction > 1')
611 call crash('fuel_left_cell_1: fuel fraction > 1')
614 !out = max(out,0.) ! make sure out is between 0 and 1
618 fire_frac_area = area
620 end subroutine fuel_left_cell_1
623 !****************************************
626 real function fuel_left_cell_2( &
627 lfn00,lfn01,lfn10,lfn11, &
628 tign00,tign01,tign10,tign11,&
629 time_now, fuel_time_cell)
630 !*** purpose: compute the fuel fraction left in one cell
633 real, intent(in)::lfn00,lfn01,lfn10,lfn11 ! level set function at 4 corners of the cell
634 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the 4 corners of the cell
635 real, intent(in)::time_now ! the time now
636 real, intent(in)::fuel_time_cell ! time to burns off to 1/e
638 ! The area burning is given by the condition L <= 0, where the function P is
639 ! interpolated from lfn(i,j)
641 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
642 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
645 ! The function computes an approxmation of the integral
650 ! fuel_frac_left = 1 - | 1 - exp(-T(x,y)/fuel_time_cell)) dxdy
657 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
658 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
659 ! so dx=1 and dy=1 assumed.
664 ! (0,1) -----O--(1,1) O = points on the fireline, T=tnow
665 ! | \ | A = the burning area for computing
666 ! | \| fuel_frac(i,j)
670 ! (0,0)---------(1,0)
673 ! Approximations allowed:
674 ! The fireline can be approximated by straight line(s).
675 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
678 ! 1. The output should be a continuous function of the arrays lfn and
679 ! tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
680 ! 2. The output should be invariant to the symmetries of the input in each cell.
681 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
682 ! 4. The result should be at least 1st order accurate in the sense that it is
683 ! exact if the time from ignition is a linear function.
685 ! If time from ignition is approximated by polynomial in the burnt
686 ! region of the cell, this is integral of polynomial times exponential
687 ! over a polygon, which can be computed exactly.
689 ! Requirement 4 is particularly important when there is a significant decrease
690 ! of the fuel fraction behind the fireline on the mesh scale, because the
691 ! rate of fuel decrease right behind the fireline is much larger
692 ! (exponential...). This will happen when
694 ! change of time from ignition within one mesh cell * fuel speed is not << 1
696 ! This is the same as
698 ! mesh cell size*fuel_speed
699 ! ------------------------- is not << 1
703 ! When X is large then the fuel burnt in one timestep in the cell is
704 ! approximately proportional to length of fireline in that cell.
706 ! When X is small then the fuel burnt in one timestep in the cell is
707 ! approximately proportional to the area of the burning region.
711 call crash('fuel_left_cell_2: not implemented, please use fire_fuel_left_method=1')
712 fuel_left_cell_2=0. ! to avoid compiler warning about value not set
713 end function fuel_left_cell_2
720 real::ps,aps,area,ta,out
721 real::t00,t01,t10,t11
722 real,parameter::safe=tiny(aps)
723 character(len=128)::msg
729 integer::mmax,nb,nmax,pmax,nin,nout
730 parameter(mmax=3,nb=64,nmax=8,pmax=8)
731 integer lda, ldb, lwork, info
732 parameter (lda=nmax, ldb=nmax, lwork=mmax+nmax+nb*(nmax+pmax))
734 real,dimension(lda,mmax):: mA
735 real,dimension(nmax):: vecD
736 real,dimension(lwork):: WORK
737 real,dimension(pmax):: vecY
738 real,dimension(mmax):: vecX
739 real,dimension(ldb,pmax)::mB
740 real,dimension(mmax)::u
745 real,dimension(8,2)::xylist,xytlist
746 real,dimension(8)::tlist,llist,xt
747 real,dimension(5)::xx,yy
748 real,dimension(5)::lfn,tign
751 real::tt,x0,y0,xts,xte,yts,yte,xt1,xt2
752 real::lfn0,lfn1,dist,nr,c,s,errQ,ae,ce,ceae,a0,a1,a2,d,cet
754 real::upper,lower,ah,ch,aa,cc,aupp,cupp,alow,clow
755 real,dimension(2,2)::mQ
756 real,dimension(2)::ut
761 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
766 double precision :: dnrm2
769 ! external subroutines
772 !executable statements
774 ! a very crude approximation - replace by a better code
775 ! call check_mesh_2dim(ids,ide+1,jds,jde+1,ims,ime,jms,jme)
777 if(lfn00>=0. .or. t00>0.)t00=0.
779 if(lfn01>=0. .or. t01>0.)t01=0.
781 if(lfn10>=0. .or. t10>0.)t10=0.
783 if(lfn11>=0. .or. t11>0.)t11=0.
785 !if (t00<0 .or. t01 <0 .or. t10<0 .or. t11<0) then
786 ! print *,'tign00=',tign00,'tign10=',tign10,&
787 ! 'tign01=',tign01,'tign11=',tign11
792 !*** case0 Do nothing
793 if ( lfn00>=0 .and. lfn10>=0 .and. lfn01>=0 .and. lfn11>=0 ) then
794 out = 1.0 ! fuel_left, no burning
795 !*** case4 all four coners are burning
796 else if (lfn00<=0 .and. lfn10<=0 .and. lfn01<=0 .and. lfn11<=0) then
797 ! least squares, A matrix for points
810 ! D vector, time from ignition
811 vecD(1)=time_now-tign00
812 vecD(2)=time_now-tign10
813 vecD(3)=time_now-tign01
814 vecD(4)=time_now-tign11
823 n=4 ! rows of matrix A and B
824 m=3 ! columns of matrix A
825 p=4 ! columns of matrix B
826 ! call least squqres in LAPACK
827 call SGGGLM(N,M,P,mA,LDA,mB,LDB,vecD,vecX,vecY, &
829 rnorm=snrm2(p,vecY,1)
831 u(1)=-vecX(1)/fuel_time_cell
832 u(2)=-vecX(2)/fuel_time_cell
833 u(3)=-vecX(3)/fuel_time_cell
834 !fuel_burn(i,j)=1-exp(u(3))*intexp(u(1)*dx)*intexp(u(2)*dy)
837 out=1-exp(u(3))*intexp(s1)*intexp(s2)
839 if ( out<0 .or. out>1.0 ) then
840 print *,'case4, out should be between 0 and 1'
844 ! set xx, yy for the coner points
845 ! move these values out of i and j loop to speed up
866 npoint = 0 ! number of points in polygon
867 !print *,'time_now=',time_now
868 !print *,'lfn00=',lfn00,'lfn10=',lfn10,&
869 ! 'lfn01=',lfn01,'lfn11=',lfn11
870 !print *,'t00=',t00,'t10=',t10,&
871 ! 't01=',t01,'t11=',t11
876 if ( lfn0 <= 0.0 ) then
878 xylist(npoint,1)=xx(k)
879 xylist(npoint,2)=yy(k)
880 tlist(npoint)=-tign(k)
883 if ( lfn0*lfn1 < 0 ) then
886 x0=xx(k)+( xx(k+1)-xx(k) )*tt
887 y0=yy(k)+( yy(k+1)-yy(k) )*tt
890 tlist(npoint)=0 ! on fireline
895 ! make the list circular
896 tlist(npoint+1)=tlist(1)
897 llist(npoint+1)=llist(1)
898 xylist(npoint+1,1)=xylist(1,1)
899 xylist(npoint+1,2)=xylist(1,2)
900 !* least squares, A matrix for points
902 mA(kk,1)=xylist(kk,1)
903 mA(kk,2)=xylist(kk,2)
905 vecD(kk)=tlist(kk) ! D vector,time from ignition
910 mB(kk,ll)=0.0 ! clear
915 mb(kk,kk) = 10 ! large enough
917 if ( kk .ne. ll ) then
918 dist = sqrt( (xylist(kk,1)-xylist(ll,1))**2+ &
919 (xylist(kk,2)-xylist(ll,2))**2 )
920 mB(kk,kk)=min( mB(kk,kk) , dist )
923 mB(kk,kk)=mB(kk,kk)+1.
926 n=npoint ! rows of matrix A and B
927 m=3 ! columns of matrix A
928 p=npoint ! columns of matrix B
929 !* call least squqres in LAPACK
930 call SGGGLM(N,M,P,mA,LDA,mB,LDB,vecD,vecX,vecY, &
932 !print *,'after LS in case3'
933 !print *,'vecX from LS',vecX
934 !print *,'tign inputed',tign00,tign10,tign11,tign01
935 rnorm=snrm2(p,vecY,1)
939 ! rotate to gradient on x only
940 nr = sqrt(u(1)**2+u(2)**2)
941 if(.not.nr.gt.eps)then
951 ! mat vec multiplication
952 call matvec(mQ,2,2,u,3,ut,2,2,2)
953 errQ = ut(2) ! should be zero
954 ae = -ut(1)/fuel_time_cell
955 ce = -u(3)/fuel_time_cell
957 call matmatp(xylist,8,2,mQ,2,2,xytlist,8,2,npoint+1,2,2)
958 call sortxt( xytlist, 8,2, xt,8,npoint )
971 if ( xt2-xt1 > eps*100 ) then
979 if ( (xts>xt1 .and. xte>xt1) .or. &
980 (xts<xt2 .and. xte<xt2) ) then
984 aa = (yts-yte)/(xts-xte)
985 cc = (xts*yte-xte*yts)/(xts-xte)
997 end if!(xts>xt1 .and. xte>xt1)
999 ce=cet !use stored ce
1000 if (ae*xt1+ce > 0 ) then
1001 ce=ce-(ae*xt1+ce)!shift small amounts exp(-**)
1003 if (ae*xt2+ce > 0) then
1009 ! integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2
1010 ! numerically sound for ae->0, ae -> infty
1011 ! this can be important for different model scales
1012 ! esp. if someone runs the model in single precision!!
1013 ! s1=int((ah*x+ch),x,xt1,xt2)
1014 s1 = (xt2-xt1)*((1./2.)*ah*(xt2+xt1)+ch)
1015 ! s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)
1017 s2 = -ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1))
1018 ! s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)
1019 ! s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)
1020 ! expand in Taylor series around ae=0
1021 ! collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)
1022 ! =(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
1023 ! + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae
1024 ! + 1/2*xt1^2-1/2*xt2^2
1026 ! coefficient at ae^2 in the expansion, after some algebra
1027 a2=(xt1-xt2)*((1./4.)*(xt1+xt2)*ceae**2+(1./3.)* &
1028 (xt1**2+xt1*xt2+xt2**2)*ceae+(1./8.)* &
1029 (xt1**3+xt1*(xt2**2)+xt1**2*xt2+xt2**3))
1032 if (abs(d)>eps) then
1033 ! since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae
1034 ! for ae, ce -> 0 rounding error approx eps/ae^2
1035 s3=( exp(ae*(xt1+ceae))*(ae*xt1-1)-&
1036 exp(ae*(xt2+ceae))*(ae*xt2-1) )/(ae**2)
1038 !we do not worry about rounding as xt1 -> xt2, then s3 -> 0
1040 ! coefficient at ae^1 in the expansion
1041 a1=(xt1-xt2)*((1./2.)*ceae*(xt1+xt2)+(1./3.)*&
1042 (xt1**2+xt1*xt2+xt2**2))
1043 ! coefficient at ae^0 in the expansion for ae->0
1044 a0=(1./2.)*(xt1-xt2)*(xt1+xt2)
1045 s3=a0+a1*ae+a2*ae**2; ! approximate the integral
1050 out=1-out !fuel_left
1051 if(out<0 .or. out>1) then
1052 print *,':fuel_fraction should be between 0 and 1'
1053 !print *, 'eps= ', eps
1054 !print *, 'xt1= ', xt1, 'xt2= ', xt2
1055 !print *,'ae= ',ae,'ce= ',ce,'ah= ',ah,'ch= ',ch
1056 !print *,'a0= ', a0,'a1= ', a1,'a2= ', a2
1057 print *,'s1= ', s1,'s2= ', s2,'s3= ', s3
1058 print *,':fuel_fraction =',out
1064 end if ! if case0, elseif case4 ,else case123
1067 if(out>1. .or. out<0.)call crash('fuel_left_cell_2: fuel fraction out of bounds [0,1]')
1068 fuel_left_cell_2 = out
1069 end function fuel_left_cell_2
1072 !****************************************
1074 real function intexp(ab)
1080 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
1082 !eps = 2.2204*(10.0**(-8))!from matlab
1083 if ( eps < abs(ab)**3/6. ) then
1084 intexp=(exp(ab)-1)/ab
1090 !****************************************
1092 subroutine sortxt(xytlist,nrow,ncolumn,xt,nxt,nvec)
1094 integer::nrow,ncolumn,nxt,nvec
1095 real,dimension(nrow,ncolumn)::xytlist
1096 real,dimension(nxt)::xt
1107 if ( xt(i) > xt(j) ) then
1115 end subroutine !sortxt
1117 !****************************************
1119 subroutine matvec(A,m,n,V,nv,out,nout,nrow,ncolumn)
1121 integer::m,n,nv,nout,nrow,ncolumn
1122 real,dimension(m,n)::A ! allocated m by n
1123 real,dimension(nv)::V ! allocated nv
1124 real,dimension(nout)::out! allocated nout
1131 out(i)=out(i)+A(i,j)*V(j)
1136 !****************************************
1138 subroutine matmatp(A,mA,nA,B,mB,nB,C,mC,nC,nrow,ncolumn,nP)
1140 integer::mA,nA,mB,nB,mC,nC,nrow,ncolumn,nP
1141 real,dimension(mA,nA)::A ! allocated m by n
1142 real,dimension(mB,nB)::B ! allocated m by n
1143 real,dimension(mC,nC)::C ! allocated m by n
1149 C(i,j)=C(i,j)+A(i,k)*B(j,k) ! B'
1156 !****************************************
1160 subroutine prop_ls( id, & ! for debug
1161 ids,ide,jds,jde, & ! domain dims
1162 ims,ime,jms,jme, & ! memory dims
1163 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1164 its,ite,jts,jte, & ! tile dims
1165 ts,dt,dx,dy, & ! scalars in
1166 tbound, & ! scalars out
1167 lfn_in,lfn_out,tign,ros, & ! arrays inout
1172 !*** purpose: advance level function in time
1174 ! Jan Mandel August 2007 - February 2008
1178 ! Propagation of closed curve by a level function method. The level function
1179 ! lfn is defined by its values at the nodes of a rectangular grid.
1180 ! The area where lfn < 0 is inside the curve. The curve is
1181 ! described implicitly by lfn=0. Points where the curve intersects gridlines
1182 ! can be found by linear interpolation from nodes.
1184 ! The level function is advanced from time ts to time ts + dt.
1186 ! The level function should be initialized to (an approximation of) the signed
1187 ! distance from the curve. If the initial curve is a circle, the initial level
1188 ! function is simply the distance from the center minus the radius.
1190 ! The curve moves outside with speed given by function speed_func.
1192 ! Method: Godunov/ENO method for the normal motion. The timestep is checked for
1193 ! CFL condition. For a straight segment in a constant field and locally linear
1194 ! level function, the method reduces to the exact normal motion. The advantage of
1195 ! the level set method is that it treats automatically special cases such as
1196 ! the curve approaching itself and merging components of the area inside the curve.
1198 ! Based on S. Osher and R. Fedkiw, Level set methods and dynamic implicit surfaces,
1199 ! Springer, 2003, Sec. 6.4, as implemented in toolboxLS for Matlab by
1200 ! I. Mitchell, A toolbox of Level Set Methods (Version 1.1), TR-2007-11,
1201 ! Dept. Computer Science, University of British Columbia, 2007
1202 ! http://www.cs.ubc.ca/\~mitchell/Toolbo\LS
1207 ! id in unique identification for prints and dumps
1208 ! ids,ide,jds,jde in domain dimensions
1209 ! ims,ime,jms,jme in memory dimensions
1210 ! its,ite,jts,jte in tile dimensions
1213 ! dx,dy in grid spacing
1214 ! tbound out bound on stable time step from CFL condition, if tbound>=dt then OK
1215 ! lfn_in,lfn_out inout,out the level set function at nodes
1216 ! tign inout the ignition time at nodes
1218 ! The dimensions are cell-based, the nodal value is associated with the south-west corner.
1219 ! The whole computation is on domain indices ids:ide+1,jds:jde+1.
1221 ! The region where new lfn and tign are computed is the tile its:ite,jts:jte
1222 ! except when the tile is at domain upper boundary, an extra band of points is added:
1223 ! if ite=ide then region goes up to ite+1, if jte=jde then region goes up to jte+1.
1225 ! The time step requires values from 2 rows of nodes beyond the region except when at the
1226 ! domain boundary one-sided derivatives are used. This is implemented by extending the input
1227 ! beyond the domain boundary so sufficient memory bounds must be allocated.
1228 ! The update on all tiles can be done in parallel. To avoid the race condition (different regions
1229 ! of the same array updated by different threads), the in and out versions of the
1230 ! arrays lft and tign are distinct. If the time step dt is larger
1231 ! that the returned tbound, the routine should be called again with timestep td<=tbound, and then
1232 ! having distinct inputs and outputs comes handy.
1239 integer,intent(in)::id,ims,ime,jms,jme,ids,ide,jds,jde,its,ite,jts,jte,ips,ipe,jps,jpe
1240 real,dimension(ims:ime,jms:jme),intent(inout)::lfn_in,tign
1241 real,dimension(ims:ime,jms:jme),intent(out)::lfn_out,ros
1242 real,intent(in)::dx,dy,ts,dt
1243 real,intent(out)::tbound
1244 type(fire_params),intent(in)::fp
1252 real,dimension(IMTS:IMTE,JMTS:JMTE):: tend, lfn1 ! region-sized with halo
1254 real::grad2,rr,tbound2,a,a1 ! a=0 euler, a=0.5 heun
1256 real::gradx,grady,aspeed,err,aerr,time_now
1257 integer::ihs,ihe,jhs,jhe
1258 integer::ihs2,ihe2,jhs2,jhe2
1259 integer::i,j,its1,ite1,jts1,jte1,k,kk,id1
1260 character(len=128)msg
1261 integer::nfirenodes,nfireline
1262 real::sum_err,min_err,max_err,sum_aerr,min_aerr,max_aerr
1265 integer,parameter :: mstep=1000, printl=1
1266 real, parameter:: zero=0.,one=1.,eps=epsilon(zero),tol=100*eps, &
1267 safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
1269 ! f90 intrinsic function
1271 intrinsic max,min,sqrt,nint,epsilon,tiny,huge
1275 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1276 write(msg,'(a8,i5,a6,i5,3(a1,i5))')'prop_ls:',id,' tile ',its,':',ite,',',jts,':',jte
1277 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1280 a=fire_back_weight ! from module_fr_sfire_util
1285 ihs2=max(its-2,ids) ! need lfn two beyond the tile but not outside the domain
1290 ihs=max(its-1,ids) ! compute tend one beyond the tile but not outside the domain
1296 call write_array_m(ihs,ihe,jhs,jhe,ims,ime,jms,jme,lfn_in,'lfn_in',id)
1299 ! check array dimensions
1300 call check_mesh_2dim(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme)
1301 call print_2d_stats(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme, &
1302 lfn_in,'prop_ls: lfn in')
1304 ! NOTE: tend_ls will extrapolate to one node strip at domain boundaries
1305 ! so that it can compute gradient at domain boundaries.
1306 ! To avoid copying, lfn_in is declared inout.
1307 ! At tile boundaries that are not domain boundaries values of lfn_in two nodes
1308 ! outside of the tile are needed.
1309 id1 = id ! for debug prints
1310 if(id1.ne.0)id1=id1+1000
1311 call tend_ls( id1, &
1312 ims,ime,jms,jme, & ! memory dims for lfn_in
1313 IMTS,IMTE,JMTS,JMTE, & ! memory dims for tend
1314 ids,ide,jds,jde, & ! domain dims - where lfn exists
1315 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1316 ihs,ihe,jhs,jhe, & ! where tend computed
1317 ims,ime,jms,jme, & ! memory dims for ros
1318 its,ite,jts,jte, & ! tile dims - where to set ros
1319 ts,dt,dx,dy, & ! scalars in
1320 lfn_in, & ! arrays in
1321 tbound, & ! scalars out
1322 tend, ros, & ! arrays out
1327 call write_array_m(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE,tend,'tend1',id)
1330 ! Euler method, the half-step, same region as ted
1333 lfn1(i,j) = lfn_in(i,j) + dt*tend(i,j)
1337 call print_2d_stats(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE, &
1338 lfn1,'prop_ls: lfn1')
1339 ! tend = F(lfn1) on the tile (not beyond)
1341 if(id1.ne.0)id1=id1+1000
1343 IMTS,IMTE,JMTS,JMTE, & ! memory dims for lfn
1344 IMTS,IMTE,JMTS,JMTE, & ! memory dims for tend
1345 ids,ide,jds,jde, & ! domain dims - where lfn exists
1346 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1347 its,ite,jts,jte, & ! tile dims - where is tend computed
1348 ims,ime,jms,jme, & ! memory dims for ros
1349 its,ite,jts,jte, & ! tile dims - where is ros set
1350 ts+dt,dt,dx,dy, & ! scalars in
1352 tbound2, & ! scalars out
1353 tend,ros, & ! arrays out
1358 call write_array_m(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'tend2',id)
1361 call print_2d_stats(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'prop_ls: tend2')
1363 tbound=min(tbound,tbound2)
1365 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1366 write(msg,'(a,f10.2,4(a,f7.2))')'prop_ls: time',ts,' dt=',dt,' bound',min(tbound,999.99), &
1368 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1371 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1372 write(msg,'(2(a,f10.2))')'prop_ls: WARNING: time step ',dt, &
1374 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1378 ! combine lfn1 and lfn_in + dt*tend -> lfn_out
1382 lfn_out(i,j) = a1*lfn1(i,j) + a*(lfn_in(i,j) + dt*tend(i,j))
1386 ! compute ignition time by interpolation
1387 ! the node was not burning at start but it is burning at end
1388 ! interpolate from the level functions at start and at end
1389 ! lfn_in is the level set function value at time ts
1390 ! lfn_out is the level set function value at time ts+dt
1391 ! 0 is the level set function value at time tign(i,j)
1392 ! thus assuming the level function is approximately linear =>
1393 ! tign(i,j)= ts + ((ts + td) - ts) * lfn_in / (lfn_in - lfn_out)
1394 ! = ts + dt * lfn_in / (lfn_in - lfn_out)
1397 time_now = time_now + abs(time_now)*epsilon(time_now)*2.
1400 ! interpolate the cross-over time
1401 if (.not. lfn_out(i,j)>0 .and. lfn_in(i,j)>0)then
1402 tign(i,j) = ts + dt * lfn_in(i,j) / (lfn_in(i,j) - lfn_out(i,j))
1404 ! set the ignition time outside of burning region
1405 if(lfn_out(i,j)>0.)tign(i,j)=time_now
1409 ! check local speed error and stats
1410 ! may not work correctly in parallel
1424 ! loop over right inside of the domain
1425 ! cannot use values outside of the domain, would have to reflect and that
1426 ! would change lfn_out
1427 ! cannot use values outside of tile, not synchronized yet
1428 ! so in parallel mode the statistics is not accurate
1431 if(lfn_out(i,j)>0.0)then ! a point out of burning region
1432 if(lfn_out(i+1,j)<=0.or.lfn_out(i,j+1)<=0.or. & ! neighbor in burning region
1433 lfn_out(i-1,j)<=0.or.lfn_out(i,j-1)<=0)then ! point next to fireline
1434 gradx=(lfn_out(i+1,j)-lfn_out(i-1,j))/(2.0*dx) ! central differences
1435 grady=(lfn_out(i,j+1)-lfn_out(i,j-1))/(2.0*dy)
1436 grad2=sqrt(gradx*gradx+grady*grady)
1437 aspeed = (lfn_in(i,j)-lfn_out(i,j))/(dt*max(grad2,rmin))
1438 rr = speed_func(gradx,grady,dx,dy,i,j,fp)
1441 min_err=min(min_err,err)
1442 max_err=max(max_err,err)
1444 sum_aerr=sum_aerr+aerr
1445 min_aerr=min(min_aerr,aerr)
1446 max_aerr=max(max_aerr,aerr)
1447 nfireline=nfireline+1
1450 nfirenodes=nfirenodes+1
1454 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1455 write(msg,'(2(a,i6,f8.4))')'prop_ls: nodes burning',nfirenodes, &
1456 (100.*nfirenodes)/((ite1-its1+1)*(jte1-jts1+1)),'% next to fireline',nfireline
1457 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1460 call print_stat_line('speed error',its1,ite1,jts1,jte1,min_err,max_err,sum_err/nfireline)
1461 call print_stat_line('abs(speed error)',its1,ite1,jts1,jte1,min_aerr,max_aerr,sum_aerr/nfireline)
1464 ! check if the fire did not get to the domain boundary
1466 ! do kk=1,(boundary_guard*(ide-ids+1))/100 ! in %
1467 do kk=1,boundary_guard ! measured in cells
1469 if(i.ge.its.and.i.le.ite)then
1471 if(lfn_out(i,j)<=0.)goto 9
1475 ! do kk=1,(boundary_guard*(jde-jds+1))/100
1476 do kk=1,boundary_guard ! measured in cells
1478 if(j.ge.jts.and.j.le.jte)then
1480 if(lfn_out(i,j)<=0.)goto 9
1487 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1488 write(msg,'(a,i2,a,2i8)')'prop_ls: fire',boundary_guard, &
1489 ' cells from domain boundary at node ',i,j
1490 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1492 call crash('prop_ls: increase the fire region')
1495 call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme, &
1496 lfn_out,'prop_ls: lfn out')
1498 end subroutine prop_ls
1501 !*****************************
1504 subroutine tend_ls( id, &
1505 lims,lime,ljms,ljme, & ! memory dims for lfn
1506 tims,time,tjms,tjme, & ! memory dims for tend
1507 ids,ide,jds,jde, & ! domain - nodes where lfn defined
1508 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1509 ints,inte,jnts,jnte, & ! region - nodes where tend computed
1510 ims,ime,jms,jme, & ! memory dims for ros
1511 its,ite,jts,jte, & ! tile dims - where is ros set
1512 t,dt,dx,dy, & ! scalars in
1514 tbound, & ! scalars out
1515 tend, ros, & ! arrays out
1521 ! compute the right hand side of the level set equation
1524 integer,intent(in)::id,lims,lime,ljms,ljme,tims,time,tjms,tjme
1525 integer,intent(in)::ims,ime,jms,jme,its,ite,jts,jte
1526 integer, intent(in)::ids,ide,jds,jde,ints,inte,jnts,jnte,ips,ipe,jps,jpe
1527 real,intent(in)::t ! time
1528 real,intent(in)::dt,dx,dy ! mesh step
1529 real,dimension(lims:lime,ljms:ljme),intent(inout)::lfn ! level set function
1530 real,intent(out)::tbound ! max allowed time step
1531 real,dimension(tims:time,tjms:tjme),intent(out)::tend ! tendency (rhs of the level set pde)
1532 real,dimension(ims:ime,jms:jme),intent(out)::ros ! rate of spread
1533 type(fire_params),intent(in)::fp
1536 real:: te,diffLx,diffLy,diffRx,diffRy, &
1537 diffCx,diffCy,diff2x,diff2y,grad,rr, &
1538 ros_back,ros_wind,ros_slope,advx,advy,scale,nvx,nvy,speed,tanphi
1539 integer::i,j,itso,iteo,jtso,jteo
1540 character(len=128)msg
1543 real, parameter:: eps=epsilon(0.0)
1545 real, parameter:: zero=0.,one=1.,tol=100*eps, &
1546 safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
1549 ! f90 intrinsic function
1551 intrinsic max,min,sqrt,nint,tiny,huge
1555 real,dimension(tims:time,tjms:tjme)::rra,grada,speeda,tanphia
1560 ! check array dimensions
1561 call check_mesh_2dim(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme)
1562 call check_mesh_2dim(ints,inte,jnts,jnte,tims,time,tjms,tjme)
1564 call continue_at_boundary(1,1,fire_lfn_ext_up, & !extend by extrapolation but never down
1565 lims,lime,ljms,ljme, & ! memory dims
1566 ids,ide,jds,jde, & ! domain - nodes where lfn defined
1567 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1568 ints,inte,jnts,jnte, & ! tile - nodes update by this thread
1569 itso,iteo,jtso,jteo, & ! where set now
1572 call print_2d_stats(itso,iteo,jtso,jteo,lims,lime,ljms,ljme, &
1573 lfn,'tend_ls: lfn cont')
1576 call write_array_m(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme,lfn,'tend_lfn_in',id)
1582 ! one sided differences
1583 diffRx = (lfn(i+1,j)-lfn(i,j))/dx
1584 diffLx = (lfn(i,j)-lfn(i-1,j))/dx
1585 diffRy = (lfn(i,j+1)-lfn(i,j))/dy
1586 diffLy = (lfn(i,j)-lfn(i,j-1))/dy
1587 diffCx = diffLx+diffRx ! TWICE CENTRAL DIFFERENCE
1588 diffCy = diffLy+diffRy
1590 !upwinding - select right or left derivative
1591 select case(fire_upwinding)
1593 grad=sqrt(diffCx**2 + diffCy**2)
1595 diff2x=select_upwind(diffLx,diffRx)
1596 diff2y=select_upwind(diffLy,diffRy)
1597 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1598 case(2) ! godunov per osher/fedkiw
1599 diff2x=select_godunov(diffLx,diffRx)
1600 diff2y=select_godunov(diffLy,diffRy)
1601 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1603 diff2x=select_eno(diffLx,diffRx)
1604 diff2y=select_eno(diffLy,diffRy)
1605 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1606 case(4) ! Sethian - twice stronger pushdown of bumps
1607 grad=sqrt(max(diffLx,0.)**2+min(diffRx,0.)**2 &
1608 + max(diffLy,0.)**2+min(diffRy,0.)**2)
1613 ! normal direction, from central differences
1614 scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps)
1618 ! wind speed in direction of spread
1619 ! speed = vx(i,j)*nvx + vy(i,j)*nvy
1622 ! get rate of spread from wind speed and slope
1624 call fire_ros(ros_back,ros_wind,ros_slope, &
1627 rr=ros_back + ros_wind + ros_slope
1628 if(fire_grows_only.gt.0)rr=max(rr,0.)
1630 ! set ros for output
1631 if(i.ge.its.and.i.le.ite.and.j.ge.jts.and.j.le.jte)ros(i,j)=rr
1633 if(fire_upwind_split.eq.0)then
1635 ! get rate of spread
1636 te = -rr*grad ! normal term
1640 ! normal direction backing rate only
1641 te = - ros_back*grad
1643 ! advection in wind direction
1644 if (abs(speed)> eps) then
1645 advx=fp%vx(i,j)*ros_wind/speed
1646 advy=fp%vy(i,j)*ros_wind/speed
1652 ! tanphi = dzdxf*nvx + dzdyf*nvy
1653 ! advection from slope direction
1654 if(abs(tanphi)>eps) then
1655 advx=advx+fp%dzdxf(i,j)*ros_slope/tanphi
1656 advy=advy+fp%dzdyf(i,j)*ros_slope/tanphi
1659 if(fire_upwind_split.eq.1)then
1661 ! one-sided upwinding
1662 te = te - max(advx,0.)*diffLx - min(advx,0.)*diffRy &
1663 - max(advy,0.)*diffLy - min(advy,0.)*diffRy
1666 elseif(fire_upwind_split.eq.2)then
1669 call crash('prop_ls: bad fire_upwind_split, Lax-Friedrichs not done yet')
1673 call crash('prop_ls: bad fire_upwind_split')
1680 tbound = max(tbound,rr*(abs(diff2x)/dx+abs(diff2y)/dy)/grad)
1683 ! add numerical viscosity
1684 te=te + fire_viscosity*abs(rr)*((diffRx-diffLx)+(diffRy-diffLy))
1693 !write(msg,*)i,j,grad,dzdx,dzdy
1696 !if(abs(te)>1e4)then
1697 ! write(msg,'(a,2i5,e12.3)')'tend_ls: tend out of bounds at ',i,j,te
1704 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,rra,'rr',id)
1705 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,grada,'grad',id)
1706 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,speeda,'speed',id)
1707 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,tanphia,'tanphi',id)
1708 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,tend,'tend',id)
1711 call print_2d_stats(ints,inte,jnts,jnte,tims,time,tjms,tjme, &
1712 tend,'tend_ls: tend out')
1714 ! the final CFL bound
1715 tbound = 1/(tbound+tol)
1717 end subroutine tend_ls
1720 !**************************
1723 real function select_upwind(diffLx,diffRx)
1725 real, intent(in):: diffLx, diffRx
1728 ! upwind differences, L or R if bith same sign, otherwise zero
1731 if (diffLx>0.and.diffRx>0.)diff2x=diffLx
1732 if (diffLx<0.and.diffRx<0.)diff2x=diffRx
1734 select_upwind=diff2x
1735 end function select_upwind
1739 !**************************
1743 real function select_godunov(diffLx,diffRx)
1745 real, intent(in):: diffLx, diffRx
1748 ! Godunov scheme: upwind differences, L or R or none
1749 ! always test on > or < never = , much faster because of IEEE
1750 ! central diff >= 0 => take left diff if >0, ortherwise 0
1751 ! central diff <= 0 => take right diff if <0, ortherwise 0
1754 diffCx=diffRx+diffLx
1755 if (diffLx>0.and..not.diffCx<0)diff2x=diffLx
1756 if (diffRx<0.and. diffCx<0)diff2x=diffRx
1758 select_godunov=diff2x
1759 end function select_godunov
1762 !**************************
1765 real function select_eno(diffLx,diffRx)
1767 real, intent(in):: diffLx, diffRx
1770 ! 1st order ENO scheme
1772 if (.not.diffLx>0 .and. .not.diffRx>0)then
1774 elseif(.not.diffLx<0 .and. .not.diffRx<0)then
1776 elseif(.not.diffLx<0 .and. .not.diffRx>0)then
1777 if(.not. abs(diffRx) < abs(diffLx))then
1787 end function select_eno
1790 !**************************
1793 real function speed_func(diffCx,diffCy,dx,dy,i,j,fp)
1795 ! the level set method speed function
1798 real, intent(in)::diffCx,diffCy ! x and y coordinates of the direction of propagation
1799 real, intent(in)::dx,dy ! x and y coordinates of the direction of propagation
1800 integer, intent(in)::i,j ! indices of the node to compute the speed at
1801 type(fire_params),intent(in)::fp
1803 real::scale,nvx,nvy,r
1804 real::ros_back , ros_wind , ros_slope
1805 real, parameter:: eps=epsilon(0.0)
1807 ! normal direction, from central differences
1808 scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps)
1812 ! get rate of spread from wind speed and slope
1814 call fire_ros(ros_back,ros_wind,ros_slope, &
1817 r=ros_back + ros_wind + ros_slope
1818 if(fire_grows_only.gt.0)r=max(r,0.)
1821 end function speed_func
1823 end module module_fr_sfire_core