| blob | 50e8c2f8e4407ba63aa7baa6268ac43ba0b8fdf0 |
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_2
8 !#define DEBUG_OUT_FUEL_LEFT
10 module module_fr_sfire_core
12 use module_fr_sfire_phys
13 use module_fr_sfire_util
15 ! The mathematical core of the fire spread model. No physical constants here.
16 !
17 ! subroutine sfire_core: only this routine should be called from the outside.
18 ! subroutine fuel_left: compute remaining fuel from time of ignition.
19 ! subroutine prop_ls: propagation of curve in normal direction.
21 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 !$OMP CRITICAL(SFIRE_CORE_CRIT)
214 write(msg,'(a,2f11.6,a,2f11.6)')'ignite_fire: from',sx,sy,' to ',ex,ey
215 call message(msg)
216 write(msg,'(a,2f11.2,a,f8.1,a)')'units ',unit_xf,unit_yf,' m max dist ',dmax,' m'
217 call message(msg)
218 write(msg,'(a,f4.1,a,f8.1,a,i10)')' radius ',r,' time',time_ign,' ignited nodes',ignited
219 call message(msg)
220 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
221 end subroutine ignite_fire
223 !
224 !**********************
225 !
227 subroutine fuel_left(&
228 ifds,ifde,jfds,jfde&
229 ims,ime,jms,jme, &
230 its,ite,jts,jte, &
231 ifs,ife,jfs,jfe, &
232 lfn, tign, fuel_time, time_now, fuel_frac, fire_area)
233 implicit none
235 !*** purpose: determine fraction of fuel remaining
236 !*** NOTE: because variables are cell centered, need halo/sync width 1 before
238 !*** Jan Mandel August 2007 email: jmandel@ucar.edu or Jan.Mandel@gmail.com
240 !*** arguments
242 integer, intent(in) ::ifds,ifde,jfds,jfde,its,ite,jts,jte,ims,ime &
243 ,jms,jme,ifs,ife,jfs,jfe
244 real, intent(in), dimension(ims:ime,jms:jme)::lfn,tign,fuel_time
245 real, intent(in):: time_now
246 real, intent(out), dimension(ifs:ife,jfs:jfe)::fuel_frac
247 real, intent(out), dimension(ims:ime,jms:jme):: fire_area
249 ! ims,ime,jms,jme in memory dimensions
250 ! its,ite,jts,jte in tile dimensions (cells where fuel_frac computed)
251 ! ifs,ife,jfs,jfe in fuel_frac memory dimensions
252 ! lfn in level function, at nodes at midpoints of cells
253 ! tign in ignition time, at nodes at nodes at midpoints of cells
254 ! fuel_time in time constant of fuel, per cell
255 ! time_now in time now
256 ! fuel_frac out fraction of fuel remaining, per cell
257 ! fire_area out fraction of cell area on fire
259 !*** local
261 integer::i,j,ir,jr,icl,jcl,isubcl,jsubcl,i2,j2,ii,jj,its1,jts1,ite1,jte1
262 real::fmax,frat,helpsum1,helpsum2,fuel_left_ff,fire_area_ff,rx,ry,tignf(2,2)
263 real,dimension(3,3)::tff,lff
264 ! help variables instead of arrays fuel_left_f and fire_area_f
265 real::lffij,lffi1j,lffij1,lffi1j1,tifij,tifi1j,tifij1,tifi1j1,tx,ty,txx,tyy
266 ! variables for calculation instead of lff(i,j) and tif(i,j)is lffij,tifij etc..#define IFCELLS (ite-its+1)*fuel_left_irl
267 #define JFCELLS (jte-jts+1)*fuel_left_jrl
268 character(len=128)::msg
269 integer::m,omp_get_thread_num
272 call check_mesh_2dim(its-1,ite+1,jts-1,jte+1,ims,ime,jms,jme)
273 call check_mesh_2dim(its,ite,jts,jte,ifs,ife,jfs,jfe)
275 ! refinement
276 ir=fuel_left_irl
277 jr=fuel_left_jrl
279 if ((ir.ne.2).or.(jr.ne.2)) then
280 call crash('fuel_left: ir.ne.2 or jr.ne.2 ')
281 endif
283 rx=1./ir
284 ry=1./jr
286 ! interpolate level set function to finer grid
287 #ifdef DEBUG_OUT_FUEL_LEFT
288 call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,lff,'lff',0)
289 call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,tif,'tif',0)
290 #endif
292 !
293 ! example for ir=2:
294 !
295 ! | coarse cell |
296 ! its-1 its ite ite+1
297 ! -------X------------|-----.-----X-----.-----|--........----|----------X----------|---------X
298 ! fine node 1 2 3 2*(ite-its+1)
299 ! fine cell 1 2 cell 2*(ite-its+1)
303 ! Loop over cells in Tile
304 ! Changes made by Volodymyr Kondratenko 09/24/2009
305 its1=max(its,ifds+1)
306 do icl=its1,ite
307 do jcl=jts,jte
308 helpsum1=0
309 helpsum2=0
311 if (icl.eq.ifs) then
312 if (lfn(icl,jcl).lt.0) then
313 call crash('tign_interp: lfn on the boundary should be not burning but &
314 lfn<0, when icl=its ')
315 elseif (jcl.eq.jfs) then
316 call tign_lfn_interpolation(time_now,icl+1,jcl+1,ims,ime,jms,jme, &
317 tign,lfn,tff,lff)
319 elseif (jcl.eq.jfe) then
320 call tign_lfn_interpolation(time_now,icl+1,jcl-1,ims,ime,jms,jme, &
321 tign,lfn,tff,lff)
322 else
324 call tign_lfn_interpolation(time_now,icl+1,jcl,ims,ime,jms,jme, &
325 tign,lfn,tff,lff)
326 endif
328 elseif (icl.eq.ife) then
329 if (lfn(icl,jcl).lt.0) then
330 call crash('tign_interp: lfn on the boundary should be not burning but &
331 lfn<0 when icl=ite')
333 elseif (jcl.eq.jfs) then
334 call tign_lfn_interpolation(time_now,icl-1,jcl+1,ims,ime,jms,jme, &
335 tign,lfn,tff,lff)
337 elseif (jcl.eq.jfe) then
338 call tign_lfn_interpolation(time_now,icl-1,jcl-1,ims,ime,jms,jme, &
339 tign,lfn,tff,lff)
340 else
342 call tign_lfn_interpolation(time_now,icl-1,jcl,ims,ime,jms,jme, &
343 tign,lfn,tff,lff)
344 endif
346 elseif (jcl.eq.jfs) then
347 if (lfn(icl,jcl).lt.0) then
348 call crash('tign_interp: lfn on the boundary should be not burning but &
349 lfn<0 when icl=ite')
351 else
352 call tign_lfn_interpolation(time_now,icl,jcl+1,ims,ime,jms,jme, &
353 tign,lfn,tff,lff)
354 endif
356 elseif (jcl.eq.jfe) then
357 if (lfn(icl,jcl).lt.0) then
358 call crash('tign_interp: lfn on the boundary should be not burning but &
359 lfn<0 when icl=ite')
361 else
362 call tign_lfn_interpolation(time_now,icl,jcl-1,ims,ime,jms,jme, &
363 tign,lfn,tff,lff)
364 endif
365 else
369 ! Loop over subcells in cell #(icl,jcl)
370 write(*,*)"ifs,ife",ifs,ife
371 write(*,*)"jfs,jfe",jfs,jfe
372 write(*,*)"its,ite,jts,jte",its,ite,jts,jte
373 write(*,*)"ims,ime,jms,jme",ims,ime,jms,jme
376 call tign_lfn_interpolation(time_now,icl,jcl,ims,ime,jms,jme, &
377 tign,lfn,tff,lff)
378 endif
380 do isubcl=1,ir
381 do jsubcl=1,jr
382 !we use 4 points of each subcell in fuel_left_cell_1
383 ! and in fuel_left_cell_2: bottome left are (1,1), (1,2), (2,1), (2,2)
384 if(fuel_left_method.eq.1)then
385 call fuel_left_cell_1( fuel_left_ff, fire_area_ff, &
386 lff(isubcl,jsubcl),lff(isubcl,jsubcl+1),lff(isubcl+1,jsubcl),lff(isubcl+1,jsubcl+1), &
387 tff(isubcl,jsubcl),tff(isubcl,jsubcl+1),tff(isubcl+1,jsubcl),tff(isubcl+1,jsubcl+1), &
388 time_now, fuel_time(icl,jcl))
389 elseif(fuel_left_method.eq.2)then
390 fire_area_ff=0 ! initialize to something until computed in fuel_left_f(i,j)
391 fuel_left_ff=fuel_left_cell_3( &
392 lff(isubcl,jsubcl),lff(isubcl,jsubcl+1),lff(isubcl+1,jsubcl),lff(isubcl+1,jsubcl+1), &
393 tff(isubcl,jsubcl),tff(isubcl,jsubcl+1),tff(isubcl+1,jsubcl),tff(isubcl+1,jsubcl+1), &
394 time_now, fuel_time(icl,jcl))
395 ! dont forget to change fire_area_ff here
396 else
397 call crash('fuel_left: unknown fuel_left_method')
398 endif
400 ! consistency check
401 if(fire_area_ff.lt.-1e-6 .or. &
402 (fire_area_ff.eq.0. .and. fuel_left_ff.lt.1.-1e-6))then
403 !$OMP CRITICAL(SFIRE_CORE_CRIT)
404 write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
405 ' of refined mesh fuel burnt',1-fuel_left_ff,' fire area',fire_area_ff
406 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
407 call crash(msg)
408 endif
410 helpsum1=helpsum1+fuel_left_ff
411 helpsum2=helpsum2+fire_area_ff
412 enddo
413 enddo
414 fuel_frac(icl,jcl)=helpsum1
415 fire_area(icl,jcl)=helpsum2
416 enddo
417 enddo
422 #ifdef DEBUG_OUT_FUEL_LEFT
423 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fire_area,'fire_area',0)
424 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fuel_frac,'fuel_frac',0)
425 call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fuel_left_f,'fuel_left_f',0)
426 call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fire_area_f,'fire_area_f',0)
427 #endif
429 ! finish the averaging
430 do j=jts,jte
431 do i=its,ite
432 fuel_frac(i,j) = fuel_frac(i,j) /(ir*jr) ! multiply by weight for averaging over coarse cell
433 fire_area(i,j) = fire_area(i,j) /(ir*jr) !
434 enddo
435 enddo
437 ! consistency check after sum
438 fmax=0
439 do j=jts,jte
440 do i=its,ite
441 if(fire_area(i,j).eq.0.)then
442 if(fuel_frac(i,j).lt.1.-1e-6)then
443 !$OMP CRITICAL(SFIRE_CORE_CRIT)
444 write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
445 ' fuel burnt',1-fuel_frac(i,j),' but fire area',fire_area(i,j)
446 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
447 call crash(msg)
448 endif
449 else
450 frat=(1-fuel_frac(i,j))/fire_area(i,j)
451 fmax=max(fmax,frat)
452 endif
453 enddo
454 enddo
455 !$OMP CRITICAL(SFIRE_CORE_CRIT)
456 write(msg,'(a,4i6,a,f10.7)')'fuel_left: tile',its,ite,jts,jte,' max fuel burnt/area',fmax
457 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
458 call message(msg)
459 return
462 end subroutine fuel_left
464 !
465 !************************
466 !
469 !
470 !*************************
471 !Subroutine that is calculating tign and lfn of four endpoints of the subcell
472 ! that is located at isubcl,jsubcl of the cell -(icl,jcl)
473 !
474 subroutine tign_lfn_interpolation(time_now,icl,jcl,ims,ime,jms,jme, &
475 tign,lfn,tff,lff)
476 real, intent(in):: time_now ! not ignited nodes will have tign set to >= time_now
477 integer, intent(in) :: icl,jcl
478 integer, intent(in) :: ims,ime,jms,jme ! memory dimensions of all arrays
479 real, intent(in), dimension(ims:ime,jms:jme)::lfn,tign
480 real, intent(out),dimension(3,3)::tff,lff
482 ! (3,1)-------------(3,2)-------------(3,3)
483 ! | | |
484 ! | (2,1) | (2,2) |
485 ! | | |
486 ! | | |
487 ! | | |
488 ! | | |
489 ! (2,1)--------node-(icl,jcl)---------(2,3)-----------(icl,jcl+1)-------------|
490 ! | sub-node (2,2) |
491 ! | | |
492 ! | (1,1) | (1,2) | each fire mesh cell is decomposed in 4
493 ! | | | tff,lff is computed at the nodes of
494 ! | | | the subcells, numbered (1,1)...(3,3)
495 ! (1,1)-------------(1,2)-------------(1,3)--
496 ! | | |
497 ! | (2,1) | (2,2) |
498 ! | | |
499 ! | | |
500 ! | node |
501 ! -------------(icl-1,jcl------------------
502 ! | | | )
503 !
505 !**********************
507 ! Direct calculation tif and lff, avoiding arrays, just for case ir=jr=2
508 ! Checking whether icl or jcl is on the boundary
511 lff(1,1)=0.25*(lfn(icl-1,jcl-1)+lfn(icl-1,jcl)+lfn(icl,jcl-1)+lfn(icl,jcl))
512 call tign_four_pnts_interp(tign(icl-1,jcl-1),tign(icl-1,jcl),tign(icl,jcl-1), &
513 tign(icl,jcl),lfn(icl-1,jcl-1),lfn(icl-1,jcl), &
514 lfn(icl,jcl-1),lfn(icl,jcl),lff(1,1),tff(1,1),time_now)
516 lff(1,2)=0.5*(lfn(icl-1,jcl)+lfn(icl,jcl))
517 call tign_line_interp(tign(icl-1,jcl),tign(icl,jcl),lfn(icl-1,jcl),lfn(icl,jcl), &
518 lff(1,2),tff(1,2),time_now)
521 lff(1,3)=0.25*(lfn(icl-1,jcl+1)+lfn(icl-1,jcl)+lfn(icl,jcl+1)+lfn(icl,jcl))
522 call tign_four_pnts_interp(tign(icl-1,jcl),tign(icl-1,jcl+1),tign(icl,jcl), &
523 tign(icl,jcl+1),lfn(icl-1,jcl),lfn(icl-1,jcl+1), &
524 lfn(icl,jcl),lfn(icl,jcl+1),lff(1,3),tff(1,3),time_now)
526 lff(2,1)=0.5*(lfn(icl,jcl-1)+lfn(icl,jcl))
527 call tign_line_interp(tign(icl,jcl-1),tign(icl,jcl),lfn(icl,jcl-1),lfn(icl,jcl), &
528 lff(2,1),tff(2,1),time_now)
529 lff(2,2)=lfn(icl,jcl)
530 tff(2,2)=tign(icl,jcl)
532 lff(2,3)=0.5*(lfn(icl,jcl+1)+lfn(icl,jcl))
533 call tign_line_interp(tign(icl,jcl+1),tign(icl,jcl),lfn(icl,jcl+1),lfn(icl,jcl), &
534 lff(2,3),tff(2,3),time_now)
536 lff(3,1)=0.25*(lfn(icl,jcl-1)+lfn(icl,jcl)+lfn(icl+1,jcl-1)+lfn(icl+1,jcl))
537 call tign_four_pnts_interp(tign(icl,jcl-1),tign(icl,jcl),tign(icl+1,jcl-1), &
538 tign(icl+1,jcl),lfn(icl,jcl-1),lfn(icl,jcl), &
539 lfn(icl+1,jcl-1),lfn(icl+1,jcl),lff(3,1),tff(3,1),time_now)
541 lff(3,2)=0.5*(lfn(icl+1,jcl)+lfn(icl,jcl))
542 call tign_line_interp(tign(icl+1,jcl),tign(icl,jcl),lfn(icl+1,jcl),lfn(icl,jcl), &
543 lff(3,2),tff(3,2),time_now)
545 lff(3,3)=0.25*(lfn(icl,jcl)+lfn(icl,jcl+1)+lfn(icl+1,jcl)+lfn(icl+1,jcl+1))
546 call tign_four_pnts_interp(tign(icl,jcl),tign(icl,jcl+1),tign(icl+1,jcl), &
547 tign(icl+1,jcl+1),lfn(icl,jcl),lfn(icl,jcl+1), &
548 lfn(icl+1,jcl),lfn(icl+1,jcl+1),lff(3,3),tff(3,3),time_now)
552 end subroutine tign_lfn_interpolation
556 subroutine tign_line_interp(tign1,tign2,lfn1,lfn2,lfn_subcl,tign_subcl,time_now)
558 real,intent(in) :: tign1,tign2,lfn1,lfn2,lfn_subcl,time_now
559 real,intent(out) :: tign_subcl
560 real :: a,b,c,err
561 err=0.1
563 ! write(*,*)"lfn1,lfn2,tign1,tign2,lfn_subcl,timenow",lfn1,lfn2,tign1,tign2,lfn_subcl,time_now
565 if((lfn1.le.0).AND.(lfn2.le.0)) then
566 tign_subcl=0.5*(tign1+tign2)
567 elseif((lfn1.gt.0).AND.(lfn2.gt.0))then
568 if ((abs(tign1-time_now).gt.err).or.(abs(tign2-time_now).gt.err)) then
569 write(*,*)"lfn1,lfn2,tign1,tign2,timenow",lfn1,lfn2,tign1,tign2,time_now
570 call crash('line_interp: when lfn1(2)>0 their tign1(2) should = time_now')
571 else
572 tign_subcl=time_now;
573 endif
574 elseif(lfn_subcl.gt.0) then
575 if ((abs(tign1-time_now).gt.err).OR.(abs(tign2-time_now).gt.err)) then
576 call crash('tign_line_interp one of tign1,2 should be equal time_now')
577 else
578 tign_subcl=time_now;
579 endif
580 else
581 !lfn_subcl<=0;
582 !case when E is on fire
583 ! tign_subcl~=c*lfn_subcl+time_now;
584 if (lfn1.le.0) then
585 a=lfn1;
586 b=tign1-time_now;
587 elseif (lfn2.le.0) then
588 a=lfn2;
589 b=tign2-time_now;
590 else
591 call crash('tign_line_interp: if E is on fire then one of lfn1 or lfn2 &
592 should be negative');
593 endif
594 if (b.gt.0) then
595 call crash('tign_ should be less than time_now');
596 else
597 c=b/a;
598 tign_subcl=c*lfn_subcl+time_now;
599 endif
600 endif
601 end subroutine tign_line_interp
602 !
603 !************************
604 !
606 subroutine tign_four_pnts_interp(tign1,tign2,tign3,tign4,lfn1,lfn2, &
607 lfn3,lfn4,lfn_subcl,tign_subcl,time_now)
609 real,intent(in) :: tign1,tign2,tign3,tign4
610 real,intent(in) :: lfn1,lfn2,lfn3,lfn4,lfn_subcl,time_now
611 real,intent(out) :: tign_subcl
613 real :: a,b,c,err
614 err=0.0001
616 if((lfn1.le.0).AND.(lfn2.le.0).AND.(lfn3.le.0).AND.(lfn4.le.0)) then
617 tign_subcl=0.25*(tign1+tign2+tign3+tign4)
618 elseif((lfn1.gt.0).AND.(lfn2.gt.0).AND.(lfn3.gt.0).AND.(lfn4.gt.0))then
619 !if ((abs(tign1-time_now).gt.err).OR.(abs(tign2-time_now).gt.err).OR.(abs(tign3-time_now).gt.err).OR.(abs(tign4-time_now).gt.err)) then
620 ! call crash('tign_four_pnts_interp: when lfn1(2,3,4)>0 their tign1(2,3,4) should = time_now')
621 !else
622 tign_subcl=time_now;
623 !endif
624 elseif(lfn_subcl.gt.0) then
625 ! if ((abs(tign1-time_now).gt.err).OR.(abs(tign2-time_now).gt.err).OR.(abs(tign3-time_now).gt.err).OR.(abs(tign4-time_now).gt.err)) then
626 ! call crash('tign_line_interp one of tign1(2,3,4) should be equal time_now')
627 ! else
628 tign_subcl=time_now;
629 ! endif
630 else
631 !lfn_subcl<=0;
632 !case when E is on fire
633 ! tign_subcl~=c*lfn_subcl+time_now;
634 a=0;
635 b=0;
636 if (lfn1.le.0) then
637 ! if (tign1.gt.time_now)
638 ! call crash('tign_four_pnts_interp tign1 should be less then time_now');
639 ! else
640 ! Can not assign to a named constant
641 a=a+lfn1*lfn1;
642 b=b+(tign1-time_now)*lfn1;
643 ! endif
644 endif
645 if (lfn2.le.0) then
646 ! if (tign2.gt.time_now)
647 ! call crash('tign_four_pnts_interp tign2 should be less then time_now');
648 ! else
649 ! Can not assign to a named constant
650 a=a+lfn2*lfn2;
651 b=b+(tign2-time_now)*lfn2;
652 ! endif
653 endif
654 if (lfn3.le.0) then
655 ! if (tign3.gt.time_now)
657 ! call crash('tign_four_pnts_interp tign3 should be less then time_now');
658 ! else
659 ! Can not assign to a named constant
660 a=a+lfn3*lfn3;
661 b=b+(tign3-time_now)*lfn3;
662 ! endif
663 endif
664 if (lfn4.le.0) then
665 ! if (tign4.gt.time_now)
666 ! call crash('tign_four_pnts_interp tign4 should be less then time_now');
667 ! else
668 ! Can not assign to a named constant
669 a=a+lfn4*lfn4;
670 b=b+(tign4-time_now)*lfn4;
671 ! endif
672 endif
673 if ((abs(a).lt.err).or.(b.lt.0)) then
674 call crash('tign_four_pnts_interp: if E is on fire then one of cells &
675 should be on fire or tign_ should be less than time_now')
676 else
677 c=b/a;
678 tign_subcl=c*lfn_subcl+time_now;
679 endif
680 endif
682 end subroutine tign_four_pnts_interp
685 !
686 !************************
687 !
690 subroutine fuel_left_cell_1( fuel_frac_left, fire_frac_area, &
691 lfn00,lfn01,lfn10,lfn11, &
692 tign00,tign01,tign10,tign11,&
693 time_now, fuel_time_cell)
694 !*** purpose: compute the fuel fraction left in one cell
695 implicit none
696 !*** arguments
697 real, intent(out):: fuel_frac_left, fire_frac_area !
698 real, intent(in)::lfn00,lfn01,lfn10,lfn11 ! level set function at 4 corners of the cell
699 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the 4 corners of the cell
700 real, intent(in)::time_now ! the time now
701 real, intent(in)::fuel_time_cell ! time to burns off to 1/e
702 !*** Description
703 ! The area burning is given by the condition L <= 0, where the function P is
704 ! interpolated from lfn(i,j)
705 !
706 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
707 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
708 ! when lfn(i,j)=0.
709 !
710 ! The function computes an approxmation of the integral
711 !
712 !
713 ! /\
714 ! |
715 ! fuel_frac_left = 1 - | 1 - exp(-T(x,y)/fuel_time_cell)) dxdy
716 ! |
717 ! \/
718 ! 0<x<1
719 ! 0<y<1
720 ! L(x,y)<=0
721 !
722 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
723 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
724 ! so dx=1 and dy=1 assumed.
725 !
726 ! Example:
727 !
728 ! lfn<0 lfn>0
729 ! (0,1) -----O--(1,1) O = points on the fireline, T=tnow
730 ! | \ | A = the burning area for computing
731 ! | \| fuel_frac(i,j)
732 ! | A O
733 ! | |
734 ! | |
735 ! (0,0)---------(1,0)
736 ! lfn<0 lfn<0
737 !
738 ! Approximations allowed:
739 ! The fireline can be approximated by straight line(s).
740 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
741 !
742 ! Requirements:
743 ! 1. The output should be a continuous function of the arrays lfn and
744 ! tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
745 ! 2. The output should be invariant to the symmetries of the input in each cell.
746 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
747 ! 4. The result should be at least 1st order accurate in the sense that it is
748 ! exact if the time from ignition is a linear function.
749 !
750 ! If time from ignition is approximated by polynomial in the burnt
751 ! region of the cell, this is integral of polynomial times exponential
752 ! over a polygon, which can be computed exactly.
753 !
754 ! Requirement 4 is particularly important when there is a significant decrease
755 ! of the fuel fraction behind the fireline on the mesh scale, because the
756 ! rate of fuel decrease right behind the fireline is much larger
757 ! (exponential...). This will happen when
758 !
759 ! change of time from ignition within one mesh cell / fuel_time_cell is not << 1
760 !
761 ! This is the same as
762 !
763 ! mesh cell size
764 ! X = ------------------------- is not << 1
765 ! fireline speed * fuel_time_cell
766 !
767 !
768 ! When X is large then the fuel burnt in one timestep in the cell is
769 ! approximately proportional to length of fireline in that cell.
770 !
771 ! When X is small then the fuel burnt in one timestep in the cell is
772 ! approximately proportional to the area of the burning region.
773 !
775 !*** calls
776 intrinsic tiny
778 !*** local
779 real::ps,aps,area,ta,out
780 real::t00,t01,t10,t11
781 real,parameter::safe=tiny(aps)
782 character(len=128)::msg
784 ! the following algorithm is a very crude approximation
786 ! minus time since ignition, 0 if no ignition yet
787 ! it is possible to have 0 in fire region when ignitin time falls in
788 ! inside the time step because lfn is updated at the beginning of the time step
790 t00=tign00-time_now
791 if(lfn00>0. .or. t00>0.)t00=0.
792 t01=tign01-time_now
793 if(lfn01>0. .or. t01>0.)t01=0.
794 t10=tign10-time_now
795 if(lfn10>0. .or. t10>0.)t10=0.
796 t11=tign11-time_now
797 if(lfn11>0. .or. t11>0.)t11=0.
799 ! approximate burning area, between 0 and 1
800 ps = lfn00+lfn01+lfn10+lfn11
801 aps = abs(lfn00)+abs(lfn01)+abs(lfn10)+abs(lfn11)
802 aps=max(aps,safe)
803 area =(-ps/aps+1.)/2.
804 area = max(area,0.) ! make sure area is between 0 and 1
805 area = min(area,1.)
807 ! average negative time since ignition
808 ta=0.25*(t00+t01+t10+t11)
810 ! exp decay in the burning area
811 out=1.
812 !if(area>0.)out=1. - area*(1. - exp(ta/fuel_time_cell))
813 if(area>0)out=area*exp(ta/fuel_time_cell) + (1. - area)
815 if(out>1.)then
816 !$OMP CRITICAL(SFIRE_CORE_CRIT)
817 write(msg,*)'out=',out,'>1 area=',area,' ta=',ta
818 call message(msg)
819 write(msg,*)'tign=',tign00,tign01,tign10,tign11,' time_now=',time_now
820 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
821 call message(msg)
822 !call message('WARNING: fuel_left_cell_1: fuel fraction > 1')
823 call crash('fuel_left_cell_1: fuel fraction > 1')
824 endif
826 !out = max(out,0.) ! make sure out is between 0 and 1
827 !out = min(out,1.)
829 fuel_frac_left = out
830 fire_frac_area = area
832 end subroutine fuel_left_cell_1
834 !
835 !****************************************
836 !
837 ! function calculation fuel_frac made by Volodymyr Kondratenko on the base of
838 ! the code made by Jan Mandel and Minjeong
840 real function fuel_left_cell_3( &
841 lfn00,lfn01,lfn10,lfn11, &
842 tign00,tign01,tign10,tign11,&
843 time_now, fuel_time_cell)
844 !*** purpose: compute the fuel fraction left in one cell
845 implicit none
846 !*** arguments
847 real, intent(in)::lfn00,lfn01,lfn10,lfn11 ! level set function at 4 corners of the cell
848 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the 4 corners of the cell
849 real, intent(in)::time_now ! the time now
850 real, intent(in)::fuel_time_cell ! time to burns off to 1/e
851 !*** Description
852 ! The area burning is given by the condition L <= 0, where the function P is
853 ! interpolated from lfn(i,j)
854 !
855 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
856 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
857 ! when lfn(i,j)=0.
858 !
859 ! The function computes an approxmation of the integral
860 !
861 !
862 ! /\
863 ! |
864 ! fuel_frac_left = 1 - | 1 - exp(-T(x,y)/fuel_time_cell)) dxdy
865 ! |
866 ! \/
867 ! 0<x<1
868 ! 0<y<1
869 ! L(x,y)<=0
870 !
871 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
872 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
873 ! so dx=1 and dy=1 assumed.
874 !
875 ! Example:
876 !
877 ! lfn<0 lfn>0
878 ! (0,1) -----O--(1,1) O = points on the fireline, T=tnow
879 ! | \ | A = the burning area for computing
880 ! | \| fuel_frac(i,j)
881 ! | A O
882 ! | |
883 ! | |
884 ! (0,0)---------(1,0)
885 ! lfn<0 lfn<0
886 !
887 ! Approximations allowed:
888 ! The fireline can be approximated by straight line(s).
889 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
890 !
891 ! Requirements:
892 ! 1. The output should be a continuous function of the arrays lfn and
893 ! tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
894 ! 2. The output should be invariant to the symmetries of the input in each cell.
895 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
896 ! 4. The result should be at least 1st order accurate in the sense that it is
897 ! exact if the time from ignition is a linear function.
898 !
899 ! If time from ignition is approximated by polynomial in the burnt
900 ! region of the cell, this is integral of polynomial times exponential
901 ! over a polygon, which can be computed exactly.
902 !
903 ! Requirement 4 is particularly important when there is a significant decrease
904 ! of the fuel fraction behind the fireline on the mesh scale, because the
905 ! rate of fuel decrease right behind the fireline is much larger
906 ! (exponential...). This will happen when
907 !
908 ! change of time from ignition within one mesh cell * fuel speed is not << 1
909 !
910 ! This is the same as
911 !
912 ! mesh cell size*fuel_speed
913 ! ------------------------- is not << 1
914 ! fireline speed
915 !
916 !
917 ! When X is large then the fuel burnt in one timestep in the cell is
918 ! approximately proportional to length of fireline in that cell.
919 !
920 ! When X is small then the fuel burnt in one timestep in the cell is
921 ! approximately proportional to the area of the burning region.
923 !#ifndef FUEL_LEFT
924 !call crash('fuel_left_cell_3: not implemented, please use fire_fuel_left_method=1')
925 !fuel_left_cell_3=0. ! to avoid compiler warning about value not set
926 !end function fuel_left_cell_3
927 !#else
929 !*** calls
930 intrinsic tiny
932 !*** local
933 real::ps,aps,area,ta,out
934 real::t00,t01,t10,t11
935 real,parameter::safe=tiny(aps)
936 character(len=128)::msg
937 real::dx,dy ! mesh sizes
938 !*** local
939 integer::i,j,k
941 ! least squares
942 integer::mmax,nb,nmax,pmax,nin,nout
943 parameter(mmax=3,nb=64,nmax=8,pmax=8)
944 integer lda, ldb, lwork, info
945 parameter (lda=nmax, ldb=nmax, lwork=mmax+nmax+nb*(nmax+pmax))
946 integer n,m,p
947 real,dimension(lda,mmax):: mA
948 real,dimension(nmax):: vecD
949 real,dimension(lwork):: WORK
950 real,dimension(pmax):: vecY
951 real,dimension(mmax):: vecX
952 real,dimension(ldb,pmax)::mB
953 real,dimension(mmax)::u
955 real::tweight,tdist
956 integer::kk,ll,ss
957 real::rnorm
958 real,dimension(8,2)::xylist,xytlist
959 real,dimension(8)::tlist,llist,xt
960 real,dimension(5)::xx,yy
961 real,dimension(5)::lfn,tign
963 integer:: npoint
964 real::tt,x0,y0,xts,xte,yts,yte,xt1,xt2
965 real::lfn0,lfn1,dist,nr,s,errQ,ae,ce,ceae,a0,a1,a2,d,cet
966 real::s1,s2,s3
967 real::upper,lower,ah,ch,aa,cc,aupp,cupp,alow,clow
968 real,dimension(2,2)::mQ
969 real,dimension(2)::ut
971 !calls
972 intrinsic epsilon
974 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
976 !!!! For finite differences by VK
977 real::tign_middle,dt_dx,dt_dy,lfn_middle,a,b,c
979 ! external functions
980 real::snrm2
981 double precision :: dnrm2
982 external dnrm2
983 external snrm2
984 ! external subroutines
985 external sggglm
986 external dggglm
987 !executable statements
989 ! a very crude approximation - replace by a better code
990 ! call check_mesh_2dim(ids,ide+1,jds,jde+1,ims,ime,jms,jme)
992 dx=1
993 dy=1
994 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
995 !comment - changed tign-time_now to time_now-tign
996 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
997 t00=time_now-tign00
998 if(lfn00>=0. .or. t00<0.)t00=0.
999 t01=time_now-tign01
1000 if(lfn01>=0. .or. t01<0.)t01=0.
1001 t10=time_now-tign10
1002 if(lfn10>=0. .or. t10<0.)t10=0.
1003 t11=time_now-tign11
1004 if(lfn11>=0. .or. t11<0.)t11=0.
1006 !*** case0 Do nothing
1007 if ( lfn00>=0 .and. lfn10>=0 .and. lfn01>=0 .and. lfn11>=0 ) then
1008 print*,"Case 0"
1009 out = 1.0 ! fuel_left, no burning
1010 !*** case4 all four coners are burning
1011 else if (lfn00<=0 .and. lfn10<=0 .and. lfn01<=0 .and. lfn11<=0) then
1012 !!!!!!!!!!!
1014 print*,"Case 4"
1015 ! All burning
1016 ! T=u(1)*x+u(2)*y+u(3)
1017 ! t(0,0)=tign(1,1)
1018 ! t(0,fd(2))=tign(1,2)
1019 ! t(fd(1),0)=tign(2,1)
1020 ! t(fd(1),fd(2))=tign(2,2)
1021 ! t(g/2,h/2)=sum(tign(i,i))/4
1022 ! dt/dx=(1/2h)*(t10-t00+t11-t01)
1023 ! dt/dy=(1/2h)*(t01-t00+t11-t10)
1024 ! approximate T(x,y)=u(1)*x+u(2)*y+u(3) Using finite differences
1025 ! t(x,y)=t(h/2,h/2)+(x-h/2)*dt/dx+(y-h/2)*dt/dy
1026 ! u(1)=dt/dx
1027 ! u(2)=dt/dy
1028 ! u(3)=t(h/2,h/2)-h/2(dt/dx+dt/dy)
1030 tign_middle=(t00+t01+t10+t11)/4
1031 write(*,*)"tign_middle",tign_middle
1032 write(*,*)"t00,t01,t10,t11",t00,t01,t10,t11
1034 ! since mesh_size is 1 we replace fd(1) and fd(2) by 1
1035 dt_dx=(t10-t00+t11-t01)/2
1036 dt_dy=(t01-t00+t11-t10)/2
1038 write(*,*)"dt_dx,dt_dy",dt_dx,dt_dy
1039 write(*,*)"dx,dy",dx,dy
1041 u(1)=dt_dx
1042 u(2)=dt_dy
1043 u(3)=tign_middle-(dt_dx+dt_dy)/2
1045 Write(*,*)"u=",u(1),u(2),u(3)
1047 ! integrate
1048 u(1)=-u(1)/fuel_time_cell
1049 u(2)=-u(2)/fuel_time_cell
1050 u(3)=-u(3)/fuel_time_cell
1051 write(*,*)"u/fuel_time_cell",u(1),u(2),u(3)
1053 write(*,*)"intexp(u(1)),intexp(u(2))",intexp(u(1)*dx),intexp(u(2)*dy)
1055 !fuel_burn(i,j)=1-exp(u(3))*intexp(u(1)*dx)*intexp(u(2)*dy)
1056 s1=u(1)
1057 s2=u(2)
1058 out=1-exp(u(3))*intexp(s1)*intexp(s2)
1059 print *,intexp(s1),intexp(s2),out
1060 if ( out<0 .or. out>1.0 ) then
1061 print *,'case4, out should be between 0 and 1'
1062 end if
1063 !*** case 1,2,3
1065 else
1066 print*,"Case 123"
1067 ! set xx, yy for the coner points
1068 ! move these values out of i and j loop to speed up
1069 xx(1) = -0.5
1070 xx(2) = 0.5
1071 xx(3) = 0.5
1072 xx(4) = -0.5
1073 xx(5) = -0.5
1074 yy(1) = -0.5
1075 yy(2) = -0.5
1076 yy(3) = 0.5
1077 yy(4) = 0.5
1078 yy(5) = -0.5
1079 lfn(1)=lfn00
1080 lfn(2)=lfn10
1081 lfn(3)=lfn11
1082 lfn(4)=lfn01
1083 lfn(5)=lfn00
1084 tign(1)=t00
1085 tign(2)=t10
1086 tign(3)=t11
1087 tign(4)=t01
1088 tign(5)=t00
1089 npoint = 0 ! number of points in polygon
1091 do k=1,4
1092 lfn0=lfn(k )
1093 lfn1=lfn(k+1)
1094 if ( lfn0 <= 0.0 ) then
1095 npoint = npoint + 1
1096 xylist(npoint,1)=xx(k)
1097 xylist(npoint,2)=yy(k)
1098 tlist(npoint)=-tign(k)
1099 llist(npoint)=lfn0
1100 end if
1101 if ( lfn0*lfn1 < 0 ) then
1102 npoint = npoint + 1
1103 tt=lfn0/(lfn0-lfn1)
1104 x0=xx(k)+( xx(k+1)-xx(k) )*tt
1105 y0=yy(k)+( yy(k+1)-yy(k) )*tt
1106 xylist(npoint,1)=x0
1107 xylist(npoint,2)=y0
1108 tlist(npoint)=0 ! on fireline
1109 llist(npoint)=0
1110 end if
1111 end do
1113 ! make the list circular
1114 tlist(npoint+1)=tlist(1)
1115 llist(npoint+1)=llist(1)
1116 xylist(npoint+1,1)=xylist(1,1)
1117 xylist(npoint+1,2)=xylist(1,2)
1121 !print *,'after LS in case3'pproximation of the plane for lfn using finite
1122 !differences
1123 ! approximate L(x,y)=u(1)*x+u(2)*y+u(3)
1124 lfn_middle=(lfn00+lfn01+lfn10+lfn11)/4
1125 print *,"lfn_middle",lfn_middle
1126 dt_dx=(lfn10-lfn00+lfn11-lfn01)/2
1127 dt_dy=(lfn01-lfn00+lfn11-lfn10)/2
1128 print *,"dt_dx,dt_dy",dt_dx,dt_dy
1129 u(1)=dt_dx
1130 u(2)=dt_dy
1131 u(3)=lfn_middle-(dt_dx+dt_dy)/2
1132 ! finding the coefficient c, reminder we work over one subcell only
1133 ! T(x,y)=c*L(x,y)+time_now
1134 write(*,*)"vector u before c",u(1),u(2),u(3)
1135 a=0
1136 b=0
1138 if (lfn00 <= 0) then
1139 a=a+lfn00*lfn00
1140 if (t00 < 0) then
1141 call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
1142 else
1143 b=b+t00*lfn00
1144 end if
1145 write(*,*)"t00,time_now,lfn00",t00,time_now,lfn00
1146 write(*,*)"a,b",a,b
1147 end if
1150 if (lfn01 <= 0) then
1151 a=a+lfn01*lfn01
1152 if (t01< 0) then
1153 call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
1154 else
1155 b=b+t01*lfn01
1156 end if
1157 write(*,*)"a,b",a,b
1158 end if
1161 if (lfn10<=0) then
1162 a=a+lfn10*lfn10
1163 if (t10<0) then
1164 call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
1165 else
1166 b=b+t10*lfn10
1167 end if
1168 write(*,*)"a,b",a,b
1169 end if
1171 if (lfn11<=0) then
1172 a=a+lfn11*lfn11
1173 if (t11<0) then
1174 call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
1175 else
1176 b=b+t11*lfn11
1177 end if
1178 write(*,*)"a,b",a,b
1179 end if
1183 if (a==0) then
1184 call crash('fuel_burnt_fd: if c is on fire then one of cells should be on fire')
1185 end if
1186 c=b/a
1187 u(1)=u(1)*c
1188 u(2)=u(2)*c
1189 u(3)=u(3)*c
1191 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1193 !print *,'vecX from LS',vecX
1194 !print *,'tign inputed',tign00,tign10,tign11,tign01
1196 write(*,*)"vector u",u(1),u(2),u(3)
1197 ! rotate to gradient on x only
1198 nr = sqrt(u(1)**2+u(2)**2)
1199 if(.not.nr.gt.eps)then
1200 out=1.
1201 goto 900
1202 endif
1203 c = u(1)/nr
1204 s = u(2)/nr
1205 mQ(1,1)=c
1206 mQ(1,2)=s
1207 mQ(2,1)=-s
1208 mQ(2,2)=c
1209 ! mat vec multiplication
1210 call matvec(mQ,2,2,u,3,ut,2,2,2)
1211 errQ = ut(2) ! should be zero
1212 ae = -ut(1)/fuel_time_cell
1213 ce = -u(3)/fuel_time_cell
1214 cet=ce!keep ce
1215 call matmatp(xylist,8,2,mQ,2,2,xytlist,8,2,npoint+1,2,2)
1216 call sortxt( xytlist, 8,2, xt,8,npoint )
1217 out=0.0
1218 aupp=0.0
1219 cupp=0.0
1220 alow=0.0
1221 clow=0.0
1222 do k=1,npoint-1
1223 xt1=xt(k)
1224 xt2=xt(k+1)
1225 upper=0
1226 lower=0
1227 ah=0
1228 ch=0
1229 if ( xt2-xt1 > eps*100 ) then
1231 do ss=1,npoint
1232 xts=xytlist(ss,1)
1233 yts=xytlist(ss,2)
1234 xte=xytlist(ss+1,1)
1235 yte=xytlist(ss+1,2)
1237 if ( (xts>xt1 .and. xte>xt1) .or. &
1238 (xts<xt2 .and. xte<xt2) ) then
1239 aa = 0 ! do nothing
1240 cc = 0
1241 else
1242 aa = (yts-yte)/(xts-xte)
1243 cc = (xts*yte-xte*yts)/(xts-xte)
1244 if (xte<xts) then
1245 aupp = aa
1246 cupp = cc
1247 ah=ah+aa
1248 ch=ch+cc
1249 upper=upper+1
1250 else
1251 alow = aa
1252 clow = cc
1253 lower=lower+1
1254 end if
1255 end if!(xts>xt1 .and. xte>xt1)
1256 end do ! ss
1257 ce=cet !use stored ce
1258 if (ae*xt1+ce > 0 ) then
1259 ce=ce-(ae*xt1+ce)!shift small amounts exp(-**)
1260 end if
1261 if (ae*xt2+ce > 0) then
1262 ce=ce-(ae*xt2+ce)
1263 end if
1265 ah = aupp-alow
1266 ch = cupp-clow
1267 ! integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2
1268 ! numerically sound for ae->0, ae -> infty
1269 ! this can be important for different model scales
1270 ! esp. if someone runs the model in single precision!!
1271 ! s1=int((ah*x+ch),x,xt1,xt2)
1272 s1 = (xt2-xt1)*((1./2.)*ah*(xt2+xt1)+ch)
1273 ! s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)
1274 ceae=ce/ae;
1275 s2 = -ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1))
1276 ! s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)
1277 ! s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)
1278 ! expand in Taylor series around ae=0
1279 ! collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)
1280 ! =(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
1281 ! + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae
1282 ! + 1/2*xt1^2-1/2*xt2^2
1283 !
1284 ! coefficient at ae^2 in the expansion, after some algebra
1285 a2=(xt1-xt2)*((1./4.)*(xt1+xt2)*ceae**2+(1./3.)* &
1286 (xt1**2+xt1*xt2+xt2**2)*ceae+(1./8.)* &
1287 (xt1**3+xt1*(xt2**2)+xt1**2*xt2+xt2**3))
1288 d=(ae**4)*a2
1290 if (abs(d)>eps) then
1291 ! since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae
1292 ! for ae, ce -> 0 rounding error approx eps/ae^2
1293 s3=( exp(ae*(xt1+ceae))*(ae*xt1-1)-&
1294 exp(ae*(xt2+ceae))*(ae*xt2-1) )/(ae**2)
1296 !we do not worry about rounding as xt1 -> xt2, then s3 -> 0
1297 else
1298 ! coefficient at ae^1 in the expansion
1299 a1=(xt1-xt2)*((1./2.)*ceae*(xt1+xt2)+(1./3.)*&
1300 (xt1**2+xt1*xt2+xt2**2))
1301 ! coefficient at ae^0 in the expansion for ae->0
1302 a0=(1./2.)*(xt1-xt2)*(xt1+xt2)
1303 s3=a0+a1*ae+a2*ae**2; ! approximate the integral
1304 end if
1306 s3=ah*s3
1307 out=out+s1+s2+s3
1308 if(out<0 .or. out>1) then
1309 print *,':fuel_fraction should be between 0 and 1'
1310 !print *, 'eps= ', eps
1311 !print *, 'xt1= ', xt1, 'xt2= ', xt2
1312 !print *,'ae= ',ae,'ce= ',ce,'ah= ',ah,'ch= ',ch
1313 !print *,'a0= ', a0,'a1= ', a1,'a2= ', a2
1314 print *,'s1= ', s1,'s2= ', s2,'s3= ', s3
1315 print *,':fuel_fraction =',out
1316 end if!print
1318 end if
1319 end do ! k
1321 out=1-out !fuel_left
1322 end if ! if case0, elseif case4 ,else case123
1324 write(*,*)"out=",out
1326 900 continue
1327 if(out>1. .or. out<0.)call crash('fuel_left_cell_3: fuel fraction out of bounds [0,1]')
1328 fuel_left_cell_3 = out
1329 end function fuel_left_cell_3
1336 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1337 real function intexp(ab)
1338 implicit none
1339 real::ab
1340 !calls
1341 intrinsic epsilon
1343 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
1345 !eps = 2.2204*(10.0**(-8))!from matlab
1346 if ( eps < abs(ab)**3/6. ) then
1347 intexp=(exp(ab)-1)/ab
1348 else
1349 intexp=1+ab/2.
1350 end if
1351 end function intexp
1352 !
1353 !****************************************
1354 !
1355 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1356 !
1357 !****************************************
1358 !
1361 !
1362 !****************************************
1366 !#ifdef FUEL_LEFT_2
1368 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1369 subroutine sortxt(xytlist,nrow,ncolumn,xt,nxt,nvec)
1370 implicit none
1371 integer::nrow,ncolumn,nxt,nvec
1372 real,dimension(nrow,ncolumn)::xytlist
1373 real,dimension(nxt)::xt
1375 integer::i,j
1376 real::temp
1378 do i=1,nvec
1379 xt(i)=xytlist(i,1)
1380 end do
1382 do i=1,nvec-1
1383 do j=i+1,nvec
1384 if ( xt(i) > xt(j) ) then
1385 temp = xt(i)
1386 xt(i)=xt(j)
1387 xt(j)=temp
1388 end if
1389 end do
1390 end do
1392 end subroutine !sortxt
1393 !
1394 !****************************************
1395 !
1396 subroutine matvec(A,m,n,V,nv,out,nout,nrow,ncolumn)
1397 implicit none
1398 integer::m,n,nv,nout,nrow,ncolumn
1399 real,dimension(m,n)::A ! allocated m by n
1400 real,dimension(nv)::V ! allocated nv
1401 real,dimension(nout)::out! allocated nout
1403 integer::i,j
1405 do i=1,nrow
1406 out(i)=0.0
1407 do j=1,ncolumn
1408 out(i)=out(i)+A(i,j)*V(j)
1409 end do
1410 end do
1411 end subroutine
1412 !
1413 !****************************************
1414 !
1415 subroutine matmatp(A,mA,nA,B,mB,nB,C,mC,nC,nrow,ncolumn,nP)
1416 implicit none
1417 integer::mA,nA,mB,nB,mC,nC,nrow,ncolumn,nP
1418 real,dimension(mA,nA)::A ! allocated m by n
1419 real,dimension(mB,nB)::B ! allocated m by n
1420 real,dimension(mC,nC)::C ! allocated m by n
1421 integer::i,j,k
1422 do i=1,nrow
1423 do j=1,ncolumn
1424 C(i,j)=0.0
1425 do k=1,nP
1426 C(i,j)=C(i,j)+A(i,k)*B(j,k) ! B'
1427 end do
1428 end do
1429 end do
1430 end subroutine
1432 !
1433 !****************************************
1434 !
1435 !#endif
1437 subroutine prop_ls( id, & ! for debug
1438 ids,ide,jds,jde, & ! domain dims
1439 ims,ime,jms,jme, & ! memory dims
1440 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1441 its,ite,jts,jte, & ! tile dims
1442 ts,dt,dx,dy, & ! scalars in
1443 tbound, & ! scalars out
1444 lfn_in,lfn_out,tign,ros, & ! arrays inout
1445 fp &
1446 )
1447 implicit none
1449 !*** purpose: advance level function in time
1451 ! Jan Mandel August 2007 - February 2008
1453 !*** description
1454 !
1455 ! Propagation of closed curve by a level function method. The level function
1456 ! lfn is defined by its values at the nodes of a rectangular grid.
1457 ! The area where lfn < 0 is inside the curve. The curve is
1458 ! described implicitly by lfn=0. Points where the curve intersects gridlines
1459 ! can be found by linear interpolation from nodes.
1460 !
1461 ! The level function is advanced from time ts to time ts + dt.
1462 !
1463 ! The level function should be initialized to (an approximation of) the signed
1464 ! distance from the curve. If the initial curve is a circle, the initial level
1465 ! function is simply the distance from the center minus the radius.
1466 !
1467 ! The curve moves outside with speed given by function speed_func.
1468 !
1469 ! Method: Godunov/ENO method for the normal motion. The timestep is checked for
1470 ! CFL condition. For a straight segment in a constant field and locally linear
1471 ! level function, the method reduces to the exact normal motion. The advantage of
1472 ! the level set method is that it treats automatically special cases such as
1473 ! the curve approaching itself and merging components of the area inside the curve.
1474 !
1475 ! Based on S. Osher and R. Fedkiw, Level set methods and dynamic implicit surfaces,
1476 ! Springer, 2003, Sec. 6.4, as implemented in toolboxLS for Matlab by
1477 ! I. Mitchell, A toolbox of Level Set Methods (Version 1.1), TR-2007-11,
1478 ! Dept. Computer Science, University of British Columbia, 2007
1479 ! http://www.cs.ubc.ca/\~mitchell/Toolbo\LS
1480 !
1482 !*** arguments
1484 ! id in unique identification for prints and dumps
1485 ! ids,ide,jds,jde in domain dimensions
1486 ! ims,ime,jms,jme in memory dimensions
1487 ! its,ite,jts,jte in tile dimensions
1488 ! ts in start time
1489 ! dt in time step
1490 ! dx,dy in grid spacing
1491 ! tbound out bound on stable time step from CFL condition, if tbound>=dt then OK
1492 ! lfn_in,lfn_out inout,out the level set function at nodes
1493 ! tign inout the ignition time at nodes
1495 ! The dimensions are cell-based, the nodal value is associated with the south-west corner.
1496 ! The whole computation is on domain indices ids:ide+1,jds:jde+1.
1497 !
1498 ! The region where new lfn and tign are computed is the tile its:ite,jts:jte
1499 ! except when the tile is at domain upper boundary, an extra band of points is added:
1500 ! if ite=ide then region goes up to ite+1, if jte=jde then region goes up to jte+1.
1502 ! The time step requires values from 2 rows of nodes beyond the region except when at the
1503 ! domain boundary one-sided derivatives are used. This is implemented by extending the input
1504 ! beyond the domain boundary so sufficient memory bounds must be allocated.
1505 ! The update on all tiles can be done in parallel. To avoid the race condition (different regions
1506 ! of the same array updated by different threads), the in and out versions of the
1507 ! arrays lft and tign are distinct. If the time step dt is larger
1508 ! that the returned tbound, the routine should be called again with timestep td<=tbound, and then
1509 ! having distinct inputs and outputs comes handy.
1511 !*** calls
1512 !
1513 ! tend_ls
1514 !
1516 integer,intent(in)::id,ims,ime,jms,jme,ids,ide,jds,jde,its,ite,jts,jte,ips,ipe,jps,jpe
1517 real,dimension(ims:ime,jms:jme),intent(inout)::lfn_in,tign
1518 real,dimension(ims:ime,jms:jme),intent(out)::lfn_out,ros
1519 real,intent(in)::dx,dy,ts,dt
1520 real,intent(out)::tbound
1521 type(fire_params),intent(in)::fp
1523 !*** local
1524 ! arrays
1525 #define IMTS its-1
1526 #define IMTE ite+1
1527 #define JMTS jts-1
1528 #define JMTE jte+1
1529 real,dimension(IMTS:IMTE,JMTS:JMTE):: tend, lfn1 ! region-sized with halo
1530 ! scalars
1531 real::grad2,rr,tbound2,a,a1 ! a=0 euler, a=0.5 heun
1533 real::gradx,grady,aspeed,err,aerr,time_now
1534 integer::ihs,ihe,jhs,jhe
1535 integer::ihs2,ihe2,jhs2,jhe2
1536 integer::i,j,its1,ite1,jts1,jte1,k,kk,id1
1537 character(len=128)msg
1538 integer::nfirenodes,nfireline
1539 real::sum_err,min_err,max_err,sum_aerr,min_aerr,max_aerr
1541 ! constants
1542 integer,parameter :: mstep=1000, printl=1
1543 real, parameter:: zero=0.,one=1.,eps=epsilon(zero),tol=100*eps, &
1544 safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
1546 ! f90 intrinsic function
1548 intrinsic max,min,sqrt,nint,epsilon,tiny,huge
1550 !*** executable
1552 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1553 write(msg,'(a8,i5,a6,i5,3(a1,i5))')'prop_ls:',id,' tile ',its,':',ite,',',jts,':',jte
1554 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1555 call message(msg)
1557 a=fire_back_weight ! from module_fr_sfire_util
1558 a1=1. - a
1560 ! tend = F(lfn)
1562 ihs2=max(its-2,ids) ! need lfn two beyond the tile but not outside the domain
1563 ihe2=min(ite+2,ide)
1564 jhs2=max(jts-2,jds)
1565 jhe2=min(jte+2,jde)
1567 ihs=max(its-1,ids) ! compute tend one beyond the tile but not outside the domain
1568 ihe=min(ite+1,ide)
1569 jhs=max(jts-1,jds)
1570 jhe=min(jte+1,jde)
1572 #ifdef DEBUG_OUT
1573 call write_array_m(ihs,ihe,jhs,jhe,ims,ime,jms,jme,lfn_in,'lfn_in',id)
1574 #endif
1576 ! check array dimensions
1577 call check_mesh_2dim(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme)
1578 call print_2d_stats(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme, &
1579 lfn_in,'prop_ls: lfn in')
1581 ! NOTE: tend_ls will extrapolate to one node strip at domain boundaries
1582 ! so that it can compute gradient at domain boundaries.
1583 ! To avoid copying, lfn_in is declared inout.
1584 ! At tile boundaries that are not domain boundaries values of lfn_in two nodes
1585 ! outside of the tile are needed.
1586 id1 = id ! for debug prints
1587 if(id1.ne.0)id1=id1+1000
1588 call tend_ls( id1, &
1589 ims,ime,jms,jme, & ! memory dims for lfn_in
1590 IMTS,IMTE,JMTS,JMTE, & ! memory dims for tend
1591 ids,ide,jds,jde, & ! domain dims - where lfn exists
1592 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1593 ihs,ihe,jhs,jhe, & ! where tend computed
1594 ims,ime,jms,jme, & ! memory dims for ros
1595 its,ite,jts,jte, & ! tile dims - where to set ros
1596 ts,dt,dx,dy, & ! scalars in
1597 lfn_in, & ! arrays in
1598 tbound, & ! scalars out
1599 tend, ros, & ! arrays out
1600 fp & ! params
1601 )
1603 #ifdef DEBUG_OUT
1604 call write_array_m(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE,tend,'tend1',id)
1605 #endif
1607 ! Euler method, the half-step, same region as ted
1608 do j=jhs,jhe
1609 do i=ihs,ihe
1610 lfn1(i,j) = lfn_in(i,j) + dt*tend(i,j)
1611 enddo
1612 enddo
1614 call print_2d_stats(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE, &
1615 lfn1,'prop_ls: lfn1')
1616 ! tend = F(lfn1) on the tile (not beyond)
1618 if(id1.ne.0)id1=id1+1000
1619 call tend_ls( id1,&
1620 IMTS,IMTE,JMTS,JMTE, & ! memory dims for lfn
1621 IMTS,IMTE,JMTS,JMTE, & ! memory dims for tend
1622 ids,ide,jds,jde, & ! domain dims - where lfn exists
1623 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1624 its,ite,jts,jte, & ! tile dims - where is tend computed
1625 ims,ime,jms,jme, & ! memory dims for ros
1626 its,ite,jts,jte, & ! tile dims - where is ros set
1627 ts+dt,dt,dx,dy, & ! scalars in
1628 lfn1, & ! arrays in
1629 tbound2, & ! scalars out
1630 tend,ros, & ! arrays out
1631 fp &
1632 )
1634 #ifdef DEBUG_OUT
1635 call write_array_m(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'tend2',id)
1636 #endif
1638 call print_2d_stats(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'prop_ls: tend2')
1640 tbound=min(tbound,tbound2)
1642 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1643 write(msg,'(a,f10.2,4(a,f7.2))')'prop_ls: time',ts,' dt=',dt,' bound',min(tbound,999.99), &
1644 ' dx=',dx,' dy=',dy
1645 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1646 call message(msg)
1647 if(dt>tbound)then
1648 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1649 write(msg,'(2(a,f10.2))')'prop_ls: WARNING: time step ',dt, &
1650 ' > bound =',tbound
1651 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1652 call message(msg)
1653 endif
1655 ! combine lfn1 and lfn_in + dt*tend -> lfn_out
1657 do j=jts,jte
1658 do i=its,ite
1659 lfn_out(i,j) = a1*lfn1(i,j) + a*(lfn_in(i,j) + dt*tend(i,j))
1660 enddo
1661 enddo
1663 ! compute ignition time by interpolation
1664 ! the node was not burning at start but it is burning at end
1665 ! interpolate from the level functions at start and at end
1666 ! lfn_in is the level set function value at time ts
1667 ! lfn_out is the level set function value at time ts+dt
1668 ! 0 is the level set function value at time tign(i,j)
1669 ! thus assuming the level function is approximately linear =>
1670 ! tign(i,j)= ts + ((ts + td) - ts) * lfn_in / (lfn_in - lfn_out)
1671 ! = ts + dt * lfn_in / (lfn_in - lfn_out)
1673 time_now=ts+dt
1674 time_now = time_now + abs(time_now)*epsilon(time_now)*2.
1675 do j=jts,jte
1676 do i=its,ite
1677 ! interpolate the cross-over time
1678 if (.not. lfn_out(i,j)>0 .and. lfn_in(i,j)>0)then
1679 tign(i,j) = ts + dt * lfn_in(i,j) / (lfn_in(i,j) - lfn_out(i,j))
1680 endif
1681 ! set the ignition time outside of burning region
1682 if(lfn_out(i,j)>0.)tign(i,j)=time_now
1683 enddo
1684 enddo
1686 ! check local speed error and stats
1687 ! may not work correctly in parallel
1688 ! init stats
1689 nfirenodes=0
1690 nfireline=0
1691 sum_err=0.
1692 min_err=rmax
1693 max_err=rmin
1694 sum_aerr=0.
1695 min_aerr=rmax
1696 max_aerr=rmin
1697 its1=its+1
1698 jts1=jts+1
1699 ite1=ite-1
1700 jte1=jte-1
1701 ! loop over right inside of the domain
1702 ! cannot use values outside of the domain, would have to reflect and that
1703 ! would change lfn_out
1704 ! cannot use values outside of tile, not synchronized yet
1705 ! so in parallel mode the statistics is not accurate
1706 do j=jts1,jte1
1707 do i=its1,ite1
1708 if(lfn_out(i,j)>0.0)then ! a point out of burning region
1709 if(lfn_out(i+1,j)<=0.or.lfn_out(i,j+1)<=0.or. & ! neighbor in burning region
1710 lfn_out(i-1,j)<=0.or.lfn_out(i,j-1)<=0)then ! point next to fireline
1711 gradx=(lfn_out(i+1,j)-lfn_out(i-1,j))/(2.0*dx) ! central differences
1712 grady=(lfn_out(i,j+1)-lfn_out(i,j-1))/(2.0*dy)
1713 grad2=sqrt(gradx*gradx+grady*grady)
1714 aspeed = (lfn_in(i,j)-lfn_out(i,j))/(dt*max(grad2,rmin))
1715 rr = speed_func(gradx,grady,dx,dy,i,j,fp)
1716 err=aspeed-rr
1717 sum_err=sum_err+err
1718 min_err=min(min_err,err)
1719 max_err=max(max_err,err)
1720 aerr=abs(err)
1721 sum_aerr=sum_aerr+aerr
1722 min_aerr=min(min_aerr,aerr)
1723 max_aerr=max(max_aerr,aerr)
1724 nfireline=nfireline+1
1725 endif
1726 else
1727 nfirenodes=nfirenodes+1
1728 endif
1729 enddo
1730 enddo
1731 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1732 write(msg,'(2(a,i6,f8.4))')'prop_ls: nodes burning',nfirenodes, &
1733 (100.*nfirenodes)/((ite1-its1+1)*(jte1-jts1+1)),'% next to fireline',nfireline
1734 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1735 call message(msg)
1736 if(nfireline>0)then
1737 call print_stat_line('speed error',its1,ite1,jts1,jte1,min_err,max_err,sum_err/nfireline)
1738 call print_stat_line('abs(speed error)',its1,ite1,jts1,jte1,min_aerr,max_aerr,sum_aerr/nfireline)
1739 endif
1741 ! check if the fire did not get to the domain boundary
1742 do k=-1,1,2
1743 ! do kk=1,(boundary_guard*(ide-ids+1))/100 ! in %
1744 do kk=1,boundary_guard ! measured in cells
1745 i=ids+k*kk
1746 if(i.ge.its.and.i.le.ite)then
1747 do j=jts,jte
1748 if(lfn_out(i,j)<=0.)goto 9
1749 enddo
1750 endif
1751 enddo
1752 ! do kk=1,(boundary_guard*(jde-jds+1))/100
1753 do kk=1,boundary_guard ! measured in cells
1754 j=jds+k*kk
1755 if(j.ge.jts.and.j.le.jte)then
1756 do i=its,ite
1757 if(lfn_out(i,j)<=0.)goto 9
1758 enddo
1759 endif
1760 enddo
1761 enddo
1762 goto 10
1763 9 continue
1764 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1765 write(msg,'(a,i2,a,2i8)')'prop_ls: fire',boundary_guard, &
1766 ' cells from domain boundary at node ',i,j
1767 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1768 call message(msg)
1769 call crash('prop_ls: increase the fire region')
1770 10 continue
1772 call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme, &
1773 lfn_out,'prop_ls: lfn out')
1775 end subroutine prop_ls
1777 !
1778 !*****************************
1779 !
1781 subroutine tend_ls( id, &
1782 lims,lime,ljms,ljme, & ! memory dims for lfn
1783 tims,time,tjms,tjme, & ! memory dims for tend
1784 ids,ide,jds,jde, & ! domain - nodes where lfn defined
1785 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1786 ints,inte,jnts,jnte, & ! region - nodes where tend computed
1787 ims,ime,jms,jme, & ! memory dims for ros
1788 its,ite,jts,jte, & ! tile dims - where is ros set
1789 t,dt,dx,dy, & ! scalars in
1790 lfn, & ! arrays in
1791 tbound, & ! scalars out
1792 tend, ros, & ! arrays out
1793 fp &
1794 )
1796 implicit none
1797 ! purpose
1798 ! compute the right hand side of the level set equation
1800 !*** arguments
1801 integer,intent(in)::id,lims,lime,ljms,ljme,tims,time,tjms,tjme
1802 integer,intent(in)::ims,ime,jms,jme,its,ite,jts,jte
1803 integer, intent(in)::ids,ide,jds,jde,ints,inte,jnts,jnte,ips,ipe,jps,jpe
1804 real,intent(in)::t ! time
1805 real,intent(in)::dt,dx,dy ! mesh step
1806 real,dimension(lims:lime,ljms:ljme),intent(inout)::lfn ! level set function
1807 real,intent(out)::tbound ! max allowed time step
1808 real,dimension(tims:time,tjms:tjme),intent(out)::tend ! tendency (rhs of the level set pde)
1809 real,dimension(ims:ime,jms:jme),intent(out)::ros ! rate of spread
1810 type(fire_params),intent(in)::fp
1812 !*** local
1813 real:: te,diffLx,diffLy,diffRx,diffRy, &
1814 diffCx,diffCy,diff2x,diff2y,grad,rr, &
1815 ros_back,ros_wind,ros_slope,advx,advy,scale,nvx,nvy,speed,tanphi
1816 integer::i,j
1817 character(len=128)msg
1819 ! constants
1820 real, parameter:: eps=epsilon(0.0)
1821 !intrinsic epsilon
1822 real, parameter:: zero=0.,one=1.,tol=100*eps, &
1823 safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
1826 ! f90 intrinsic function
1828 intrinsic max,min,sqrt,nint,tiny,huge
1831 #ifdef DEBUG_OUT
1832 real,dimension(tims:time,tjms:tjme)::rra,grada,speeda,tanphia
1833 #endif
1835 !*** executable
1837 ! check array dimensions
1838 call check_mesh_2dim(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme)
1839 call check_mesh_2dim(ints,inte,jnts,jnte,tims,time,tjms,tjme)
1841 call continue_at_boundary(1,1,fire_lfn_ext_up, & !extend by extrapolation but never down
1842 lims,lime,ljms,ljme, & ! memory dims
1843 ids,ide,jds,jde, & ! domain - nodes where lfn defined
1844 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1845 ints,inte,jnts,jnte, & ! tile - nodes update by this thread
1846 lfn) ! array
1848 call print_2d_stats(ints-1,inte+1,jnts,jnte,lims,lime,ljms,ljme, &
1849 lfn,'tend_ls: lfn cont dir x')
1850 call print_2d_stats(ints,inte,jnts-1,jnte+1,lims,lime,ljms,ljme, &
1851 lfn,'tend_ls: lfn cont dir y')
1853 #ifdef DEBUG_OUT
1854 call write_array_m(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme,lfn,'tend_lfn_in',id)
1855 #endif
1857 tbound=0
1858 do j=jnts,jnte
1859 do i=ints,inte
1860 ! one sided differences
1861 diffRx = (lfn(i+1,j)-lfn(i,j))/dx
1862 diffLx = (lfn(i,j)-lfn(i-1,j))/dx
1863 diffRy = (lfn(i,j+1)-lfn(i,j))/dy
1864 diffLy = (lfn(i,j)-lfn(i,j-1))/dy
1865 diffCx = diffLx+diffRx ! TWICE CENTRAL DIFFERENCE
1866 diffCy = diffLy+diffRy
1868 !upwinding - select right or left derivative
1869 select case(fire_upwinding)
1870 case(0) ! none
1871 grad=sqrt(diffCx**2 + diffCy**2)
1872 case(1) ! standard
1873 diff2x=select_upwind(diffLx,diffRx)
1874 diff2y=select_upwind(diffLy,diffRy)
1875 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1876 case(2) ! godunov per osher/fedkiw
1877 diff2x=select_godunov(diffLx,diffRx)
1878 diff2y=select_godunov(diffLy,diffRy)
1879 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1880 case(3) ! eno
1881 diff2x=select_eno(diffLx,diffRx)
1882 diff2y=select_eno(diffLy,diffRy)
1883 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
1884 case(4) ! Sethian - twice stronger pushdown of bumps
1885 grad=sqrt(max(diffLx,0.)**2+min(diffRx,0.)**2 &
1886 + max(diffLy,0.)**2+min(diffRy,0.)**2)
1887 case default
1888 grad=0.
1889 end select
1891 ! normal direction, from central differences
1892 scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps)
1893 nvx=diffCx/scale
1894 nvy=diffCy/scale
1896 ! wind speed in direction of spread
1897 ! speed = vx(i,j)*nvx + vy(i,j)*nvy
1900 ! get rate of spread from wind speed and slope
1902 call fire_ros(ros_back,ros_wind,ros_slope, &
1903 nvx,nvy,i,j,fp)
1905 rr=ros_back + ros_wind + ros_slope
1906 if(fire_grows_only.gt.0)rr=max(rr,0.)
1908 ! set ros for output
1909 if(i.ge.its.and.i.le.ite.and.j.ge.jts.and.j.le.jte)ros(i,j)=rr
1911 if(fire_upwind_split.eq.0)then
1913 ! get rate of spread
1914 te = -rr*grad ! normal term
1916 else
1918 ! normal direction backing rate only
1919 te = - ros_back*grad
1921 ! advection in wind direction
1922 if (abs(speed)> eps) then
1923 advx=fp%vx(i,j)*ros_wind/speed
1924 advy=fp%vy(i,j)*ros_wind/speed
1925 else
1926 advx=0
1927 advy=0
1928 endif
1930 ! tanphi = dzdxf*nvx + dzdyf*nvy
1931 ! advection from slope direction
1932 if(abs(tanphi)>eps) then
1933 advx=advx+fp%dzdxf(i,j)*ros_slope/tanphi
1934 advy=advy+fp%dzdyf(i,j)*ros_slope/tanphi
1935 endif
1937 if(fire_upwind_split.eq.1)then
1939 ! one-sided upwinding
1940 te = te - max(advx,0.)*diffLx - min(advx,0.)*diffRy &
1941 - max(advy,0.)*diffLy - min(advy,0.)*diffRy
1944 elseif(fire_upwind_split.eq.2)then
1946 ! Lax-Friedrichs
1947 call crash('prop_ls: bad fire_upwind_split, Lax-Friedrichs not done yet')
1949 else
1951 call crash('prop_ls: bad fire_upwind_split')
1953 endif
1954 endif
1956 ! cfl condition
1957 if (grad > 0.) then
1958 tbound = max(tbound,rr*(abs(diff2x)/dx+abs(diff2y)/dy)/grad)
1959 endif
1961 ! add numerical viscosity
1962 te=te + fire_viscosity*abs(rr)*((diffRx-diffLx)+(diffRy-diffLy))
1964 tend(i,j)=te
1965 #ifdef DEBUG_OUT
1966 rra(i,j)=rr
1967 grada(i,j)=grad
1968 speeda(i,j)=speed
1969 tanphia(i,j)=tanphi
1970 #endif
1971 !write(msg,*)i,j,grad,dzdx,dzdy
1972 !call message(msg)
1974 !if(abs(te)>1e4)then
1975 ! write(msg,'(a,2i5,e12.3)')'tend_ls: tend out of bounds at ',i,j,te
1976 ! call crash(msg)
1977 !endif
1978 enddo
1979 enddo
1981 #ifdef DEBUG_OUT
1982 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,rra,'rr',id)
1983 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,grada,'grad',id)
1984 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,speeda,'speed',id)
1985 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,tanphia,'tanphi',id)
1986 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,tend,'tend',id)
1987 #endif
1989 call print_2d_stats(ints,inte,jnts,jnte,tims,time,tjms,tjme, &
1990 tend,'tend_ls: tend out')
1992 ! the final CFL bound
1993 tbound = 1/(tbound+tol)
1995 end subroutine tend_ls
1997 !
1998 !**************************
1999 !
2001 real function select_upwind(diffLx,diffRx)
2002 implicit none
2003 real, intent(in):: diffLx, diffRx
2004 real diff2x
2006 ! upwind differences, L or R if bith same sign, otherwise zero
2008 diff2x=0
2009 if (diffLx>0.and.diffRx>0.)diff2x=diffLx
2010 if (diffLx<0.and.diffRx<0.)diff2x=diffRx
2012 select_upwind=diff2x
2013 end function select_upwind
2016 !
2017 !**************************
2018 !
2021 real function select_godunov(diffLx,diffRx)
2022 implicit none
2023 real, intent(in):: diffLx, diffRx
2024 real diff2x,diffCx
2026 ! Godunov scheme: upwind differences, L or R or none
2027 ! always test on > or < never = , much faster because of IEEE
2028 ! central diff >= 0 => take left diff if >0, ortherwise 0
2029 ! central diff <= 0 => take right diff if <0, ortherwise 0
2031 diff2x=0
2032 diffCx=diffRx+diffLx
2033 if (diffLx>0.and..not.diffCx<0)diff2x=diffLx
2034 if (diffRx<0.and. diffCx<0)diff2x=diffRx
2036 select_godunov=diff2x
2037 end function select_godunov
2039 !
2040 !**************************
2041 !
2043 real function select_eno(diffLx,diffRx)
2044 implicit none
2045 real, intent(in):: diffLx, diffRx
2046 real diff2x
2048 ! 1st order ENO scheme
2050 if (.not.diffLx>0 .and. .not.diffRx>0)then
2051 diff2x=diffRx
2052 elseif(.not.diffLx<0 .and. .not.diffRx<0)then
2053 diff2x=diffLx
2054 elseif(.not.diffLx<0 .and. .not.diffRx>0)then
2055 if(.not. abs(diffRx) < abs(diffLx))then
2056 diff2x=diffRx
2057 else
2058 diff2x=diffLx
2059 endif
2060 else
2061 diff2x=0.
2062 endif
2064 select_eno=diff2x
2065 end function select_eno
2067 !
2068 !**************************
2069 !
2071 real function speed_func(diffCx,diffCy,dx,dy,i,j,fp)
2072 !*** purpose
2073 ! the level set method speed function
2074 implicit none
2075 !*** arguments
2076 real, intent(in)::diffCx,diffCy ! x and y coordinates of the direction of propagation
2077 real, intent(in)::dx,dy ! x and y coordinates of the direction of propagation
2078 integer, intent(in)::i,j ! indices of the node to compute the speed at
2079 type(fire_params),intent(in)::fp
2080 !*** local
2081 real::scale,nvx,nvy,r
2082 real::ros_back , ros_wind , ros_slope
2083 real, parameter:: eps=epsilon(0.0)
2084 !*** executable
2085 ! normal direction, from central differences
2086 scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps)
2087 nvx=diffCx/scale
2088 nvy=diffCy/scale
2090 ! get rate of spread from wind speed and slope
2092 call fire_ros(ros_back,ros_wind,ros_slope, &
2093 nvx,nvy,i,j,fp)
2095 r=ros_back + ros_wind + ros_slope
2096 if(fire_grows_only.gt.0)r=max(r,0.)
2097 speed_func=r
2099 end function speed_func
2101 end module module_fr_sfire_core