source: lmdz_wrf/trunk/WRFV3/phys/module_fr_sfire_model.F @ 1404

Last change on this file since 1404 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 17.4 KB
Line 
1!
2#define DEBUG_OUT
3
4module module_fr_sfire_model
5
6use module_fr_sfire_core
7use module_fr_sfire_util
8use module_fr_sfire_phys
9
10contains
11
12subroutine sfire_model (                    &
13    id,                                     & ! unique number for prints and debug
14    ifun,                                   & ! what to do see below
15    restart,                                &
16    need_lfn_update,                          & ! if lfn needs to be synced between tiles
17    num_ignitions,                          & ! number of ignitions before advancing
18    ifuelread,nfuel_cat0,                   & ! initialize fuel categories
19    ifds,ifde,jfds,jfde,                    & ! fire domain dims - the whole domain
20    ifms,ifme,jfms,jfme,                    & ! fire memory dims - how declared
21    ifps,ifpe,jfps,jfpe,                    & ! patch - nodes owned by this process
22    ifts,ifte,jfts,jfte,                    & ! fire tile dims  - this thread
23    time_start,dt,                          & ! time and increment
24    fdx,fdy,                                & ! fire mesh spacing,
25    ignition_line,                          & ! small array of ignition line descriptions
26    ignitions_done,ignited_tile,            &
27    coord_xf,coord_yf,unit_xf,unit_yf,      & ! fire mesh coordinates
28    lfn,lfn_out,tign,fuel_frac,fire_area,   & ! state: level function, ign time, fuel left, area burning
29    grnhfx,grnqfx,                          & ! output: heat fluxes
30    ros,                                    & ! output: rate of spread
31    nfuel_cat,                              & ! fuel data per point
32    fuel_time,                              & ! save derived internal data
33    fp &
34)
35
36! This subroutine implements the fire spread model.
37! All quantities are on the fire grid. It inputs
38! winds given on the nodes of the fire grid
39! and outputs the heat fluxes on the cells of the fire grid.
40! This subroutine has no knowledge of any atmospheric model.
41! This code was written to conform with the WRF parallelism model, however it
42! does not depend on it. It can be called with domain equal to tile.
43! Wind and height must be given on 1 more node beyond the domain bounds.
44! The subroutine changes only array entries of the arguments in the tile.
45! Upon exit with ifun=2 (time step), lfn_out is to be copied into lfn by the caller.
46! When this subroutine is used on separate tiles that make a domain the value, the
47! it uses lfn on a strip of width 2 from neighboring tiles.
48!
49! All computation is done on one tile.
50!
51! This subroutine is intended to be called in a loop like
52!
53!
54! do ifun=1,6 (if initizalize run, otherwise 3,6)
55!   start parallel loop over tiles
56!       if ifun=1, set z and fuel data
57!       if ifun=3, set the wind arrays
58!       call sfire_model(....)
59!   end parallel loop over tiles
60!
61!   if need_lfn_update, halo exchange on lfn width 2
62!   
63!   if ifun=0
64!       halo exchange on z width 2
65!       halo exchange on fuel data width 1
66!   endif
67!   
68!   if ifun=3, halo exchange on winds width 2
69!   
70! enddo
71
72implicit none
73
74!*** arguments
75
76! control switches
77integer, intent(in) :: id
78integer, intent(in) :: ifun                 ! 1 = initialize run pass 1
79                                            ! 2 = initialize run pass 2
80                                            ! 3 = initialize timestep
81                                            ! 4 = do one timestep
82                                            ! 5 = copy timestep output to input
83                                            ! 6 = compute output fluxes
84logical, intent(in):: restart               ! if true, use existing state
85logical, intent(out)::need_lfn_update       ! if true, halo update on lfn afterwards
86! scalar data
87integer, intent(in) :: num_ignitions        ! number of ignition lines
88integer, intent(in) :: ifuelread,nfuel_cat0 ! for set_fire_params
89integer, intent(in) :: ifds,ifde,jfds,jfde,&  ! fire domain bounds
90        ifps,ifpe,jfps,jfpe                ! patch - nodes owned by this process
91integer, intent(in) :: ifts,ifte,jfts,jfte  ! fire tile bounds
92integer, intent(in) :: ifms,ifme,jfms,jfme  ! fire memory array bounds
93REAL,INTENT(in) :: time_start,dt            ! starting time, time step
94REAL,INTENT(in) :: fdx,fdy                  ! spacing of the fire mesh
95! array data
96type(ignition_line_type), dimension (num_ignitions), intent(in):: ignition_line ! descriptions of ignition lines
97integer, intent(out):: ignited_tile(num_ignitions),ignitions_done
98real, dimension(ifms:ifme, jfms:jfme), intent(in):: &
99    coord_xf,coord_yf                       !  node coordinates 
100real, intent(in):: unit_xf,unit_yf          !  coordinate units in m
101   
102! state
103REAL, INTENT(inout), dimension(ifms:ifme,jfms:jfme):: &
104    lfn   , &                               ! level function: fire is where lfn<0 (node)
105    tign  , &                               ! absolute time of ignition (node)
106    fuel_frac                               ! fuel fraction (node), currently redundant
107
108REAL, INTENT(out), dimension(ifms:ifme,jfms:jfme):: &
109    fire_area                               ! fraction of each cell burning
110   
111! output
112REAL, INTENT(out), dimension(ifms:ifme,jfms:jfme):: &
113    lfn_out, &                              !                             
114    grnhfx,grnqfx, &                        ! heat fluxes J/m^2/s  (cell)             
115    ros                                     ! output: rate of spread
116 
117! constant arrays - set at initialization
118real, intent(inout), dimension(ifms:ifme, jfms:jfme)::nfuel_cat ! cell based, data, constant
119real,intent(inout),dimension(ifms:ifme,jfms:jfme):: fuel_time
120type(fire_params),intent(inout)::fp
121
122!*** local
123
124integer :: xifms,xifme,xjfms,xjfme  ! memory bounds for pass-through arguments to normal spread
125real, dimension(ifts:ifte,jfts:jfte)::fuel_frac_burnt,fuel_frac_end
126integer::ignited,ig,i,j,itso,iteo,jtso,jteo
127real::tbound,err,erri,errj,maxgrad,grad,tfa,thf,mhf,tqf,mqf,aw,mw
128character(len=128)::msg
129logical:: freeze_fire
130integer:: stat_lev=1
131
132!*** executable
133
134call check_mesh_2dim(ifts-1,ifte+1,jfts-1,jfte+1,ifms,ifme,jfms,jfme)
135
136xifms=ifms  ! dimensions for the include file
137xifme=ifme
138xjfms=jfms
139xjfme=jfme
140
141
142! init flags
143need_lfn_update=.false.
144ignitions_done=0
145freeze_fire = fire_const_time > 0. .and. time_start < fire_const_time
146
147if(ifun.eq.1)then       ! do nothing, init pass 1 is outside only
148elseif(ifun.eq.2)then   
149        ! initialize all arrays that the model will not change later
150
151        ! assuming halo on zsf done
152        ! extrapolate on 1 row of cells beyond the domain boundary
153        ! including on the halo regions
154
155        call continue_at_boundary(1,1,0., & ! do x direction or y direction
156            ifms,ifme,jfms,jfme,           &                ! memory dims
157            ifds,ifde,jfds,jfde, &                     ! domain dims
158            ifps,ifpe,jfps,jfpe, &            ! patch dims - winds defined up to +1
159            ifts,ifte,jfts,jfte, &                ! tile dims
160            itso,iteo,jtso,jteo, &              ! where set now
161            fp%zsf)                               ! array
162
163!       compute the gradients once for all
164        err=0.
165        maxgrad=0.
166        do j=jfts,jfte
167            do i=ifts,ifte
168                erri = fp%dzdxf(i,j) - (fp%zsf(i+1,j)-fp%zsf(i-1,j))/(2.*fdx)
169                errj = fp%dzdyf(i,j) - (fp%zsf(i,j+1)-fp%zsf(i,j-1))/(2.*fdy)
170                err=max(err,abs(erri),abs(errj))
171                grad=sqrt(fp%dzdxf(i,j)**2+fp%dzdyf(i,j)**2)
172                maxgrad=max(maxgrad,grad)
173            enddo
174        enddo
175!$OMP CRITICAL(SFIRE_MODEL_CRIT)
176        write(msg,*)'max gradient ',maxgrad,' max error against zsf',err
177!$OMP END CRITICAL(SFIRE_MODEL_CRIT)
178        call message(msg)
179
180        if(.not.restart)call set_nfuel_cat( &
181            ifms,ifme,jfms,jfme, &
182            ifts,ifte,jfts,jfte, &
183            ifuelread,nfuel_cat0,&
184            fp%zsf,nfuel_cat)            ! better not use the extrapolated zsf!!
185
186        ! uses nfuel_cat to set the other fuel data arrays
187        ! needs zsf on halo width 1 to compute the terrain gradient
188        if(.not.restart)call set_fire_params(   &
189            ifds,ifde,jfds,jfde, &
190            ifms,ifme,jfms,jfme, &
191            ifts,ifte,jfts,jfte, &
192            fdx,fdy,nfuel_cat0,  &
193            nfuel_cat,fuel_time, &
194            fp  &
195)
196
197        ! initialize model state to no fire
198        if(.not.restart)then
199            call init_no_fire  ( &
200            ifds,ifde,jfds,jfde, &
201            ifms,ifme,jfms,jfme, &
202            ifts,ifte,jfts,jfte, &
203            fdx,fdy,time_start,  &
204            fuel_frac,fire_area,lfn,tign)
205           
206            need_lfn_update=.true. ! because we have set lfn
207
208        endif
209
210elseif(ifun.eq.3)then   ! ignition if so specified
211
212   
213elseif (ifun.eq.4) then  ! do the timestep
214
215    if(fire_print_msg.ge.stat_lev)then
216      aw=fun_real(RNRM_SUM,  &
217        ifms,ifme,1,1,jfms,jfme, &                ! memory dims
218        ifds,ifde,1,1,jfds,jfde, &                ! domain dims
219        ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
220        0,0,0,       &                            ! staggering
221        fp%vx,fp%vy)/((ifde-ifds+1)*(jfde-jfds+1))
222      mw=fun_real(RNRM_MAX,  &
223        ifms,ifme,1,1,jfms,jfme, &                ! memory dims
224        ifds,ifde,1,1,jfds,jfde, &                ! domain dims
225        ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
226        0,0,0,       &                            ! staggering
227        fp%vx,fp%vy)
228!$OMP MASTER
229      write(msg,91)time_start,'Average wind        ',aw,'m/s'
230      call message(msg,stat_lev)
231      write(msg,91)time_start,'Maximum wind        ',mw,'m/s'
232      call message(msg,stat_lev)
233!$OMP END MASTER
234    endif
235
236!   compute fuel fraction at start
237!    call fuel_left( &
238!        ifms,ifme,jfms,jfme, &
239!        ifts,ifte,jfts,jfte, &
240!        ifms,ifme,jfms,jfme, &
241!        lfn,tign,fuel_time,time_start,fuel_frac,fire_area) ! fuel frac is shared
242
243    call print_2d_stats(ifts,ifte,jfts,jfte, &
244                   ifms,ifme,jfms,jfme, &
245                   fuel_frac,'model: fuel_frac start')
246
247    ! advance the model from time_start to time_start+dt
248    ! return the fuel fraction burnt this call in each fire cell
249    ! will call module_fr_sfire_speed::normal_spread for propagation speed
250    ! We cannot simply compute the spread rate here because that will change with the
251    ! angle of the wind and the direction of propagation, thus it is done in subroutine
252    ! normal_spread at each fire time step. Instead, we pass arguments that
253    ! the speed function may use as fp.
254
255!   propagate level set function in time
256!   set lfn_out tign
257!   lfn does not change, tign has no halos
258
259    if(.not. freeze_fire)then
260
261    call prop_ls(id,     &
262        ifds,ifde,jfds,jfde,                      & ! fire domain dims - the whole domain
263        ifms,ifme,jfms,jfme,                      &
264        ifps,ifpe,jfps,jfpe, &                ! patch - nodes owned by this process
265        ifts,ifte,jfts,jfte,                      &
266        time_start,dt,fdx,fdy,tbound,  &
267        lfn,lfn_out,tign,ros, fp &
268    )
269    else
270        call message('sfire_model: EXPERIMENTAL: skipping fireline propagation')
271
272    endif
273   
274elseif (ifun.eq.5) then ! copy the result of timestep back to input
275    ! this cannot be done in the time step itself because of race condition
276    ! some thread may still be using lfn as input in their tile halo
277
278    if(.not. freeze_fire)then
279   
280    do j=jfts,jfte
281        do i=ifts,ifte
282            lfn(i,j)=lfn_out(i,j)
283            ! if want to try timestep again treat tign the same way here
284            ! even if tign does not need a halo
285        enddo
286    enddo
287
288    endif
289
290    ! check for ignitions
291    do ig = 1,num_ignitions
292   
293!  for now, check for ignition every time step...
294!        if(ignition_line(ig)%end_time>=time_start.and.ignition_line(ig)%start_time<time_start+dt)then
295            call ignite_fire(                             &
296                ifds,ifde,jfds,jfde,                      & ! fire domain dims - the whole domain
297                ifms,ifme,jfms,jfme,                      &
298                ifts,ifte,jfts,jfte,                      &
299                ignition_line(ig),                        &
300                time_start,time_start+dt,                 &
301                coord_xf,coord_yf,unit_xf,unit_yf,        &
302                lfn,tign,ignited)
303
304            ignitions_done=ignitions_done+1
305            ignited_tile(ignitions_done)=ignited
306               
307!            need_lfn_update=.true. ! if ignition, lfn changed
308#ifdef DEBUG_OUT   
309            call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,lfn,'lfn_ig',id)
310            call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,coord_xf,'coord_xf_ig',id)
311            call write_array_m(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme,coord_yf,'coord_yf_ig',id)
312#endif
313!        endif
314       
315    enddo
316           
317    call print_2d_stats(ifts,ifte,jfts,jfte,ifms,ifme,jfms,jfme, &
318                   lfn,'sfire_model: lfn out')
319
320   
321    need_lfn_update=.true. ! duh
322
323elseif (ifun.eq.6) then ! timestep postprocessing
324
325  if(.not. freeze_fire)then
326
327    ! compute the heat fluxes from the fuel burned
328    ! needs lfn and tign from neighbors so halo must be updated before
329    call fuel_left(&
330        ifms,ifme,jfms,jfme, &
331        ifts,ifte,jfts,jfte, &
332        ifts,ifte,jfts,jfte, &
333        lfn,tign,fuel_time,time_start+dt,fuel_frac_end,fire_area) !fuel_frac_end is private and tile based
334
335    call print_2d_stats(ifts,ifte,jfts,jfte, &
336                   ifts,ifte,jfts,jfte, &
337                   fuel_frac_end,'model: fuel_frac end')
338   
339    do j=jfts,jfte
340        do i=ifts,ifte
341            fuel_frac_burnt(i,j)=fuel_frac(i,j)-fuel_frac_end(i,j) ! fuel lost this timestep
342            fuel_frac(i,j)=fuel_frac_end(i,j) ! copy new value to state array
343        enddo
344    enddo
345
346    call print_2d_stats(ifts,ifte,jfts,jfte, &
347                   ifts,ifte,jfts,jfte, &
348                   fuel_frac_burnt,'model: fuel_frac burned')
349       
350    call heat_fluxes(dt,                          &
351        ifms,ifme,jfms,jfme,                      &
352        ifts,ifte,jfts,jfte,                      &
353        ifts,ifte,jfts,jfte,                      &  ! fuel_frac_burned is tile dimensioned
354        fp%fgip,                                     &
355        fuel_frac_burnt,                          & !
356        grnhfx,grnqfx)                              !out
357
358    if(fire_print_msg.ge.stat_lev)then
359      tfa=fun_real(REAL_SUM,  &
360        ifms,ifme,1,1,jfms,jfme, &                ! memory dims
361        ifds,ifde,1,1,jfds,jfde, &                ! domain dims
362        ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
363        0,0,0,       &                            ! staggering
364        fire_area,fire_area) * fdx * fdy
365      thf=fun_real(REAL_SUM,  &
366        ifms,ifme,1,1,jfms,jfme, &                ! memory dims
367        ifds,ifde,1,1,jfds,jfde, &                ! domain dims
368        ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
369        0,0,0,       &                            ! staggering
370        grnhfx,grnhfx) * fdx * fdy
371      mhf=fun_real(REAL_MAX,  &
372        ifms,ifme,1,1,jfms,jfme, &                ! memory dims
373        ifds,ifde,1,1,jfds,jfde, &                ! domain dims
374        ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
375        0,0,0,       &                            ! staggering
376        grnhfx,grnhfx)
377      tqf=fun_real(REAL_SUM,  &
378        ifms,ifme,1,1,jfms,jfme, &                ! memory dims
379        ifds,ifde,1,1,jfds,jfde, &                ! domain dims
380        ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
381        0,0,0,       &                            ! staggering
382        grnqfx,grnqfx) * fdx * fdy
383      mqf=fun_real(REAL_MAX,  &
384        ifms,ifme,1,1,jfms,jfme, &                ! memory dims
385        ifds,ifde,1,1,jfds,jfde, &                ! domain dims
386        ifts,ifte,1,1,jfts,jfte, &                ! patch or tile dims
387        0,0,0,       &                            ! staggering
388        grnqfx,grnqfx)
389!$OMP MASTER
390      write(msg,91)time_start,'Fire area           ',tfa,'m^2'
391      call message(msg,stat_lev)
392      write(msg,91)time_start,'Heat output         ',thf,'W'
393      call message(msg,stat_lev)
394      write(msg,91)time_start,'Max heat flux       ',mhf,'W/m^2'
395      call message(msg,stat_lev)
396      write(msg,91)time_start,'Latent heat output  ',tqf,'W'
397      call message(msg,stat_lev)
398      write(msg,91)time_start,'Max latent heat flux',mqf,'W/m^2'
399      call message(msg,stat_lev)
400!$OMP END MASTER
40191  format('Time ',f11.3,' s ',a,e12.3,1x,a)
402    endif
403       
404
405  else
406     call message('sfire_model: EXPERIMENTAL: skipping fuel burnt computation')
407
408     if (fire_const_grnhfx >= 0. .and. fire_const_grnqfx >= 0.) then
409
410!$OMP CRITICAL(SFIRE_MODEL_CRIT)
411        write(msg,'(a,2e12.3,a)')'sfire_model: EXPERIMENTAL output constant heat flux', &
412           fire_const_grnhfx, fire_const_grnqfx, ' W/s'
413!$OMP END CRITICAL(SFIRE_MODEL_CRIT)
414        call message(msg)
415       
416        do j=jfts,jfte
417            do i=ifts,ifte
418                grnhfx(i,j)=fire_const_grnhfx
419                grnqfx(i,j)=fire_const_grnqfx
420            enddo
421        enddo
422
423      endif
424
425   endif
426
427    call print_2d_stats(ifts,ifte,jfts,jfte, &
428                   ifms,ifme,jfms,jfme, &
429                   grnhfx,'model: heat flux(J/m^2/s)')
430
431else
432!$OMP CRITICAL(SFIRE_MODEL_CRIT)
433    write(msg,*)'sfire_model: bad ifun=',ifun
434!$OMP END CRITICAL(SFIRE_MODEL_CRIT)
435    call crash(msg)
436endif
437
438end subroutine sfire_model
439
440!
441!*****************
442!
443           
444end module module_fr_sfire_model
Note: See TracBrowser for help on using the repository browser.