| blob | 474ba806770e0d0117c88d60c54af2938b08a55d |
1 !
2 #define DEBUG_OUT
3 #define DEBUG_PRINT
4 !#define FUEL_LEFT
5 !#define DEBUG_OUT_FUEL_LEFT
7 module module_fr_fire_core
9 use module_fr_fire_phys, only: fire_params , fire_ros
10 use module_fr_fire_util
12 ! The mathematical core of the fire spread model. No physical constants here.
13 !
14 ! subroutine fire_core: only this routine should be called from the outside.
15 ! subroutine fuel_left: compute remaining fuel from time of ignition.
16 ! subroutine prop_ls: propagation of curve in normal direction.
18 ! describe one ignition line
19 type ignition_line_type
20 REAL ros, & ! subscale rate of spread during the ignition process
21 stop_time, & ! when the ignition process stops from ignition start (s)
22 wind_red, & ! wind reduction factor at the ignition line
23 wrdist, & ! distance from the ignition line when the wind reduction stops
24 wrupwind, & ! use distance interpolated between 0. = nearest 1. = upwind
25 start_x, & ! x coordinate of the ignition line start point (m, or long/lat)
26 start_y, & ! y coordinate of the ignition line start point
27 end_x, & ! x coordinate of the ignition line end point
28 end_y, & ! y coordinate of the ignition line end point
29 start_time, & ! ignition time for the start point from simulation start (s)
30 end_time, & ! ignition time for the end poin from simulation start (s)
31 radius ! all within this radius ignites immediately
32 end type ignition_line_type
34 contains
36 !
37 !****************************************
38 !
40 subroutine init_no_fire(&
41 ifds,ifde,jfds,jfde, &
42 ifms,ifme,jfms,jfme, &
43 ifts,ifte,jfts,jfte, &
44 fdx,fdy,time_now, & ! scalars in
45 fuel_frac,fire_area,lfn,tign) ! arrays out
46 implicit none
48 !*** purpose: initialize model to no fire
50 !*** arguments
51 integer, intent(in):: ifds,ifde,jfds,jfde ! fire domain bounds
52 integer, intent(in):: ifts,ifte,jfts,jfte ! fire tile bounds
53 integer, intent(in):: ifms,ifme,jfms,jfme ! array bounds
54 real, intent(in) :: fdx,fdy,time_now ! mesh spacing, time
55 real, intent(out), dimension (ifms:ifme,jfms:jfme) :: &
56 fuel_frac,fire_area,lfn,tign ! model state
58 !*** calls
59 intrinsic epsilon
61 !*** local
62 integer:: i,j
63 real lfn_init,time_init
65 lfn_init = 2*max((ifde-ifds+1)*fdx,(jfde-jfds+1)*fdy) ! more than domain diameter
66 time_init=time_now + max(time_now,1.0)*epsilon(time_now) ! a bit in future
68 do j=jfts,jfte
69 do i=ifts,ifte
70 fuel_frac(i,j)=1. ! fuel at start is 1 by definition
71 fire_area(i,j)=0. ! nothing burning
72 tign(i,j) = time_init ! ignition in future
73 lfn(i,j) = lfn_init ! no fire
74 enddo
75 enddo
76 call message('init_no_fire: state set to no fire')
78 end subroutine init_no_fire
80 !
81 !******************
82 !
85 subroutine ignite_fire( ifds,ifde,jfds,jfde, & ! fire domain dims - the whole domain
86 ifms,ifme,jfms,jfme, &
87 ifts,ifte,jfts,jfte, &
88 ignition_line, &
89 start_ts,end_ts, &
90 coord_xf,coord_yf, &
91 unit_xf,unit_yf, &
92 lfn,tign,ignited)
93 implicit none
95 !*** purpose: ignite a circular/line/prescribed fire
97 !*** description
98 ! ignite fire in the region within radius r from the line (sx,sy) to (ex,ey).
99 ! the coordinates of nodes are given as the arrays coord_xf and coord_yf
100 ! r is given in m
101 ! one unit of coord_xf is unit_xf m
102 ! one unit of coord_yf is unit_yf m
103 ! so a node (i,j) will be ignited iff for some (x,y) on the line
104 ! || ( (coord_xf(i,j) - x)*unit_xf , (coord_yf(i,j) - y)*unit_yf ) || <= r
107 !*** arguments
108 integer, intent(in):: ifds,ifde,jfds,jfde ! fire domain bounds
109 integer, intent(in):: ifts,ifte,jfts,jfte ! fire tile bounds
110 integer, intent(in):: ifms,ifme,jfms,jfme ! array bounds
111 type(ignition_line_type), intent(in):: ignition_line ! description of the ignition line
112 real, intent(in):: start_ts,end_ts ! the time step start and end
113 real, dimension(ifms:ifme, jfms:jfme), intent(in):: &
114 coord_xf,coord_yf ! node coordinates
115 real, intent(in):: unit_xf,unit_yf ! coordinate units in m
116 real, intent(inout), dimension (ifms:ifme,jfms:jfme) :: &
117 lfn, tign ! level function, ignition time (state)
118 integer, intent(out):: ignited ! number of nodes newly ignited
120 !*** local
121 integer:: i,j
122 real::lfn_new,time_ign,ax,ay,rels,rele,d
123 real:: sx,sy ! start of ignition line, from lower left corner
124 real:: ex,ey ! end of ignition line, or zero
125 real:: st,et ! start and end of time of the ignition line
126 character(len=128):: msg
127 real::cx2,cy2,dmax,axmin,axmax,aymin,aymax,dmin
128 real:: start_x,start_y ! start of ignition line, from lower left corner
129 real:: end_x,end_y ! end of ignition line, or zero
130 real:: radius ! all within the radius of the line will ignite
131 real:: start_time,end_time ! the ignition time for the start and the end of the line
132 real:: ros,tos ! ignition rate and time of spread
134 !*** executable
136 ! copy ignition line fields to local variables
137 start_x = ignition_line%start_x ! x coordinate of the ignition line start point (m, or long/lat)
138 start_y = ignition_line%start_y ! y coordinate of the ignition line start point
139 end_x = ignition_line%end_x ! x coordinate of the ignition line end point
140 end_y = ignition_line%end_y ! y coordinate of the ignition line end point
141 start_time = ignition_line%start_time ! ignition time for the start point from simulation start (s)
142 end_time = ignition_line%end_time! ignition time for the end poin from simulation start (s)
143 radius = ignition_line%radius ! all within this radius ignites immediately
144 ros = ignition_line%ros ! rate of spread
145 tos = radius/ros ! time of spread to the given radius
146 st = start_time ! the start time of ignition considered in this time step
147 et = min(end_ts,end_time) ! the end time of the ignition segment in this time step
149 ! this should be called whenever (start_ts, end_ts) \subset (start_time, end_time + tos)
150 if(start_ts>et+tos .or. end_ts<st)return ! too late or too early, nothing to do
151 !print *, 'IGNITING'
152 if(start_time < end_time)then ! we really want to test start_time .ne. end_time, but avoiding test of floats on equality
153 ! segment of nonzero length
154 !rels = (st - start_time) / (end_time - start_time) ! relative position of st in the segment (start,end)
155 !sx = start_x + rels * (end_x - start_x)
156 !sy = start_y + rels * (end_y - start_y)
157 !paj rels = 0.
158 sx = start_x
159 sy = start_y
160 rele = (et - start_time) / (end_time - start_time) ! relative position of et in the segment (start,end)
161 ex = start_x + rele * (end_x - start_x)
162 ey = start_y + rele * (end_y - start_y)
163 else
164 ! just a point
165 !paj rels = 0.
166 !paj rele = 1.
167 sx = start_x
168 sy = start_y
169 ex = end_x
170 ey = end_y
171 endif
174 cx2=unit_xf*unit_xf
175 cy2=unit_yf*unit_yf
177 axmin=coord_xf(ifts,jfts)
178 aymin=coord_yf(ifts,jfts)
179 axmax=coord_xf(ifte,jfte)
180 aymax=coord_yf(ifte,jfte)
181 !$OMP CRITICAL(FIRE_CORE_CRIT)
182 write(msg,'(a,2g13.6,a,2g13.6)')'IGN from ',sx,sy,' to ',ex,ey
183 call message(msg)
184 write(msg,'(a,2f10.2,a,2f10.2,a)')'IGN timestep [',start_ts,end_ts,'] in [',start_time,end_time,']'
185 call message(msg)
186 write(msg,'(a,2g13.6,a,2g13.6)')'IGN tile coord from ',axmin,aymin,' to ',axmax,aymax
187 call message(msg)
188 !$OMP END CRITICAL(FIRE_CORE_CRIT)
189 ignited=0
190 dmax=0
191 dmin=huge(dmax)
192 11 format('IGN ',6(a,g17.7,1x))
193 12 format('IGN ',4(a,2g13.7,1x))
194 do j=jfts,jfte
195 do i=ifts,ifte
196 ax=coord_xf(i,j)
197 ay=coord_yf(i,j)
199 ! get d= distance from the nearest point on the ignition segment
200 ! and time_ign = the ignition time there
201 call nearest(d,time_ign,ax,ay,sx,sy,st,ex,ey,et,cx2,cy2)
202 dmax=max(d,dmax)
203 dmin=min(d,dmin)
205 lfn_new=d - min( radius, ros*(end_ts - time_ign) ) ! lft at end_ts
206 if(fire_print_msg.ge.3)then
207 !$OMP CRITICAL(FIRE_CORE_CRIT)
208 write(msg,*)'IGN1 i,j=',i,j,' lfn(i,j)=',lfn(i,j),' tign(i,j)=',tign(i,j)
209 call message(msg)
210 write(msg,*)'IGN2 i,j=',i,j,' lfn_new= ',lfn_new, ' time_ign= ',time_ign,' d=',d
211 call message(msg)
212 !$OMP END CRITICAL(FIRE_CORE_CRIT)
213 endif
214 if(.not.lfn_new>0.) then
215 ignited=ignited+1 ! count
216 endif
217 if(lfn(i,j)>0. .and. .not. lfn_new > 0.) then
218 tign(i,j)=time_ign + d/ros ! newly ignited now
219 if(fire_print_msg.ge.3)then
220 !$OMP CRITICAL(FIRE_CORE_CRIT)
221 write(msg,'(a,2i6,a,2g13.6,a,f10.2,a,2f10.2,a)')'IGN ignited cell ',i,j,' at',ax,ay, &
222 ' time',tign(i,j),' in [',start_ts,end_ts,']'
223 call message(msg)
224 write(msg,'(a,g10.3,a,f10.2,a,2f10.2,a)')'IGN distance',d,' from ignition line at',time_ign,' in [',st,et,']'
225 call message(msg)
226 !$OMP END CRITICAL(FIRE_CORE_CRIT)
227 endif
228 if(tign(i,j) < start_ts .or. tign(i,j) > end_ts )then
229 !$OMP CRITICAL(FIRE_CORE_CRIT)
230 write(msg,'(a,2i6,a,f11.6,a,2f11.6,a)')'WARNING ',i,j, &
231 ' fixing ignition time ',tign(i,j),' outside of the time step [',start_ts,end_ts,']'
232 call message (msg)
233 !$OMP END CRITICAL(FIRE_CORE_CRIT)
234 tign(i,j) = min(max(tign(i,j),start_ts),end_ts)
235 endif
236 endif
237 lfn(i,j)=min(lfn(i,j),lfn_new) ! update the level set function
238 if(fire_print_msg.ge.3)then
239 !$OMP CRITICAL(FIRE_CORE_CRIT)
240 write(msg,*)'IGN3 i,j=',i,j,' lfn(i,j)=',lfn(i,j),' tign(i,j)=',tign(i,j)
241 call message(msg)
242 !$OMP END CRITICAL(FIRE_CORE_CRIT)
243 endif
244 enddo
245 enddo
246 !$OMP CRITICAL(FIRE_CORE_CRIT)
247 write(msg,'(a,2g13.2,a,g10.2,a,g10.2)')'IGN units ',unit_xf,unit_yf,' m max dist ',dmax,' min',dmin
248 call message(msg)
249 write(msg,'(a,f6.1,a,f8.1,a,i10)')'IGN radius ',radius,' time of spread',tos,' ignited nodes',ignited
250 call message(msg)
251 !$OMP END CRITICAL(FIRE_CORE_CRIT)
253 return
254 99 continue
256 end subroutine ignite_fire
258 ! called from the inside of a loop, inline if possible
259 !DEC$ ATTRIBUTES FORCEINLINE
260 SUBROUTINE nearest(d,t,ax,ay,sx,sy,st,ex,ey,et,cx2,cy2)
261 implicit none
262 !*** arguments
263 real, intent(out):: d,t
264 real, intent(in):: ax,ay,sx,sy,st,ex,ey,et,cx2,cy2
265 ! input:
266 ! ax, ay coordinates of point a
267 ! sx,sy,ex,ey coordinates of endpoints of segment [x,y]
268 ! st,et values at the endpoints x,y
269 ! cx2,cy2 x and y unit squared for computing distance
271 ! output
272 ! d the distance of a and the nearest point z on the segment [x,y]
273 ! t linear interpolation from the values st,et to the point z
274 !
275 ! method: compute d as the distance (ax,ay) from the midpoint (mx,my)
276 ! minus a correction (because of rounding errors): |a-c|^2 = |a-m|^2 - |m-c|^2
277 ! when |m-c| >= |s-e|/2 the nearest point is one of the endpoints
278 ! the computation work also for the case when s=e exactly or approximately
279 !
280 !
281 ! a
282 ! /| \
283 ! s---m-c--e
284 !
285 ! |m-c| = |a-m| cos (a-m,e-s)
286 ! = |a-m| (a-m).(e-s))/(|a-m|*|e-s|)
287 !*** local
288 real:: mx,my,dam2,dames,am_es,cos2,dmc2,mcrel,mid_t,dif_t,des2,cx,cy
289 character(len=128):: msg
290 !*** executable
292 11 format('IGN ',6(a,g17.7,1x))
293 12 format('IGN ',4(a,2g13.7,1x))
295 ! midpoint m = (mx,my)
296 mx = (sx + ex)*0.5
297 my = (sy + ey)*0.5
298 dam2=(ax-mx)*(ax-mx)*cx2+(ay-my)*(ay-my)*cy2 ! |a-m|^2
299 des2 = (ex-sx)*(ex-sx)*cx2+(ey-sy)*(ey-sy)*cy2 ! des2 = |e-s|^2
300 dames = dam2*des2
301 am_es=(ax-mx)*(ex-sx)*cx2+(ay-my)*(ey-sy)*cy2 ! am_es = (a-m).(e-s)
302 if(dames>0)then
303 cos2 = (am_es*am_es)/dames ! cos2 = cos^2 (a-m,e-s)
304 else ! point a already is the midpoint
305 cos2 = 0.
306 endif
307 dmc2 = dam2*cos2 ! dmc2 = |m-c|^2
308 if(4.*dmc2 < des2)then ! if |m-c|<=|e-s|/2
309 ! d = sqrt(max(dam2 - dmc2,0.)) ! d=|a-m|^2 - |m-c|^2, guard rounding
310 mcrel = sign(sqrt(4.*dmc2/des2),am_es) ! relative distance of c from m
311 elseif(am_es>0)then ! if cos > 0, closest is e
312 mcrel = 1.0
313 else ! closest is s
314 mcrel = -1.0
315 endif
316 cx = (ex + sx)*0.5 + mcrel*(ex - sx)*0.5 ! interpolate to c by going from m
317 cy = (ey + sy)*0.5 + mcrel*(ey - sy)*0.5 ! interpolate to c by going from m
318 d=sqrt((ax-cx)*(ax-cx)*cx2+(ay-cy)*(ay-cy)*cy2) ! |a-c|^2
319 t = (et + st)*0.5 + mcrel*(et - st)*0.5 ! interpolate to c by going from m
320 if(fire_print_msg.ge.3)then
321 !$OMP CRITICAL(FIRE_CORE_CRIT)
322 write(msg,12)'find nearest to [',ax,ay,'] from [',sx,sy,'] [',ex,ey,']' ! DEB
323 call message(msg)
324 write(msg,12)'end times',st,et,' scale squared',cx2,cy2 ! DEB
325 call message(msg)
326 write(msg,11)'nearest at mcrel=',mcrel,'from the midpoint, t=',t ! DEB
327 call message(msg)
328 write(msg,12)'nearest is [',cx,cy,'] d=',d ! DEB
329 call message(msg)
330 write(msg,11)'dam2=',dam2,'des2=',des2,'dames=',dames
331 call message(msg)
332 write(msg,11)'am_es=',am_es,'cos2=',cos2,'dmc2=',dmc2 ! DEB
333 call message(msg)
334 !$OMP END CRITICAL(FIRE_CORE_CRIT)
335 endif
336 END SUBROUTINE nearest
339 !
340 !**********************
341 !
343 subroutine fuel_left(&
344 ims,ime,jms,jme, &
345 its,ite,jts,jte, &
346 ifs,ife,jfs,jfe, &
347 lfn, tign, fuel_time, time_now, fuel_frac, fire_area)
348 implicit none
350 !*** purpose: determine fraction of fuel remaining
351 !*** NOTE: because variables are cell centered, need halo/sync width 1 before
353 !*** arguments
355 integer, intent(in) :: its,ite,jts,jte,ims,ime,jms,jme,ifs,ife,jfs,jfe
356 real, intent(in), dimension(ims:ime,jms:jme)::lfn,tign,fuel_time
357 real, intent(in):: time_now
358 real, intent(out), dimension(ifs:ife,jfs:jfe)::fuel_frac
359 real, intent(out), dimension(ims:ime,jms:jme):: fire_area
361 ! ims,ime,jms,jme in memory dimensions
362 ! its,ite,jts,jte in tile dimensions (cells where fuel_frac computed)
363 ! ifs,ife,jfs,jfe in fuel_frac memory dimensions
364 ! lfn in level function, at nodes at midpoints of cells
365 ! tign in ignition time, at nodes at nodes at midpoints of cells
366 ! fuel_time in time constant of fuel, per cell
367 ! time_now in time now
368 ! fuel_frac out fraction of fuel remaining, per cell
369 ! fire_area out fraction of cell area on fire
371 !*** local
373 integer::i,j,ir,jr,icl,jcl,isubcl,jsubcl,i2,j2,ii,jj
374 real::fmax,frat,helpsum1,helpsum2,fuel_left_ff,fire_area_ff,rx,ry,tignf(2,2)
375 ! help variables instead of arrays fuel_left_f and fire_area_f
376 real::lffij,lffi1j,lffij1,lffi1j1,tifij,tifi1j,tifij1,tifi1j1,tx,ty,txx,tyy
377 ! variables for calculation instead of lff(i,j) and tif(i,j)is lffij,tifij etc..#define IFCELLS (ite-its+1)*fuel_left_irl
378 #define JFCELLS (jte-jts+1)*fuel_left_jrl
379 character(len=128)::msg
380 integer::m,omp_get_thread_num
383 ! refinement
384 ir=fuel_left_irl
385 jr=fuel_left_jrl
387 if ((ir.ne.2).or.(jr.ne.2)) then
388 call crash('fuel_left: ir.ne.2 or jr.ne.2 ')
389 endif
391 rx=1./ir
392 ry=1./jr
394 ! interpolate level set function to finer grid
395 #ifdef DEBUG_OUT_FUEL_LEFT
396 call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,lff,'lff',0)
397 call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,tif,'tif',0)
398 #endif
400 !
401 ! example for ir=2:
402 !
403 ! | coarse cell |
404 ! its-1 its ite ite+1
405 ! -------X------------|-----.-----X-----.-----|--........----|----------X----------|---------X
406 ! fine node 1 2 3 2*(ite-its+1)
407 ! fine cell 1 2 cell 2*(ite-its+1)
411 ! Loop over cells in Tile
412 ! Changes made by Volodymyr Kondratenko 09/24/2009
413 do icl=its,ite
414 do jcl=jts,jte
415 helpsum1=0
416 helpsum2=0
417 ! Loop over subcells in cell #(icl,jcl)
418 do isubcl=1,ir
419 do jsubcl=1,jr
420 i=(icl-its)*ir+isubcl
421 j=(jcl-jts)*jr+jsubcl
423 ! Direct calculation tif and lff, avoiding arrays, just for case ir=jr=2
424 if ((isubcl.eq.1).and.(jsubcl.eq.1)) then
425 i2=icl-1
426 j2=jcl-1
427 ty=0.5
428 tx=0.5
429 tyy=1.0
430 txx=1.0
431 else if ((isubcl.eq.2).and.(jsubcl.eq.1)) then
432 i2=icl
433 j2=jcl-1
434 ty=0.5
435 tx=0
436 tyy=1.0
437 txx=0.5
438 else if ((isubcl.eq.1).and.(jsubcl.eq.2)) then
439 i2=icl-1
440 j2=jcl
441 tx=0.5
442 ty=0
443 txx=1.0
444 tyy=0.5
445 else if ((isubcl.eq.2).and.(jsubcl.eq.2)) then
446 i2=icl
447 j2=jcl
448 tx=0
449 ty=0
450 txx=0.5
451 tyy=0.5
452 else
453 call crash('fuel_left: isubcl,jsubcl should be only 1 or 2')
454 endif
456 ! calculation lff and tif in 4 endpoints of subcell
457 lffij= &
458 (1-tx)*(1-ty)*lfn(i2,j2) &
459 + (1-tx)*ty *lfn(i2,j2+1) &
460 + tx*(1-ty)*lfn(i2+1,j2) &
461 + tx*ty *lfn(i2+1,j2+1)
462 lffi1j= &
463 (1-txx)*(1-ty)*lfn(i2,j2) &
464 + (1-txx)*ty *lfn(i2,j2+1) &
465 + (txx)*(1-ty)*lfn(i2+1,j2) &
466 + (txx)*ty *lfn(i2+1,j2+1)
467 lffij1= &
468 (1-tx)*(1-tyy)*lfn(i2,j2) &
469 + (1-tx)*(tyy) *lfn(i2,j2+1) &
470 + tx*(1-tyy)*lfn(i2+1,j2) &
471 + tx*(tyy) *lfn(i2+1,j2+1)
472 lffi1j1 = &
473 (1-txx)*(1-tyy)*lfn(i2,j2) &
474 + (1-txx)*(tyy) *lfn(i2,j2+1) &
475 + (txx)*(1-tyy)*lfn(i2+1,j2) &
476 + (txx)*(tyy) *lfn(i2+1,j2+1)
478 ! get ready to fix up tign values to be interpolated
479 do ii=1,2
480 do jj=1,2
481 tignf(ii,jj)=tign(i2+ii-1,j2+jj-1)
482 enddo
483 enddo
484 tifij= &
485 (1-tx)*(1-ty)*tignf(1,1) &
486 + (1-tx)*ty*tignf(1,1+1) &
487 + tx*(1-ty)*tignf(1+1,1) &
488 + tx*ty*tignf(1+1,1+1)
489 tifi1j= &
490 (1-txx)*(1-ty)*tignf(1,1) &
491 + (1-txx)*ty*tignf(1,1+1) &
492 + (txx)*(1-ty)*tignf(1+1,1) &
493 + (txx)*(ty)*tignf(1+1,1+1)
494 tifij1= &
495 (1-tx)*(1-tyy)*tignf(1,1) &
496 + (1-tx)*(tyy)*tignf(1,1+1) &
497 + tx*(1-tyy)*tignf(1+1,1) &
498 + tx*(tyy)*tignf(1+1,1+1)
499 tifi1j1= &
500 (1-txx)*(1-tyy)*tignf(1,1) &
501 + (1-txx)*(tyy)*tignf(1,1+1) &
502 + (txx)*(1-tyy)*tignf(1+1,1) &
503 + (txx)*(tyy)*tignf(1+1,1+1)
506 if(fuel_left_method.eq.1)then
507 call fuel_left_cell_1( fuel_left_ff, fire_area_ff, &
508 lffij,lffij1,lffi1j,lffi1j1,&
509 tifij,tifij1,tifi1j,tifi1j1,&
510 time_now, fuel_time(icl,jcl))
511 elseif(fuel_left_method.eq.2)then
512 fire_area_ff=0 ! initialize to something until computed in fuel_left_f(i,j)
513 fuel_left_ff=fuel_left_cell_2( &
514 lffij,lffij1,lffi1j,lffi1j1,&
515 tifij,tifij1,tifi1j,tifi1j1,&
516 time_now, fuel_time(icl,jcl))
517 else
518 call crash('fuel_left: unknown fuel_left_method')
519 endif
521 ! consistency check
522 if(fire_area_ff.lt.-1e-6 .or. &
523 (fire_area_ff.eq.0. .and. fuel_left_ff.lt.1.-1e-6))then
524 !$OMP CRITICAL(FIRE_CORE_CRIT)
525 write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
526 ' of refined mesh fuel burnt',1-fuel_left_ff,' fire area',fire_area_ff
527 !$OMP END CRITICAL(FIRE_CORE_CRIT)
528 call crash(msg)
529 endif
531 helpsum1=helpsum1+fuel_left_ff
532 helpsum2=helpsum2+fire_area_ff
533 enddo
534 enddo
535 fuel_frac(icl,jcl)=helpsum1
536 fire_area(icl,jcl)=helpsum2
537 enddo
538 enddo
543 #ifdef DEBUG_OUT_FUEL_LEFT
544 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fire_area,'fire_area',0)
545 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fuel_frac,'fuel_frac',0)
546 call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fuel_left_f,'fuel_left_f',0)
547 call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fire_area_f,'fire_area_f',0)
548 #endif
550 ! finish the averaging
551 do j=jts,jte
552 do i=its,ite
553 fuel_frac(i,j) = fuel_frac(i,j) /(ir*jr) ! multiply by weight for averaging over coarse cell
554 fire_area(i,j) = fire_area(i,j) /(ir*jr) !
555 enddo
556 enddo
558 ! consistency check after sum
559 fmax=0
560 do j=jts,jte
561 do i=its,ite
562 if(fire_area(i,j).eq.0.)then
563 if(fuel_frac(i,j).lt.1.-1e-6)then
564 !$OMP CRITICAL(FIRE_CORE_CRIT)
565 write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
566 ' fuel burnt',1-fuel_frac(i,j),' but fire area',fire_area(i,j)
567 !$OMP END CRITICAL(FIRE_CORE_CRIT)
568 call crash(msg)
569 endif
570 else
571 frat=(1-fuel_frac(i,j))/fire_area(i,j)
572 fmax=max(fmax,frat)
573 endif
574 enddo
575 enddo
576 !$OMP CRITICAL(FIRE_CORE_CRIT)
577 write(msg,'(a,4i6,a,f10.7)')'fuel_left: tile',its,ite,jts,jte,' max fuel burnt/area',fmax
578 !$OMP END CRITICAL(FIRE_CORE_CRIT)
579 call message(msg)
580 return
583 end subroutine fuel_left
585 !
586 !************************
587 !
589 subroutine fuel_left_cell_1( fuel_frac_left, fire_frac_area, &
590 lfn00,lfn01,lfn10,lfn11, &
591 tign00,tign01,tign10,tign11,&
592 time_now, fuel_time_cell)
593 !*** purpose: compute the fuel fraction left in one cell
594 implicit none
595 !*** arguments
596 real, intent(out):: fuel_frac_left, fire_frac_area !
597 real, intent(in)::lfn00,lfn01,lfn10,lfn11 ! level set function at 4 corners of the cell
598 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the 4 corners of the cell
599 real, intent(in)::time_now ! the time now
600 real, intent(in)::fuel_time_cell ! time to burns off to 1/e
601 !*** Description
602 ! The area burning is given by the condition L <= 0, where the function P is
603 ! interpolated from lfn(i,j)
604 !
605 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
606 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
607 ! when lfn(i,j)=0.
608 !
609 ! The function computes an approxmation of the integral
610 !
611 !
612 ! /\
613 ! |
614 ! fuel_frac_left = 1 - | 1 - exp(-T(x,y)/fuel_time_cell)) dxdy
615 ! |
616 ! \/
617 ! 0<x<1
618 ! 0<y<1
619 ! L(x,y)<=0
620 !
621 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
622 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
623 ! so dx=1 and dy=1 assumed.
624 !
625 ! Example:
626 !
627 ! lfn<0 lfn>0
628 ! (0,1) -----O--(1,1) O = points on the fireline, T=tnow
629 ! | \ | A = the burning area for computing
630 ! | \| fuel_frac(i,j)
631 ! | A O
632 ! | |
633 ! | |
634 ! (0,0)---------(1,0)
635 ! lfn<0 lfn<0
636 !
637 ! Approximations allowed:
638 ! The fireline can be approximated by straight line(s).
639 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
640 !
641 ! Requirements:
642 ! 1. The output should be a continuous function of the arrays lfn and
643 ! tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
644 ! 2. The output should be invariant to the symmetries of the input in each cell.
645 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
646 ! 4. The result should be at least 1st order accurate in the sense that it is
647 ! exact if the time from ignition is a linear function.
648 !
649 ! If time from ignition is approximated by polynomial in the burnt
650 ! region of the cell, this is integral of polynomial times exponential
651 ! over a polygon, which can be computed exactly.
652 !
653 ! Requirement 4 is particularly important when there is a significant decrease
654 ! of the fuel fraction behind the fireline on the mesh scale, because the
655 ! rate of fuel decrease right behind the fireline is much larger
656 ! (exponential...). This will happen when
657 !
658 ! change of time from ignition within one mesh cell / fuel_time_cell is not << 1
659 !
660 ! This is the same as
661 !
662 ! mesh cell size
663 ! X = ------------------------- is not << 1
664 ! fireline speed * fuel_time_cell
665 !
666 !
667 ! When X is large then the fuel burnt in one timestep in the cell is
668 ! approximately proportional to length of fireline in that cell.
669 !
670 ! When X is small then the fuel burnt in one timestep in the cell is
671 ! approximately proportional to the area of the burning region.
672 !
674 !*** calls
675 intrinsic tiny
677 !*** local
678 real::ps,aps,area,ta,out
679 real::t00,t01,t10,t11
680 real,parameter::safe=tiny(aps)
681 character(len=128)::msg
683 ! the following algorithm is a very crude approximation
685 ! minus time since ignition, 0 if no ignition yet
686 ! it is possible to have 0 in fire region when ignitin time falls in
687 ! inside the time step because lfn is updated at the beginning of the time step
689 t00=tign00-time_now
690 if(lfn00>0. .or. t00>0.)t00=0.
691 t01=tign01-time_now
692 if(lfn01>0. .or. t01>0.)t01=0.
693 t10=tign10-time_now
694 if(lfn10>0. .or. t10>0.)t10=0.
695 t11=tign11-time_now
696 if(lfn11>0. .or. t11>0.)t11=0.
698 ! approximate burning area, between 0 and 1
699 ps = lfn00+lfn01+lfn10+lfn11
700 aps = abs(lfn00)+abs(lfn01)+abs(lfn10)+abs(lfn11)
701 aps=max(aps,safe)
702 area =(-ps/aps+1.)/2.
703 area = max(area,0.) ! make sure area is between 0 and 1
704 area = min(area,1.)
706 ! average negative time since ignition
707 ta=0.25*(t00+t01+t10+t11)
709 ! exp decay in the burning area
710 out=1.
711 !if(area>0.)out=1. - area*(1. - exp(ta/fuel_time_cell))
712 if(area>0)out=area*exp(ta/fuel_time_cell) + (1. - area)
714 if(out>1.)then
715 !$OMP CRITICAL(FIRE_CORE_CRIT)
716 write(msg,*)'out=',out,'>1 area=',area,' ta=',ta
717 call message(msg)
718 write(msg,*)'tign=',tign00,tign01,tign10,tign11,' time_now=',time_now
719 !$OMP END CRITICAL(FIRE_CORE_CRIT)
720 call message(msg)
721 !call message('WARNING: fuel_left_cell_1: fuel fraction > 1')
722 call crash('fuel_left_cell_1: fuel fraction > 1')
723 endif
725 !out = max(out,0.) ! make sure out is between 0 and 1
726 !out = min(out,1.)
728 fuel_frac_left = out
729 fire_frac_area = area
731 end subroutine fuel_left_cell_1
733 !
734 !****************************************
735 !
737 real function fuel_left_cell_2( &
738 lfn00,lfn01,lfn10,lfn11, &
739 tign00,tign01,tign10,tign11,&
740 time_now, fuel_time_cell)
741 !*** purpose: compute the fuel fraction left in one cell
742 implicit none
743 !*** arguments
744 real, intent(in)::lfn00,lfn01,lfn10,lfn11 ! level set function at 4 corners of the cell
745 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the 4 corners of the cell
746 real, intent(in)::time_now ! the time now
747 real, intent(in)::fuel_time_cell ! time to burns off to 1/e
748 !*** Description
749 ! The area burning is given by the condition L <= 0, where the function P is
750 ! interpolated from lfn(i,j)
751 !
752 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
753 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
754 ! when lfn(i,j)=0.
755 !
756 ! The function computes an approxmation of the integral
757 !
758 !
759 ! /\
760 ! |
761 ! fuel_frac_left = 1 - | 1 - exp(-T(x,y)/fuel_time_cell)) dxdy
762 ! |
763 ! \/
764 ! 0<x<1
765 ! 0<y<1
766 ! L(x,y)<=0
767 !
768 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
769 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
770 ! so dx=1 and dy=1 assumed.
771 !
772 ! Example:
773 !
774 ! lfn<0 lfn>0
775 ! (0,1) -----O--(1,1) O = points on the fireline, T=tnow
776 ! | \ | A = the burning area for computing
777 ! | \| fuel_frac(i,j)
778 ! | A O
779 ! | |
780 ! | |
781 ! (0,0)---------(1,0)
782 ! lfn<0 lfn<0
783 !
784 ! Approximations allowed:
785 ! The fireline can be approximated by straight line(s).
786 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
787 !
788 ! Requirements:
789 ! 1. The output should be a continuous function of the arrays lfn and
790 ! tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
791 ! 2. The output should be invariant to the symmetries of the input in each cell.
792 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
793 ! 4. The result should be at least 1st order accurate in the sense that it is
794 ! exact if the time from ignition is a linear function.
795 !
796 ! If time from ignition is approximated by polynomial in the burnt
797 ! region of the cell, this is integral of polynomial times exponential
798 ! over a polygon, which can be computed exactly.
799 !
800 ! Requirement 4 is particularly important when there is a significant decrease
801 ! of the fuel fraction behind the fireline on the mesh scale, because the
802 ! rate of fuel decrease right behind the fireline is much larger
803 ! (exponential...). This will happen when
804 !
805 ! change of time from ignition within one mesh cell * fuel speed is not << 1
806 !
807 ! This is the same as
808 !
809 ! mesh cell size*fuel_speed
810 ! ------------------------- is not << 1
811 ! fireline speed
812 !
813 !
814 ! When X is large then the fuel burnt in one timestep in the cell is
815 ! approximately proportional to length of fireline in that cell.
816 !
817 ! When X is small then the fuel burnt in one timestep in the cell is
818 ! approximately proportional to the area of the burning region.
819 !
821 #ifndef FUEL_LEFT
822 call crash('fuel_left_cell_2: not implemented, please use fire_fuel_left_method=1')
823 fuel_left_cell_2=0. ! to avoid compiler warning about value not set
824 end function fuel_left_cell_2
825 #else
827 !*** calls
828 intrinsic tiny
830 !*** local
831 real::ps,aps,area,ta,out
832 real::t00,t01,t10,t11
833 real,parameter::safe=tiny(aps)
834 character(len=128)::msg
836 !*** local
837 integer::i,j,k
839 ! least squares
840 integer::mmax,nb,nmax,pmax,nin,nout
841 parameter(mmax=3,nb=64,nmax=8,pmax=8)
842 integer lda, ldb, lwork, info
843 parameter (lda=nmax, ldb=nmax, lwork=mmax+nmax+nb*(nmax+pmax))
844 integer n,m,p
845 real,dimension(lda,mmax):: mA
846 real,dimension(nmax):: vecD
847 real,dimension(lwork):: WORK
848 real,dimension(pmax):: vecY
849 real,dimension(mmax):: vecX
850 real,dimension(ldb,pmax)::mB
851 real,dimension(mmax)::u
853 real::tweight,tdist
854 integer::kk,ll,ss
855 real::rnorm
856 real,dimension(8,2)::xylist,xytlist
857 real,dimension(8)::tlist,llist,xt
858 real,dimension(5)::xx,yy
859 real,dimension(5)::lfn,tign
861 integer:: npoint
862 real::tt,x0,y0,xts,xte,yts,yte,xt1,xt2
863 real::lfn0,lfn1,dist,nr,c,s,errQ,ae,ce,ceae,a0,a1,a2,d,cet
864 real::s1,s2,s3
865 real::upper,lower,ah,ch,aa,cc,aupp,cupp,alow,clow
866 real,dimension(2,2)::mQ
867 real,dimension(2)::ut
869 !calls
870 intrinsic epsilon
872 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
875 ! external functions
876 real::snrm2
877 double precision :: dnrm2
878 external dnrm2
879 external snrm2
880 ! external subroutines
881 external sggglm
882 external dggglm
883 !executable statements
885 ! a very crude approximation - replace by a better code
886 ! call check_mesh_2dim(ids,ide+1,jds,jde+1,ims,ime,jms,jme)
887 t00=tign00-time_now
888 if(lfn00>=0. .or. t00>0.)t00=0.
889 t01=tign01-time_now
890 if(lfn01>=0. .or. t01>0.)t01=0.
891 t10=tign10-time_now
892 if(lfn10>=0. .or. t10>0.)t10=0.
893 t11=tign11-time_now
894 if(lfn11>=0. .or. t11>0.)t11=0.
896 !if (t00<0 .or. t01 <0 .or. t10<0 .or. t11<0) then
897 ! print *,'tign00=',tign00,'tign10=',tign10,&
898 ! 'tign01=',tign01,'tign11=',tign11
899 !end if
903 !*** case0 Do nothing
904 if ( lfn00>=0 .and. lfn10>=0 .and. lfn01>=0 .and. lfn11>=0 ) then
905 out = 1.0 ! fuel_left, no burning
906 !*** case4 all four coners are burning
907 else if (lfn00<=0 .and. lfn10<=0 .and. lfn01<=0 .and. lfn11<=0) then
908 ! least squares, A matrix for points
909 mA(1,1)=0.0
910 mA(2,1)=1.0
911 mA(3,1)=0.0
912 mA(4,1)=1.0
913 mA(1,2)=0.0
914 mA(2,2)=0.0
915 mA(3,2)=1.0
916 mA(4,2)=1.0
917 mA(1,3)=1.0
918 mA(2,3)=1.0
919 mA(3,3)=1.0
920 mA(4,3)=1.0
921 ! D vector, time from ignition
922 vecD(1)=time_now-tign00
923 vecD(2)=time_now-tign10
924 vecD(3)=time_now-tign01
925 vecD(4)=time_now-tign11
926 ! B matrix, weights
927 do kk=1,4
928 do ll=1,4
929 mB(kk,ll)=0.0
930 end do
931 mB(kk,kk)=2.0
932 end do
933 ! set the m,n,p
934 n=4 ! rows of matrix A and B
935 m=3 ! columns of matrix A
936 p=4 ! columns of matrix B
937 ! call least squqres in LAPACK
938 call SGGGLM(N,M,P,mA,LDA,mB,LDB,vecD,vecX,vecY, &
939 WORK,LWORK,INFO)
940 rnorm=snrm2(p,vecY,1)
941 ! integrate
942 u(1)=-vecX(1)/fuel_time_cell
943 u(2)=-vecX(2)/fuel_time_cell
944 u(3)=-vecX(3)/fuel_time_cell
945 !fuel_burn(i,j)=1-exp(u(3))*intexp(u(1)*dx)*intexp(u(2)*dy)
946 s1=u(1)
947 s2=u(2)
948 out=1-exp(u(3))*intexp(s1)*intexp(s2)
949 !print *,'intexp
950 if ( out<0 .or. out>1.0 ) then
951 print *,'case4, out should be between 0 and 1'
952 end if
953 !*** case 1,2,3
954 else
955 ! set xx, yy for the coner points
956 ! move these values out of i and j loop to speed up
957 xx(1) = -0.5
958 xx(2) = 0.5
959 xx(3) = 0.5
960 xx(4) = -0.5
961 xx(5) = -0.5
962 yy(1) = -0.5
963 yy(2) = -0.5
964 yy(3) = 0.5
965 yy(4) = 0.5
966 yy(5) = -0.5
967 lfn(1)=lfn00
968 lfn(2)=lfn10
969 lfn(3)=lfn11
970 lfn(4)=lfn01
971 lfn(5)=lfn00
972 tign(1)=t00
973 tign(2)=t10
974 tign(3)=t11
975 tign(4)=t01
976 tign(5)=t00
977 npoint = 0 ! number of points in polygon
978 !print *,'time_now=',time_now
979 !print *,'lfn00=',lfn00,'lfn10=',lfn10,&
980 ! 'lfn01=',lfn01,'lfn11=',lfn11
981 !print *,'t00=',t00,'t10=',t10,&
982 ! 't01=',t01,'t11=',t11
984 do k=1,4
985 lfn0=lfn(k )
986 lfn1=lfn(k+1)
987 if ( lfn0 <= 0.0 ) then
988 npoint = npoint + 1
989 xylist(npoint,1)=xx(k)
990 xylist(npoint,2)=yy(k)
991 tlist(npoint)=-tign(k)
992 llist(npoint)=lfn0
993 end if
994 if ( lfn0*lfn1 < 0 ) then
995 npoint = npoint + 1
996 tt=lfn0/(lfn0-lfn1)
997 x0=xx(k)+( xx(k+1)-xx(k) )*tt
998 y0=yy(k)+( yy(k+1)-yy(k) )*tt
999 xylist(npoint,1)=x0
1000 xylist(npoint,2)=y0
1001 tlist(npoint)=0 ! on fireline
1002 llist(npoint)=0
1003 end if
1004 end do
1006 ! make the list circular
1007 tlist(npoint+1)=tlist(1)
1008 llist(npoint+1)=llist(1)
1009 xylist(npoint+1,1)=xylist(1,1)
1010 xylist(npoint+1,2)=xylist(1,2)
1011 !* least squares, A matrix for points
1012 do kk=1,npoint
1013 mA(kk,1)=xylist(kk,1)
1014 mA(kk,2)=xylist(kk,2)
1015 mA(kk,3)=1.0
1016 vecD(kk)=tlist(kk) ! D vector,time from ignition
1017 end do
1018 ! B matrix, weights
1019 do kk=1,ldb
1020 do ll=1,pmax
1021 mB(kk,ll)=0.0 ! clear
1022 end do
1023 end do
1025 do kk=1,npoint
1026 mb(kk,kk) = 10 ! large enough
1027 do ll=1,npoint
1028 if ( kk .ne. ll ) then
1029 dist = sqrt( (xylist(kk,1)-xylist(ll,1))**2+ &
1030 (xylist(kk,2)-xylist(ll,2))**2 )
1031 mB(kk,kk)=min( mB(kk,kk) , dist )
1032 end if
1033 end do !ll
1034 mB(kk,kk)=mB(kk,kk)+1.
1035 end do ! kk
1036 ! set the m,n,p
1037 n=npoint ! rows of matrix A and B
1038 m=3 ! columns of matrix A
1039 p=npoint ! columns of matrix B
1040 !* call least squqres in LAPACK
1041 call SGGGLM(N,M,P,mA,LDA,mB,LDB,vecD,vecX,vecY, &
1042 WORK,LWORK,INFO)
1043 !print *,'after LS in case3'
1044 !print *,'vecX from LS',vecX
1045 !print *,'tign inputed',tign00,tign10,tign11,tign01
1046 rnorm=snrm2(p,vecY,1)
1047 u(1)=vecX(1)
1048 u(2)=vecX(2)
1049 u(3)=vecX(3)
1050 ! rotate to gradient on x only
1051 nr = sqrt(u(1)**2+u(2)**2)
1052 if(.not.nr.gt.eps)then
1053 out=1.
1054 goto 900
1055 endif
1056 c = u(1)/nr
1057 s = u(2)/nr
1058 mQ(1,1)=c
1059 mQ(1,2)=s
1060 mQ(2,1)=-s
1061 mQ(2,2)=c
1062 ! mat vec multiplication
1063 call matvec(mQ,2,2,u,3,ut,2,2,2)
1064 errQ = ut(2) ! should be zero
1065 ae = -ut(1)/fuel_time_cell
1066 ce = -u(3)/fuel_time_cell
1067 cet=ce!keep ce
1068 call matmatp(xylist,8,2,mQ,2,2,xytlist,8,2,npoint+1,2,2)
1069 call sortxt( xytlist, 8,2, xt,8,npoint )
1070 out=0.0
1071 aupp=0.0
1072 cupp=0.0
1073 alow=0.0
1074 clow=0.0
1075 do k=1,npoint-1
1076 xt1=xt(k)
1077 xt2=xt(k+1)
1078 upper=0
1079 lower=0
1080 ah=0
1081 ch=0
1082 if ( xt2-xt1 > eps*100 ) then
1084 do ss=1,npoint
1085 xts=xytlist(ss,1)
1086 yts=xytlist(ss,2)
1087 xte=xytlist(ss+1,1)
1088 yte=xytlist(ss+1,2)
1090 if ( (xts>xt1 .and. xte>xt1) .or. &
1091 (xts<xt2 .and. xte<xt2) ) then
1092 aa = 0 ! do nothing
1093 cc = 0
1094 else
1095 aa = (yts-yte)/(xts-xte)
1096 cc = (xts*yte-xte*yts)/(xts-xte)
1097 if (xte<xts) then
1098 aupp = aa
1099 cupp = cc
1100 ah=ah+aa
1101 ch=ch+cc
1102 upper=upper+1
1103 else
1104 alow = aa
1105 clow = cc
1106 lower=lower+1
1107 end if
1108 end if!(xts>xt1 .and. xte>xt1)
1109 end do ! ss
1110 ce=cet !use stored ce
1111 if (ae*xt1+ce > 0 ) then
1112 ce=ce-(ae*xt1+ce)!shift small amounts exp(-**)
1113 end if
1114 if (ae*xt2+ce > 0) then
1115 ce=ce-(ae*xt2+ce)
1116 end if
1118 ah = aupp-alow
1119 ch = cupp-clow
1120 ! integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2
1121 ! numerically sound for ae->0, ae -> infty
1122 ! this can be important for different model scales
1123 ! esp. if someone runs the model in single precision!!
1124 ! s1=int((ah*x+ch),x,xt1,xt2)
1125 s1 = (xt2-xt1)*((1./2.)*ah*(xt2+xt1)+ch)
1126 ! s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)
1127 ceae=ce/ae;
1128 s2 = -ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1))
1129 ! s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)
1130 ! s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)
1131 ! expand in Taylor series around ae=0
1132 ! collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)
1133 ! =(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
1134 ! + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae
1135 ! + 1/2*xt1^2-1/2*xt2^2
1136 !
1137 ! coefficient at ae^2 in the expansion, after some algebra
1138 a2=(xt1-xt2)*((1./4.)*(xt1+xt2)*ceae**2+(1./3.)* &
1139 (xt1**2+xt1*xt2+xt2**2)*ceae+(1./8.)* &
1140 (xt1**3+xt1*(xt2**2)+xt1**2*xt2+xt2**3))
1141 d=(ae**4)*a2
1143 if (abs(d)>eps) then
1144 ! since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae
1145 ! for ae, ce -> 0 rounding error approx eps/ae^2
1146 s3=( exp(ae*(xt1+ceae))*(ae*xt1-1)-&
1147 exp(ae*(xt2+ceae))*(ae*xt2-1) )/(ae**2)
1149 !we do not worry about rounding as xt1 -> xt2, then s3 -> 0
1150 else
1151 ! coefficient at ae^1 in the expansion
1152 a1=(xt1-xt2)*((1./2.)*ceae*(xt1+xt2)+(1./3.)*&
1153 (xt1**2+xt1*xt2+xt2**2))
1154 ! coefficient at ae^0 in the expansion for ae->0
1155 a0=(1./2.)*(xt1-xt2)*(xt1+xt2)
1156 s3=a0+a1*ae+a2*ae**2; ! approximate the integral
1157 end if
1159 s3=ah*s3
1160 out=out+s1+s2+s3
1161 out=1-out !fuel_left
1162 if(out<0 .or. out>1) then
1163 print *,':fuel_fraction should be between 0 and 1'
1164 !print *, 'eps= ', eps
1165 !print *, 'xt1= ', xt1, 'xt2= ', xt2
1166 !print *,'ae= ',ae,'ce= ',ce,'ah= ',ah,'ch= ',ch
1167 !print *,'a0= ', a0,'a1= ', a1,'a2= ', a2
1168 print *,'s1= ', s1,'s2= ', s2,'s3= ', s3
1169 print *,':fuel_fraction =',out
1170 end if!print
1172 end if
1173 end do ! k
1175 end if ! if case0, elseif case4 ,else case123
1177 900 continue
1178 if(out>1. .or. out<0.)call crash('fuel_left_cell_2: fuel fraction out of bounds [0,1]')
1179 fuel_left_cell_2 = out
1180 end function fuel_left_cell_2
1182 !
1183 !****************************************
1184 !
1185 real function intexp(ab)
1186 implicit none
1187 real::ab
1188 !calls
1189 intrinsic epsilon
1191 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
1193 !eps = 2.2204*(10.0**(-8))!from matlab
1194 if ( eps < abs(ab)**3/6. ) then
1195 intexp=(exp(ab)-1)/ab
1196 else
1197 intexp=1+ab/2.
1198 end if
1199 end function
1200 !
1201 !****************************************
1202 !
1203 subroutine sortxt(xytlist,nrow,ncolumn,xt,nxt,nvec)
1204 implicit none
1205 integer::nrow,ncolumn,nxt,nvec
1206 real,dimension(nrow,ncolumn)::xytlist
1207 real,dimension(nxt)::xt
1209 integer::i,j
1210 real::temp
1212 do i=1,nvec
1213 xt(i)=xytlist(i,1)
1214 end do
1216 do i=1,nvec-1
1217 do j=i+1,nvec
1218 if ( xt(i) > xt(j) ) then
1219 temp = xt(i)
1220 xt(i)=xt(j)
1221 xt(j)=temp
1222 end if
1223 end do
1224 end do
1226 end subroutine !sortxt
1227 !
1228 !****************************************
1229 !
1230 subroutine matvec(A,m,n,V,nv,out,nout,nrow,ncolumn)
1231 implicit none
1232 integer::m,n,nv,nout,nrow,ncolumn
1233 real,dimension(m,n)::A ! allocated m by n
1234 real,dimension(nv)::V ! allocated nv
1235 real,dimension(nout)::out! allocated nout
1237 integer::i,j
1239 do i=1,nrow
1240 out(i)=0.0
1241 do j=1,ncolumn
1242 out(i)=out(i)+A(i,j)*V(j)
1243 end do
1244 end do
1245 end subroutine
1246 !
1247 !****************************************
1248 !
1249 subroutine matmatp(A,mA,nA,B,mB,nB,C,mC,nC,nrow,ncolumn,nP)
1250 implicit none
1251 integer::mA,nA,mB,nB,mC,nC,nrow,ncolumn,nP
1252 real,dimension(mA,nA)::A ! allocated m by n
1253 real,dimension(mB,nB)::B ! allocated m by n
1254 real,dimension(mC,nC)::C ! allocated m by n
1255 integer::i,j,k
1256 do i=1,nrow
1257 do j=1,ncolumn
1258 C(i,j)=0.0
1259 do k=1,nP
1260 C(i,j)=C(i,j)+A(i,k)*B(j,k) ! B'
1261 end do
1262 end do
1263 end do
1264 end subroutine
1266 !
1267 !****************************************
1268 !
1269 #endif
1271 subroutine prop_ls_rk3(id, &
1272 ifds,ifde,jfds,jfde, &
1273 ifms,ifme,jfms,jfme, &
1274 ifps,ifpe,jfps,jfpe, &
1275 ifts,ifte,jfts,jfte, &
1276 ts,dt,dx,dy, &
1277 tbound, &
1278 lfn_in, &
1279 lfn_0,lfn_1,lfn_2, &
1280 lfn_out,tign,ros, &
1281 fp, &
1282 grid, &
1283 ids,ide,jds,jde,kds,kde, &
1284 ims,ime,jms,jme,kms,kme, &
1285 ips,ipe,jps,jpe,kps,kpe &
1286 )
1289 #ifdef DM_PARALLEL
1290 USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks
1291 USE module_comm_dm , ONLY : halo_fire_lfn_1_sub, halo_fire_lfn_2_sub, halo_fire_lfn_0_sub
1292 #endif
1293 USE module_domain , only: domain
1295 implicit none
1297 !*** purpose: advance level function in time
1299 !***************************************************************************************!
1300 !*** Level-set method recoded for true mpi communications and added ***!
1301 !*** 3rd-order Runge-Kutta time integration and 3rd/5th-order WENO advection schemes ***!
1302 !*** Implemented by: Domingo Munoz-Esparza (NCAR/RAL, April 2016) ***!
1303 !*** Reference: D. Munoz-Esparza, B. Kosovic, P. Jimenez, J. Coen: "An accurate ***!
1304 !*** fire-spread algorithm in the Weather Research and Forecasting model using the ***!
1305 !*** level-set method", Journal of Advances in Modeling Earth Systems, 2018 ***!
1306 !*** https://doi.org/10.1002/2017MS001108 ***!
1307 !***************************************************************************************!
1309 !*** description
1310 !
1311 ! Propagation of closed curve by a level function method. The level function
1312 ! lfn is defined by its values at the nodes of a rectangular grid.
1313 ! The area where lfn < 0 is inside the curve. The curve is
1314 ! described implicitly by lfn=0. Points where the curve intersects gridlines
1315 ! can be found by linear interpolation from nodes.
1316 !
1317 ! The level function is advanced from time ts to time ts + dt.
1318 !
1319 ! The level function should be initialized to (an approximation of) the signed
1320 ! distance from the curve. If the initial curve is a circle, the initial level
1321 ! function is simply the distance from the center minus the radius.
1322 !
1323 ! The curve moves outside with speed given by function speed_func.
1324 !
1325 ! Method: Godunov/ENO method for the normal motion. The timestep is checked for
1326 ! CFL condition. For a straight segment in a constant field and locally linear
1327 ! level function, the method reduces to the exact normal motion. The advantage of
1328 ! the level set method is that it treats automatically special cases such as
1329 ! the curve approaching itself and merging components of the area inside the curve.
1330 !
1331 ! Based on S. Osher and R. Fedkiw, Level set methods and dynamic implicit surfaces,
1332 ! Springer, 2003, Sec. 6.4, as implemented in toolboxLS for Matlab by
1333 ! I. Mitchell, A toolbox of Level Set Methods (Version 1.1), TR-2007-11,
1334 ! Dept. Computer Science, University of British Columbia, 2007
1335 ! http://www.cs.ubc.ca/\~mitchell/Toolbo\LS
1336 !
1338 !*** arguments
1340 ! id in unique identification for prints and dumps
1341 ! ids,ide,jds,jde in domain dimensions
1342 ! ims,ime,jms,jme in memory dimensions
1343 ! its,ite,jts,jte in tile dimensions
1344 ! ts in start time
1345 ! dt in time step
1346 ! dx,dy in grid spacing
1347 ! tbound out bound on stable time step from CFL condition, if tbound>=dt then OK
1348 ! lfn_in,lfn_out inout,out the level set function at nodes
1349 ! tign inout the ignition time at nodes
1351 !*** calls
1352 !
1353 ! DME added for halo update
1354 type(domain) , target :: grid
1355 integer, intent(in):: ids,ide,jds,jde,kds,kde, &
1356 ims,ime,jms,jme,kms,kme, &
1357 ips,ipe,jps,jpe,kps,kpe
1358 !
1360 integer,intent(in)::id,ifms,ifme,jfms,jfme,ifds,ifde,jfds,jfde,ifts,ifte,jfts,jfte,ifps,ifpe,jfps,jfpe
1361 real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_in,tign
1362 real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_1,lfn_2,lfn_0
1363 real,dimension(ifms:ifme,jfms:jfme),intent(out)::lfn_out,ros
1364 real,intent(in)::dx,dy,ts,dt
1365 real,intent(out)::tbound
1366 type(fire_params),intent(in)::fp
1368 real,dimension(ifms:ifme,jfms:jfme):: tend
1369 ! scalars
1370 real::grad2,rr,tbound2,tbound3
1372 real::gradx,grady,aspeed,err,aerr,time_now
1373 integer::i,j,its1,ite1,jts1,jte1,k,kk,id1
1374 character(len=128)::msg
1375 integer::nfirenodes,nfireline
1376 real::sum_err,min_err,max_err,sum_aerr,min_aerr,max_aerr
1378 ! constants
1379 integer,parameter :: mstep=1000, printl=1
1380 real, parameter:: zero=0.,one=1.,eps=epsilon(zero),tol=100*eps, &
1381 safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
1383 ! f90 intrinsic function
1385 intrinsic max,min,sqrt,nint,epsilon,tiny,huge
1387 !*** executable
1389 !$OMP CRITICAL(FIRE_CORE_CRIT)
1390 write(msg,'(a8,i5,a6,i5,3(a1,i5))')'prop_ls:',id,' tile ',ifts,':',ifte,',',jfts,':',jfte
1391 !$OMP END CRITICAL(FIRE_CORE_CRIT)
1392 call message(msg)
1394 !!!!!!!!!!!!!!!!!!!!!!
1395 ! Runge-Kutta step 1 !
1396 !!!!!!!!!!!!!!!!!!!!!!
1398 do j=jfts,jfte
1399 do i=ifts,ifte
1400 lfn_0(i,j) = lfn_in(i,j)
1401 enddo
1402 enddo
1404 #ifdef DM_PARALLEL
1405 # include "HALO_FIRE_LFN_0.inc"
1406 #endif
1408 id1 = id ! for debug prints
1409 if(id1.ne.0)id1=id1+1000
1410 call tend_ls(id1, &
1411 ifds,ifde,jfds,jfde, & ! domain dims
1412 ifps,ifpe,jfps,jfpe, & ! patch dims
1413 ifts,ifte,jfts,jfte, & ! tile dims
1414 ifms,ifme,jfms,jfme, & ! memory dims
1415 ts,dt,dx,dy, & ! scalars in
1416 lfn_0, & ! arrays in
1417 tbound, & ! scalars out
1418 tend, ros, & ! arrays out
1419 fp & ! params
1420 )
1422 do j=jfts,jfte
1423 do i=ifts,ifte
1424 lfn_1(i,j) = lfn_0(i,j) + (dt/3.0)*tend(i,j)
1425 enddo
1426 enddo
1428 #ifdef DM_PARALLEL
1429 # include "HALO_FIRE_LFN_1.inc"
1430 #endif
1432 !!!!!!!!!!!!!!!!!!!!!!
1433 ! Runge-Kutta step 2 !
1434 !!!!!!!!!!!!!!!!!!!!!!
1436 if(id1.ne.0)id1=id1+1000
1437 call tend_ls(id1, &
1438 ifds,ifde,jfds,jfde, &
1439 ifps,ifpe,jfps,jfpe, &
1440 ifts,ifte,jfts,jfte, &
1441 ifms,ifme,jfms,jfme, &
1442 ts+dt,dt,dx,dy, &
1443 lfn_1, &
1444 tbound2, &
1445 tend,ros, &
1446 fp &
1447 )
1449 do j=jfts,jfte
1450 do i=ifts,ifte
1451 lfn_2(i,j) = lfn_0(i,j) + (dt/2.0)*tend(i,j)
1452 enddo
1453 enddo
1455 #ifdef DM_PARALLEL
1456 # include "HALO_FIRE_LFN_2.inc"
1457 #endif
1459 !!!!!!!!!!!!!!!!!!!!!!
1460 ! Runge-Kutta step 3 !
1461 !!!!!!!!!!!!!!!!!!!!!!
1463 if(id1.ne.0)id1=id1+1000
1464 call tend_ls(id1, &
1465 ifds,ifde,jfds,jfde, &
1466 ifps,ifpe,jfps,jfpe, &
1467 ifts,ifte,jfts,jfte, &
1468 ifms,ifme,jfms,jfme, &
1469 ts+dt,dt,dx,dy, &
1470 lfn_2, &
1471 tbound3, &
1472 tend,ros, &
1473 fp &
1474 )
1476 do j=jfts,jfte
1477 do i=ifts,ifte
1478 lfn_out(i,j) = lfn_0(i,j) + dt*tend(i,j)
1479 lfn_2(i,j) = lfn_out(i,j) ! lfn_2=lfn_out (needed for reinitialization purposes)
1480 enddo
1481 enddo
1483 #ifdef DM_PARALLEL
1484 # include "HALO_FIRE_LFN_2.inc"
1485 #endif
1487 ! CFL check
1489 tbound=min(tbound,tbound2,tbound3)
1491 !$OMP CRITICAL(FIRE_CORE_CRIT)
1492 write(msg,'(a,f10.2,4(a,f7.2))')'prop_ls: time',ts,' dt=',dt,' bound',min(tbound,999.99), &
1493 ' dx=',dx,' dy=',dy
1494 !$OMP END CRITICAL(FIRE_CORE_CRIT)
1495 call message(msg)
1496 if(dt>tbound)then
1497 !$OMP CRITICAL(FIRE_CORE_CRIT)
1498 write(msg,'(2(a,f10.2))')'prop_ls: WARNING: time step ',dt, &
1499 ' > bound =',tbound
1500 !$OMP END CRITICAL(FIRE_CORE_CRIT)
1501 call message(msg)
1502 endif
1504 end subroutine prop_ls_rk3
1506 !
1507 !*****************************
1508 !
1510 subroutine reinit_ls_rk3(id, &
1511 ifts,ifte,jfts,jfte, &
1512 ifms,ifme,jfms,jfme, &
1513 ifds,ifde,jfds,jfde, &
1514 ifps,ifpe,jfps,jfpe, &
1515 ts,dt,dx,dy, &
1516 lfn_in, &
1517 lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3, &
1518 lfn_out,tign, &
1519 grid, &
1520 ids,ide,jds,jde,kds,kde, &
1521 ims,ime,jms,jme,kms,kme, &
1522 ips,ipe,jps,jpe,kps,kpe &
1523 )
1526 #ifdef DM_PARALLEL
1527 USE module_dm , ONLY : ntasks_x,ntasks_y,local_communicator,mytask,ntasks
1528 USE module_comm_dm , ONLY : halo_fire_lfn_s1_sub,halo_fire_lfn_s2_sub,halo_fire_lfn_s3_sub
1529 #endif
1530 USE module_domain , only: domain
1531 USE module_configure, only: grid_config_rec_type
1533 implicit none
1535 !*************************************************************************************!
1536 !*** Level-set function reinitialization following: ***!
1537 !*** Sussman, Smereka, Osher. Journal of Computational Physics 114, 146-159 (1994) ***!
1538 !*** implemented by Domingo Munoz-Esparza (NCAR/RAL, April 2016) ***!
1539 !*** Reference: D. Munoz-Esparza, B. Kosovic, P. Jimenez, J. Coen: "An accurate ***!
1540 !*** fire-spread algorithm in the Weather Research and Forecasting model using the ***!
1541 !*** level-set method", Journal of Advances in Modeling Earth Systems, 2018 ***!
1542 !*** https://doi.org/10.1002/2017MS001108 ***!
1543 !*************************************************************************************!
1545 !*** purpose: reinitialize the level-set function
1547 type(domain) , target :: grid
1548 integer, intent(in):: ids,ide,jds,jde,kds,kde, &
1549 ims,ime,jms,jme,kms,kme, &
1550 ips,ipe,jps,jpe,kps,kpe
1552 integer,intent(in)::ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,id
1553 integer,intent(in)::ifds,ifde,jfds,jfde,ifps,ifpe,jfps,jfpe
1554 real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_in,tign
1555 real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_2,lfn_s0,lfn_s1,lfn_s2,lfn_s3
1556 real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_out
1557 real,intent(in)::dx,dy,ts,dt
1558 real::dt_s
1559 real,dimension(ifts:ifte,jfts:jfte):: tend_1,tend_2,tend_3
1560 real::diffLx,diffLy,diffRx,diffRy,diff2x,diff2y,grad,time_now
1561 integer::nts,i,j,k,kk
1562 intrinsic epsilon
1563 character(len=128)::msg
1564 integer::itso,iteo,jtso,jteo
1565 real::threshold_HLl,threshold_HLu
1566 integer,parameter::bdy_eno1=10
1567 !*** executable
1569 threshold_HLl=-fire_lsm_band_ngp*dx
1570 threshold_HLu=fire_lsm_band_ngp*dx
1572 ! Define S0 based on current lfn values
1574 do j=jfts,jfte
1575 do i=ifts,ifte
1576 lfn_s0(i,j) = lfn_2(i,j)/sqrt(lfn_2(i,j)**2.0+dx**2.0)
1577 lfn_s3(i,j) = lfn_2(i,j)
1578 enddo
1579 enddo
1580 #ifdef DM_PARALLEL
1581 # include "HALO_FIRE_LFN_S3.inc"
1582 #endif
1584 call continue_at_boundary(1,1,fire_lfn_ext_up, & !extend by extrapolation but never down
1585 ifms,ifme,jfms,jfme, & ! memory dims
1586 ifds,ifde,jfds,jfde, & ! domain
1587 ifps,ifpe,jfps,jfpe, & ! patch
1588 ifts,ifte,jfts,jfte, & ! tile
1589 itso,iteo,jtso,jteo, & ! where set now
1590 lfn_s3)
1593 dt_s=0.01*dx
1594 do nts=1,fire_lsm_reinit_iter ! iterate to solve to steady state (1 iter each time step is enoguh)
1595 ! Runge-Kutta step 1 !
1597 call advance_ls_reinit( &
1598 ifms,ifme,jfms,jfme, &
1599 ifds,ifde,jfds,jfde, &
1600 ifts,ifte,jfts,jfte, &
1601 dx,dy,dt_s,bdy_eno1,threshold_HLu, &
1602 lfn_s0,lfn_s3,lfn_s3,lfn_s1,1.0/3.0) ! sign funcition, initial ls, current stage ls, next stage advanced ls, RK coefficient
1604 #ifdef DM_PARALLEL
1605 # include "HALO_FIRE_LFN_S1.inc"
1606 #endif
1608 call continue_at_boundary(1,1,fire_lfn_ext_up, &
1609 ifms,ifme,jfms,jfme, &
1610 ifds,ifde,jfds,jfde, &
1611 ifps,ifpe,jfps,jfpe, &
1612 ifts,ifte,jfts,jfte, &
1613 itso,iteo,jtso,jteo, &
1614 lfn_s1)
1616 ! Runge-Kutta step 2 !
1618 call advance_ls_reinit( &
1619 ifms,ifme,jfms,jfme, &
1620 ifds,ifde,jfds,jfde, &
1621 ifts,ifte,jfts,jfte, &
1622 dx,dy,dt_s,bdy_eno1,threshold_HLu, &
1623 lfn_s0,lfn_s3,lfn_s1,lfn_s2,1.0/2.0)
1625 #ifdef DM_PARALLEL
1626 # include "HALO_FIRE_LFN_S2.inc"
1627 #endif
1629 call continue_at_boundary(1,1,fire_lfn_ext_up, &
1630 ifms,ifme,jfms,jfme, &
1631 ifds,ifde,jfds,jfde, &
1632 ifps,ifpe,jfps,jfpe, &
1633 ifts,ifte,jfts,jfte, &
1634 itso,iteo,jtso,jteo, &
1635 lfn_s2)
1637 ! Runge-Kutta step 3 !
1639 call advance_ls_reinit( &
1640 ifms,ifme,jfms,jfme, &
1641 ifds,ifde,jfds,jfde, &
1642 ifts,ifte,jfts,jfte, &
1643 dx,dy,dt_s,bdy_eno1,threshold_HLu, &
1644 lfn_s0,lfn_s3,lfn_s2,lfn_s3,1.0)
1646 #ifdef DM_PARALLEL
1647 # include "HALO_FIRE_LFN_S3.inc"
1648 #endif
1650 call continue_at_boundary(1,1,fire_lfn_ext_up, &
1651 ifms,ifme,jfms,jfme, &
1652 ifds,ifde,jfds,jfde, &
1653 ifps,ifpe,jfps,jfpe, &
1654 ifts,ifte,jfts,jfte, &
1655 itso,iteo,jtso,jteo, &
1656 lfn_s3)
1658 enddo ! end iterations for steady state solution of reinitialization PDE
1660 do j=jfts,jfte
1661 do i=ifts,ifte
1662 lfn_out(i,j)=lfn_s3(i,j) ! assing to lfn_out the reinitialized level-set function
1663 lfn_out(i,j)=min(lfn_out(i,j),lfn_in(i,j)) ! fire area can only increase
1664 enddo
1665 enddo
1667 end subroutine reinit_ls_rk3
1669 !
1670 !*****************************
1671 !
1673 subroutine advance_ls_reinit(ifms,ifme,jfms,jfme, &
1674 ifds,ifde,jfds,jfde, &
1675 ifts,ifte,jfts,jfte, &
1676 dx,dy,dt_s,bdy_eno1,threshold_HLu, &
1677 lfn_s0,lfn_ini,lfn_curr,lfn_fin,rk_coeff)
1679 ! Calculates right-hand-side forcing and advances a RK-stage the level-set reinitialization PDE
1680 ! Domingo Munoz-Esparza, NCAR/RAL
1681 ! Reference: D. Munoz-Esparza, B. Kosovic, P. Jimenez, J. Coen: "An accurate
1682 ! fire-spread algorithm in the Weather Research and Forecasting model using the
1683 ! level-set method", Journal of Advances in Modeling Earth Systems, 2018
1684 ! https://doi.org/10.1002/2017MS001108
1687 implicit none
1689 integer,intent(in)::ifms,ifme,jfms,jfme,ifts,ifte,jfts,jfte,ifds,ifde,jfds,jfde,bdy_eno1
1690 real,dimension(ifms:ifme,jfms:jfme),intent(in)::lfn_s0,lfn_ini,lfn_curr
1691 real,dimension(ifms:ifme,jfms:jfme),intent(inout)::lfn_fin
1692 real,intent(in)::dx,dy,dt_s,threshold_HLu,rk_coeff
1693 !
1694 integer::i,j
1695 real::diffLx,diffLy,diffRx,diffRy,diff2x,diff2y,grad,tend_r
1696 !
1698 do j=jfts,jfte
1699 do i=ifts,ifte
1700 if (i.lt.ifds+bdy_eno1 .OR. i.gt.ifde-bdy_eno1 .OR. j.lt.jfds+bdy_eno1 .OR. j.gt.jfde-bdy_eno1) then !
1701 diffLx=(lfn_curr(i,j)-lfn_curr(i-1,j))/dx
1702 diffLy=(lfn_curr(i,j)-lfn_curr(i,j-1))/dy
1703 diffRx=(lfn_curr(i+1,j)-lfn_curr(i,j))/dx
1704 diffRy=(lfn_curr(i,j+1)-lfn_curr(i,j))/dy
1705 diff2x=select_eno(diffLx,diffRx)
1706 diff2y=select_eno(diffLy,diffRy)
1707 else
1708 select case(fire_upwinding_reinit)
1709 case(1)
1710 diff2x=select_4th(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j))
1711 diff2y=select_4th(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2))
1712 diff2x=select_weno3(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j),lfn_s0(i,j)*diff2x)
1713 diff2y=select_weno3(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2),lfn_s0(i,j)*diff2y)
1714 case(2)
1715 diff2x=select_4th(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j))
1716 diff2y=select_4th(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2))
1717 diff2x=select_weno5(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i-3,j),lfn_curr(i+1,j),lfn_curr(i+2,j),lfn_curr(i+3,j),lfn_s0(i,j)*diff2x)
1718 diff2y=select_weno5(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j-3),lfn_curr(i,j+1),lfn_curr(i,j+2),lfn_curr(i,j+3),lfn_s0(i,j)*diff2y)
1719 case(3)
1720 if (lfn_curr(i,j).lt.threshold_HLu) then
1721 diff2x=select_4th(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j))
1722 diff2y=select_4th(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2))
1723 diff2x=select_weno3(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j),lfn_s0(i,j)*diff2x)
1724 diff2y=select_weno3(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2),lfn_s0(i,j)*diff2y)
1725 else
1726 diffLx=(lfn_curr(i,j)-lfn_curr(i-1,j))/dx
1727 diffLy=(lfn_curr(i,j)-lfn_curr(i,j-1))/dy
1728 diffRx=(lfn_curr(i+1,j)-lfn_curr(i,j))/dx
1729 diffRy=(lfn_curr(i,j+1)-lfn_curr(i,j))/dy
1730 diff2x=select_eno(diffLx,diffRx)
1731 diff2y=select_eno(diffLy,diffRy)
1732 endif
1733 case(4)
1734 if (lfn_curr(i,j).lt.threshold_HLu) then
1735 diff2x=select_4th(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j))
1736 diff2y=select_4th(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2))
1737 diff2x=select_weno5(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i-3,j),lfn_curr(i+1,j),lfn_curr(i+2,j),lfn_curr(i+3,j),lfn_s0(i,j)*diff2x)
1738 diff2y=select_weno5(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j-3),lfn_curr(i,j+1),lfn_curr(i,j+2),lfn_curr(i,j+3),lfn_s0(i,j)*diff2y)
1739 else
1740 diffLx=(lfn_curr(i,j)-lfn_curr(i-1,j))/dx
1741 diffLy=(lfn_curr(i,j)-lfn_curr(i,j-1))/dy
1742 diffRx=(lfn_curr(i+1,j)-lfn_curr(i,j))/dx
1743 diffRy=(lfn_curr(i,j+1)-lfn_curr(i,j))/dy
1744 diff2x=select_eno(diffLx,diffRx)
1745 diff2y=select_eno(diffLy,diffRy)
1746 endif
1747 case default
1748 if (lfn_curr(i,j).lt.threshold_HLu) then
1749 diff2x=select_4th(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i+1,j),lfn_curr(i+2,j))
1750 diff2y=select_4th(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j+1),lfn_curr(i,j+2))
1751 diff2x=select_weno5(dx,lfn_curr(i,j),lfn_curr(i-1,j),lfn_curr(i-2,j),lfn_curr(i-3,j),lfn_curr(i+1,j),lfn_curr(i+2,j),lfn_curr(i+3,j),lfn_s0(i,j)*diff2x)
1752 diff2y=select_weno5(dy,lfn_curr(i,j),lfn_curr(i,j-1),lfn_curr(i,j-2),lfn_curr(i,j-3),lfn_curr(i,j+1),lfn_curr(i,j+2),lfn_curr(i,j+3),lfn_s0(i,j)*diff2y)
1753 else
1754 diffLx=(lfn_curr(i,j)-lfn_curr(i-1,j))/dx
1755 diffLy=(lfn_curr(i,j)-lfn_curr(i,j-1))/dy
1756 diffRx=(lfn_curr(i+1,j)-lfn_curr(i,j))/dx
1757 diffRy=(lfn_curr(i,j+1)-lfn_curr(i,j))/dy
1758 diff2x=select_eno(diffLx,diffRx)
1759 diff2y=select_eno(diffLy,diffRy)
1760 endif
1761 end select
1762 endif
1763 grad=sqrt(diff2x*diff2x+diff2y*diff2y)
1764 tend_r=lfn_s0(i,j)*(1.0-grad)
1765 lfn_fin(i,j)=lfn_ini(i,j)+(dt_s*rk_coeff)*tend_r
1766 enddo
1767 enddo
1769 end subroutine advance_ls_reinit
1771 !
1772 !*****************************
1773 !
1775 subroutine tign_update(ifts,ifte,jfts,jfte, & ! tile dims
1776 ifms,ifme,jfms,jfme, & ! memory dims
1777 ifds,jfds,ifde,jfde, & ! domain dims
1778 ts,dt, & ! scalars in
1779 lfn_in,lfn_out,tign & ! arrays inout
1780 )
1782 implicit none
1784 integer,intent(in)::ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,ifds,jfds,ifde,jfde
1785 real,dimension(ifms:ifme,jfms:jfme),intent(inout)::tign
1786 real,dimension(ifms:ifme,jfms:jfme),intent(in)::lfn_in,lfn_out
1787 real,intent(in)::ts,dt
1788 real::time_now
1789 integer::i,j,k,kk
1790 intrinsic epsilon
1791 character(len=128)::msg
1792 !*** executable
1794 ! compute ignition time by interpolation
1795 ! the node was not burning at start but it is burning at end
1796 ! interpolate from the level functions at start and at end
1797 ! lfn_in==lfn is the level set function value at time ts
1798 ! lfn_out is the level set function value at time ts+dt (after reinitialization)
1799 ! 0 is the level set function value at time tign(i,j)
1800 ! thus assuming the level function is approximately linear =>
1801 ! tign(i,j)= ts + ((ts + td) - ts) * lfn_in / (lfn_in - lfn_out)
1802 ! = ts + dt * lfn_in / (lfn_in - lfn_out)
1804 time_now=ts+dt
1805 time_now = time_now + abs(time_now)*epsilon(time_now)*2.
1806 do j=jfts,jfte
1807 do i=ifts,ifte
1808 ! interpolate the cross-over time
1809 if (.not. lfn_out(i,j)>0 .and. lfn_in(i,j)>0)then
1810 tign(i,j) = ts + dt * lfn_in(i,j) / (lfn_in(i,j) - lfn_out(i,j))
1811 endif
1812 ! set the ignition time outside of burning region
1813 if(lfn_out(i,j)>0.)tign(i,j)=time_now
1814 enddo
1815 enddo
1816 ! stop simulation if fire is within boundary_guard grid points from the domain boundaries
1817 do j=jfts,jfte
1818 if (j.le.boundary_guard .or. j.gt.(jfde-boundary_guard)) then
1819 do i=ifts,ifte
1820 if (lfn_out(i,j).lt.0.) then
1821 WRITE(msg,*)'j-boundary reached'
1822 CALL wrf_debug(0,msg)
1823 WRITE(msg,*)'i,j,lfn_out=',i,j,lfn_out(i,j)
1824 CALL wrf_debug(0,msg)
1825 call crash('wrf: SUCCESS COMPLETE WRF. Fire has reached domain boundary.')
1826 endif
1827 enddo
1828 endif
1829 enddo
1830 do i=ifts,ifte
1831 if (i.le.boundary_guard .or. i.gt.(ifde-boundary_guard)) then
1832 do j=jfts,jfte
1833 if (lfn_out(i,j).lt.0.) then
1834 WRITE(msg,*)'i-boundary reached'
1835 CALL wrf_debug(0,msg)
1836 WRITE(msg,*)'i,j,lfn_out=',i,j,lfn_out(i,j)
1837 CALL wrf_debug(0,msg)
1838 call crash('wrf: SUCCESS COMPLETE WRF. Fire has reached domain boundary.')
1839 endif
1840 enddo
1841 endif
1842 enddo
1844 call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme, &
1845 lfn_out,'prop_ls: lfn out')
1847 end subroutine tign_update
1849 !
1850 !*****************************
1851 !
1853 subroutine calc_flame_length(ifts,ifte,jfts,jfte, &
1854 ifms,ifme,jfms,jfme, &
1855 ros,iboros,flame_length,ros_fl,fire_area)
1857 ! Diagnoses flame length
1858 ! Domingo Munoz-Esparza, NCAR/RAL
1860 implicit none
1862 integer,intent(in)::ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme
1863 real,dimension(ifms:ifme,jfms:jfme),intent(in)::ros,iboros,fire_area
1864 real,dimension(ifms:ifme,jfms:jfme),intent(out)::flame_length,ros_fl
1865 ! local
1866 integer::i,j
1868 do j=jfts,jfte
1869 do i=ifts,ifte
1870 if (fire_area(i,j).gt.0.0 .AND. fire_area(i,j).lt.1.0) then
1871 flame_length(i,j)=0.0775*(iboros(i,j)*ros(i,j))**0.46 ! flame length according to Byram (1959)
1872 ros_fl(i,j)=ros(i,j)
1873 else
1874 flame_length(i,j)=0.0
1875 ros_fl(i,j)=0.0
1876 endif
1877 enddo
1878 enddo
1880 end subroutine calc_flame_length
1882 !
1883 !*****************************
1884 !
1886 subroutine tend_ls( id, &
1887 ids,ide,jds,jde, & ! domain - nodes where lfn defined
1888 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1889 its,ite,jts,jte, & ! region - nodes where tend computed
1890 ims,ime,jms,jme, & ! memory dims for ros
1891 t,dt,dx,dy, & ! scalars in
1892 lfn, & ! arrays in
1893 tbound, & ! scalars out
1894 tend, ros, & ! arrays out
1895 fp &
1896 )
1898 implicit none
1899 ! purpose
1900 ! compute the right hand side of the level set equation
1902 !*** arguments
1903 integer,intent(in)::id
1904 integer,intent(in)::ims,ime,jms,jme,its,ite,jts,jte
1905 integer, intent(in)::ids,ide,jds,jde,ips,ipe,jps,jpe
1906 real,intent(in)::t ! time
1907 real,intent(in)::dt,dx,dy ! mesh step
1908 real,dimension(ims:ime,jms:jme),intent(inout)::lfn ! level set function
1909 real,dimension(ims:ime,jms:jme),intent(out)::tend ! tendency (rhs of the level set pde)
1910 real,dimension(ims:ime,jms:jme),intent(out)::ros ! rate of spread
1911 real,intent(out)::tbound ! max allowed time step
1912 type(fire_params),intent(in)::fp
1914 !*** local
1915 real:: te,diffLx,diffLy,diffRx,diffRy, &
1916 diffCx,diffCy,diff2x,diff2y,grad,rr, &
1917 ros_base,ros_wind,ros_slope,ros_back,advx,advy,scale,nvx,nvy, &
1918 speed,tanphi
1919 integer::i,j,itso,iteo,jtso,jteo
1920 character(len=128)msg
1922 ! constants
1923 real, parameter:: eps=epsilon(0.0)
1924 !intrinsic epsilon
1925 real, parameter:: zero=0.,one=1.,tol=100*eps, &
1926 safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
1927 !
1928 real :: diff2xn, diff2yn
1929 real :: a_valor, signo_x, signo_y
1930 real :: threshold_HLl,threshold_HLu
1931 real :: threshold_AV,fire_viscosity_var
1932 integer,parameter :: bdy_eno1=10
1934 ! f90 intrinsic function
1936 intrinsic max,min,sqrt,nint,tiny,huge
1939 !*** executable
1941 threshold_HLl=-fire_lsm_band_ngp*dx
1942 threshold_HLu=fire_lsm_band_ngp*dx
1943 threshold_AV=fire_viscosity_ngp*dx
1945 call continue_at_boundary(1,1,fire_lfn_ext_up, & !extend by extrapolation but never down
1946 ims,ime,jms,jme, & ! memory dims
1947 ids,ide,jds,jde, & ! domain
1948 ips,ipe,jps,jpe, & ! patch
1949 its,ite,jts,jte, & ! tile
1950 itso,iteo,jtso,jteo, & ! where set now
1951 lfn) ! field
1953 tbound=0
1954 do j=jts,jte
1955 do i=its,ite
1956 ! one sided differences
1957 diffRx = (lfn(i+1,j)-lfn(i,j))/dx
1958 diffLx = (lfn(i,j)-lfn(i-1,j))/dx
1959 diffRy = (lfn(i,j+1)-lfn(i,j))/dy
1960 diffLy = (lfn(i,j)-lfn(i,j-1))/dy
1961 diffCx = diffLx+diffRx ! TWICE CENTRAL DIFFERENCE
1962 diffCy = diffLy+diffRy
1964 if (i.lt.ids+bdy_eno1 .OR. i.gt.ide-bdy_eno1 .OR. j.lt.jds+bdy_eno1 .OR. j.gt.jde-bdy_eno1) then ! DME: use eno1 near domain boundaries
1965 diff2x=select_eno(diffLx,diffRx)
1966 diff2y=select_eno(diffLy,diffRy)
1967 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1968 else
1969 select case(fire_upwinding)
1970 case(0) ! none
1971 grad=sqrt(diffCx**2 + diffCy**2)
1972 case(1) ! standard
1973 diff2x=select_upwind(diffLx,diffRx)
1974 diff2y=select_upwind(diffLy,diffRy)
1975 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1976 case(2) ! godunov per osher/fedkiw
1977 diff2x=select_godunov(diffLx,diffRx)
1978 diff2y=select_godunov(diffLy,diffRy)
1979 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1980 case(3) ! eno
1981 diff2x=select_eno(diffLx,diffRx)
1982 diff2y=select_eno(diffLy,diffRy)
1983 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1984 case(4) ! Sethian - twice stronger pushdown of bumps
1985 grad=sqrt(max(diffLx,0.)**2+min(diffRx,0.)**2 &
1986 + max(diffLy,0.)**2+min(diffRy,0.)**2)
1987 case(5) ! 2nd order (DME)
1988 diff2x=select_2nd(rr,dx,lfn(i,j),lfn(i-1,j),lfn(i+1,j))
1989 diff2y=select_2nd(rr,dy,lfn(i,j),lfn(i,j-1),lfn(i,j+1))
1990 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1991 case(6) ! 3rd-order WENO (DME)
1992 a_valor=select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))*fp%vx(i,j)+ &
1993 select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))*fp%vy(i,j)
1994 signo_x=a_valor*select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))
1995 signo_y=a_valor*select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))
1996 diff2x=select_weno3(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j),signo_x)
1997 diff2y=select_weno3(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2),signo_y)
1998 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1999 case(7) ! 5th-order WENO (DME)
2000 a_valor=select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))*fp%vx(i,j)+ &
2001 select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))*fp%vy(i,j)
2002 signo_x=a_valor*select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))
2003 signo_y=a_valor*select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))
2004 diff2x=select_weno5(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i-3,j),lfn(i+1,j),lfn(i+2,j),lfn(i+3,j),signo_x)
2005 diff2y=select_weno5(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j-3),lfn(i,j+1),lfn(i,j+2),lfn(i,j+3),signo_y)
2006 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
2007 case(8) ! high-order local band method WENO3/ENO1 (DME)
2008 if (abs(lfn(i,j)).lt.threshold_HLu) then
2009 a_valor=select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))*fp%vx(i,j)+ &
2010 select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))*fp%vy(i,j)
2011 signo_x=a_valor*select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))
2012 signo_y=a_valor*select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))
2013 diff2x=select_weno3(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j),signo_x)
2014 diff2y=select_weno3(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2),signo_y)
2015 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
2016 else
2017 diff2x=select_eno(diffLx,diffRx)
2018 diff2y=select_eno(diffLy,diffRy)
2019 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
2020 endif
2021 case(9) ! high-order local band method WENO5/ENO1 (DME)
2022 if (abs(lfn(i,j)).lt.threshold_HLu) then
2023 a_valor=select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))*fp%vx(i,j)+ &
2024 select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))*fp%vy(i,j)
2025 signo_x=a_valor*select_4th(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i+1,j),lfn(i+2,j))
2026 signo_y=a_valor*select_4th(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j+1),lfn(i,j+2))
2027 diff2x=select_weno5(dx,lfn(i,j),lfn(i-1,j),lfn(i-2,j),lfn(i-3,j),lfn(i+1,j),lfn(i+2,j),lfn(i+3,j),signo_x)
2028 diff2y=select_weno5(dy,lfn(i,j),lfn(i,j-1),lfn(i,j-2),lfn(i,j-3),lfn(i,j+1),lfn(i,j+2),lfn(i,j+3),signo_y)
2029 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
2030 else
2031 diff2x=select_eno(diffLx,diffRx)
2032 diff2y=select_eno(diffLy,diffRy)
2033 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
2034 endif
2035 case default
2036 grad=0.
2037 end select
2038 endif
2039 ! DME use the same discretization scheme from fire_upwinding option to calculate the normals
2040 scale=sqrt(grad**2.0+eps)
2041 nvx=diff2x/scale
2042 nvy=diff2y/scale
2044 ! wind speed in direction of spread
2045 speed = fp%vx(i,j)*nvx + fp%vy(i,j)*nvy
2046 ! get rate of spread from wind speed and slope
2047 call fire_ros(ros_base,ros_wind,ros_slope, &
2048 nvx,nvy,i,j,fp)
2049 rr=ros_base + ros_wind + fire_slope_factor*ros_slope
2051 if(fire_grows_only.gt.0)rr=max(rr,0.)
2052 ! set ros for output
2053 if(i.ge.its.and.i.le.ite.and.j.ge.jts.and.j.le.jte)ros(i,j)=rr
2055 if(fire_upwind_split.eq.0)then
2057 ! get rate of spread
2058 te = -rr*grad ! normal term
2060 else
2062 ! normal direction backing rate only
2063 te = - ros_base*grad
2065 ! advection in wind direction
2066 if (abs(speed)> eps) then
2067 advx=fp%vx(i,j)*ros_wind/speed
2068 advy=fp%vy(i,j)*ros_wind/speed
2069 else
2070 advx=0
2071 advy=0
2072 endif
2074 tanphi = fp%dzdxf(i,j)*nvx + fp%dzdyf(i,j)*nvy
2075 ! advection from slope direction
2076 if(abs(tanphi)>eps) then
2077 advx=advx+fp%dzdxf(i,j)*ros_slope/tanphi
2078 advy=advy+fp%dzdyf(i,j)*ros_slope/tanphi
2079 endif
2081 if(fire_upwind_split.eq.1)then
2083 ! one-sided upwinding
2084 te = te - max(advx,0.)*diffLx - min(advx,0.)*diffRy &
2085 - max(advy,0.)*diffLy - min(advy,0.)*diffRy
2088 elseif(fire_upwind_split.eq.2)then
2090 ! Lax-Friedrichs
2091 call crash('prop_ls: bad fire_upwind_split, Lax-Friedrichs not done yet')
2093 else
2095 call crash('prop_ls: bad fire_upwind_split')
2097 endif
2098 endif
2100 ! cfl condition
2101 if (grad > 0.) then
2102 tbound = max(tbound,rr*(abs(diff2x)/dx+abs(diff2y)/dy)/grad)
2103 endif
2105 ! DME: adds spatially variable artificial viscosity
2106 if (abs(lfn(i,j)).lt.threshold_AV .AND. (i.gt.ids+bdy_eno1 .and. i.lt.ide-bdy_eno1) .AND. (j.gt.jds+bdy_eno1 .and. j.lt.jde-bdy_eno1)) then
2107 fire_viscosity_var=fire_viscosity_bg
2108 elseif (abs(lfn(i,j)).ge.threshold_AV .AND. abs(lfn(i,j)).lt.threshold_AV*(1.0+fire_viscosity_band) .AND. (i.gt.ids+10 .and. i.lt.ide-10) .AND. (j.gt.jds+10 .and. j.lt.jde-10)) then
2109 fire_viscosity_var=min(fire_viscosity_bg+(fire_viscosity-fire_viscosity_bg)*(abs(lfn(i,j))-threshold_AV)/(fire_viscosity_band*threshold_AV),fire_viscosity)
2110 else
2111 fire_viscosity_var=fire_viscosity
2112 endif
2113 te=te + fire_viscosity_var*abs(rr)*((diffRx-diffLx)+(diffRy-diffLy))
2114 tend(i,j)=te
2115 enddo
2116 enddo
2118 call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme, &
2119 tend,'tend_ls: tend out')
2121 ! the final CFL bound
2122 tbound = 1/(tbound+tol)
2124 end subroutine tend_ls
2126 !
2127 !**************************
2128 !
2129 !
2131 real function select_upwind(diffLx,diffRx)
2132 implicit none
2133 real, intent(in):: diffLx, diffRx
2134 real diff2x
2136 ! upwind differences, L or R if bith same sign, otherwise zero
2138 diff2x=0
2139 if (diffLx>0.and.diffRx>0.)diff2x=diffLx
2140 if (diffLx<0.and.diffRx<0.)diff2x=diffRx
2142 select_upwind=diff2x
2143 end function select_upwind
2146 !
2147 !**************************
2148 !
2150 real function select_godunov(diffLx,diffRx)
2151 implicit none
2152 real, intent(in):: diffLx, diffRx
2153 real diff2x,diffCx
2155 ! Godunov scheme: upwind differences, L or R or none
2156 ! always test on > or < never = , much faster because of IEEE
2157 ! central diff >= 0 => take left diff if >0, ortherwise 0
2158 ! central diff <= 0 => take right diff if <0, ortherwise 0
2160 diff2x=0
2161 diffCx=diffRx+diffLx
2162 if (diffLx>0.and..not.diffCx<0)diff2x=diffLx
2163 if (diffRx<0.and. diffCx<0)diff2x=diffRx
2165 select_godunov=diff2x
2166 end function select_godunov
2168 !
2169 !**************************
2170 !
2172 real function select_eno(diffLx,diffRx)
2173 implicit none
2174 real, intent(in):: diffLx, diffRx
2175 real diff2x
2177 ! 1st order ENO scheme
2179 if (.not.diffLx>0 .and. .not.diffRx>0)then
2180 diff2x=diffRx
2181 elseif(.not.diffLx<0 .and. .not.diffRx<0)then
2182 diff2x=diffLx
2183 elseif(.not.diffLx<0 .and. .not.diffRx>0)then
2184 if(.not. abs(diffRx) < abs(diffLx))then
2185 diff2x=diffRx
2186 else
2187 diff2x=diffLx
2188 endif
2189 else
2190 diff2x=0.
2191 endif
2193 select_eno=diff2x
2194 end function select_eno
2196 !
2197 !**************************
2198 !
2201 real function select_2nd(ros,dx,lfn_i,lfn_im1,lfn_ip1)
2202 implicit none
2203 real, intent(in):: lfn_i, lfn_im1, lfn_ip1
2204 real, intent(in):: ros, dx
2205 real:: diff2x_p, diff2x_m
2207 ! 2nd-order advection scheme in the x,y-direction (DME)
2209 diff2x_p=0.
2210 diff2x_m=0.
2211 diff2x_p=(lfn_ip1+lfn_i)/(2.*dx)
2212 diff2x_m=(lfn_i+lfn_im1)/(2.*dx)
2213 select_2nd=diff2x_p-diff2x_m
2214 end function select_2nd
2216 !
2217 !**************************
2218 !
2221 real function select_4th(dx,lfn_i,lfn_im1,lfn_im2,lfn_ip1,lfn_ip2)
2222 implicit none
2223 real, intent(in):: lfn_i,lfn_im1,lfn_im2,lfn_ip1,lfn_ip2
2224 real, intent(in):: dx
2225 real:: diff2x_p,diff2x_m
2227 ! 4th-order advection scheme in the x,y-direction (DME)
2229 diff2x_p=0.
2230 diff2x_m=0.
2231 diff2x_p=(7.0*lfn_ip1+7.0*lfn_i-lfn_ip2-lfn_im1)/(12.*dx)
2232 diff2x_m=(7.0*lfn_i+7.0*lfn_im1-lfn_ip1-lfn_im2)/(12.*dx)
2233 select_4th=diff2x_p-diff2x_m
2234 end function select_4th
2236 !
2237 !**************************
2238 !
2241 real function select_weno3(dx,lfn_it,lfn_im1t,lfn_im2t,lfn_ip1t,lfn_ip2t,uf)
2242 implicit none
2243 real, intent(in)::lfn_it,lfn_im1t,lfn_im2t,lfn_ip1t,lfn_ip2t,uf
2244 real:: lfn_i,lfn_im1,lfn_im2,lfn_ip1,lfn_ip2
2245 real, intent(in):: dx
2246 real:: flux_p,flux_m
2247 real, parameter:: gamma1=1./3.,gamma2=2./3.,tol=1e-6
2248 real:: w1,w2,w1t,w2t,beta1,beta2
2249 real:: fh_1,fh_2
2251 ! 3rd-order advection WENO scheme in the x,y-direction (DME)
2253 flux_p=0.
2254 flux_m=0.
2256 if (uf .ge. 0.0 ) then
2257 lfn_i=lfn_it
2258 lfn_im1=lfn_im1t
2259 lfn_im2=lfn_im2t
2260 lfn_ip1=lfn_ip1t
2261 else
2262 lfn_i=lfn_it
2263 lfn_im1=lfn_ip1t
2264 lfn_im2=lfn_ip2t
2265 lfn_ip1=lfn_im1t
2266 endif
2268 ! numerical flux at i,j+1/2 face
2269 fh_1=-0.5*lfn_im1+1.5*lfn_i
2270 fh_2=0.5*lfn_i+0.5*lfn_ip1
2271 beta1=(lfn_i-lfn_im1)**2
2272 beta2=(lfn_ip1-lfn_i)**2
2273 w1t=gamma1/(beta1+tol)**2
2274 w2t=gamma2/(beta2+tol)**2
2275 w1=w1t/(w1t+w2t)
2276 w2=w2t/(w1t+w2t)
2277 flux_p=w1*fh_1+w2*fh_2
2278 ! numerical flux at i,j-1/2 face
2279 fh_1=-0.5*lfn_im2+1.5*lfn_im1
2280 fh_2=0.5*lfn_im1+0.5*lfn_i
2281 beta1=(lfn_im1-lfn_im2)**2
2282 beta2=(lfn_i-lfn_im1)**2
2283 w1t=gamma1/(beta1+tol)**2
2284 w2t=gamma2/(beta2+tol)**2
2285 w1=w1t/(w1t+w2t)
2286 w2=w2t/(w1t+w2t)
2287 flux_m=w1*fh_1+w2*fh_2
2289 if (uf .ge. 0.0 ) then
2290 select_weno3=(flux_p-flux_m)/dx
2291 else
2292 select_weno3=(flux_m-flux_p)/dx
2293 endif
2294 end function select_weno3
2296 !
2297 !**************************
2298 !
2301 real function select_weno5(dx,lfn_it,lfn_im1t,lfn_im2t,lfn_im3t,lfn_ip1t,lfn_ip2t,lfn_ip3t,uf)
2302 implicit none
2303 real, intent(in):: lfn_it,lfn_im1t,lfn_im2t,lfn_im3t,lfn_ip1t,lfn_ip2t,lfn_ip3t,uf
2304 real:: lfn_i,lfn_im1,lfn_im2,lfn_im3,lfn_ip1,lfn_ip2,lfn_ip3
2305 real, intent(in):: dx
2306 real:: flux_p,flux_m
2307 real, parameter:: gamma1=1./10.,gamma2=3./5.,gamma3=3./10.,tol=1e-6
2308 real:: w1,w2,w3,w1t,w2t,w3t,beta1,beta2,beta3
2309 real:: fh_1,fh_2,fh_3
2311 ! 5th-order advection WENO scheme in the x,y-direction (DME)
2313 flux_p=0.
2314 flux_m=0.
2316 if (uf .ge. 0.0 ) then
2317 lfn_i=lfn_it
2318 lfn_im1=lfn_im1t
2319 lfn_im2=lfn_im2t
2320 lfn_im3=lfn_im3t
2321 lfn_ip1=lfn_ip1t
2322 lfn_ip2=lfn_ip2t
2323 else
2324 lfn_i=lfn_it
2325 lfn_im1=lfn_ip1t
2326 lfn_im2=lfn_ip2t
2327 lfn_im3=lfn_ip3t
2328 lfn_ip1=lfn_im1t
2329 lfn_ip2=lfn_im2t
2330 endif
2332 ! numerical flux at i,j+1/2 face
2333 fh_1=(2.0*lfn_im2-7.0*lfn_im1+11.0*lfn_i)/6.0
2334 fh_2=(-1.0*lfn_im1+5.0*lfn_i+2.0*lfn_ip1)/6.0
2335 fh_3=(2.0*lfn_i+5.0*lfn_ip1-1.0*lfn_ip2)/6.0
2336 beta1=(13.0/12.0)*(lfn_im2-2.0*lfn_im1+lfn_i)**2+0.25*(lfn_im2-4.0*lfn_im1+3.0*lfn_i)**2
2337 beta2=(13.0/12.0)*(lfn_im1-2.0*lfn_i+lfn_ip1)**2+0.25*(lfn_im1-lfn_ip1)**2
2338 beta3=(13.0/12.0)*(lfn_i-2.0*lfn_ip1+lfn_ip2)**2+0.25*(3.0*lfn_i-4.0*lfn_ip1+lfn_ip2)**2
2339 w1t=gamma1/(beta1+tol)**2
2340 w2t=gamma2/(beta2+tol)**2
2341 w3t=gamma3/(beta3+tol)**2
2342 w1=w1t/(w1t+w2t+w3t)
2343 w2=w2t/(w1t+w2t+w3t)
2344 w3=w3t/(w1t+w2t+w3t)
2345 flux_p=w1*fh_1+w2*fh_2+w3*fh_3
2346 ! numerical flux at i,j-1/2 face
2347 fh_1=(2.0*lfn_im3-7.0*lfn_im2+11.0*lfn_im1)/6.0
2348 fh_2=(-1.0*lfn_im2+5.0*lfn_im1+2.0*lfn_i)/6.0
2349 fh_3=(2.0*lfn_im1+5.0*lfn_i-1.0*lfn_ip1)/6.0
2350 beta1=(13.0/12.0)*(lfn_im3-2.0*lfn_im2+lfn_im1)**2+0.25*(lfn_im3-4.0*lfn_im2+3.0*lfn_im1)**2
2351 beta2=(13.0/12.0)*(lfn_im2-2.0*lfn_im1+lfn_i)**2+0.25*(lfn_im2-lfn_i)**2
2352 beta3=(13.0/12.0)*(lfn_im1-2.0*lfn_i+lfn_ip1)**2+0.25*(3.0*lfn_im1-4.0*lfn_i+lfn_ip1)**2
2353 w1t=gamma1/(beta1+tol)**2
2354 w2t=gamma2/(beta2+tol)**2
2355 w3t=gamma3/(beta3+tol)**2
2356 w1=w1t/(w1t+w2t+w3t)
2357 w2=w2t/(w1t+w2t+w3t)
2358 w3=w3t/(w1t+w2t+w3t)
2359 flux_m=w1*fh_1+w2*fh_2+w3*fh_3
2361 if (uf .ge. 0.0 ) then
2362 select_weno5=(flux_p-flux_m)/dx
2363 else
2364 select_weno5=(flux_m-flux_p)/dx
2365 endif
2366 end function select_weno5
2368 !
2369 !**************************
2370 !
2372 real function speed_func(diffCx,diffCy,dx,dy,i,j,fp)
2373 !*** purpose
2374 ! the level set method speed function
2375 implicit none
2376 !*** arguments
2377 real, intent(in)::diffCx,diffCy ! x and y coordinates of the direction of propagation
2378 real, intent(in)::dx,dy ! x and y coordinates of the direction of propagation
2379 integer, intent(in)::i,j ! indices of the node to compute the speed at
2380 type(fire_params),intent(in)::fp
2381 !*** local
2382 real::scale,nvx,nvy,r
2383 real::ros_base , ros_wind , ros_slope
2384 real, parameter:: eps=epsilon(0.0)
2385 !*** executable
2386 ! normal direction, from central differences
2387 scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps)
2388 nvx=diffCx/scale
2389 nvy=diffCy/scale
2391 ! get rate of spread from wind speed and slope
2393 call fire_ros(ros_base,ros_wind,ros_slope, &
2394 nvx,nvy,i,j,fp)
2396 r=ros_base + ros_wind + ros_slope
2397 if(fire_grows_only.gt.0)r=max(r,0.)
2398 speed_func=r
2400 end function speed_func
2402 end module module_fr_fire_core