| blob | 8401a6994c0a39d60bbc56c89ebc0e5f9fe9215c |
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 !$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 ims,ime,jms,jme, &
229 its,ite,jts,jte, &
230 ifs,ife,jfs,jfe, &
231 lfn, tign, fuel_time, time_now, fuel_frac, fire_area)
232 implicit none
234 !*** purpose: determine fraction of fuel remaining
235 !*** NOTE: because variables are cell centered, need halo/sync width 1 before
237 !*** Jan Mandel August 2007 email: jmandel@ucar.edu or Jan.Mandel@gmail.com
239 !*** arguments
241 integer, intent(in) :: its,ite,jts,jte,ims,ime,jms,jme,ifs,ife,jfs,jfe
242 real, intent(in), dimension(ims:ime,jms:jme)::lfn,tign,fuel_time
243 real, intent(in):: time_now
244 real, intent(out), dimension(ifs:ife,jfs:jfe)::fuel_frac
245 real, intent(out), dimension(ims:ime,jms:jme):: fire_area
247 ! ims,ime,jms,jme in memory dimensions
248 ! its,ite,jts,jte in tile dimensions (cells where fuel_frac computed)
249 ! ifs,ife,jfs,jfe in fuel_frac memory dimensions
250 ! lfn in level function, at nodes at midpoints of cells
251 ! tign in ignition time, at nodes at nodes at midpoints of cells
252 ! fuel_time in time constant of fuel, per cell
253 ! time_now in time now
254 ! fuel_frac out fraction of fuel remaining, per cell
255 ! fire_area out fraction of cell area on fire
257 !*** local
259 integer::i,j,ir,jr,icl,jcl,isubcl,jsubcl,i2,j2,ii,jj
260 real::fmax,frat,helpsum1,helpsum2,fuel_left_ff,fire_area_ff,rx,ry,tignf(2,2)
261 real,dimension(3,3)::tff,lff
262 ! help variables instead of arrays fuel_left_f and fire_area_f
263 real::lffij,lffi1j,lffij1,lffi1j1,tifij,tifi1j,tifij1,tifi1j1,tx,ty,txx,tyy
264 ! variables for calculation instead of lff(i,j) and tif(i,j)is lffij,tifij etc..#define IFCELLS (ite-its+1)*fuel_left_irl
265 #define JFCELLS (jte-jts+1)*fuel_left_jrl
266 character(len=128)::msg
267 integer::m,omp_get_thread_num
270 call check_mesh_2dim(its-1,ite+1,jts-1,jte+1,ims,ime,jms,jme)
271 call check_mesh_2dim(its,ite,jts,jte,ifs,ife,jfs,jfe)
273 ! refinement
274 ir=fuel_left_irl
275 jr=fuel_left_jrl
277 if ((ir.ne.2).or.(jr.ne.2)) then
278 call crash('fuel_left: ir.ne.2 or jr.ne.2 ')
279 endif
281 rx=1./ir
282 ry=1./jr
284 ! interpolate level set function to finer grid
285 #ifdef DEBUG_OUT_FUEL_LEFT
286 call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,lff,'lff',0)
287 call write_array_m(1,IFCELLS+1,1,JFCELLS+1,1,IFCELLS+1,1,JFCELLS+1,tif,'tif',0)
288 #endif
290 !
291 ! example for ir=2:
292 !
293 ! | coarse cell |
294 ! its-1 its ite ite+1
295 ! -------X------------|-----.-----X-----.-----|--........----|----------X----------|---------X
296 ! fine node 1 2 3 2*(ite-its+1)
297 ! fine cell 1 2 cell 2*(ite-its+1)
301 ! Loop over cells in Tile
302 ! Changes made by Volodymyr Kondratenko 09/24/2009
303 do icl=its,ite
304 do jcl=jts,jte
305 helpsum1=0
306 helpsum2=0
307 ! Loop over subcells in cell #(icl,jcl)
308 call tign_lfn_interpolation(time_now,icl,jcl,ims,ime,jms,jme, &
309 tign,lfn,tff,lff)
310 do isubcl=1,ir
311 do jsubcl=1,jr
312 !we use 4 points of each subcell in fuel_left_cell_1
313 ! and in fuel_left_cell_2: bottome left are (1,1), (1,2), (2,1), (2,2)
314 if(fuel_left_method.eq.1)then
315 call fuel_left_cell_1( fuel_left_ff, fire_area_ff, &
316 lff(isubcl,jsubcl),lff(isubcl,jsubcl+1),lff(isubcl+1,jsubcl),lff(isubcl+1,jsubcl+1), &
317 tff(isubcl,jsubcl),tff(isubcl,jsubcl+1),tff(isubcl+1,jsubcl),tff(isubcl+1,jsubcl+1), &
318 time_now, fuel_time(icl,jcl))
319 elseif(fuel_left_method.eq.2)then
320 fire_area_ff=0 ! initialize to something until computed in fuel_left_f(i,j)
321 fuel_left_ff=fuel_left_cell_2( &
322 lff(isubcl,jsubcl),lff(isubcl,jsubcl+1),lff(isubcl+1,jsubcl),lff(isubcl+1,jsubcl+1), &
323 tff(isubcl,jsubcl),tff(isubcl,jsubcl+1),tff(isubcl+1,jsubcl),tff(isubcl+1,jsubcl+1), &
324 time_now, fuel_time(icl,jcl))
325 else
326 call crash('fuel_left: unknown fuel_left_method')
327 endif
329 ! consistency check
330 if(fire_area_ff.lt.-1e-6 .or. &
331 (fire_area_ff.eq.0. .and. fuel_left_ff.lt.1.-1e-6))then
332 !$OMP CRITICAL(SFIRE_CORE_CRIT)
333 write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
334 ' of refined mesh fuel burnt',1-fuel_left_ff,' fire area',fire_area_ff
335 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
336 call crash(msg)
337 endif
339 helpsum1=helpsum1+fuel_left_ff
340 helpsum2=helpsum2+fire_area_ff
341 enddo
342 enddo
343 fuel_frac(icl,jcl)=helpsum1
344 fire_area(icl,jcl)=helpsum2
345 enddo
346 enddo
351 #ifdef DEBUG_OUT_FUEL_LEFT
352 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fire_area,'fire_area',0)
353 call write_array_m(its,ite,jts,jte,ims,ime,jms,jme,fuel_frac,'fuel_frac',0)
354 call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fuel_left_f,'fuel_left_f',0)
355 call write_array_m(1,IFCELLS,1,JFCELLS,1,IFCELLS,1,JFCELLS,fire_area_f,'fire_area_f',0)
356 #endif
358 ! finish the averaging
359 do j=jts,jte
360 do i=its,ite
361 fuel_frac(i,j) = fuel_frac(i,j) /(ir*jr) ! multiply by weight for averaging over coarse cell
362 fire_area(i,j) = fire_area(i,j) /(ir*jr) !
363 enddo
364 enddo
366 ! consistency check after sum
367 fmax=0
368 do j=jts,jte
369 do i=its,ite
370 if(fire_area(i,j).eq.0.)then
371 if(fuel_frac(i,j).lt.1.-1e-6)then
372 !$OMP CRITICAL(SFIRE_CORE_CRIT)
373 write(msg,'(a,2i6,2(a,f11.8))')'fuel_left: at node',i,j, &
374 ' fuel burnt',1-fuel_frac(i,j),' but fire area',fire_area(i,j)
375 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
376 call crash(msg)
377 endif
378 else
379 frat=(1-fuel_frac(i,j))/fire_area(i,j)
380 fmax=max(fmax,frat)
381 endif
382 enddo
383 enddo
384 !$OMP CRITICAL(SFIRE_CORE_CRIT)
385 write(msg,'(a,4i6,a,f10.7)')'fuel_left: tile',its,ite,jts,jte,' max fuel burnt/area',fmax
386 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
387 call message(msg)
388 return
391 end subroutine fuel_left
393 !
394 !************************
395 !
398 !
399 !*************************
400 !Subroutine that is calculating tign and lfn of four endpoints of the subcell
401 ! that is located at isubcl,jsubcl of the cell -(icl,jcl)
402 !
403 subroutine tign_lfn_interpolation(time_now,icl,jcl,ims,ime,jms,jme, &
404 tign,lfn,tff,lff)
405 real, intent(in):: time_now ! not ignited nodes will have tign set to >= time_now
406 integer, intent(in) :: icl,jcl
407 !integer, intent(in) :: isubcl,jsubcl
408 integer, intent(in) :: ims,ime,jms,jme ! memory dimensions of all arrays
409 real, intent(in), dimension(ims:ime,jms:jme)::lfn,tign
410 real, intent(out),dimension(3,3)::tff,lff
412 ! (3,1)-------------(3,2)-------------(3,3)
413 ! | | |
414 ! | (2,1) | (2,2) |
415 ! | | |
416 ! | | |
417 ! | | |
418 ! | | |
419 ! (2,1)--------node-(icl,jcl)---------(2,3)-----------(icl,jcl+1)-------------|
420 ! | sub-node (2,2) |
421 ! | | |
422 ! | (1,1) | (1,2) | each fire mesh cell is decomposed in 4
423 ! | | | tff,lff is computed at the nodes of
424 ! | | | the subcells, numbered (1,1)...(3,3)
425 ! (1,1)-------------(1,2)-------------(1,3)--
426 ! | | |
427 ! | (2,1) | (2,2) |
428 ! | | |
429 ! | | |
430 ! | node |
431 ! -------------(icl-1,jcl------------------
432 ! | | | )
433 !
435 !**********************
437 ! Direct calculation tif and lff, avoiding arrays, just for case ir=jr=2
439 lff(1,1)=0.25*(lfn(icl-1,jcl-1)+lfn(icl-1,jcl)+lfn(icl,jcl-1)+lfn(icl,jcl))
440 call tign_four_pnts_interp(tign(icl-1,jcl-1),tign(icl-1,jcl),tign(icl,jcl-1), &
441 tign(icl,jcl),lfn(icl-1,jcl-1),lfn(icl-1,jcl), &
442 lfn(icl,jcl-1),lfn(icl,jcl),lff(1,1),tff(1,1),time_now)
444 lff(1,2)=0.5*(lfn(icl-1,jcl)+lfn(icl,jcl))
445 call tign_line_interp(tign(icl-1,jcl),tign(icl,jcl),lfn(icl-1,jcl),lfn(icl,jcl), &
446 lff(1,2),tff(1,2),time_now)
448 lff(1,3)=0.25*(lfn(icl-1,jcl+1)+lfn(icl-1,jcl)+lfn(icl,jcl+1)+lfn(icl,jcl))
449 call tign_four_pnts_interp(tign(icl-1,jcl),tign(icl-1,jcl+1),tign(icl,jcl), &
450 tign(icl,jcl+1),lfn(icl-1,jcl),lfn(icl-1,jcl+1), &
451 lfn(icl,jcl),lfn(icl,jcl+1),lff(1,3),tff(1,3),time_now)
453 lff(2,1)=0.5*(lfn(icl,jcl-1)+lfn(icl,jcl))
454 call tign_line_interp(tign(icl,jcl-1),tign(icl,jcl),lfn(icl,jcl-1),lfn(icl,jcl), &
455 lff(2,1),tff(2,1),time_now)
457 lff(2,2)=lfn(icl,jcl)
458 tff(2,2)=tign(icl,jcl)
460 lff(2,3)=0.5*(lfn(icl,jcl+1)+lfn(icl,jcl))
461 call tign_line_interp(tign(icl,jcl+1),tign(icl,jcl),lfn(icl,jcl+11),lfn(icl,jcl), &
462 lff(2,3),tff(2,3),time_now)
464 lff(3,1)=0.25*(lfn(icl,jcl-1)+lfn(icl,jcl)+lfn(icl+1,jcl-1)+lfn(icl+1,jcl))
465 call tign_four_pnts_interp(tign(icl,jcl-1),tign(icl,jcl),tign(icl+1,jcl-1), &
466 tign(icl+1,jcl),lfn(icl,jcl-1),lfn(icl,jcl), &
467 lfn(icl+1,jcl-1),lfn(icl+1,jcl),lff(3,1),tff(3,1),time_now)
469 lff(3,2)=0.5*(lfn(icl+1,jcl)+lfn(icl,jcl))
470 call tign_line_interp(tign(icl+1,jcl),tign(icl,jcl),lfn(icl+1,jcl),lfn(icl,jcl), &
471 lff(3,2),tff(3,2),time_now)
473 lff(3,3)=0.25*(lfn(icl,jcl)+lfn(icl,jcl+1)+lfn(icl+1,jcl)+lfn(icl+1,jcl+1))
474 call tign_four_pnts_interp(tign(icl,jcl),tign(icl,jcl+1),tign(icl+1,jcl), &
475 tign(icl+1,jcl+1),lfn(icl,jcl),lfn(icl,jcl+1), &
476 lfn(icl+1,jcl),lfn(icl+1,jcl+1),lff(3,3),tff(3,3),time_now)
481 end subroutine tign_lfn_interpolation
485 subroutine tign_line_interp(tign1,tign2,lfn1,lfn2,lfn_subcl,tign_subcl,time_now)
487 real,intent(in) :: tign1,tign2,lfn1,lfn2,lfn_subcl,time_now
488 real,intent(out) :: tign_subcl
490 real :: a,b,c,err
491 err=0.0001
492 if((lfn1.le.0).AND.(lfn2.le.0)) then
493 tign_subcl=0.5*(tign1+tign2)
494 elseif((lfn1.gt.0).AND.(lfn2.gt.0))then
495 ! if ((tign1.ne.time_now).OR.(tign.ne.time_now))then
496 if ((abs(tign1-time_now).gt.err).or.(abs(tign2-time_now).gt.err)) then
497 call crash('line_interp: when lfn1(2)>0 their tign1(2) should = time_now')
498 else
499 tign_subcl=time_now;
500 endif
501 elseif(lfn_subcl.gt.0) then
502 if ((abs(tign1-time_now).gt.err).OR.(abs(tign2-time_now).gt.err)) then
503 call crash('tign_line_interp one of tign1,2 should be equal time_now')
504 else
505 tign_subcl=time_now;
506 endif
507 else
508 !lfn_subcl<=0;
509 !case when E is on fire
510 ! tign_subcl~=c*lfn_subcl+time_now;
511 if (lfn1.le.0) then
512 a=lfn1;
513 b=tign1-time_now;
514 elseif (lfn2.le.0) then
515 a=lfn2;
516 b=tign2-time_now;
517 else
518 call crash('tign_line_interp: if E is on fire then one of lfn1 or lfn2 &
519 should be negative');
520 endif
521 if (b.gt.0) then
522 call crash('tign_ should be less than time_now');
523 else
524 c=b/a;
525 tign_subcl=c*lfn_subcl+time_now;
526 endif
527 endif
528 end subroutine tign_line_interp
529 !
530 !************************
531 !
533 subroutine tign_four_pnts_interp(tign1,tign2,tign3,tign4,lfn1,lfn2, &
534 lfn3,lfn4,lfn_subcl,tign_subcl,time_now)
536 real,intent(in) :: tign1,tign2,tign3,tign4
537 real,intent(in) :: lfn1,lfn2,lfn3,lfn4,lfn_subcl,time_now
538 real,intent(out) :: tign_subcl
540 real :: a,b,c,err
541 err=0.0001
542 ! to check later why does it appear to be Dummy variable and where
543 !tign_subcl=time_now;
544 !!!!!! will remove later
546 if((lfn1.le.0).AND.(lfn2.le.0).AND.(lfn3.le.0).AND.(lfn4.le.0)) then
547 tign_subcl=0.25*(tign1+tign2+tign3+tign4)
548 elseif((lfn1.gt.0).AND.(lfn2.gt.0).AND.(lfn3.gt.0).AND.(lfn4.gt.0))then
549 !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
550 ! call crash('tign_four_pnts_interp: when lfn1(2,3,4)>0 their tign1(2,3,4) should = time_now')
551 !else
552 tign_subcl=time_now;
553 !endif
554 elseif(lfn_subcl.gt.0) then
555 ! 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
556 ! call crash('tign_line_interp one of tign1(2,3,4) should be equal time_now')
557 ! else
558 tign_subcl=time_now;
559 ! endif
560 else
561 !lfn_subcl<=0;
562 !case when E is on fire
563 ! tign_subcl~=c*lfn_subcl+time_now;
564 a=0;
565 b=0;
566 if (lfn1.le.0) then
567 ! if (tign1.gt.time_now)
568 ! call crash('tign_four_pnts_interp tign1 should be less then time_now');
569 ! else
570 ! Can not assign to a named constant
571 a=a+lfn1*lfn1;
572 b=b+(tign1-time_now)*lfn1;
573 ! endif
574 endif
575 if (lfn2.le.0) then
576 ! if (tign2.gt.time_now)
577 ! call crash('tign_four_pnts_interp tign2 should be less then time_now');
578 ! else
579 ! Can not assign to a named constant
580 a=a+lfn2*lfn2;
581 b=b+(tign2-time_now)*lfn2;
582 ! endif
583 endif
584 if (lfn3.le.0) then
585 ! if (tign3.gt.time_now)
586 ! call crash('tign_four_pnts_interp tign3 should be less then time_now');
587 ! else
588 ! Can not assign to a named constant
589 a=a+lfn3*lfn3;
590 b=b+(tign3-time_now)*lfn3;
591 ! endif
592 endif
593 if (lfn4.le.0) then
594 ! if (tign4.gt.time_now)
595 ! call crash('tign_four_pnts_interp tign4 should be less then time_now');
596 ! else
597 ! Can not assign to a named constant
598 a=a+lfn4*lfn4;
599 b=b+(tign4-time_now)*lfn4;
600 ! endif
601 endif
602 if ((abs(a).lt.err).or.(b.lt.0)) then
603 call crash('tign_four_pnts_interp: if E is on fire then one of cells &
604 should be on fire or tign_ should be less than time_now')
605 else
606 c=b/a;
607 tign_subcl=c*lfn_subcl+time_now;
608 endif
609 endif
611 end subroutine tign_four_pnts_interp
614 !
615 !************************
616 !
619 subroutine fuel_left_cell_1( fuel_frac_left, fire_frac_area, &
620 lfn00,lfn01,lfn10,lfn11, &
621 tign00,tign01,tign10,tign11,&
622 time_now, fuel_time_cell)
623 !*** purpose: compute the fuel fraction left in one cell
624 implicit none
625 !*** arguments
626 real, intent(out):: fuel_frac_left, fire_frac_area !
627 real, intent(in)::lfn00,lfn01,lfn10,lfn11 ! level set function at 4 corners of the cell
628 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the 4 corners of the cell
629 real, intent(in)::time_now ! the time now
630 real, intent(in)::fuel_time_cell ! time to burns off to 1/e
631 !*** Description
632 ! The area burning is given by the condition L <= 0, where the function P is
633 ! interpolated from lfn(i,j)
634 !
635 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
636 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
637 ! when lfn(i,j)=0.
638 !
639 ! The function computes an approxmation of the integral
640 !
641 !
642 ! /\
643 ! |
644 ! fuel_frac_left = 1 - | 1 - exp(-T(x,y)/fuel_time_cell)) dxdy
645 ! |
646 ! \/
647 ! 0<x<1
648 ! 0<y<1
649 ! L(x,y)<=0
650 !
651 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
652 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
653 ! so dx=1 and dy=1 assumed.
654 !
655 ! Example:
656 !
657 ! lfn<0 lfn>0
658 ! (0,1) -----O--(1,1) O = points on the fireline, T=tnow
659 ! | \ | A = the burning area for computing
660 ! | \| fuel_frac(i,j)
661 ! | A O
662 ! | |
663 ! | |
664 ! (0,0)---------(1,0)
665 ! lfn<0 lfn<0
666 !
667 ! Approximations allowed:
668 ! The fireline can be approximated by straight line(s).
669 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
670 !
671 ! Requirements:
672 ! 1. The output should be a continuous function of the arrays lfn and
673 ! tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
674 ! 2. The output should be invariant to the symmetries of the input in each cell.
675 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
676 ! 4. The result should be at least 1st order accurate in the sense that it is
677 ! exact if the time from ignition is a linear function.
678 !
679 ! If time from ignition is approximated by polynomial in the burnt
680 ! region of the cell, this is integral of polynomial times exponential
681 ! over a polygon, which can be computed exactly.
682 !
683 ! Requirement 4 is particularly important when there is a significant decrease
684 ! of the fuel fraction behind the fireline on the mesh scale, because the
685 ! rate of fuel decrease right behind the fireline is much larger
686 ! (exponential...). This will happen when
687 !
688 ! change of time from ignition within one mesh cell / fuel_time_cell is not << 1
689 !
690 ! This is the same as
691 !
692 ! mesh cell size
693 ! X = ------------------------- is not << 1
694 ! fireline speed * fuel_time_cell
695 !
696 !
697 ! When X is large then the fuel burnt in one timestep in the cell is
698 ! approximately proportional to length of fireline in that cell.
699 !
700 ! When X is small then the fuel burnt in one timestep in the cell is
701 ! approximately proportional to the area of the burning region.
702 !
704 !*** calls
705 intrinsic tiny
707 !*** local
708 real::ps,aps,area,ta,out
709 real::t00,t01,t10,t11
710 real,parameter::safe=tiny(aps)
711 character(len=128)::msg
713 ! the following algorithm is a very crude approximation
715 ! minus time since ignition, 0 if no ignition yet
716 ! it is possible to have 0 in fire region when ignitin time falls in
717 ! inside the time step because lfn is updated at the beginning of the time step
719 t00=tign00-time_now
720 if(lfn00>0. .or. t00>0.)t00=0.
721 t01=tign01-time_now
722 if(lfn01>0. .or. t01>0.)t01=0.
723 t10=tign10-time_now
724 if(lfn10>0. .or. t10>0.)t10=0.
725 t11=tign11-time_now
726 if(lfn11>0. .or. t11>0.)t11=0.
728 ! approximate burning area, between 0 and 1
729 ps = lfn00+lfn01+lfn10+lfn11
730 aps = abs(lfn00)+abs(lfn01)+abs(lfn10)+abs(lfn11)
731 aps=max(aps,safe)
732 area =(-ps/aps+1.)/2.
733 area = max(area,0.) ! make sure area is between 0 and 1
734 area = min(area,1.)
736 ! average negative time since ignition
737 ta=0.25*(t00+t01+t10+t11)
739 ! exp decay in the burning area
740 out=1.
741 !if(area>0.)out=1. - area*(1. - exp(ta/fuel_time_cell))
742 if(area>0)out=area*exp(ta/fuel_time_cell) + (1. - area)
744 if(out>1.)then
745 !$OMP CRITICAL(SFIRE_CORE_CRIT)
746 write(msg,*)'out=',out,'>1 area=',area,' ta=',ta
747 call message(msg)
748 write(msg,*)'tign=',tign00,tign01,tign10,tign11,' time_now=',time_now
749 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
750 call message(msg)
751 !call message('WARNING: fuel_left_cell_1: fuel fraction > 1')
752 call crash('fuel_left_cell_1: fuel fraction > 1')
753 endif
755 !out = max(out,0.) ! make sure out is between 0 and 1
756 !out = min(out,1.)
758 fuel_frac_left = out
759 fire_frac_area = area
761 end subroutine fuel_left_cell_1
763 !
764 !****************************************
765 !
767 real function fuel_left_cell_2( &
768 lfn00,lfn01,lfn10,lfn11, &
769 tign00,tign01,tign10,tign11,&
770 time_now, fuel_time_cell)
771 !*** purpose: compute the fuel fraction left in one cell
772 implicit none
773 !*** arguments
774 real, intent(in)::lfn00,lfn01,lfn10,lfn11 ! level set function at 4 corners of the cell
775 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the 4 corners of the cell
776 real, intent(in)::time_now ! the time now
777 real, intent(in)::fuel_time_cell ! time to burns off to 1/e
778 !*** Description
779 ! The area burning is given by the condition L <= 0, where the function P is
780 ! interpolated from lfn(i,j)
781 !
782 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
783 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
784 ! when lfn(i,j)=0.
785 !
786 ! The function computes an approxmation of the integral
787 !
788 !
789 ! /\
790 ! |
791 ! fuel_frac_left = 1 - | 1 - exp(-T(x,y)*fuel_speed)) dxdy
792 ! |
793 ! \/
794 ! 0<x<1
795 ! 0<y<1
796 ! L(x,y)<=0
797 !
798 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
799 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
800 ! so dx=1 and dy=1 assumed.
801 !
802 ! Example:
803 !
804 ! lfn<0 lfn>0
805 ! (0,1) -----O--(1,1) O = points on the fireline, T=tnow
806 ! | \ | A = the burning area for computing
807 ! | \| fuel_frac(i,j)
808 ! | A O
809 ! | |
810 ! | |
811 ! (0,0)---------(1,0)
812 ! lfn<0 lfn<0
813 !
814 ! Approximations allowed:
815 ! The fireline can be approximated by straight line(s).
816 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
817 !
818 ! Requirements:
819 ! 1. The output should be a continuous function of the arrays lfn and
820 ! tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
821 ! 2. The output should be invariant to the symmetries of the input in each cell.
822 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
823 ! 4. The result should be at least 1st order accurate in the sense that it is
824 ! exact if the time from ignition is a linear function.
825 !
826 ! If time from ignition is approximated by polynomial in the burnt
827 ! region of the cell, this is integral of polynomial times exponential
828 ! over a polygon, which can be computed exactly.
829 !
830 ! Requirement 4 is particularly important when there is a significant decrease
831 ! of the fuel fraction behind the fireline on the mesh scale, because the
832 ! rate of fuel decrease right behind the fireline is much larger
833 ! (exponential...). This will happen when
834 !
835 ! change of time from ignition within one mesh cell * fuel speed is not << 1
836 !
837 ! This is the same as
838 !
839 ! mesh cell size*fuel_speed
840 ! ------------------------- is not << 1
841 ! fireline speed
842 !
843 !
844 ! When X is large then the fuel burnt in one timestep in the cell is
845 ! approximately proportional to length of fireline in that cell.
846 !
847 ! When X is small then the fuel burnt in one timestep in the cell is
848 ! approximately proportional to the area of the burning region.
849 !
851 #ifndef FUEL_LEFT
852 call crash('fuel_left_cell_2: not implemented, please use fire_fuel_left_method=1')
853 fuel_left_cell_2=0. ! to avoid compiler warning about value not set
854 end function fuel_left_cell_2
855 #else
857 !*** calls
858 intrinsic tiny
860 !*** local
861 real::ps,aps,area,ta,out
862 real::t00,t01,t10,t11
863 real,parameter::safe=tiny(aps)
864 character(len=128)::msg
866 !*** local
867 integer::i,j,k
869 ! least squares
870 integer::mmax,nb,nmax,pmax,nin,nout
871 parameter(mmax=3,nb=64,nmax=8,pmax=8)
872 integer lda, ldb, lwork, info
873 parameter (lda=nmax, ldb=nmax, lwork=mmax+nmax+nb*(nmax+pmax))
874 integer n,m,p
875 real,dimension(lda,mmax):: mA
876 real,dimension(nmax):: vecD
877 real,dimension(lwork):: WORK
878 real,dimension(pmax):: vecY
879 real,dimension(mmax):: vecX
880 real,dimension(ldb,pmax)::mB
881 real,dimension(mmax)::u
883 real::tweight,tdist
884 integer::kk,ll,ss
885 real::rnorm
886 real,dimension(8,2)::xylist,xytlist
887 real,dimension(8)::tlist,llist,xt
888 real,dimension(5)::xx,yy
889 real,dimension(5)::lfn,tign
891 integer:: npoint
892 real::tt,x0,y0,xts,xte,yts,yte,xt1,xt2
893 real::lfn0,lfn1,dist,nr,c,s,errQ,ae,ce,ceae,a0,a1,a2,d,cet
894 real::s1,s2,s3
895 real::upper,lower,ah,ch,aa,cc,aupp,cupp,alow,clow
896 real,dimension(2,2)::mQ
897 real,dimension(2)::ut
899 !calls
900 intrinsic epsilon
902 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
905 ! external functions
906 real::snrm2
907 double precision :: dnrm2
908 external dnrm2
909 external snrm2
910 ! external subroutines
911 external sggglm
912 external dggglm
913 !executable statements
915 ! a very crude approximation - replace by a better code
916 ! call check_mesh_2dim(ids,ide+1,jds,jde+1,ims,ime,jms,jme)
917 t00=tign00-time_now
918 if(lfn00>=0. .or. t00>0.)t00=0.
919 t01=tign01-time_now
920 if(lfn01>=0. .or. t01>0.)t01=0.
921 t10=tign10-time_now
922 if(lfn10>=0. .or. t10>0.)t10=0.
923 t11=tign11-time_now
924 if(lfn11>=0. .or. t11>0.)t11=0.
926 !if (t00<0 .or. t01 <0 .or. t10<0 .or. t11<0) then
927 ! print *,'tign00=',tign00,'tign10=',tign10,&
928 ! 'tign01=',tign01,'tign11=',tign11
929 !end if
933 !*** case0 Do nothing
934 if ( lfn00>=0 .and. lfn10>=0 .and. lfn01>=0 .and. lfn11>=0 ) then
935 out = 1.0 ! fuel_left, no burning
936 !*** case4 all four coners are burning
937 else if (lfn00<=0 .and. lfn10<=0 .and. lfn01<=0 .and. lfn11<=0) then
938 ! least squares, A matrix for points
939 mA(1,1)=0.0
940 mA(2,1)=1.0
941 mA(3,1)=0.0
942 mA(4,1)=1.0
943 mA(1,2)=0.0
944 mA(2,2)=0.0
945 mA(3,2)=1.0
946 mA(4,2)=1.0
947 mA(1,3)=1.0
948 mA(2,3)=1.0
949 mA(3,3)=1.0
950 mA(4,3)=1.0
951 ! D vector, time from ignition
952 vecD(1)=time_now-tign00
953 vecD(2)=time_now-tign10
954 vecD(3)=time_now-tign01
955 vecD(4)=time_now-tign11
956 ! B matrix, weights
957 do kk=1,4
958 do ll=1,4
959 mB(kk,ll)=0.0
960 end do
961 mB(kk,kk)=2.0
962 end do
963 ! set the m,n,p
964 n=4 ! rows of matrix A and B
965 m=3 ! columns of matrix A
966 p=4 ! columns of matrix B
967 ! call least squqres in LAPACK
968 call SGGGLM(N,M,P,mA,LDA,mB,LDB,vecD,vecX,vecY, &
969 WORK,LWORK,INFO)
970 rnorm=snrm2(p,vecY,1)
971 ! integrate
972 u(1)=-vecX(1)/fuel_time_cell
973 u(2)=-vecX(2)/fuel_time_cell
974 u(3)=-vecX(3)/fuel_time_cell
975 !fuel_burn(i,j)=1-exp(u(3))*intexp(u(1)*dx)*intexp(u(2)*dy)
976 s1=u(1)
977 s2=u(2)
978 out=1-exp(u(3))*intexp(s1)*intexp(s2)
979 !print *,'intexp
980 if ( out<0 .or. out>1.0 ) then
981 print *,'case4, out should be between 0 and 1'
982 end if
983 !*** case 1,2,3
984 else
985 ! set xx, yy for the coner points
986 ! move these values out of i and j loop to speed up
987 xx(1) = -0.5
988 xx(2) = 0.5
989 xx(3) = 0.5
990 xx(4) = -0.5
991 xx(5) = -0.5
992 yy(1) = -0.5
993 yy(2) = -0.5
994 yy(3) = 0.5
995 yy(4) = 0.5
996 yy(5) = -0.5
997 lfn(1)=lfn00
998 lfn(2)=lfn10
999 lfn(3)=lfn11
1000 lfn(4)=lfn01
1001 lfn(5)=lfn00
1002 tign(1)=t00
1003 tign(2)=t10
1004 tign(3)=t11
1005 tign(4)=t01
1006 tign(5)=t00
1007 npoint = 0 ! number of points in polygon
1008 !print *,'time_now=',time_now
1009 !print *,'lfn00=',lfn00,'lfn10=',lfn10,&
1010 ! 'lfn01=',lfn01,'lfn11=',lfn11
1011 !print *,'t00=',t00,'t10=',t10,&
1012 ! 't01=',t01,'t11=',t11
1014 do k=1,4
1015 lfn0=lfn(k )
1016 lfn1=lfn(k+1)
1017 if ( lfn0 <= 0.0 ) then
1018 npoint = npoint + 1
1019 xylist(npoint,1)=xx(k)
1020 xylist(npoint,2)=yy(k)
1021 tlist(npoint)=-tign(k)
1022 llist(npoint)=lfn0
1023 end if
1024 if ( lfn0*lfn1 < 0 ) then
1025 npoint = npoint + 1
1026 tt=lfn0/(lfn0-lfn1)
1027 x0=xx(k)+( xx(k+1)-xx(k) )*tt
1028 y0=yy(k)+( yy(k+1)-yy(k) )*tt
1029 xylist(npoint,1)=x0
1030 xylist(npoint,2)=y0
1031 tlist(npoint)=0 ! on fireline
1032 llist(npoint)=0
1033 end if
1034 end do
1036 ! make the list circular
1037 tlist(npoint+1)=tlist(1)
1038 llist(npoint+1)=llist(1)
1039 xylist(npoint+1,1)=xylist(1,1)
1040 xylist(npoint+1,2)=xylist(1,2)
1041 !* least squares, A matrix for points
1042 do kk=1,npoint
1043 mA(kk,1)=xylist(kk,1)
1044 mA(kk,2)=xylist(kk,2)
1045 mA(kk,3)=1.0
1046 vecD(kk)=tlist(kk) ! D vector,time from ignition
1047 end do
1048 ! B matrix, weights
1049 do kk=1,ldb
1050 do ll=1,pmax
1051 mB(kk,ll)=0.0 ! clear
1052 end do
1053 end do
1055 do kk=1,npoint
1056 mb(kk,kk) = 10 ! large enough
1057 do ll=1,npoint
1058 if ( kk .ne. ll ) then
1059 dist = sqrt( (xylist(kk,1)-xylist(ll,1))**2+ &
1060 (xylist(kk,2)-xylist(ll,2))**2 )
1061 mB(kk,kk)=min( mB(kk,kk) , dist )
1062 end if
1063 end do !ll
1064 mB(kk,kk)=mB(kk,kk)+1.
1065 end do ! kk
1066 ! set the m,n,p
1067 n=npoint ! rows of matrix A and B
1068 m=3 ! columns of matrix A
1069 p=npoint ! columns of matrix B
1070 !* call least squqres in LAPACK
1071 call SGGGLM(N,M,P,mA,LDA,mB,LDB,vecD,vecX,vecY, &
1072 WORK,LWORK,INFO)
1073 !print *,'after LS in case3'
1074 !print *,'vecX from LS',vecX
1075 !print *,'tign inputed',tign00,tign10,tign11,tign01
1076 rnorm=snrm2(p,vecY,1)
1077 u(1)=vecX(1)
1078 u(2)=vecX(2)
1079 u(3)=vecX(3)
1080 ! rotate to gradient on x only
1081 nr = sqrt(u(1)**2+u(2)**2)
1082 if(.not.nr.gt.eps)then
1083 out=1.
1084 goto 900
1085 endif
1086 c = u(1)/nr
1087 s = u(2)/nr
1088 mQ(1,1)=c
1089 mQ(1,2)=s
1090 mQ(2,1)=-s
1091 mQ(2,2)=c
1092 ! mat vec multiplication
1093 call matvec(mQ,2,2,u,3,ut,2,2,2)
1094 errQ = ut(2) ! should be zero
1095 ae = -ut(1)/fuel_time_cell
1096 ce = -u(3)/fuel_time_cell
1097 cet=ce!keep ce
1098 call matmatp(xylist,8,2,mQ,2,2,xytlist,8,2,npoint+1,2,2)
1099 call sortxt( xytlist, 8,2, xt,8,npoint )
1100 out=0.0
1101 aupp=0.0
1102 cupp=0.0
1103 alow=0.0
1104 clow=0.0
1105 do k=1,npoint-1
1106 xt1=xt(k)
1107 xt2=xt(k+1)
1108 upper=0
1109 lower=0
1110 ah=0
1111 ch=0
1112 if ( xt2-xt1 > eps*100 ) then
1114 do ss=1,npoint
1115 xts=xytlist(ss,1)
1116 yts=xytlist(ss,2)
1117 xte=xytlist(ss+1,1)
1118 yte=xytlist(ss+1,2)
1120 if ( (xts>xt1 .and. xte>xt1) .or. &
1121 (xts<xt2 .and. xte<xt2) ) then
1122 aa = 0 ! do nothing
1123 cc = 0
1124 else
1125 aa = (yts-yte)/(xts-xte)
1126 cc = (xts*yte-xte*yts)/(xts-xte)
1127 if (xte<xts) then
1128 aupp = aa
1129 cupp = cc
1130 ah=ah+aa
1131 ch=ch+cc
1132 upper=upper+1
1133 else
1134 alow = aa
1135 clow = cc
1136 lower=lower+1
1137 end if
1138 end if!(xts>xt1 .and. xte>xt1)
1139 end do ! ss
1140 ce=cet !use stored ce
1141 if (ae*xt1+ce > 0 ) then
1142 ce=ce-(ae*xt1+ce)!shift small amounts exp(-**)
1143 end if
1144 if (ae*xt2+ce > 0) then
1145 ce=ce-(ae*xt2+ce)
1146 end if
1148 ah = aupp-alow
1149 ch = cupp-clow
1150 ! integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2
1151 ! numerically sound for ae->0, ae -> infty
1152 ! this can be important for different model scales
1153 ! esp. if someone runs the model in single precision!!
1154 ! s1=int((ah*x+ch),x,xt1,xt2)
1155 s1 = (xt2-xt1)*((1./2.)*ah*(xt2+xt1)+ch)
1156 ! s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)
1157 ceae=ce/ae;
1158 s2 = -ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1))
1159 ! s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)
1160 ! s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)
1161 ! expand in Taylor series around ae=0
1162 ! collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)
1163 ! =(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
1164 ! + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae
1165 ! + 1/2*xt1^2-1/2*xt2^2
1166 !
1167 ! coefficient at ae^2 in the expansion, after some algebra
1168 a2=(xt1-xt2)*((1./4.)*(xt1+xt2)*ceae**2+(1./3.)* &
1169 (xt1**2+xt1*xt2+xt2**2)*ceae+(1./8.)* &
1170 (xt1**3+xt1*(xt2**2)+xt1**2*xt2+xt2**3))
1171 d=(ae**4)*a2
1173 if (abs(d)>eps) then
1174 ! since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae
1175 ! for ae, ce -> 0 rounding error approx eps/ae^2
1176 s3=( exp(ae*(xt1+ceae))*(ae*xt1-1)-&
1177 exp(ae*(xt2+ceae))*(ae*xt2-1) )/(ae**2)
1179 !we do not worry about rounding as xt1 -> xt2, then s3 -> 0
1180 else
1181 ! coefficient at ae^1 in the expansion
1182 a1=(xt1-xt2)*((1./2.)*ceae*(xt1+xt2)+(1./3.)*&
1183 (xt1**2+xt1*xt2+xt2**2))
1184 ! coefficient at ae^0 in the expansion for ae->0
1185 a0=(1./2.)*(xt1-xt2)*(xt1+xt2)
1186 s3=a0+a1*ae+a2*ae**2; ! approximate the integral
1187 end if
1189 s3=ah*s3
1190 out=out+s1+s2+s3
1191 out=1-out !fuel_left
1192 if(out<0 .or. out>1) then
1193 print *,':fuel_fraction should be between 0 and 1'
1194 !print *, 'eps= ', eps
1195 !print *, 'xt1= ', xt1, 'xt2= ', xt2
1196 !print *,'ae= ',ae,'ce= ',ce,'ah= ',ah,'ch= ',ch
1197 !print *,'a0= ', a0,'a1= ', a1,'a2= ', a2
1198 print *,'s1= ', s1,'s2= ', s2,'s3= ', s3
1199 print *,':fuel_fraction =',out
1200 end if!print
1202 end if
1203 end do ! k
1205 end if ! if case0, elseif case4 ,else case123
1207 900 continue
1208 if(out>1. .or. out<0.)call crash('fuel_left_cell_2: fuel fraction out of bounds [0,1]')
1209 fuel_left_cell_2 = out
1210 end function fuel_left_cell_2
1211 #endif
1212 !
1213 !****************************************
1214 !
1216 real function fuel_left_cell_3( &
1217 lfn00,lfn01,lfn10,lfn11, &
1218 tign00,tign01,tign10,tign11,&
1219 time_now, fuel_time_cell)
1220 !*** purpose: compute the fuel fraction left in one cell
1221 implicit none
1222 !*** arguments
1223 real, intent(in)::lfn00,lfn01,lfn10,lfn11 ! level set function at 4 corners of the cell
1224 real, intent(in)::tign00,tign01,tign10,tign11! ignition time at the 4 corners of the cell
1225 real, intent(in)::time_now ! the time now
1226 real, intent(in)::fuel_time_cell ! time to burns off to 1/e
1227 !*** Description
1228 ! The area burning is given by the condition L <= 0, where the function P is
1229 ! interpolated from lfn(i,j)
1230 !
1231 ! The time since ignition is the function T, interpolated in from tign(i,j)-time_now.
1232 ! The values of tign(i,j) where lfn(i,j)>=0 are ignored, tign(i,j)=0 is taken
1233 ! when lfn(i,j)=0.
1234 !
1235 ! The function computes an approxmation of the integral
1236 !
1237 !
1238 ! /\
1239 ! |
1240 ! fuel_frac_left = 1 - | 1 - exp(-T(x,y)*fuel_speed)) dxdy
1241 ! |
1242 ! \/
1243 ! 0<x<1
1244 ! 0<y<1
1245 ! L(x,y)<=0
1246 !
1247 ! When the cell is not burning at all (all lfn>=0), then fuel_frac(i,j)=1.
1248 ! Because of symmetries, the result should not depend on the mesh spacing dx dy
1249 ! so dx=1 and dy=1 assumed.
1250 !
1251 ! Example:
1252 !
1253 ! lfn<0 lfn>0
1254 ! (0,1) -----O--(1,1) O = points on the fireline, T=tnow
1255 ! | \ | A = the burning area for computing
1256 ! | \| fuel_frac(i,j)
1257 ! | A O
1258 ! | |
1259 ! | |
1260 ! (0,0)---------(1,0)
1261 ! lfn<0 lfn<0
1262 !
1263 ! Approximations allowed:
1264 ! The fireline can be approximated by straight line(s).
1265 ! When all cell is burning, approximation by 1 point Gaussian quadrature is OK.
1266 !
1267 ! Requirements:
1268 ! 1. The output should be a continuous function of the arrays lfn and
1269 ! tign whenever lfn(i,j)=0 implies tign(i,j)=tnow.
1270 ! 2. The output should be invariant to the symmetries of the input in each cell.
1271 ! 3. Arbitrary combinations of the signs of lfn(i,j) should work.
1272 ! 4. The result should be at least 1st order accurate in the sense that it is
1273 ! exact if the time from ignition is a linear function.
1274 !
1275 ! If time from ignition is approximated by polynomial in the burnt
1276 ! region of the cell, this is integral of polynomial times exponential
1277 ! over a polygon, which can be computed exactly.
1278 !
1279 ! Requirement 4 is particularly important when there is a significant decrease
1280 ! of the fuel fraction behind the fireline on the mesh scale, because the
1281 ! rate of fuel decrease right behind the fireline is much larger
1282 ! (exponential...). This will happen when
1283 !
1284 ! change of time from ignition within one mesh cell * fuel speed is not << 1
1285 !
1286 ! This is the same as
1287 !
1288 ! mesh cell size*fuel_speed
1289 ! ------------------------- is not << 1
1290 ! fireline speed
1291 !
1292 !
1293 ! When X is large then the fuel burnt in one timestep in the cell is
1294 ! approximately proportional to length of fireline in that cell.
1295 !
1296 ! When X is small then the fuel burnt in one timestep in the cell is
1297 ! approximately proportional to the area of the burning region.
1299 #ifndef FUEL_LEFT
1300 call crash('fuel_left_cell_3: not implemented, please use fire_fuel_left_method=1')
1301 fuel_left_cell_3=0. ! to avoid compiler warning about value not set
1302 end function fuel_left_cell_3
1303 #else
1305 !*** calls
1306 intrinsic tiny
1308 !*** local
1309 real::ps,aps,area,ta,out
1310 real::t00,t01,t10,t11
1311 real,parameter::safe=tiny(aps)
1312 character(len=128)::msg
1314 !*** local
1315 integer::i,j,k
1317 ! least squares
1318 integer::mmax,nb,nmax,pmax,nin,nout
1319 parameter(mmax=3,nb=64,nmax=8,pmax=8)
1320 integer lda, ldb, lwork, info
1321 parameter (lda=nmax, ldb=nmax, lwork=mmax+nmax+nb*(nmax+pmax))
1322 integer n,m,p
1323 real,dimension(lda,mmax):: mA
1324 real,dimension(nmax):: vecD
1325 real,dimension(lwork):: WORK
1326 real,dimension(pmax):: vecY
1327 real,dimension(mmax):: vecX
1328 real,dimension(ldb,pmax)::mB
1329 real,dimension(mmax)::u
1331 real::tweight,tdist
1332 integer::kk,ll,ss
1333 real::rnorm
1334 real,dimension(8,2)::xylist,xytlist
1335 real,dimension(8)::tlist,llist,xt
1336 real,dimension(5)::xx,yy
1337 real,dimension(5)::lfn,tign
1339 integer:: npoint
1340 real::tt,x0,y0,xts,xte,yts,yte,xt1,xt2
1341 real::lfn0,lfn1,dist,nr,c,s,errQ,ae,ce,ceae,a0,a1,a2,d,cet
1342 real::s1,s2,s3
1343 real::upper,lower,ah,ch,aa,cc,aupp,cupp,alow,clow
1344 real,dimension(2,2)::mQ
1345 real,dimension(2)::ut
1347 !calls
1348 intrinsic epsilon
1350 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
1352 !!!! For finite differences by VK
1353 real::tign_middle,dt_dx,dt_dy,lfn_middle,a,b,c
1355 ! external functions
1356 real::snrm2
1357 double precision :: dnrm2
1358 external dnrm2
1359 external snrm2
1360 ! external subroutines
1361 external sggglm
1362 external dggglm
1363 !executable statements
1365 ! a very crude approximation - replace by a better code
1366 ! call check_mesh_2dim(ids,ide+1,jds,jde+1,ims,ime,jms,jme)
1367 t00=tign00-time_now
1368 if(lfn00>=0. .or. t00>0.)t00=0.
1369 t01=tign01-time_now
1370 if(lfn01>=0. .or. t01>0.)t01=0.
1371 t10=tign10-time_now
1372 if(lfn10>=0. .or. t10>0.)t10=0.
1373 t11=tign11-time_now
1374 if(lfn11>=0. .or. t11>0.)t11=0.
1376 !if (t00<0 .or. t01 <0 .or. t10<0 .or. t11<0) then
1377 ! print *,'tign00=',tign00,'tign10=',tign10,&
1378 ! 'tign01=',tign01,'tign11=',tign11
1379 !end if
1383 !*** case0 Do nothing
1384 if ( lfn00>=0 .and. lfn10>=0 .and. lfn01>=0 .and. lfn11>=0 ) then
1385 out = 1.0 ! fuel_left, no burning
1386 !*** case4 all four coners are burning
1387 else if (lfn00<=0 .and. lfn10<=0 .and. lfn01<=0 .and. lfn11<=0) then
1388 !!!!!!!!!!!
1390 ! All burning
1391 ! T=u(1)*x+u(2)*y+u(3)
1392 ! t(0,0)=tign(1,1)
1393 ! t(0,fd(2))=tign(1,2)
1394 ! t(fd(1),0)=tign(2,1)
1395 ! t(fd(1),fd(2))=tign(2,2)
1396 ! t(g/2,h/2)=sum(tign(i,i))/4
1397 ! dt/dx=(1/2h)*(t10-t00+t11-t01)
1398 ! dt/dy=(1/2h)*(t01-t00+t11-t10)
1399 ! approximate T(x,y)=u(1)*x+u(2)*y+u(3) Using finite differences
1400 ! t(x,y)=t(h/2,h/2)+(x-h/2)*dt/dx+(y-h/2)*dt/dy
1401 ! u(1)=dt/dx
1402 ! u(2)=dt/dy
1403 ! u(3)=t(h/2,h/2)-h/2(dt/dx+dt/dy)
1405 tign_middle=(t00+t01+t10+t11)/4
1406 ! since mesh_size is 1 we replace fd(1) and fd(2) by 1
1407 dt_dx=(t00-t10+t01-t11)/2
1408 dt_dy=(t00-t01+t10-t11)/2
1409 u(1)=dt_dx
1410 u(2)=dt_dy
1411 u(3)=tign_middle-(dt_dx+dt_dy)/2
1413 ! integrate
1414 u(1)=-u(1)/fuel_time_cell
1415 u(2)=-u(2)/fuel_time_cell
1416 u(3)=-u(3)/fuel_time_cell
1418 !fuel_burn(i,j)=1-exp(u(3))*intexp(u(1)*dx)*intexp(u(2)*dy)
1419 s1=u(1)
1420 s2=u(2)
1421 out=1-exp(u(3))*intexp(s1)*intexp(s2)
1422 !print *,'intexp
1423 if ( out<0 .or. out>1.0 ) then
1424 print *,'case4, out should be between 0 and 1'
1425 end if
1426 !*** case 1,2,3
1427 else
1428 ! set xx, yy for the coner points
1429 ! move these values out of i and j loop to speed up
1430 xx(1) = -0.5
1431 xx(2) = 0.5
1432 xx(3) = 0.5
1433 xx(4) = -0.5
1434 xx(5) = -0.5
1435 yy(1) = -0.5
1436 yy(2) = -0.5
1437 yy(3) = 0.5
1438 yy(4) = 0.5
1439 yy(5) = -0.5
1440 lfn(1)=lfn00
1441 lfn(2)=lfn10
1442 lfn(3)=lfn11
1443 lfn(4)=lfn01
1444 lfn(5)=lfn00
1445 tign(1)=t00
1446 tign(2)=t10
1447 tign(3)=t11
1448 tign(4)=t01
1449 tign(5)=t00
1450 npoint = 0 ! number of points in polygon
1451 !print *,'time_now=',time_now
1452 !print *,'lfn00=',lfn00,'lfn10=',lfn10,&
1453 ! 'lfn01=',lfn01,'lfn11=',lfn11
1454 !print *,'t00=',t00,'t10=',t10,&
1455 ! 't01=',t01,'t11=',t11
1457 do k=1,4
1458 lfn0=lfn(k )
1459 lfn1=lfn(k+1)
1460 if ( lfn0 <= 0.0 ) then
1461 npoint = npoint + 1
1462 xylist(npoint,1)=xx(k)
1463 xylist(npoint,2)=yy(k)
1464 tlist(npoint)=-tign(k)
1465 llist(npoint)=lfn0
1466 end if
1467 if ( lfn0*lfn1 < 0 ) then
1468 npoint = npoint + 1
1469 tt=lfn0/(lfn0-lfn1)
1470 x0=xx(k)+( xx(k+1)-xx(k) )*tt
1471 y0=yy(k)+( yy(k+1)-yy(k) )*tt
1472 xylist(npoint,1)=x0
1473 xylist(npoint,2)=y0
1474 tlist(npoint)=0 ! on fireline
1475 llist(npoint)=0
1476 end if
1477 end do
1479 ! make the list circular
1480 tlist(npoint+1)=tlist(1)
1481 llist(npoint+1)=llist(1)
1482 xylist(npoint+1,1)=xylist(1,1)
1483 xylist(npoint+1,2)=xylist(1,2)
1487 !print *,'after LS in case3'pproximation of the plane for lfn using finite
1488 !differences
1489 ! approximate L(x,y)=u(1)*x+u(2)*y+u(3)
1490 lfn_middle=(lfn00+lfn01+lfn10+lfn11)/4
1491 dt_dx=(lfn00-lfn10+lfn01-lfn11)/2
1492 dt_dy=(lfn00-lfn01+lfn10-lfn11)/2
1493 u(1)=dt_dx
1494 u(2)=dt_dy
1495 u(3)=lfn_middle-(dt_dx+dt_dy)/2
1496 % finding the coefficient c, reminder we work over one subcell only
1497 % T(x,y)=c*L(x,y)+time_now
1498 a=0
1499 b=0
1501 if (lfn00 <= 0) then
1502 a=a+lfn00*lfn00
1503 if (t00 > time_now) then
1504 call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
1505 else
1506 b=b+(t00-time_now)*lfn00
1507 end if
1508 end if
1511 if (lfn01 <= 0) then
1512 a=a+lfn01*lfn01
1513 if (t01>time_now) then
1514 call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
1515 else
1516 b=b+(t01-time_now)*lfn01
1517 end if
1518 end if
1521 if (lfn10<=0) then
1522 a=a+lfn10*lfn10
1523 if (t10>time_now) then
1524 call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
1525 else
1526 b=b+(t10-time_now)*lfn10
1527 end if
1528 end if
1530 if (lfn11<=0) then
1531 a=a+lfn11*lfn11
1532 if (t11>time_now) then
1533 call crash('fuel_burnt_fd: tign(i1) should be less then time_now')
1534 else
1535 b=b+(t11-time_now)*lfn11
1536 end if
1537 end if
1541 if (a==0) then
1542 call crash('fuel_burnt_fd: if c is on fire then one of cells should be on fire')
1543 end if
1544 c=b/a
1545 u(1)=u(1)*c
1546 u(2)=u(2)*c
1547 u(3)=u(3)*c
1549 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1551 !print *,'vecX from LS',vecX
1552 !print *,'tign inputed',tign00,tign10,tign11,tign01
1553 rnorm=snrm2(p,vecY,1)
1554 u(1)=vecX(1)
1555 u(2)=vecX(2)
1556 u(3)=vecX(3)
1557 ! rotate to gradient on x only
1558 nr = sqrt(u(1)**2+u(2)**2)
1559 if(.not.nr.gt.eps)then
1560 out=1.
1561 goto 900
1562 endif
1563 c = u(1)/nr
1564 s = u(2)/nr
1565 mQ(1,1)=c
1566 mQ(1,2)=s
1567 mQ(2,1)=-s
1568 mQ(2,2)=c
1569 ! mat vec multiplication
1570 call matvec(mQ,2,2,u,3,ut,2,2,2)
1571 errQ = ut(2) ! should be zero
1572 ae = -ut(1)/fuel_time_cell
1573 ce = -u(3)/fuel_time_cell
1574 cet=ce!keep ce
1575 call matmatp(xylist,8,2,mQ,2,2,xytlist,8,2,npoint+1,2,2)
1576 call sortxt( xytlist, 8,2, xt,8,npoint )
1577 out=0.0
1578 aupp=0.0
1579 cupp=0.0
1580 alow=0.0
1581 clow=0.0
1582 do k=1,npoint-1
1583 xt1=xt(k)
1584 xt2=xt(k+1)
1585 upper=0
1586 lower=0
1587 ah=0
1588 ch=0
1589 if ( xt2-xt1 > eps*100 ) then
1591 do ss=1,npoint
1592 xts=xytlist(ss,1)
1593 yts=xytlist(ss,2)
1594 xte=xytlist(ss+1,1)
1595 yte=xytlist(ss+1,2)
1597 if ( (xts>xt1 .and. xte>xt1) .or. &
1598 (xts<xt2 .and. xte<xt2) ) then
1599 aa = 0 ! do nothing
1600 cc = 0
1601 else
1602 aa = (yts-yte)/(xts-xte)
1603 cc = (xts*yte-xte*yts)/(xts-xte)
1604 if (xte<xts) then
1605 aupp = aa
1606 cupp = cc
1607 ah=ah+aa
1608 ch=ch+cc
1609 upper=upper+1
1610 else
1611 alow = aa
1612 clow = cc
1613 lower=lower+1
1614 end if
1615 end if!(xts>xt1 .and. xte>xt1)
1616 end do ! ss
1617 ce=cet !use stored ce
1618 if (ae*xt1+ce > 0 ) then
1619 ce=ce-(ae*xt1+ce)!shift small amounts exp(-**)
1620 end if
1621 if (ae*xt2+ce > 0) then
1622 ce=ce-(ae*xt2+ce)
1623 end if
1625 ah = aupp-alow
1626 ch = cupp-clow
1627 ! integrate (ah*x+ch)*(1-exp(ae*x+ce) from xt1 to xt2
1628 ! numerically sound for ae->0, ae -> infty
1629 ! this can be important for different model scales
1630 ! esp. if someone runs the model in single precision!!
1631 ! s1=int((ah*x+ch),x,xt1,xt2)
1632 s1 = (xt2-xt1)*((1./2.)*ah*(xt2+xt1)+ch)
1633 ! s2=int((ch)*(-exp(ae*x+ce)),x,xt1,xt2)
1634 ceae=ce/ae;
1635 s2 = -ch*exp(ae*(xt1+ceae))*(xt2-xt1)*intexp(ae*(xt2-xt1))
1636 ! s3=int((ah*x)*(-exp(ae*x+ce)),x,xt1,xt2)
1637 ! s3=int((ah*x)*(-exp(ae*(x+ceae))),x,xt1,xt2)
1638 ! expand in Taylor series around ae=0
1639 ! collect(expand(taylor(int(x*(-exp(ae*(x+ceae))),x,xt1,xt2)*ae^2,ae,4)/ae^2),ae)
1640 ! =(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
1641 ! + (-1/3*xt2^3-1/2*xt2^2*ceae+1/3*xt1^3+1/2*xt1^2*ceae)*ae
1642 ! + 1/2*xt1^2-1/2*xt2^2
1643 !
1644 ! coefficient at ae^2 in the expansion, after some algebra
1645 a2=(xt1-xt2)*((1./4.)*(xt1+xt2)*ceae**2+(1./3.)* &
1646 (xt1**2+xt1*xt2+xt2**2)*ceae+(1./8.)* &
1647 (xt1**3+xt1*(xt2**2)+xt1**2*xt2+xt2**3))
1648 d=(ae**4)*a2
1650 if (abs(d)>eps) then
1651 ! since ae*xt1+ce<=0 ae*xt2+ce<=0 all fine for large ae
1652 ! for ae, ce -> 0 rounding error approx eps/ae^2
1653 s3=( exp(ae*(xt1+ceae))*(ae*xt1-1)-&
1654 exp(ae*(xt2+ceae))*(ae*xt2-1) )/(ae**2)
1656 !we do not worry about rounding as xt1 -> xt2, then s3 -> 0
1657 else
1658 ! coefficient at ae^1 in the expansion
1659 a1=(xt1-xt2)*((1./2.)*ceae*(xt1+xt2)+(1./3.)*&
1660 (xt1**2+xt1*xt2+xt2**2))
1661 ! coefficient at ae^0 in the expansion for ae->0
1662 a0=(1./2.)*(xt1-xt2)*(xt1+xt2)
1663 s3=a0+a1*ae+a2*ae**2; ! approximate the integral
1664 end if
1665 ldskjflksdjflksdjf
1667 s3=ah*s3
1668 out=out+s1+s2+s3
1669 out=1-out !fuel_left
1670 if(out<0 .or. out>1) then
1671 print *,':fuel_fraction should be between 0 and 1'
1672 !print *, 'eps= ', eps
1673 !print *, 'xt1= ', xt1, 'xt2= ', xt2
1674 !print *,'ae= ',ae,'ce= ',ce,'ah= ',ah,'ch= ',ch
1675 !print *,'a0= ', a0,'a1= ', a1,'a2= ', a2
1676 print *,'s1= ', s1,'s2= ', s2,'s3= ', s3
1677 print *,':fuel_fraction =',out
1678 end if!print
1680 end if
1681 end do ! k
1683 end if ! if case0, elseif case4 ,else case123
1685 900 continue
1686 if(out>1. .or. out<0.)call crash('fuel_left_cell_3: fuel fraction out of bounds [0,1]')
1687 fuel_left_cell_3 = out
1688 end function fuel_left_cell_3
1695 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1696 real function intexp(ab)
1697 implicit none
1698 real::ab
1699 !calls
1700 intrinsic epsilon
1702 real, parameter:: zero=0.,one=1.,eps=epsilon(zero)
1704 !eps = 2.2204*(10.0**(-8))!from matlab
1705 if ( eps < abs(ab)**3/6. ) then
1706 intexp=(exp(ab)-1)/ab
1707 else
1708 intexp=1+ab/2.
1709 end if
1710 end function intexp
1711 !
1712 !****************************************
1713 !
1714 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1715 !
1716 !****************************************
1717 !
1720 !
1721 !****************************************
1727 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1728 subroutine sortxt(xytlist,nrow,ncolumn,xt,nxt,nvec)
1729 implicit none
1730 integer::nrow,ncolumn,nxt,nvec
1731 real,dimension(nrow,ncolumn)::xytlist
1732 real,dimension(nxt)::xt
1734 integer::i,j
1735 real::temp
1737 do i=1,nvec
1738 xt(i)=xytlist(i,1)
1739 end do
1741 do i=1,nvec-1
1742 do j=i+1,nvec
1743 if ( xt(i) > xt(j) ) then
1744 temp = xt(i)
1745 xt(i)=xt(j)
1746 xt(j)=temp
1747 end if
1748 end do
1749 end do
1751 end subroutine !sortxt
1752 !
1753 !****************************************
1754 !
1755 subroutine matvec(A,m,n,V,nv,out,nout,nrow,ncolumn)
1756 implicit none
1757 integer::m,n,nv,nout,nrow,ncolumn
1758 real,dimension(m,n)::A ! allocated m by n
1759 real,dimension(nv)::V ! allocated nv
1760 real,dimension(nout)::out! allocated nout
1762 integer::i,j
1764 do i=1,nrow
1765 out(i)=0.0
1766 do j=1,ncolumn
1767 out(i)=out(i)+A(i,j)*V(j)
1768 end do
1769 end do
1770 end subroutine
1771 !
1772 !****************************************
1773 !
1774 subroutine matmatp(A,mA,nA,B,mB,nB,C,mC,nC,nrow,ncolumn,nP)
1775 implicit none
1776 integer::mA,nA,mB,nB,mC,nC,nrow,ncolumn,nP
1777 real,dimension(mA,nA)::A ! allocated m by n
1778 real,dimension(mB,nB)::B ! allocated m by n
1779 real,dimension(mC,nC)::C ! allocated m by n
1780 integer::i,j,k
1781 do i=1,nrow
1782 do j=1,ncolumn
1783 C(i,j)=0.0
1784 do k=1,nP
1785 C(i,j)=C(i,j)+A(i,k)*B(j,k) ! B'
1786 end do
1787 end do
1788 end do
1789 end subroutine
1791 !
1792 !****************************************
1793 !
1794 #endif
1796 subroutine prop_ls( id, & ! for debug
1797 ids,ide,jds,jde, & ! domain dims
1798 ims,ime,jms,jme, & ! memory dims
1799 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1800 its,ite,jts,jte, & ! tile dims
1801 ts,dt,dx,dy, & ! scalars in
1802 tbound, & ! scalars out
1803 lfn_in,lfn_out,tign,ros, & ! arrays inout
1804 fp &
1805 )
1806 implicit none
1808 !*** purpose: advance level function in time
1810 ! Jan Mandel August 2007 - February 2008
1812 !*** description
1813 !
1814 ! Propagation of closed curve by a level function method. The level function
1815 ! lfn is defined by its values at the nodes of a rectangular grid.
1816 ! The area where lfn < 0 is inside the curve. The curve is
1817 ! described implicitly by lfn=0. Points where the curve intersects gridlines
1818 ! can be found by linear interpolation from nodes.
1819 !
1820 ! The level function is advanced from time ts to time ts + dt.
1821 !
1822 ! The level function should be initialized to (an approximation of) the signed
1823 ! distance from the curve. If the initial curve is a circle, the initial level
1824 ! function is simply the distance from the center minus the radius.
1825 !
1826 ! The curve moves outside with speed given by function speed_func.
1827 !
1828 ! Method: Godunov/ENO method for the normal motion. The timestep is checked for
1829 ! CFL condition. For a straight segment in a constant field and locally linear
1830 ! level function, the method reduces to the exact normal motion. The advantage of
1831 ! the level set method is that it treats automatically special cases such as
1832 ! the curve approaching itself and merging components of the area inside the curve.
1833 !
1834 ! Based on S. Osher and R. Fedkiw, Level set methods and dynamic implicit surfaces,
1835 ! Springer, 2003, Sec. 6.4, as implemented in toolboxLS for Matlab by
1836 ! I. Mitchell, A toolbox of Level Set Methods (Version 1.1), TR-2007-11,
1837 ! Dept. Computer Science, University of British Columbia, 2007
1838 ! http://www.cs.ubc.ca/\~mitchell/Toolbo\LS
1839 !
1841 !*** arguments
1843 ! id in unique identification for prints and dumps
1844 ! ids,ide,jds,jde in domain dimensions
1845 ! ims,ime,jms,jme in memory dimensions
1846 ! its,ite,jts,jte in tile dimensions
1847 ! ts in start time
1848 ! dt in time step
1849 ! dx,dy in grid spacing
1850 ! tbound out bound on stable time step from CFL condition, if tbound>=dt then OK
1851 ! lfn_in,lfn_out inout,out the level set function at nodes
1852 ! tign inout the ignition time at nodes
1854 ! The dimensions are cell-based, the nodal value is associated with the south-west corner.
1855 ! The whole computation is on domain indices ids:ide+1,jds:jde+1.
1856 !
1857 ! The region where new lfn and tign are computed is the tile its:ite,jts:jte
1858 ! except when the tile is at domain upper boundary, an extra band of points is added:
1859 ! if ite=ide then region goes up to ite+1, if jte=jde then region goes up to jte+1.
1861 ! The time step requires values from 2 rows of nodes beyond the region except when at the
1862 ! domain boundary one-sided derivatives are used. This is implemented by extending the input
1863 ! beyond the domain boundary so sufficient memory bounds must be allocated.
1864 ! The update on all tiles can be done in parallel. To avoid the race condition (different regions
1865 ! of the same array updated by different threads), the in and out versions of the
1866 ! arrays lft and tign are distinct. If the time step dt is larger
1867 ! that the returned tbound, the routine should be called again with timestep td<=tbound, and then
1868 ! having distinct inputs and outputs comes handy.
1870 !*** calls
1871 !
1872 ! tend_ls
1873 !
1875 integer,intent(in)::id,ims,ime,jms,jme,ids,ide,jds,jde,its,ite,jts,jte,ips,ipe,jps,jpe
1876 real,dimension(ims:ime,jms:jme),intent(inout)::lfn_in,tign
1877 real,dimension(ims:ime,jms:jme),intent(out)::lfn_out,ros
1878 real,intent(in)::dx,dy,ts,dt
1879 real,intent(out)::tbound
1880 type(fire_params),intent(in)::fp
1882 !*** local
1883 ! arrays
1884 #define IMTS its-1
1885 #define IMTE ite+1
1886 #define JMTS jts-1
1887 #define JMTE jte+1
1888 real,dimension(IMTS:IMTE,JMTS:JMTE):: tend, lfn1 ! region-sized with halo
1889 ! scalars
1890 real::grad2,rr,tbound2,a,a1 ! a=0 euler, a=0.5 heun
1892 real::gradx,grady,aspeed,err,aerr,time_now
1893 integer::ihs,ihe,jhs,jhe
1894 integer::ihs2,ihe2,jhs2,jhe2
1895 integer::i,j,its1,ite1,jts1,jte1,k,kk,id1
1896 character(len=128)msg
1897 integer::nfirenodes,nfireline
1898 real::sum_err,min_err,max_err,sum_aerr,min_aerr,max_aerr
1900 ! constants
1901 integer,parameter :: mstep=1000, printl=1
1902 real, parameter:: zero=0.,one=1.,eps=epsilon(zero),tol=100*eps, &
1903 safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
1905 ! f90 intrinsic function
1907 intrinsic max,min,sqrt,nint,epsilon,tiny,huge
1909 !*** executable
1911 !$OMP CRITICAL(SFIRE_CORE_CRIT)
1912 write(msg,'(a8,i5,a6,i5,3(a1,i5))')'prop_ls:',id,' tile ',its,':',ite,',',jts,':',jte
1913 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
1914 call message(msg)
1916 a=fire_back_weight ! from module_fr_sfire_util
1917 a1=1. - a
1919 ! tend = F(lfn)
1921 ihs2=max(its-2,ids) ! need lfn two beyond the tile but not outside the domain
1922 ihe2=min(ite+2,ide)
1923 jhs2=max(jts-2,jds)
1924 jhe2=min(jte+2,jde)
1926 ihs=max(its-1,ids) ! compute tend one beyond the tile but not outside the domain
1927 ihe=min(ite+1,ide)
1928 jhs=max(jts-1,jds)
1929 jhe=min(jte+1,jde)
1931 #ifdef DEBUG_OUT
1932 call write_array_m(ihs,ihe,jhs,jhe,ims,ime,jms,jme,lfn_in,'lfn_in',id)
1933 #endif
1935 ! check array dimensions
1936 call check_mesh_2dim(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme)
1937 call print_2d_stats(ihs2,ihe2,jhs2,jhe2,ims,ime,jms,jme, &
1938 lfn_in,'prop_ls: lfn in')
1940 ! NOTE: tend_ls will extrapolate to one node strip at domain boundaries
1941 ! so that it can compute gradient at domain boundaries.
1942 ! To avoid copying, lfn_in is declared inout.
1943 ! At tile boundaries that are not domain boundaries values of lfn_in two nodes
1944 ! outside of the tile are needed.
1945 id1 = id ! for debug prints
1946 if(id1.ne.0)id1=id1+1000
1947 call tend_ls( id1, &
1948 ims,ime,jms,jme, & ! memory dims for lfn_in
1949 IMTS,IMTE,JMTS,JMTE, & ! memory dims for tend
1950 ids,ide,jds,jde, & ! domain dims - where lfn exists
1951 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1952 ihs,ihe,jhs,jhe, & ! where tend computed
1953 ims,ime,jms,jme, & ! memory dims for ros
1954 its,ite,jts,jte, & ! tile dims - where to set ros
1955 ts,dt,dx,dy, & ! scalars in
1956 lfn_in, & ! arrays in
1957 tbound, & ! scalars out
1958 tend, ros, & ! arrays out
1959 fp & ! params
1960 )
1962 #ifdef DEBUG_OUT
1963 call write_array_m(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE,tend,'tend1',id)
1964 #endif
1966 ! Euler method, the half-step, same region as ted
1967 do j=jhs,jhe
1968 do i=ihs,ihe
1969 lfn1(i,j) = lfn_in(i,j) + dt*tend(i,j)
1970 enddo
1971 enddo
1973 call print_2d_stats(ihs,ihe,jhs,jhe,IMTS,IMTE,JMTS,JMTE, &
1974 lfn1,'prop_ls: lfn1')
1975 ! tend = F(lfn1) on the tile (not beyond)
1977 if(id1.ne.0)id1=id1+1000
1978 call tend_ls( id1,&
1979 IMTS,IMTE,JMTS,JMTE, & ! memory dims for lfn
1980 IMTS,IMTE,JMTS,JMTE, & ! memory dims for tend
1981 ids,ide,jds,jde, & ! domain dims - where lfn exists
1982 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
1983 its,ite,jts,jte, & ! tile dims - where is tend computed
1984 ims,ime,jms,jme, & ! memory dims for ros
1985 its,ite,jts,jte, & ! tile dims - where is ros set
1986 ts+dt,dt,dx,dy, & ! scalars in
1987 lfn1, & ! arrays in
1988 tbound2, & ! scalars out
1989 tend,ros, & ! arrays out
1990 fp &
1991 )
1993 #ifdef DEBUG_OUT
1994 call write_array_m(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'tend2',id)
1995 #endif
1997 call print_2d_stats(its,ite,jts,jte,IMTS,IMTE,JMTS,JMTE,tend,'prop_ls: tend2')
1999 tbound=min(tbound,tbound2)
2001 !$OMP CRITICAL(SFIRE_CORE_CRIT)
2002 write(msg,'(a,f10.2,4(a,f7.2))')'prop_ls: time',ts,' dt=',dt,' bound',min(tbound,999.99), &
2003 ' dx=',dx,' dy=',dy
2004 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
2005 call message(msg)
2006 if(dt>tbound)then
2007 !$OMP CRITICAL(SFIRE_CORE_CRIT)
2008 write(msg,'(2(a,f10.2))')'prop_ls: WARNING: time step ',dt, &
2009 ' > bound =',tbound
2010 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
2011 call message(msg)
2012 endif
2014 ! combine lfn1 and lfn_in + dt*tend -> lfn_out
2016 do j=jts,jte
2017 do i=its,ite
2018 lfn_out(i,j) = a1*lfn1(i,j) + a*(lfn_in(i,j) + dt*tend(i,j))
2019 enddo
2020 enddo
2022 ! compute ignition time by interpolation
2023 ! the node was not burning at start but it is burning at end
2024 ! interpolate from the level functions at start and at end
2025 ! lfn_in is the level set function value at time ts
2026 ! lfn_out is the level set function value at time ts+dt
2027 ! 0 is the level set function value at time tign(i,j)
2028 ! thus assuming the level function is approximately linear =>
2029 ! tign(i,j)= ts + ((ts + td) - ts) * lfn_in / (lfn_in - lfn_out)
2030 ! = ts + dt * lfn_in / (lfn_in - lfn_out)
2032 time_now=ts+dt
2033 time_now = time_now + abs(time_now)*epsilon(time_now)*2.
2034 do j=jts,jte
2035 do i=its,ite
2036 ! interpolate the cross-over time
2037 if (.not. lfn_out(i,j)>0 .and. lfn_in(i,j)>0)then
2038 tign(i,j) = ts + dt * lfn_in(i,j) / (lfn_in(i,j) - lfn_out(i,j))
2039 endif
2040 ! set the ignition time outside of burning region
2041 if(lfn_out(i,j)>0.)tign(i,j)=time_now
2042 enddo
2043 enddo
2045 ! check local speed error and stats
2046 ! may not work correctly in parallel
2047 ! init stats
2048 nfirenodes=0
2049 nfireline=0
2050 sum_err=0.
2051 min_err=rmax
2052 max_err=rmin
2053 sum_aerr=0.
2054 min_aerr=rmax
2055 max_aerr=rmin
2056 its1=its+1
2057 jts1=jts+1
2058 ite1=ite-1
2059 jte1=jte-1
2060 ! loop over right inside of the domain
2061 ! cannot use values outside of the domain, would have to reflect and that
2062 ! would change lfn_out
2063 ! cannot use values outside of tile, not synchronized yet
2064 ! so in parallel mode the statistics is not accurate
2065 do j=jts1,jte1
2066 do i=its1,ite1
2067 if(lfn_out(i,j)>0.0)then ! a point out of burning region
2068 if(lfn_out(i+1,j)<=0.or.lfn_out(i,j+1)<=0.or. & ! neighbor in burning region
2069 lfn_out(i-1,j)<=0.or.lfn_out(i,j-1)<=0)then ! point next to fireline
2070 gradx=(lfn_out(i+1,j)-lfn_out(i-1,j))/(2.0*dx) ! central differences
2071 grady=(lfn_out(i,j+1)-lfn_out(i,j-1))/(2.0*dy)
2072 grad2=sqrt(gradx*gradx+grady*grady)
2073 aspeed = (lfn_in(i,j)-lfn_out(i,j))/(dt*max(grad2,rmin))
2074 rr = speed_func(gradx,grady,dx,dy,i,j,fp)
2075 err=aspeed-rr
2076 sum_err=sum_err+err
2077 min_err=min(min_err,err)
2078 max_err=max(max_err,err)
2079 aerr=abs(err)
2080 sum_aerr=sum_aerr+aerr
2081 min_aerr=min(min_aerr,aerr)
2082 max_aerr=max(max_aerr,aerr)
2083 nfireline=nfireline+1
2084 endif
2085 else
2086 nfirenodes=nfirenodes+1
2087 endif
2088 enddo
2089 enddo
2090 !$OMP CRITICAL(SFIRE_CORE_CRIT)
2091 write(msg,'(2(a,i6,f8.4))')'prop_ls: nodes burning',nfirenodes, &
2092 (100.*nfirenodes)/((ite1-its1+1)*(jte1-jts1+1)),'% next to fireline',nfireline
2093 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
2094 call message(msg)
2095 if(nfireline>0)then
2096 call print_stat_line('speed error',its1,ite1,jts1,jte1,min_err,max_err,sum_err/nfireline)
2097 call print_stat_line('abs(speed error)',its1,ite1,jts1,jte1,min_aerr,max_aerr,sum_aerr/nfireline)
2098 endif
2100 ! check if the fire did not get to the domain boundary
2101 do k=-1,1,2
2102 ! do kk=1,(boundary_guard*(ide-ids+1))/100 ! in %
2103 do kk=1,boundary_guard ! measured in cells
2104 i=ids+k*kk
2105 if(i.ge.its.and.i.le.ite)then
2106 do j=jts,jte
2107 if(lfn_out(i,j)<=0.)goto 9
2108 enddo
2109 endif
2110 enddo
2111 ! do kk=1,(boundary_guard*(jde-jds+1))/100
2112 do kk=1,boundary_guard ! measured in cells
2113 j=jds+k*kk
2114 if(j.ge.jts.and.j.le.jte)then
2115 do i=its,ite
2116 if(lfn_out(i,j)<=0.)goto 9
2117 enddo
2118 endif
2119 enddo
2120 enddo
2121 goto 10
2122 9 continue
2123 !$OMP CRITICAL(SFIRE_CORE_CRIT)
2124 write(msg,'(a,i2,a,2i8)')'prop_ls: fire',boundary_guard, &
2125 ' cells from domain boundary at node ',i,j
2126 !$OMP END CRITICAL(SFIRE_CORE_CRIT)
2127 call message(msg)
2128 call crash('prop_ls: increase the fire region')
2129 10 continue
2131 call print_2d_stats(its,ite,jts,jte,ims,ime,jms,jme, &
2132 lfn_out,'prop_ls: lfn out')
2134 end subroutine prop_ls
2136 !
2137 !*****************************
2138 !
2140 subroutine tend_ls( id, &
2141 lims,lime,ljms,ljme, & ! memory dims for lfn
2142 tims,time,tjms,tjme, & ! memory dims for tend
2143 ids,ide,jds,jde, & ! domain - nodes where lfn defined
2144 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
2145 ints,inte,jnts,jnte, & ! region - nodes where tend computed
2146 ims,ime,jms,jme, & ! memory dims for ros
2147 its,ite,jts,jte, & ! tile dims - where is ros set
2148 t,dt,dx,dy, & ! scalars in
2149 lfn, & ! arrays in
2150 tbound, & ! scalars out
2151 tend, ros, & ! arrays out
2152 fp &
2153 )
2155 implicit none
2156 ! purpose
2157 ! compute the right hand side of the level set equation
2159 !*** arguments
2160 integer,intent(in)::id,lims,lime,ljms,ljme,tims,time,tjms,tjme
2161 integer,intent(in)::ims,ime,jms,jme,its,ite,jts,jte
2162 integer, intent(in)::ids,ide,jds,jde,ints,inte,jnts,jnte,ips,ipe,jps,jpe
2163 real,intent(in)::t ! time
2164 real,intent(in)::dt,dx,dy ! mesh step
2165 real,dimension(lims:lime,ljms:ljme),intent(inout)::lfn ! level set function
2166 real,intent(out)::tbound ! max allowed time step
2167 real,dimension(tims:time,tjms:tjme),intent(out)::tend ! tendency (rhs of the level set pde)
2168 real,dimension(ims:ime,jms:jme),intent(out)::ros ! rate of spread
2169 type(fire_params),intent(in)::fp
2171 !*** local
2172 real:: te,diffLx,diffLy,diffRx,diffRy, &
2173 diffCx,diffCy,diff2x,diff2y,grad,rr, &
2174 ros_back,ros_wind,ros_slope,advx,advy,scale,nvx,nvy,speed,tanphi
2175 integer::i,j
2176 character(len=128)msg
2178 ! constants
2179 real, parameter:: eps=epsilon(0.0)
2180 !intrinsic epsilon
2181 real, parameter:: zero=0.,one=1.,tol=100*eps, &
2182 safe=2.,rmin=safe*tiny(zero),rmax=huge(zero)/safe
2185 ! f90 intrinsic function
2187 intrinsic max,min,sqrt,nint,tiny,huge
2190 #ifdef DEBUG_OUT
2191 real,dimension(tims:time,tjms:tjme)::rra,grada,speeda,tanphia
2192 #endif
2194 !*** executable
2196 ! check array dimensions
2197 call check_mesh_2dim(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme)
2198 call check_mesh_2dim(ints,inte,jnts,jnte,tims,time,tjms,tjme)
2200 call continue_at_boundary(1,1,fire_lfn_ext_up, & !extend by extrapolation but never down
2201 lims,lime,ljms,ljme, & ! memory dims
2202 ids,ide,jds,jde, & ! domain - nodes where lfn defined
2203 ips,ipe,jps,jpe, & ! patch - nodes owned by this process
2204 ints,inte,jnts,jnte, & ! tile - nodes update by this thread
2205 lfn) ! array
2207 call print_2d_stats(ints-1,inte+1,jnts,jnte,lims,lime,ljms,ljme, &
2208 lfn,'tend_ls: lfn cont dir x')
2209 call print_2d_stats(ints,inte,jnts-1,jnte+1,lims,lime,ljms,ljme, &
2210 lfn,'tend_ls: lfn cont dir y')
2212 #ifdef DEBUG_OUT
2213 call write_array_m(ints-1,inte+1,jnts-1,jnte+1,lims,lime,ljms,ljme,lfn,'tend_lfn_in',id)
2214 #endif
2216 tbound=0
2217 do j=jnts,jnte
2218 do i=ints,inte
2219 ! one sided differences
2220 diffRx = (lfn(i+1,j)-lfn(i,j))/dx
2221 diffLx = (lfn(i,j)-lfn(i-1,j))/dx
2222 diffRy = (lfn(i,j+1)-lfn(i,j))/dy
2223 diffLy = (lfn(i,j)-lfn(i,j-1))/dy
2224 diffCx = diffLx+diffRx ! TWICE CENTRAL DIFFERENCE
2225 diffCy = diffLy+diffRy
2227 !upwinding - select right or left derivative
2228 select case(fire_upwinding)
2229 case(0) ! none
2230 grad=sqrt(diffCx**2 + diffCy**2)
2231 case(1) ! standard
2232 diff2x=select_upwind(diffLx,diffRx)
2233 diff2y=select_upwind(diffLy,diffRy)
2234 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
2235 case(2) ! godunov per osher/fedkiw
2236 diff2x=select_godunov(diffLx,diffRx)
2237 diff2y=select_godunov(diffLy,diffRy)
2238 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
2239 case(3) ! eno
2240 diff2x=select_eno(diffLx,diffRx)
2241 diff2y=select_eno(diffLy,diffRy)
2242 grad=sqrt(diff2x*diff2x + diff2y*diff2y)
2243 case(4) ! Sethian - twice stronger pushdown of bumps
2244 grad=sqrt(max(diffLx,0.)**2+min(diffRx,0.)**2 &
2245 + max(diffLy,0.)**2+min(diffRy,0.)**2)
2246 case default
2247 grad=0.
2248 end select
2250 ! normal direction, from central differences
2251 scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps)
2252 nvx=diffCx/scale
2253 nvy=diffCy/scale
2255 ! wind speed in direction of spread
2256 ! speed = vx(i,j)*nvx + vy(i,j)*nvy
2259 ! get rate of spread from wind speed and slope
2261 call fire_ros(ros_back,ros_wind,ros_slope, &
2262 nvx,nvy,i,j,fp)
2264 rr=ros_back + ros_wind + ros_slope
2265 if(fire_grows_only.gt.0)rr=max(rr,0.)
2267 ! set ros for output
2268 if(i.ge.its.and.i.le.ite.and.j.ge.jts.and.j.le.jte)ros(i,j)=rr
2270 if(fire_upwind_split.eq.0)then
2272 ! get rate of spread
2273 te = -rr*grad ! normal term
2275 else
2277 ! normal direction backing rate only
2278 te = - ros_back*grad
2280 ! advection in wind direction
2281 if (abs(speed)> eps) then
2282 advx=fp%vx(i,j)*ros_wind/speed
2283 advy=fp%vy(i,j)*ros_wind/speed
2284 else
2285 advx=0
2286 advy=0
2287 endif
2289 ! tanphi = dzdxf*nvx + dzdyf*nvy
2290 ! advection from slope direction
2291 if(abs(tanphi)>eps) then
2292 advx=advx+fp%dzdxf(i,j)*ros_slope/tanphi
2293 advy=advy+fp%dzdyf(i,j)*ros_slope/tanphi
2294 endif
2296 if(fire_upwind_split.eq.1)then
2298 ! one-sided upwinding
2299 te = te - max(advx,0.)*diffLx - min(advx,0.)*diffRy &
2300 - max(advy,0.)*diffLy - min(advy,0.)*diffRy
2303 elseif(fire_upwind_split.eq.2)then
2305 ! Lax-Friedrichs
2306 call crash('prop_ls: bad fire_upwind_split, Lax-Friedrichs not done yet')
2308 else
2310 call crash('prop_ls: bad fire_upwind_split')
2312 endif
2313 endif
2315 ! cfl condition
2316 if (grad > 0.) then
2317 tbound = max(tbound,rr*(abs(diff2x)/dx+abs(diff2y)/dy)/grad)
2318 endif
2320 ! add numerical viscosity
2321 te=te + fire_viscosity*abs(rr)*((diffRx-diffLx)+(diffRy-diffLy))
2323 tend(i,j)=te
2324 #ifdef DEBUG_OUT
2325 rra(i,j)=rr
2326 grada(i,j)=grad
2327 speeda(i,j)=speed
2328 tanphia(i,j)=tanphi
2329 #endif
2330 !write(msg,*)i,j,grad,dzdx,dzdy
2331 !call message(msg)
2333 !if(abs(te)>1e4)then
2334 ! write(msg,'(a,2i5,e12.3)')'tend_ls: tend out of bounds at ',i,j,te
2335 ! call crash(msg)
2336 !endif
2337 enddo
2338 enddo
2340 #ifdef DEBUG_OUT
2341 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,rra,'rr',id)
2342 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,grada,'grad',id)
2343 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,speeda,'speed',id)
2344 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,tanphia,'tanphi',id)
2345 call write_array_m(ints,inte,jnts,jnte,tims,time,tjms,tjme,tend,'tend',id)
2346 #endif
2348 call print_2d_stats(ints,inte,jnts,jnte,tims,time,tjms,tjme, &
2349 tend,'tend_ls: tend out')
2351 ! the final CFL bound
2352 tbound = 1/(tbound+tol)
2354 end subroutine tend_ls
2356 !
2357 !**************************
2358 !
2360 real function select_upwind(diffLx,diffRx)
2361 implicit none
2362 real, intent(in):: diffLx, diffRx
2363 real diff2x
2365 ! upwind differences, L or R if bith same sign, otherwise zero
2367 diff2x=0
2368 if (diffLx>0.and.diffRx>0.)diff2x=diffLx
2369 if (diffLx<0.and.diffRx<0.)diff2x=diffRx
2371 select_upwind=diff2x
2372 end function select_upwind
2375 !
2376 !**************************
2377 !
2380 real function select_godunov(diffLx,diffRx)
2381 implicit none
2382 real, intent(in):: diffLx, diffRx
2383 real diff2x,diffCx
2385 ! Godunov scheme: upwind differences, L or R or none
2386 ! always test on > or < never = , much faster because of IEEE
2387 ! central diff >= 0 => take left diff if >0, ortherwise 0
2388 ! central diff <= 0 => take right diff if <0, ortherwise 0
2390 diff2x=0
2391 diffCx=diffRx+diffLx
2392 if (diffLx>0.and..not.diffCx<0)diff2x=diffLx
2393 if (diffRx<0.and. diffCx<0)diff2x=diffRx
2395 select_godunov=diff2x
2396 end function select_godunov
2398 !
2399 !**************************
2400 !
2402 real function select_eno(diffLx,diffRx)
2403 implicit none
2404 real, intent(in):: diffLx, diffRx
2405 real diff2x
2407 ! 1st order ENO scheme
2409 if (.not.diffLx>0 .and. .not.diffRx>0)then
2410 diff2x=diffRx
2411 elseif(.not.diffLx<0 .and. .not.diffRx<0)then
2412 diff2x=diffLx
2413 elseif(.not.diffLx<0 .and. .not.diffRx>0)then
2414 if(.not. abs(diffRx) < abs(diffLx))then
2415 diff2x=diffRx
2416 else
2417 diff2x=diffLx
2418 endif
2419 else
2420 diff2x=0.
2421 endif
2423 select_eno=diff2x
2424 end function select_eno
2426 !
2427 !**************************
2428 !
2430 real function speed_func(diffCx,diffCy,dx,dy,i,j,fp)
2431 !*** purpose
2432 ! the level set method speed function
2433 implicit none
2434 !*** arguments
2435 real, intent(in)::diffCx,diffCy ! x and y coordinates of the direction of propagation
2436 real, intent(in)::dx,dy ! x and y coordinates of the direction of propagation
2437 integer, intent(in)::i,j ! indices of the node to compute the speed at
2438 type(fire_params),intent(in)::fp
2439 !*** local
2440 real::scale,nvx,nvy,r
2441 real::ros_back , ros_wind , ros_slope
2442 real, parameter:: eps=epsilon(0.0)
2443 !*** executable
2444 ! normal direction, from central differences
2445 scale=sqrt(diffCx*diffCx+diffCy*diffCy+eps)
2446 nvx=diffCx/scale
2447 nvy=diffCy/scale
2449 ! get rate of spread from wind speed and slope
2451 call fire_ros(ros_back,ros_wind,ros_slope, &
2452 nvx,nvy,i,j,fp)
2454 r=ros_back + ros_wind + ros_slope
2455 if(fire_grows_only.gt.0)r=max(r,0.)
2456 speed_func=r
2458 end function speed_func
2460 end module module_fr_sfire_core