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