| blob | aab8d25a1d1a8db7be739b6a37bcf096516fc949 |
1 !-------------------------------------------------------------------------------
2 ! This module was developed in RAL-NCAR/UCAR, 2019-2020 by:
3 ! *** M. Frediani and T. Juliano ***
4 ! This research was supervised by J. Knievel and B. Kosovic.
5 !-------------------------------------------------------------------------------
7 !=============================================================
8 !=============================================================
10 MODULE module_firebrand_spotting
12 USE, INTRINSIC :: IEEE_ARITHMETIC
13 USE module_domain, ONLY : get_ijk_from_grid, domain ! grid
14 USE module_configure, ONLY : grid_config_rec_type ! config_flags
15 USE module_symbols_util, ONLY : WRFU_TimeInterval, WRFU_TimeIntervalGet, WRFU_TimeIntervalSet
16 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
17 USE MPI
18 #endif
20 IMPLICIT NONE
22 PRIVATE
23 PUBLIC firebrand_spotting_em_init, firebrand_spotting_em_driver, get_local_ijk
25 ! THESE VARIABLES ARE IN MODULE SCOPE ! Careful with reassignments - don't reassign
26 ! SAVE attribute is default
28 !-------------------------------------------------------------------------------
29 ! variables in module scope: private, only available module-wide (host association)
30 !-------------------------------------------------------------------------------
31 ! They should not be declared again in suboutines (may not compile)
32 ! and must not be among routine dummy arguments. Consequently, cannot be IO variables
33 !
34 ! Runtime variables are not available at module level (e.g., namelist, tile dependent variables).
35 ! Include here only what can be set during compilation:
36 ! fixed parameters, allocatables, declarions (without initialization)
38 !-------------------------------------------------------------------------------
39 ! Fixed parameters ***MODULE SCOPE***
40 !-------------------------------------------------------------------------------
42 INTEGER, PARAMETER :: dp = KIND(0.d0) ! double precision
44 INTEGER, PARAMETER :: wrfdbg = 10
45 ! Fuel parameters must match values from module_fr_phys.F: SUBROUTINE set_fuel_params(nfuel_cat0, nfuel_cat)
46 INTEGER, PARAMETER :: no_fuel_cat = 14
47 INTEGER, PARAMETER :: nf_sb = 204 ! maximum category on
48 INTEGER, PARAMETER :: nfuelcats = 53 ! number of fuel categories that are specified
50 REAL, PARAMETER :: grav = 9.80616_dp ! gravity (m/s2)
51 REAL, PARAMETER :: rdry = 287.04_dp ! dry air (J/Kg-K)
52 REAL, PARAMETER :: p2jm = 100.0_dp ! mb to j/m3
53 REAL, PARAMETER :: rcd = 0.45_dp ! drag constant
54 REAL, PARAMETER :: rd = 287.15_dp ! gas constant dry air (J / kg K)
55 REAL, PARAMETER :: rv = 461.6_dp
56 REAL, PARAMETER :: t0 = 300.0_dp
57 REAL, PARAMETER :: cp = 7.0_dp*rd/2.0_dp
58 REAL, PARAMETER :: rcp = rd/cp
59 REAL, PARAMETER :: p1000mb= 100000.0_dp
60 REAL, PARAMETER :: r1o3 = 1.0_dp/3.0_dp ! ratio 1 over 3
61 REAL, PARAMETER :: sboltz = 5.67E-5_dp ! stefan-boltzmann (g / s3-K4)
62 REAL, PARAMETER :: emiss = 0.9_dp ! emissivity
63 REAL, PARAMETER :: cpw = 1466.0_dp ! specific heat wood (J / kg-K)
64 REAL, PARAMETER :: cpc = 712.0_dp ! specific heat char (J / kg-K)
65 REAL, PARAMETER :: beta0 = 4.8E-7_dp ! burning rate constant (m2 / s)
66 REAL, PARAMETER :: s_coeff= 110.4_dp ! Sutherland's law coefficient (K)
67 REAL, PARAMETER :: b_coeff= 1.458E-3_dp ! Sutherland's law coefficient [g / m-s-K-(1/2)]
68 REAL, PARAMETER :: shmt = 0.7_dp ! schmidt number
69 REAL, PARAMETER :: thcon = 27.0_dp ! thermal conductivity air
71 !USE module_state_description, ONLY: p_qv
72 !INTEGER, PARAMETER :: p_qv = 1
73 REAL, PARAMETER :: pr = 0.7_dp ! Prandtl
75 REAL, PARAMETER :: NEGLIGIBLE = 10*EPSILON(1.0)
76 REAL, PARAMETER :: ZERO_dp = 0.0_dp ! this is a real type variable, not a double precision type
77 REAL, PARAMETER :: dp05 = 0.5_dp
78 REAL, PARAMETER :: dp1 = 1.0_dp
81 !-------------------------------------------------------------------------------
82 ! ***MODULE SCOPE***
83 !-------------------------------------------------------------------------------
85 ! Mass threshold to consider burnout (g)
86 REAL, PARAMETER :: br_min_mass = EPSILON(dp1) ! max precision for real(4)
88 ! Diameter threshold to consider burnout (mm)
89 REAL, PARAMETER :: br_min_diam = 0.0000001_dp
91 !-------------------------------------------------------------------------------
92 ! Generic variables for multiple use within module ***MODULE SCOPE***
93 !-------------------------------------------------------------------------------
95 CHARACTER (LEN=200), SAVE :: msg
96 CHARACTER (LEN=256), SAVE :: fmt
97 CHARACTER (LEN=200), DIMENSION(10) :: amsg
98 INTEGER, SAVE :: imsg ! loop counters
100 !-------------------------------------------------------------------------------
101 ! variables from namelist ***MODULE SCOPE***
102 !-------------------------------------------------------------------------------
104 ! These will be available to all subroutines and will be set in init routine
105 INTEGER, SAVE :: fs_array_maxsize ! maximum number of particles carried in simulation
106 INTEGER, SAVE :: fs_gen_levels ! one per grid pt (may change if releasing at multiple levels)
107 INTEGER, SAVE :: fs_gen_lim
108 INTEGER, SAVE :: fs_gen_dt
110 REAL, SAVE :: firebrand_dens ! 513000.0_dp ! density of firebrand (g / m3)
111 REAL, SAVE :: firebrand_dens_char ! 299000.0_dp ! density of char (g m-3) [gupta et al 2002; fuel]
114 !-------------------------------------------------------------------------------
115 ! Fixed indices, ranks ***MODULE SCOPE***
116 !-------------------------------------------------------------------------------
118 LOGICAL :: this_is_ideal = .FALSE.
120 TYPE p_properties ! fb_prop
121 REAL :: p_mass ! fbmass
122 REAL :: p_diam ! fbdiam
123 REAL :: p_effd ! fbediam
124 REAL :: p_temp ! fbtemp
125 REAL :: p_tvel ! fbvelo
126 END TYPE p_properties
128 !-------------------------------------------------------------------------------
129 ! grid and cf are not here because dimensions are given at runtime (derived types)
130 ! Also, grid values change with timestep, so they need to be passed as dummy arguments
131 !-------------------------------------------------------------------------------
133 !-------------------------------------------------------------------------------
134 ! Variable bounds - Initialized in init, used in dummy arguments in driver
135 ! ***MODULE SCOPE***
136 !-------------------------------------------------------------------------------
137 INTEGER, SAVE :: ids, jds, ide, jde, kde ! domain bounds
138 INTEGER, SAVE :: ims, jms, ime, jme, kms, kme ! memory bounds
139 INTEGER, SAVE :: is, ie, js, je, ks, ke ! patch start/end
140 INTEGER, SAVE :: ifps, jfps, ifpe, jfpe ! refined fire grid bounds
142 CONTAINS
144 !******************************************************************
145 !******************************************************************
146 !* *
147 !* Module Routines *
148 !* *
149 !******************************************************************
150 !******************************************************************
154 !=============================================================
155 SUBROUTINE firebrand_spotting_em_init ( &
156 !=============================================================
157 grid, &
158 cf, &
159 fs_p_id, &
160 fs_p_src, &
161 fs_p_dt, &
162 fs_p_x, &
163 fs_p_y, &
164 fs_p_z, &
165 fs_last_gen_dt,&
166 fs_gen_idmax, &
167 fs_count_reset,&
168 fs_p_mass, &
169 fs_p_diam, &
170 fs_p_effd,&
171 fs_p_temp, &
172 fs_p_tvel, &
173 fs_count_landed_all,&
174 fs_count_landed_hist,&
175 fs_landing_mask,&
176 fs_spotting_lkhd,&
177 fs_frac_landed)
179 !=============================================================
180 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
181 USE module_firebrand_spotting_mpi, ONLY: fs_mpi_init
182 #endif
183 !-------------------------------------------------------------------------------
184 ! Initialize all I/O and variables declared in the registry
185 !-------------------------------------------------------------------------------
187 IMPLICIT NONE
189 !-------------------------------------------------------------------------------
190 ! Arguments
191 !-------------------------------------------------------------------------------
193 TYPE(domain), INTENT(IN) :: grid !
194 TYPE(grid_config_rec_type), INTENT(IN) :: cf ! run-time configuration (namelist) for domain
196 INTEGER, INTENT(INOUT) :: fs_last_gen_dt, fs_gen_idmax
197 LOGICAL, INTENT(INOUT) :: fs_count_reset
199 ! Firebrand Spotting Particle Properties (fs_p_*)
200 INTEGER, INTENT(INOUT), DIMENSION(:) :: fs_p_id, fs_p_dt, fs_p_src
201 REAL, INTENT(INOUT), DIMENSION(:) :: fs_p_x, fs_p_y, fs_p_z ! positions are relative to grid edges: particle at grid center point (x,y,z) = (1.5, 1.5, 1.5)
202 REAL, INTENT(INOUT), DIMENSION(:) :: fs_p_mass, fs_p_diam, fs_p_effd
203 REAL, INTENT(INOUT), DIMENSION(:) :: fs_p_temp, fs_p_tvel
204 REAL, INTENT(INOUT), DIMENSION(ims:,jms:) :: fs_count_landed_all, fs_count_landed_hist, fs_spotting_lkhd, fs_frac_landed
205 INTEGER, INTENT(INOUT), DIMENSION(ims:,jms:) :: fs_landing_mask
207 !-------------------------------------------------------------------------------
209 !-------------------------------------------------------------------------------
210 ! Module variables from namelist/default
212 fs_array_maxsize = cf%fs_array_maxsize
213 fs_gen_levels = cf%fs_firebrand_gen_levels
214 fs_gen_lim = cf%fs_firebrand_gen_lim
215 fs_gen_dt = cf%fs_firebrand_gen_dt
217 fmt = '(A,1x,I6)'
218 WRITE (amsg(1),*) 'SPFire_init: fspotting_em_init'
219 WRITE (amsg(2),fmt) 'SPFire_init: firebrand limit =', fs_gen_lim
220 WRITE (amsg(3),fmt) 'SPFire_init: dt =', fs_gen_dt
221 WRITE (amsg(4),fmt) 'SPFire_init: fs_array_maxsize =', fs_array_maxsize
222 WRITE (amsg(5),fmt) 'SPFire_init: fs_gen_levels =', fs_gen_levels
223 DO imsg=1,5
224 CALL wrf_debug (wrfdbg, TRIM(amsg(imsg)) )
225 ENDDO
227 !-------------------------------------------------------------------------------
228 ! Set bounds to be used as dummy arguments in driver
229 ! ***variables declared in MODULE SCOPE***
230 !-------------------------------------------------------------------------------
232 CALL get_local_ijk(grid, &
233 ifps=ifps, jfps=jfps, ifpe=ifpe, jfpe=jfpe, &
234 ids=ids, jds=jds, ide=ide, jde=jde, kde=kde, &
235 ims=ims, jms=jms, ime=ime, jme=jme, kms=kms, kme=kme, &
236 ips=is, jps=js, ipe=ie, jpe=je, kps=ks, kpe=ke)
238 WRITE (msg,'(6(i6,1x))') is, ie, js, je, ks, ke
239 CALL wrf_debug (wrfdbg, 'SPFire_init tile bounds: '//msg)
241 WRITE (msg,'(6(i6,1x))') ims, ime, jms, jme, kms, kme
242 CALL wrf_debug (wrfdbg, 'SPFire_init memory bounds: '//msg)
244 WRITE (msg,'(4(i6,1x))') ifps, ifpe, jfps, jfpe
245 CALL wrf_debug (wrfdbg, 'SPFire_init fire refined bounds: '//msg)
248 !-------------------------------------------------------------------------------
249 ! Initialize registry arrays
250 !-------------------------------------------------------------------------------
251 fs_count_reset=.FALSE.
253 fs_last_gen_dt = 0
254 fs_gen_idmax = 0
255 fs_p_id = 0
256 fs_p_dt = 0
257 fs_p_src = 0
258 fs_p_x = 0.0_dp
259 fs_p_y = 0.0_dp
260 fs_p_z = 0.0_dp
261 fs_p_mass = 0.0_dp
262 fs_p_diam = 0.0_dp
263 fs_p_effd = 0.0_dp
264 fs_p_temp = 0.0_dp
265 fs_p_tvel = 0.0_dp
267 fs_landing_mask(:,:) = 0
268 fs_count_landed_all (:,:) = 0.0_dp
269 fs_count_landed_hist(:,:) = 0.0_dp
270 fs_spotting_lkhd(:,:) = 0.0_dp
271 fs_frac_landed(:,:) = 0.0_dp
273 !fs_fire_ROSdt(:,:) = 0.0_dp
274 !fs_gen_inst(:,:) = 0
276 !-------------------------------------------------------------------------------
278 IF ( grid%this_is_an_ideal_run ) THEN
280 this_is_ideal = .TRUE.
281 CALL wrf_debug (wrfdbg, 'SPFire_init: Ideal Run detected' )
283 ELSE
285 this_is_ideal = .FALSE.
286 CALL wrf_debug (wrfdbg, 'SPFire_init: Not an Ideal Run' )
288 ENDIF
291 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
292 CALL fs_mpi_init(grid)
293 #endif
295 END SUBROUTINE firebrand_spotting_em_init
296 !=============================================================
297 !=============================================================
300 !=============================================================
301 PURE &
302 FUNCTION order_val(arr, ord)
303 !=============================================================
305 IMPLICIT NONE
307 REAL, INTENT(IN), DIMENSION(:) :: arr
308 INTEGER, INTENT(IN) :: ord
309 REAL :: order_val
310 INTEGER :: diff, besti, iord, i, error
312 !-----------------------------------------------------------
313 ! Return the value in the array at the given order (allowing an error of +-2 positions)
314 ! E.g.:
315 ! For an arr with dimension = 100
316 ! ord = 50 will return one of the values near the 50% percentile
317 ! i.e., the value in one of the positions between 48:52 if the array was ranked
318 !
319 ! The algorithm scans the array and
320 ! counts the number of elements with value below the scanned position
321 ! It returns when it finds a value where
322 ! (order - error) <= count <= (order + error)
323 !-----------------------------------------------------------
325 error = 2
327 diff = SIZE(arr) -1
328 besti = 1
330 DO i=1,SIZE(arr)
332 iord = COUNT(arr > arr(i))
333 IF (ABS(ord - iord) <= diff) besti = i
334 diff = MIN(diff, ABS(ord - iord))
335 IF (diff <= error) EXIT
336 ENDDO
338 ! WRITE(msg, '(2(i6,1x))') iord, diff
339 ! CALL wrf_debug (wrfdbg, 'SPFire_rel_order_val: '//msg)
341 order_val = arr(besti)
343 END FUNCTION order_val
344 !=============================================================
345 !=============================================================
349 !=============================================================
350 PURE &
351 FUNCTION fire2tile(fr_arr, dsx,dsy) RESULT(new_arr)
352 !=============================================================
354 IMPLICIT NONE
356 REAL, INTENT(IN), DIMENSION(ifps:ifpe,jfps:jfpe) :: fr_arr
357 INTEGER, INTENT(IN) :: dsx, dsy
358 REAL, ALLOCATABLE,DIMENSION(:,:) :: new_arr
359 INTEGER :: i, j
360 INTEGER, DIMENSION(2) :: fshp
362 !-----------------------------------------------------------
363 ! Converts a 2-D array from the fire refined grid to tile
364 ! Tile values are the sum of corresponding refined gridpoints
365 !-----------------------------------------------------------
367 fshp = SHAPE(fr_arr)
369 ALLOCATE(new_arr(1:fshp(1)/dsx, 1:fshp(2)/dsy))
370 new_arr(:,:) = ZERO_dp
372 new_arr(1:fshp(1)/dsx, 1:fshp(2)/dsy) = &
373 RESHAPE([((&
374 SUM(fr_arr(i:i+dsx-1,j:j+dsy-1)), &
375 i=ifps,ifpe-1,dsx), &
376 j=jfps,jfpe-1,dsy)], &
377 [fshp(1)/dsx,fshp(2)/dsy])
379 END FUNCTION fire2tile
380 !=============================================================
381 !=============================================================
385 !=============================================================
386 ELEMENTAL FUNCTION fuel_spotting_risk(fuel_fgi, factor, fuel_mcg)!, fuel_mcg)
387 !=============================================================
389 ! Returns the refined grid array with the fuel_spotting_risk for spotfires
390 !-------------------------------------------------------------------------------
392 ! fuel_name(k) : fuel model name
393 ! windrf(k) : wind reduction factor from 20ft to midflame height
394 ! fgi(k) : initial total mass of surface fuel (kg/m**2)
395 ! fgi_lh(k) : initial total mass of surface live fuel [sb: 1-h] (kg/m**2)
396 ! fgi_all(k) : initial total mass of surface live fuel [sb: 1-h,10-h,100-h] (kg/m**2)
397 ! fueldepthm(k) : fuel depth (m)
398 ! savr(k) : fuel particle surface-area-to-volume ratio, 1/ft
399 ! fuelmce(k) : moisture content of extinction
400 ! fueldens(k) : ovendry particle density, lb/ft^3
401 ! st(k) : fuel particle total mineral content
402 ! se(k) : fuel particle effective mineral content
403 ! weight(k) : weighting parameter that determines the slope of the mass loss curve
404 ! fci_d(k) : initial dry mass of canopy fuel
405 ! fct(k) : burn out time for canopy fuel, after dry (s)
406 ! ichap(k) : 1 if chaparral, 0 if not
407 ! fci(k) : initial total mass of canopy fuel
408 ! fcbr(k) : fuel canopy burn rate (kg/m**2/s)
409 ! hfgl : surface fire heat flux threshold to ignite canopy (w/m^2)
410 ! cmbcnst : joules per kg of dry fuel
411 ! fuelheat : fuel particle low heat content, btu/lb
412 ! fuelmc_g : fuel particle (surface) moisture content
413 ! fuelmc_c : fuel particle (canopy) moisture content
415 IMPLICIT NONE
417 REAL :: fuel_spotting_risk
418 REAL, INTENT(IN) :: factor
419 REAL, INTENT(IN) :: fuel_fgi, fuel_mcg ! fmcg2df ground dead fuel moisture content (fire grid)
421 fuel_spotting_risk = factor * (dp1 - MIN(dp1, fuel_mcg/fuel_fgi)) ! *(fuel_fgi-fuel_mcg)
423 END FUNCTION fuel_spotting_risk
424 !=============================================================
425 !=============================================================
429 !=============================================================
430 ELEMENTAL FUNCTION spotting_threat_factor(fcat)
431 !=============================================================
432 IMPLICIT NONE
434 INTEGER, INTENT(IN) :: fcat
435 REAL, DIMENSION(nfuelcats+1) :: factor
436 REAL :: spotting_threat_factor
438 ! **************** read this from external input ****************
439 ! **************** call it in firebrand_spotting_driver once, when grid%itimestep =1 ****************
440 ! **************** and load a registry array ****************
442 ! User defined coefficient to increase spotfire likelihood;
443 ! threat = 0 yields zero likelihood
445 factor(1) = 1.0_dp ! 1: 'Short grass (1 ft)',
446 factor(2:13) = 1.0_dp
447 factor(14) = ZERO_dp ! 14: 'no fuel',
448 factor(15:nfuelcats+1) = 1.0_dp
450 spotting_threat_factor = factor(fcat)
452 !CALL read_table_static
454 END FUNCTION spotting_threat_factor
456 !=============================================================
457 !=============================================================
460 ! !=============================================================
461 ! SUBROUTINE read_table_static ! Not debuggeed - from module_mp_thompson
462 ! !=============================================================
464 ! USE module_domain
465 ! USE module_dm
466 ! IMPLICIT NONE
468 ! LOGICAL, EXTERNAL:: wrf_dm_on_monitor
469 ! INTEGER:: iunit_fs1, i
470 ! LOGICAL:: opened
471 ! CHARACTER (LEN=64) :: filename
473 ! iunit_fs1 = -1
474 ! IF ( wrf_dm_on_monitor() ) THEN
475 ! DO i = 20,99
476 ! INQUIRE ( i , OPENED = opened )
477 ! IF ( .NOT. opened ) THEN
478 ! iunit_fs1 = i
479 ! EXIT
480 ! ENDIF
481 ! ENDDO
482 ! ENDIF
484 ! #if defined(DM_PARALLEL) && !defined(STUBMPI)
485 ! CALL wrf_dm_bcast_bytes ( iunit_fs1 , IWORDSIZE )
486 ! #endif
487 ! IF ( iunit_fs1 < 0 ) &
488 ! CALL wrf_error_fatal('SPFire_read_table_static: Cannot find a fortran unit to read file')
490 ! IF ( wrf_dm_on_monitor() ) THEN
492 ! WRITE(msg, '(i2,1x,a20)') iunit_fs1, filename
493 ! CALL wrf_debug(wrfdbg,'SPFire_read_table_static: opening table on unit ' // msg)
495 ! OPEN(iunit_fs1,FILE=tablename, FORM='UNFORMATTED',STATUS='OLD',ERR=9009)
497 ! ENDIF
499 ! #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes(A, size(A)*R4SIZE)
501 ! IF ( wrf_dm_on_monitor() ) READ(iunit_fs1,ERR=9010) tableval
502 ! #if defined(DM_PARALLEL) && !defined(STUBMPI)
503 ! DM_BCAST_MACRO(tableval)
504 ! #endif
506 ! RETURN
508 ! 9009 CONTINUE
509 ! CALL wrf_error_fatal('SPFire_read_table_static: error openinig static table')
510 ! RETURN
512 ! 9010 CONTINUE
513 ! CALL wrf_error_fatal('SPFire_read_table_static: error reading static table')
516 ! END SUBROUTINE read_table_static
518 ! !=============================================================
519 ! !=============================================================
522 !=============================================================
523 ELEMENTAL FUNCTION get_fuel_cat(crosswalk, cat)
524 !=============================================================
526 ! Use with Scott & Burgen fuel model (WRF CO-FSP)
528 IMPLICIT NONE
530 INTEGER :: get_fuel_cat
531 LOGICAL, INTENT(IN):: crosswalk
532 REAL, INTENT(IN):: cat
534 INTEGER, DIMENSION(1:nf_sb) :: ksb ! Anderson82 + S&B2005 fuel categories array
535 INTEGER :: i, icat
536 !-------------------------------------------------------------------------------
538 icat = INT(cat)
540 ksb = no_fuel_cat ! 14 ! Initialize all values to no-fuel
542 ! Anderson 1982
543 ksb(1:13) = [(i, i=1, 13)]
545 IF ( crosswalk ) THEN ! Scott & Burgan crosswalks to Anderson
547 ksb([101, 104, 107])=1 ! Short grass -- 1
548 ksb(102)=2 ! Timber grass and understory -- 2
549 ksb(121:124) = 2 ! Timber grass and understory -- 2
550 ksb([103,105,106,108,109])=3 ! Tall grass -- 3
551 ksb([145, 147])=4 ! Chaparral -- 4
552 ksb(142)=5 ! Brush -- 5
553 ksb([141, 146])=6 ! Dormant Brushi -- 6
554 ksb([143, 144, 148, 149])=7 ! Southern Rough -- 7
555 ksb([181, 183, 184, 187])=8 ! Compact Timber Litter -- 8
556 ksb([182, 186, 188, 189])=9 ! Hardwood Litter -- 9
557 ksb(161: 165)=10 ! Timber (understory) -- 10
558 ksb([185, 201])=11 ! Light Logging Slash -- 11
559 ksb(202)=12 ! Medium Logging Slash -- 12
560 ksb([203, 204])=13 ! Heavy Logging Slash -- 13
562 ELSE ! full Scott and Burgan (2005)
564 ksb(101:109) = [(i, i=15,23)] ! Grass (GR)
565 ksb(121:124) = [(i, i=24,27)] ! Grass-Shrub (GS)
566 ksb(141:149) = [(i, i=28,36)] ! Shrub (SH)
567 ksb(161:165) = [(i, i=37,41)] ! Timber-Understory (TU)
568 ksb(181:189) = [(i, i=42,50)] ! Timber litter (TL)
569 ksb(201:204) = [(i, i=51,54)] ! Slash-Blowdown (SB)
571 ENDIF
573 get_fuel_cat = ksb(icat)
575 END FUNCTION get_fuel_cat
576 !=============================================================
577 !=============================================================
580 !=============================================================
581 ELEMENTAL FUNCTION firebrand_gen_factor(fcat)
582 !=============================================================
583 IMPLICIT NONE
585 REAL :: firebrand_gen_factor
586 INTEGER, INTENT(IN) :: fcat
587 REAL, DIMENSION(nfuelcats+1) :: factor
590 ! **************** read this from external input ****************
591 ! **************** call it in firebrand_spotting_driver once, when grid%itimestep =1 ****************
592 ! **************** and load a registry array ****************
594 ! User defined coefficient to set release priority
595 ! factor = 0 yields zero releases from gridpoint
597 factor(1:nfuelcats+1)= 1.0_dp
598 firebrand_gen_factor = factor(fcat)
600 END FUNCTION firebrand_gen_factor
602 !=============================================================
603 !=============================================================
607 ! !=============================================================
608 ! ELEMENTAL FUNCTION firebrand_gen_maxhgt(fcat)
609 ! !=============================================================
610 ! IMPLICIT NONE
612 ! REAL :: firebrand_gen_maxhgt
613 ! INTEGER, INTENT(IN) :: fcat
614 ! REAL, DIMENSION(nfuelcats+1) :: maxhgt
617 ! ! **************** read this from external input ****************
618 ! ! **************** call it in firebrand_spotting_driver once, when grid%itimestep =1 ****************
619 ! ! **************** and load a registry array ****************
621 ! ! User defined max height to release firebrand
623 ! maxhgt(1:nfuelcats+1)= 50.0_dp
624 ! firebrand_gen_maxhgt = maxhgt(fcat)
626 ! END FUNCTION firebrand_gen_maxhgt
628 !=============================================================
629 !=============================================================
633 !=============================================================
634 ELEMENTAL FUNCTION firebrand_gen_potential(fuel_fgi, factor, fire_rate, fuel_mcg)!, fuel_mcg)
635 !=============================================================
637 ! Returns the refined grid array with the potential for firebrand generation
638 !-------------------------------------------------------------------------------
640 IMPLICIT NONE
642 REAL :: firebrand_gen_potential
643 REAL, INTENT(IN) :: factor, fire_rate
644 REAL, INTENT(IN) :: fuel_fgi, fuel_mcg ! fmcg2df ground dead fuel moisture content (fire grid)
646 firebrand_gen_potential = factor * fire_rate * fuel_fgi !*(dp1 - fuel_mcg/(dp1+fuel_mcg)) ! fuelloadm; dry load from fr_fire_phys
648 END FUNCTION firebrand_gen_potential
649 !=============================================================
650 !=============================================================
654 !=============================================================
655 ! PURE & ! can't call wrf_fatal in pure routines
656 SUBROUTINE generate_firebrands( fs_p_id, fs_p_src, fs_p_dt, fs_p_z, fs_p_x, fs_p_y, &
657 release_i, release_j, release_k, &
658 release_dt, release_src, &
659 active_br, fs_gen_idmax, myprocid, &
660 release_prop, fs_p_prop)
661 !=============================================================
663 IMPLICIT NONE
665 TYPE(p_properties), INTENT(IN), DIMENSION(:):: release_prop
666 TYPE(p_properties), INTENT(OUT),DIMENSION(:):: fs_p_prop
668 INTEGER, INTENT(INOUT), DIMENSION(:) :: fs_p_id, fs_p_src, fs_p_dt
669 REAL, INTENT(INOUT), DIMENSION(:) :: fs_p_z, fs_p_x, fs_p_y
670 REAL, INTENT(INOUT), DIMENSION(:) :: release_i, release_j, release_k
671 INTEGER, OPTIONAL, INTENT(IN), DIMENSION(:) :: release_dt, release_src
672 INTEGER, INTENT(INOUT) :: fs_gen_idmax
673 INTEGER, INTENT(IN) :: myprocid
674 INTEGER, INTENT(OUT) :: active_br
676 LOGICAL :: release_true
677 INTEGER :: br, ii
678 REAL :: rx, ry, rz
679 INTEGER :: new_release
680 LOGICAL, DIMENSION(fs_array_maxsize) :: flag_true
682 !-------------------------------------------------------------------------------
683 release_true = .FALSE. ! scalar: true when brands are ready to be released
684 flag_true = .FALSE. ! logical array: flag to find active elements in array
686 active_br = 0 ! number of active brands
687 new_release = 0 ! number of brands to release
689 ! when all brand IDs of array elements are zero: nothing active
690 active_br = COUNT( fs_p_id > 0 )
691 new_release = COUNT( (INT(release_i) > 0) .AND. (INT(release_j) > 0) )
693 IF (new_release > 0 ) release_true = .TRUE. ! to-do: include dependency on release_dt
695 IF ( .NOT. release_true) THEN
696 RETURN !CALL wrf_debug (wrfdbg, 'SPFire_generate_firebrands: Found nothing to release')
697 ENDIF
699 !-------------------------------------------------------------------------------
700 ! New brands fit in array? Release particles: set lat/lon/lev of new brands
701 !-------------------------------------------------------------------------------
703 IF (active_br + new_release <= fs_array_maxsize) THEN
705 !-------------------------------------------------------------------------------
706 ! Add new grid positions to array: create a subroutine for this: brand_release
707 !-------------------------------------------------------------------------------
709 ! Flag is True where positions are available to be filled
710 WHERE(fs_p_id == 0) flag_true = .TRUE. ! id_indx holds brand ID (unique integer)
712 ii = active_br + 1
713 DO br=1, new_release
715 ! Find an ii position available to store a new release
716 IF ( .NOT. flag_true(ii) ) THEN ! WRITE (msg,'(3(i6,1x))') ii, active_br, fs_p_id(ii) ! CALL wrf_debug (wrfdbg, 'SPFire_release flag false >>> '// msg)
717 CALL wrf_error_fatal('SPFire_generate_firebrands: Did not find free index to release brands! Did you pack fs_p_x, fs_p_y, fs_p_id, etc?')
718 ENDIF
720 ! Sanity check
721 IF (INT(release_i(br)) == 0 .OR. &
722 INT(release_j(br)) == 0 ) &
723 CALL wrf_error_fatal('SPFire_generate_firebrands: release_ijk is zero! Positions cannot be zero.')
725 IF (fs_gen_idmax + 10 >= HUGE(1)) fs_gen_idmax = 0
727 ! Assign values and convert height to z index
728 fs_p_x(ii) = release_i(br)
729 fs_p_y(ii) = release_j(br)
730 fs_p_z(ii) = release_k(br)
731 fs_p_id(ii) = fs_gen_idmax + 1 ! Assign new IDmax to particle
733 IF (.NOT. PRESENT(release_dt)) fs_p_dt(ii) = 0 ! firebrand has not been advected yet so dt is zero
734 IF (.NOT. PRESENT(release_src)) fs_p_src(ii) = myprocid + fs_p_id(ii)
736 IF (PRESENT(release_dt)) fs_p_dt(ii) = release_dt(br)
737 IF (PRESENT(release_src)) fs_p_src(ii) = release_src(br)
739 fs_p_prop(ii)%p_mass = release_prop(br)%p_mass
740 fs_p_prop(ii)%p_diam = release_prop(br)%p_diam
741 fs_p_prop(ii)%p_effd = release_prop(br)%p_effd
742 fs_p_prop(ii)%p_temp = release_prop(br)%p_temp
743 fs_p_prop(ii)%p_tvel = release_prop(br)%p_tvel
745 fs_gen_idmax = fs_p_id(ii)
747 flag_true(ii) = .FALSE.
748 ii = ii + 1
750 ENDDO
752 active_br = active_br + new_release
753 release_i = ZERO_dp
754 release_j = ZERO_dp
755 release_k = ZERO_dp
757 ! WRITE (msg,'(2(i6,1x))') new_release, fs_gen_idmax
758 ! CALL wrf_debug (wrfdbg, 'SPFire_release: Released new brands fs_gen_idmax: '//msg)
760 ! ELSE
762 ! WRITE (msg,*) active_br + new_release, fs_array_maxsize
763 ! CALL wrf_debug (wrfdbg, 'SPFire_release: brand array is full, cannot release new brands'// msg)
765 ENDIF
768 END SUBROUTINE generate_firebrands
770 !=============================================================
771 !=============================================================
776 !=============================================================
777 PURE &
778 FUNCTION hgt2k(xp, yp, hgt, z_at_w, znw)
779 !=============================================================
781 !-------------------------------------------------------------------------------
782 ! Converts height in meters to vertical level
783 !-------------------------------------------------------------------------------
785 IMPLICIT NONE
786 !INTEGER, INTENT(IN) :: ims, jms
787 REAL, INTENT(IN) :: hgt, xp, yp
788 REAL, INTENT(IN), DIMENSION(ims:,:,jms:) :: z_at_w
789 REAL, INTENT(IN), DIMENSION(:) :: znw
791 REAL :: zdz, z_lower_at_p, z_upper_at_p, x0,y0
792 REAL :: hgt2k
793 INTEGER :: k, i, j
795 ! An (xp,yp)=(1.5, 1.5) places the particle at the gridpoint center
796 ! but z_at_w is horizontally unstaggered,
797 ! so z_at_w(i,j) is at the horizontal grid center when (i,j)=(1,1)
798 !
799 ! Hence, we adjust the particle position to fetch the horizontally adjacent heights for the linear interp.
800 ! We subtract the deltax between
801 ! position at grid center and position at grid edge from the particle position xp, yp, such that
802 ! z_at_w(i,1,j) for (i,j)=(1,1) corresponds to the height at (xp,yp)=(1.5,1.5):
803 ! (i = xp - 0.5, j = yp - 0.5) = (i=1, j=1), for (xp=1.5, yp=1.5)
805 x0 = xp - dp05
806 y0 = yp - dp05
808 i = FLOOR(x0) ! shift indices by 0.5 because z_at_w is given at horizontal grid center
809 j = FLOOR(y0)
811 ! Minloc where hgt - z_k is positive: level right below brand release
812 k = MINLOC(hgt - z_at_w(i,:,j), dim=1, &
813 mask=( hgt - z_at_w(i,:,j) >= 0.0_dp ))
815 z_lower_at_p = u_2d_interp(xp=x0, yp=y0, u_i0j0=z_at_w(i, k,j), u_i0j1=z_at_w(i, k,j+1),&
816 u_i1j0=z_at_w(i+1,k,j), u_i1j1=z_at_w(i+1,k,j+1))
818 z_upper_at_p = u_2d_interp(xp=x0, yp=y0, u_i0j0=z_at_w(i, k+1,j), u_i0j1=z_at_w(i, k+1,j+1),&
819 u_i1j0=z_at_w(i+1,k+1,j), u_i1j1=z_at_w(i+1,k+1,j+1))
821 zdz = (hgt - z_lower_at_p) / (z_upper_at_p-z_lower_at_p)
822 hgt2k = k + zdz
824 END FUNCTION hgt2k
825 !=============================================================
826 !=============================================================
831 !=============================================================
832 PURE &
833 SUBROUTINE releaseijk2atm (nij_2d, sr_x, sr_y, fcat, maxhgt_usr, pi, pj, pk, nij)
834 !=============================================================
836 IMPLICIT NONE
838 INTEGER, INTENT(IN), DIMENSION(:,:) :: nij_2d ! do i need ifps here? No
839 INTEGER, INTENT(IN), DIMENSION(:,:) :: fcat ! do i need ifps here? No
840 INTEGER, INTENT(IN) :: sr_x, sr_y
841 REAL, INTENT(IN) :: maxhgt_usr
843 REAL, INTENT(INOUT), DIMENSION(:) :: pj, pi, pk
844 INTEGER, INTENT(INOUT), DIMENSION(:) :: nij
846 INTEGER, ALLOCATABLE, DIMENSION(:) :: fpi,fpj
847 INTEGER :: i, j, cnt
849 !-------------------------------------------------------------------------------
850 ! Returns a 1-D array with the release positions (i, j, k[m])
851 ! and number of particles to release (nij)
852 ! converted from fire refine grid to wrf grid
853 ! (i,j k are type float and represent the original refined grid position)
854 !-------------------------------------------------------------------------------
856 cnt = COUNT(nij_2d > 0)
858 pi(:) = ZERO_dp
859 pj(:) = ZERO_dp
860 pk(:) = ZERO_dp
861 nij(:) = 0
863 ALLOCATE(fpi(cnt), fpj(cnt))
864 fpi(:) = 0
865 fpj(:) = 0
868 ! fpi = RESHAPE([( [(i, i=ifps, ifpe)], j=jfps, jfpe)], SHAPE(fr_arr))
869 ! fpj = RESHAPE([( [(j, i=ifps, ifpe)], j=jfps, jfpe)], SHAPE(fr_arr))
871 fpi = PACK(RESHAPE([( [(i, i=ifps, ifpe)], j=jfps, jfpe)], SHAPE(nij_2d)), mask=(nij_2d > 0))
872 fpj = PACK(RESHAPE([( [(j, i=ifps, ifpe)], j=jfps, jfpe)], SHAPE(nij_2d)), mask=(nij_2d > 0))
873 nij = PACK(nij_2d, mask=(nij_2d > 0))
875 !-------------------------------------------------------------------------------
876 ! convert from refined fire grid to atm grid
877 !-------------------------------------------------------------------------------
878 pi = fire2atm_ij(fpi, sr_x )
879 pj = fire2atm_ij(fpj, sr_y )
881 ! Release height and firebrand properties can depend on fuel category (in future developments)
882 ! Because fcat is on the refined grid, hgt and fprop must be specified here
884 ! pk = firebrand_gen_maxhgt(PACK(fcat, mask=(nij_2d > 0))) ! Use this to specify hgts for fuel types
885 pk = maxhgt_usr
888 END SUBROUTINE releaseijk2atm
889 !=============================================================
890 !=============================================================
894 !=============================================================
895 ELEMENTAL FUNCTION fire2atm_ij (fp_ij, sr)
896 !=============================================================
897 ! Convert refined subgrid to default atm grid
898 ! Use function separately for i and j:
899 ! pi = fire2atm_ij(fpi, grid%sr_x )
900 ! pj = fire2atm_ij(fpj, grid%sr_y )
902 IMPLICIT NONE
904 REAL :: fire2atm_ij
905 INTEGER, INTENT(IN) :: fp_ij
906 INTEGER, INTENT(IN) :: sr
908 fire2atm_ij = ( dp1+ REAL(fp_ij - dp1)/REAL(sr) )
910 ! fpij = 298, sr = 3: 297/3+1 = 100
911 ! fpij = 299, sr = 3: 298/3+1 = 100.333
912 ! fpij = 300, sr = 3: 299/3+1 = 100.666
913 ! fpij = 301, sr = 3: 300/3+1 = 101
915 ! fpij = 297, sr = 4: 296/4+1 = 75.00
916 ! fpij = 298, sr = 4: 297/4+1 = 75.25
917 ! fpij = 299, sr = 4: 298/4+1 = 75.50
918 ! fpij = 300, sr = 4: 299/4+1 = 75.75
919 ! fpij = 301, sr = 4: 299/4+1 = 76.00
921 END FUNCTION fire2atm_ij
922 !=============================================================
923 !=============================================================
927 !=============================================================
928 ELEMENTAL FUNCTION atm2fire_ij (pij, sr)
929 !=============================================================
930 ! Convert default atm grid to lower left index of refined grid
931 ! (atm_i, atm_j) = (f_i : f_i + srx , f_j : f_j + sry)
933 IMPLICIT NONE
935 INTEGER :: atm2fire_ij
936 INTEGER, INTENT(IN) :: pij
937 INTEGER, INTENT(IN) :: sr
939 atm2fire_ij = (pij - 1) * sr + 1
941 ! pij = 100, sr = 3: 99*3+1 = 298
942 ! pij = 101, sr = 3: 100*3+1 = 301
944 END FUNCTION atm2fire_ij
945 !=============================================================
946 !=============================================================
950 ! !=============================================================
951 ! SUBROUTINE subgrid_fire_property(fire_ij, fuel_frac, grid)
952 ! !=============================================================
954 ! ! Not used - may be delted from final code version
956 ! IMPLICIT NONE
958 ! TYPE(domain), INTENT(IN) :: grid ! input data **Note: Intent IN**
959 ! REAL, INTENT(OUT), DIMENSION(:) :: fuel_frac
960 ! INTEGER, INTENT(IN), DIMENSION(:,:) :: fire_ij
962 ! INTEGER, ALLOCATABLE, DIMENSION(:) :: fpi, fpj
963 ! INTEGER :: k, cnt, sri, srj, srx, sry
965 ! srx = grid%sr_x
966 ! sry = grid%sr_y
968 ! cnt = SIZE(fire_ij(:,1))
969 ! fpi = atm2fire_ij( fire_ij(:,1), srx)
970 ! fpj = atm2fire_ij( fire_ij(:,2), sry)
971 ! !ALLOCATE( flame(cnt), fuel(cnt))
973 ! ! Get flame and fuel_frac maxvals from subgrids
974 ! DO k=1,SIZE(fire_ij(:,1))
976 ! sri = fpi(k)
977 ! srj = fpj(k)
979 ! fuel_frac(k) = MAXVAL(grid%fuel_frac( sri: sri+srx-1, srj : srj+sry-1))
981 ! !WRITE (msg,*) MAXVAL(grid%zsf(sri:sri+dsx,srj:srj+dsy)) always zero??
983 ! ENDDO
985 ! END SUBROUTINE subgrid_fire_property
986 ! !=============================================================
987 ! !=============================================================
991 ! !=============================================================
992 ! PURE FUNCTION uniq_ij (fpi, fpj)
993 ! !=============================================================
994 ! ! remove duplicated grid points (after converting from refined fire grid to grid)
996 ! ! Not used but quite handy - do not delete it.
998 ! IMPLICIT NONE
1000 ! INTEGER, INTENT(IN), DIMENSION(:) :: fpi, fpj
1001 ! INTEGER, ALLOCATABLE, DIMENSION(:,:) :: uniq_ij
1002 ! LOGICAL, ALLOCATABLE, DIMENSION(:) :: mask
1003 ! INTEGER :: cnt, i
1005 ! cnt = SIZE(fpi)
1007 ! ALLOCATE(mask(cnt))
1008 ! mask = .TRUE.
1010 ! DO i = cnt, 2, -1
1011 ! mask(i) = .NOT.( ANY(fpi(:i-1) == fpi(i) .AND. &
1012 ! fpj(:i-1) == fpj(i)))
1013 ! END DO
1015 ! ALLOCATE( uniq_ij(COUNT(mask),2) )
1016 ! uniq_ij(:,1) = REAL(PACK(fpi, mask))
1017 ! uniq_ij(:,2) = REAL(PACK(fpj, mask))
1019 ! END FUNCTION uniq_ij
1021 ! !=============================================================
1022 ! !=============================================================
1026 !=============================================================
1027 SUBROUTINE prep_release_hgt(release_i, release_j, release_k, release_n, release_prop, levrand_seed, landhgt)
1028 !=============================================================
1030 !----------------------------------------------------------
1031 ! Set release locations for new particles
1032 ! and distribute particles over fs_gen_levels heights
1033 !----------------------------------------------------------
1035 IMPLICIT NONE
1037 TYPE(p_properties), INTENT(INOUT), DIMENSION(:) :: release_prop
1038 REAL, INTENT(INOUT), DIMENSION(:) :: release_i, release_j
1039 REAL, INTENT(INOUT), DIMENSION(:) :: release_k
1040 INTEGER, INTENT(IN), DIMENSION(:) :: release_n
1041 INTEGER, INTENT(IN) :: levrand_seed
1042 REAL, INTENT(IN) :: landhgt
1044 INTEGER :: ii, kk, cnt ! counters
1045 INTEGER :: nseeds
1046 REAL :: minhgt
1047 REAL, ALLOCATABLE, DIMENSION(:) :: rand_arr, frachgt
1048 INTEGER, ALLOCATABLE, DIMENSION(:) :: seeds
1050 minhgt = dp1 + landhgt
1052 !-------------------------------------------------------------------------------
1053 ! Calculate number of releases and heights from this fire
1054 !-------------------------------------------------------------------------------
1056 ALLOCATE(frachgt(fs_gen_levels))
1057 ALLOCATE(rand_arr(SIZE(release_n)*fs_gen_levels))
1059 IF (fs_gen_levels == 1) THEN
1060 frachgt(:) = dp1
1061 ELSE
1062 frachgt(:) = [( REAL(kk-1)*(dp1/REAL(fs_gen_levels-1)), kk=fs_gen_levels, 1, -1)]
1063 ENDIF
1065 IF (levrand_seed > 0) THEN
1067 nseeds = SIZE(rand_arr)
1068 CALL random_seed(size=nseeds)
1069 ALLOCATE(seeds(nseeds))
1070 seeds = [(( release_i(ii) * release_j(ii) * levrand_seed * kk, &
1071 kk=1, fs_gen_levels), &
1072 ii=1, SIZE(release_n))] ! force seed to vary across space
1073 CALL random_seed(put = seeds)
1074 DEALLOCATE(seeds)
1075 CALL random_number(rand_arr)
1077 ENDIF
1080 ii = SIZE(release_n)+1 ! array index for release at various height (append to end of the array)
1082 DO kk=2, fs_gen_levels ! release at the remaining levels below maxhgt (k=1 is already filled with maxhgt)
1083 DO cnt=1,SIZE(release_n) ! loop over release gridpoints and set a new release at the same ij for 1 < hgt < release_k
1085 !-------------------------------------------------------------------------------
1086 ! Fill release array
1087 !-------------------------------------------------------------------------------
1089 release_i(ii) = release_i(cnt) ! ii increments inside this loop
1090 release_j(ii) = release_j(cnt)
1091 IF (release_k(cnt) >= minhgt) THEN
1092 release_k(ii) = ((release_k(cnt)-minhgt) * frachgt(kk)) + minhgt ! keep a minimum height of 1-m above landhgt
1093 ! if we want to perturb the max_release_hgt in releases_k(cnt), we can change its value now
1094 IF (levrand_seed > 0) &
1095 release_k(ii) = ((release_k(cnt)-minhgt) * rand_arr(ii)) + minhgt
1096 ELSE
1097 release_k(ii) = ((release_k(cnt)-dp1) * frachgt(kk)) + dp1 ! keep a minimum height of 1-m above landhgt
1098 IF (levrand_seed > 0) &
1099 release_k(ii) = ((release_k(cnt)-dp1) * rand_arr(ii)) + dp1
1100 ENDIF
1103 release_prop(ii) = release_prop(cnt)
1105 ! IF (cnt < 3) THEN
1106 ! WRITE (msg,'(2(i4,1x),i8,1x,6(f8.3,1x))') cnt, kk, ii, &
1107 ! release_i(ii), release_j(ii), release_k(cnt), release_k(ii), frachgt(kk), rand_arr(ii)
1108 ! CALL wrf_debug (0, 'SPFire_prep_release c,k,i, ijk: '//msg)
1109 ! ENDIF
1111 ii = ii + 1 ! increment release array index
1113 ENDDO
1114 ENDDO
1116 ! Add another loop here to perturb release_k from ii=,SIZE(release_n) -> currently set to maxhgt
1118 END SUBROUTINE prep_release_hgt
1119 !=============================================================
1120 !=============================================================
1124 !=============================================================
1125 FUNCTION firebrand_property(arr) RESULT(prop)
1126 !=============================================================
1128 IMPLICIT NONE
1130 TYPE (p_properties) :: prop
1131 REAL, INTENT(IN), DIMENSION(:) :: arr
1132 REAL, PARAMETER :: pi = 4.0_dp * ATAN (1.0_dp)
1134 prop%p_diam = arr(1) ! 10.0_dp ! mm
1135 prop%p_effd = arr(2) ! 10.0_dp ! mm
1136 prop%p_temp = arr(3) ! 900.0_dp ! K
1137 ! termvel should start with 0m/s because firebrand is released at top of trajectory
1138 prop%p_tvel = arr(4) ! 0.0_dp ! m/s
1139 prop%p_mass = (firebrand_dens * pi * (prop%p_effd/1000.0_dp)**2)/6.0_dp ! g
1141 END FUNCTION firebrand_property
1143 !=============================================================
1144 !=============================================================
1148 !=============================================================
1149 FUNCTION idx_packed_1d(mask) RESULT(mask_idx)
1150 !=============================================================
1152 ! Return an array of indices where mask is true
1153 IMPLICIT NONE
1154 LOGICAL, INTENT(IN), DIMENSION(:) :: mask
1155 !INTEGER, INTENT(IN) :: active_br
1156 INTEGER, ALLOCATABLE, DIMENSION(:) :: mask_idx
1157 INTEGER :: nresets, ii
1159 nresets = COUNT(mask)
1160 ALLOCATE(mask_idx(nresets))
1161 mask_idx = PACK([(ii, ii=1,SIZE(mask))], mask) ! Get flag indices
1163 END FUNCTION idx_packed_1d
1164 !=============================================================
1165 !=============================================================
1170 !=============================================================
1171 PURE &
1172 SUBROUTINE firebrand_physics(dt, hgt, loc_p, loc_t, loc_d, loc_w, fbprop)
1173 !=============================================================
1175 IMPLICIT NONE
1177 TYPE(p_properties), INTENT(INOUT) :: fbprop
1178 REAL, INTENT(IN) :: loc_p, loc_t, loc_d, loc_w ! local pressure, temp, density, w
1179 REAL, INTENT(IN) :: dt
1180 REAL, INTENT(INOUT) :: hgt
1181 !REAL, INTENT(INOUT) :: p_mass, p_diam, p_effd, p_temp, p_tvel
1183 !-------------------------------------------------------------------------------
1184 ! firebrand burnout
1185 !-------------------------------------------------------------------------------
1187 ! WRITE (msg,'(6(f12.8,1x))') fbprop%p_mass, fbprop%p_diam, &
1188 ! fbprop%p_effd, fbprop%p_temp, fbprop%p_tvel, hgt
1189 ! CALL wrf_debug (wrfdbg, 'SPFire br physics1 m,d,e,t,v,h : '//msg)
1191 CALL burnout(p_mass = fbprop%p_mass, & ! firebrand mass (g)
1192 p_diam = fbprop%p_diam, & ! firebrand diameter (mm)
1193 p_effd = fbprop%p_effd, & ! firebrand effective diameter (mm)
1194 p_temp = fbprop%p_temp, & ! firebrand temperature (K)
1195 p_tvel = fbprop%p_tvel, & ! firebrand terminal velocity (m/s)
1196 aird = loc_d, & ! g/m3
1197 pres = loc_p, & ! Pa
1198 temp = loc_t, & ! K
1199 loc_w = loc_w, & ! m/s
1200 dt = dt)
1202 ! WRITE (msg,'(4(f12.8,1x))') fbprop%p_mass, fbprop%p_diam, &
1203 ! fbprop%p_effd, fbprop%p_temp
1204 ! CALL wrf_debug (wrfdbg, 'SPFire br physics2 m, d, e, t: '//msg)
1206 !-------------------------------------------------------------------------------
1207 ! firebrand terminal velocity
1208 !-------------------------------------------------------------------------------
1210 CALL termvel(p_diam = fbprop%p_diam, &
1211 p_effd = fbprop%p_effd,&
1212 p_temp = fbprop%p_temp, &
1213 p_tvel = fbprop%p_tvel, &
1214 aird = loc_d, & ! g/m3
1215 pres = loc_p, & ! Pa
1216 temp = loc_t, & ! K
1217 dt = dt, & ! sec
1218 hgt = hgt) ! m
1220 ! WRITE (msg,'(6(f12.8,1x))') fbprop%p_mass, fbprop%p_diam, &
1221 ! fbprop%p_effd, fbprop%p_temp, fbprop%p_tvel, hgt
1222 ! CALL wrf_debug (wrfdbg, 'SPFire br physics2 m,d,e,t,v,h : '//msg)
1224 ! **********
1225 ! NOTE:
1226 ! Should we calculate in double precision when diam is in meters?
1227 ! double: real(8), 15 precision digits
1228 ! single: real(4), 7 precision digits
1230 END SUBROUTINE firebrand_physics
1232 !=============================================================
1233 !=============================================================
1237 !=============================================================
1238 ELEMENTAL SUBROUTINE burnout(pres, aird, temp, loc_w, dt, p_mass, p_diam, p_effd, p_temp, p_tvel)
1239 !=============================================================
1241 !-------------------------------------------------------------------------------
1242 ! This subroutine was developed in RAL-NCAR/UCAR, 2019-2020 by:
1243 ! *** T. W. Juliano ***
1244 ! based on
1245 ! Tse & Fernandez-Pello 1998 (DOI: 10.1016/ S0379-7112(97)00050-7)
1246 ! Bhutia, Jenkins & Sun 2010 (DOI:10.3894/JAMES.2010.2.4)
1247 !-------------------------------------------------------------------------------
1249 IMPLICIT NONE
1251 REAL, INTENT(IN) :: pres, aird, temp, loc_w, dt
1252 REAL, INTENT(INOUT) :: p_mass ! firebrand mass (g)
1253 REAL, INTENT(INOUT) :: p_diam ! firebrand diameter (mm)
1254 REAL, INTENT(INOUT) :: p_effd ! firebrand effective diameter (mm)
1255 REAL, INTENT(INOUT) :: p_temp ! firebrand temperature (K)
1256 REAL, INTENT(IN) :: p_tvel ! new firebrand terminal velocity (m/s)
1258 !-------------------------------------------------------------------------------
1259 ! Constants are at module top
1260 !-------------------------------------------------------------------------------
1262 REAL, PARAMETER :: pi = 4.0_dp * ATAN (1.0_dp)
1264 REAL :: aird2, wind, gama
1265 REAL :: pdia, pedia
1266 REAL :: p_effd2, pdia_new4, reyn, beta, deff, dmvc, dmvc2
1267 REAL :: parta, partb, dtemp, nuss, hbar, fbvol, qcon, qrad, partc
1268 REAL :: pratio, cpmix
1270 p_mass = MAX(p_mass, EPSILON(dp1))
1271 !-------------------------------------------------------------------------------
1273 ! air density surrounding firebrand
1274 aird2 = 1000.0_dp * pres / ( rd * p_temp )
1276 ! Assume vel diff between firebrand and wind is due only to
1277 ! vertical; that is, assume particle moves with horiz wind
1278 wind = ABS(loc_w - p_tvel)
1280 ! particle diameter (mm) to (m)
1281 pdia = p_diam / 1000.0_dp
1283 ! particle effective diameter (mm) to (m)
1284 pedia = p_effd / 1000.0_dp
1286 ! local atmospheric dynamic viscosity according to Sutherland's law
1287 dmvc = (b_coeff * temp**1.5_dp) / (temp + s_coeff)
1289 ! dynamic viscosity surrounding firebrand according to Sutherland's law
1290 dmvc2 = (b_coeff * p_temp**1.5_dp) / (p_temp + s_coeff)
1292 ! kinematic viscosity average for atmosphere and firebrand
1293 gama = ((dmvc / aird) + (dmvc2 / aird2))/2.0_dp
1295 ! reynolds number
1296 reyn = wind * pdia / gama
1298 ! nusselt number
1299 nuss = 2.0_dp + 0.6_dp * reyn**(dp05) * pr**(r1o3)
1301 ! convection heat transfer coefficient
1302 hbar = nuss * thcon / pdia
1304 ! burning rate
1305 beta = beta0 + beta0 * (0.276_dp * reyn**(dp05) * shmt**(r1o3))
1307 ! WRITE (msg,'(9(f12.8,1x),f12.3,1x,f12.8)') aird2, wind, pdia, &
1308 ! pedia, dmvc, dmvc2, &
1309 ! gama, reyn, nuss, hbar, beta
1310 ! CALL wrf_debug (wrfdbg, 'SPFire burnout1: '//msg)
1312 ! change in firebrand mass diameter (m)
1313 parta = pdia**4
1314 partb = SQRT(3.0_dp) * (beta**2) * (dt**2)
1315 pdia_new4 = parta - partb
1316 p_diam = pdia_new4**(0.25_dp)
1318 ! change in firebrand effective mass diameter (m)
1319 p_effd2 = (pedia**2) - beta * dt
1320 p_effd = SQRT(p_effd2)
1322 ! new firebrand mass (g)
1323 p_mass = (firebrand_dens * pi * p_effd2 * p_effd)/6.0_dp
1325 ! new firebrand volume (m3)
1326 fbvol = p_mass / firebrand_dens
1328 ! new firebrand temp (K)
1329 qcon = hbar * (p_temp - temp)
1330 qrad = sboltz * emiss * (p_temp**4 - temp**4)
1331 pratio = p_effd / p_diam
1332 cpmix = (pratio * cpw) + ((1.0_dp - pratio) * cpc)
1333 partc = 6.0_dp / (p_mass * cpmix * p_diam)
1334 dtemp = dt * fbvol * partc * (qcon + qrad)
1335 p_temp = p_temp - dtemp
1337 ! WRITE (msg,'(8(f14.12,1x))') parta, partb, pdia_new4, p_diam, p_effd2, p_effd, p_mass, fbvol
1338 ! CALL wrf_debug (wrfdbg, 'SPFire burnout2: '//msg)
1339 ! WRITE (msg,'(2(f12.3,1x),f12.8,1x,2(f12.3,1x),2(f12.8,1x))') qcon, qrad, pratio, cpmix, partc, dtemp, p_temp
1340 ! CALL wrf_debug (wrfdbg, 'SPFire burnout3: '//msg)
1342 ! firebrand diameter from (m) to (mm)
1343 p_diam = 1000.0_dp * p_diam
1344 p_effd = 1000.0_dp * p_effd
1346 END SUBROUTINE burnout
1348 !=============================================================
1349 !=============================================================
1353 !=============================================================
1354 ELEMENTAL &
1355 SUBROUTINE termvel(p_diam, p_effd, p_tvel, p_temp, hgt, dt, pres, aird, temp)
1356 !=============================================================
1358 !-------------------------------------------------------------------------------
1359 ! This subroutine was developed in RAL-NCAR/UCAR, 2019-2020 by:
1360 ! *** T. W. Juliano ***
1361 ! based on Bhutia, Jenkins & Sun 2010 (DOI:10.3894/JAMES.2010.2.4)
1362 !-------------------------------------------------------------------------------
1364 ! Constants are at module top
1366 IMPLICIT NONE
1368 REAL, INTENT(IN) :: pres, aird, temp, dt
1369 REAL, INTENT(IN) :: p_diam ! firebrand diameter (mm)
1370 REAL, INTENT(IN) :: p_effd ! firebrand effective diameter (mm)
1371 REAL, INTENT(IN) :: p_temp ! firebrand temperature (K)
1372 REAL, INTENT(INOUT):: hgt ! particle position at height agl (m)
1373 REAL, INTENT(INOUT):: p_tvel ! new firebrand terminal velocity (m/s)
1375 REAL :: pbot, ptop, aird2, pdia, parta, partb
1376 REAL :: drop, vt, pratio ! fbtype,
1378 !-------------------------------------------------------------------------------
1380 ! air density around particle
1381 aird2 = 1000.0_dp * pres / (rd * p_temp)
1383 ! particle diameter (mm) to (m)
1384 pdia = p_diam /1000.0_dp
1386 ! terminal velocity
1387 pratio = p_effd / p_diam
1388 parta = ( (pratio*firebrand_dens) + ( (dp1-pratio) * firebrand_dens_char) ) *pdia*grav
1389 partb = 3.0_dp * ((aird + aird2) / 2.0_dp)*rcd
1391 pratio = MAX( MIN(HUGE(1.0),(parta / partb)), TINY(1.0)) ! Too large (small) values for the float type size leads to seg fault in debug mode
1393 ! WRITE (msg,'(4(f12.6,1x))') parta, partb, pratio, SQRT(pratio)
1394 ! CALL wrf_debug (wrfdbg, 'SPFire termvel : '//msg)
1396 vt = SQRT(pratio)
1398 ! change in vertical position due to settling
1399 drop = vt * ABS(dt)
1400 pbot = MAX(0.0_dp, hgt-drop)
1401 hgt = pbot
1403 p_tvel=vt
1405 END SUBROUTINE termvel
1407 !=============================================================
1408 !=============================================================
1411 !=============================================================
1412 PURE &
1413 SUBROUTINE get_local_met(xp, yp, zp, loc_p, loc_d, loc_t, p, t, d, ihs, jhs, ihe, jhe)
1414 !=============================================================
1416 !-------------------------------------------------------------------------------
1417 !
1418 ! Calculate met from variables updated in HALO after the call to microphysics:
1419 !
1420 ! HALO_EM_E_5:48: u_1,u_2,v_1,v_2,w_1,w_2,t_1,t_2,ph_1,ph_2,tke_1,tke_2,;4:mu_1,mu_2
1424 IMPLICIT NONE
1426 !TYPE(domain), INTENT(IN) :: grid
1427 INTEGER, INTENT(IN) :: ihs, jhs, ihe, jhe ! tile +- halo array bounds
1428 REAL, INTENT(IN) :: xp, yp, zp ! particle position
1429 REAL, INTENT(OUT) :: loc_p, loc_t, loc_d ! local met properties
1430 REAL, INTENT(IN), DIMENSION(ims:,kms:,jms:):: p,t,d
1432 INTEGER :: i, j, k , kk
1433 REAL :: tmp1, tmp2, zph ! corresponding zp at half levels
1435 k = FLOOR(zp)
1436 j = MIN(MAX(FLOOR(yp-dp05), jhs), jhe)
1437 i = MIN(MAX(FLOOR(xp-dp05), ihs), ihe)
1439 ! physics variables with (i,j) = (1,1) corresponds to (x,y) = (1.5,1.5)
1440 ! T(i=1): x=1.5
1441 ! T(i=2): x=2.5
1442 ! To interpolate T(i=1:2), x=[1.5:2.5[
1443 ! See comment on function hgt2k
1445 ! boundaries ims, ..., ke are declared in module scope
1447 ! p_hydrostatic is on full levels: no vertical (de)staggering needed
1448 loc_p = u_3d_interp( x=xp-dp05, y=yp-dp05, z=zp, &
1449 u_i0j0_bot = p(i, k, j), u_i0j1_bot = p(i, k, j+1),&
1450 u_i1j0_bot = p(i+1,k, j), u_i1j1_bot = p(i+1,k, j+1),&
1451 u_i0j0_top = p(i, k+1, j), u_i0j1_top = p(i, k+1,j+1),&
1452 u_i1j0_top = p(i+1,k+1, j), u_i1j1_top = p(i+1,k+1,j+1))
1454 ! theta moist is on full levels: no vertical (de)staggering needed: th8w
1455 loc_t = u_3d_interp( x=xp-dp05, y=yp-dp05, z=zp, &
1456 u_i0j0_bot = t(i, k, j), u_i0j1_bot = t(i, k, j+1),&
1457 u_i1j0_bot = t(i+1,k, j), u_i1j1_bot = t(i+1,k, j+1),&
1458 u_i0j0_top = t(i, k+1, j), u_i0j1_top = t(i, k+1,j+1),&
1459 u_i1j0_top = t(i+1,k+1, j), u_i1j1_top = t(i+1,k+1,j+1))
1461 ! Variables on half levels: provide the k index corresponding to particle bottom & top levels
1462 ! Is a linear interpolation ok?
1464 ! when zp = 1.0, particle is at the surface on the stag grid,
1465 ! and at k=0.5 on half levels (below the first level)
1466 ! when zp = 1.5, particle is in the middle of the stag grid, and at k=1.0 on half level
1467 ! This is the location where the unstaggered variable value is valid,
1468 ! so the distance between the particle and this point must be 0.
1470 ! Physics variables are all unstaggered in x and y.
1472 zph = zp - dp05
1473 k = MAX(1, FLOOR(zph))
1475 loc_d = u_3d_interp( x=xp-dp05, y=yp-dp05, z=zph, &
1476 u_i0j0_bot = d(i, k, j), u_i0j1_bot = d(i, k, j+1),&
1477 u_i1j0_bot = d(i+1,k, j), u_i1j1_bot = d(i+1,k, j+1),&
1478 u_i0j0_top = d(i, k+1, j), u_i0j1_top = d(i, k+1,j+1),&
1479 u_i1j0_top = d(i+1,k+1, j), u_i1j1_top = d(i+1,k+1,j+1))
1481 loc_d = loc_d * 1000.0_dp ! Convert from kg/m3 to g/m3
1484 END SUBROUTINE get_local_met
1486 !=============================================================
1488 !=============================================================
1492 !=============================================================
1493 SUBROUTINE firebrand_spotting_em_driver( &
1494 !=============================================================
1495 cf, &
1496 grid, &
1497 fs_p_id, &
1498 fs_p_src, &
1499 fs_p_dt, &
1500 fs_p_x, &
1501 fs_p_y, &
1502 fs_p_z, &
1503 fs_gen_inst, &
1504 fs_p_mass, &
1505 fs_p_diam, &
1506 fs_p_effd, &
1507 fs_p_temp, &
1508 fs_p_tvel, &
1509 fs_last_gen_dt, &
1510 fs_gen_idmax, &
1511 fs_fire_ROSdt, &
1512 fs_fire_area, &
1513 fs_count_landed_all, &
1514 fs_count_landed_hist, &
1515 fs_landing_mask, &
1516 fs_spotting_lkhd, &
1517 fs_frac_landed, &
1518 fs_fuel_spotting_risk, &
1519 fs_count_reset)
1522 !-------------------------------------------------------------------------------
1523 ! Modules
1524 !-------------------------------------------------------------------------------
1526 USE module_domain, ONLY: domain_get_time_since_sim_start, &
1527 domain_clock_get, is_alarm_tstep
1528 USE module_domain_type, ONLY: HISTORY_ALARM, restart_alarm, AUXHIST23_ALARM
1529 USE module_state_description, ONLY: p_qv, num_moist, param_first_scalar
1530 USE module_utility, ONLY: WRFU_Alarm !WRFU_Clock,
1531 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
1532 USE module_firebrand_spotting_mpi, ONLY: fs_mpi_init, &
1533 fs_mpi_sendbuff1_int, &
1534 fs_mpi_sendbuff1_real, &
1535 fs_mpi_recvbuff1_int, &
1536 fs_mpi_recvbuff1_real, &
1537 fs_mpi_nothing2send, &
1538 fs_mpi_send2neighbors, &
1539 fs_mpi_sendbuffsize, &
1540 fs_mpi_sendbuff_int, &
1541 fs_mpi_sendbuff_real, &
1542 fs_mpi_recvbuffsize, &
1543 fs_mpi_recvfrompatch_int, &
1544 fs_mpi_recvfrompatch_real, &
1545 fs_mpi_sendbuff_real, &
1546 fs_mpi_checkreceive, &
1547 fs_mpi_recv
1548 USE module_firebrand_spotting_mpi, ONLY: neighbors, my_id, task_id, mpiprocs, &
1549 left_id, right_id, up_id, down_id, &
1550 upleft_id, upright_id, downleft_id, downright_id
1551 #endif
1554 IMPLICIT NONE
1555 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
1556 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
1557 #endif
1558 !-------------------------------------------------------------------------------
1559 ! WRF variables
1560 !-------------------------------------------------------------------------------
1562 TYPE(domain), INTENT(IN) :: grid ! input data **Note: Intent IN**
1563 TYPE (GRID_config_rec_type),INTENT(IN) :: cf ! type available by module host association
1565 !-------------------------------------------------------------------------------
1566 ! Dummy Arguments
1567 !-------------------------------------------------------------------------------
1569 LOGICAL, INTENT(INOUT) :: fs_count_reset
1570 INTEGER, INTENT(INOUT) :: fs_last_gen_dt, fs_gen_idmax
1572 ! Firebrand Spotting Particle Properties (fs_p_*)
1573 INTEGER, INTENT(INOUT), DIMENSION(:) :: fs_p_id, fs_p_src, fs_p_dt
1574 REAL, INTENT(INOUT), DIMENSION(:) :: fs_p_x, fs_p_y, fs_p_z ! positions are relative to grid edges: particle at grid center point (x,y,z) = (1.5, 1.5, 1.5)
1575 REAL, INTENT(INOUT), DIMENSION(:) :: fs_p_mass, fs_p_diam, fs_p_effd
1576 REAL, INTENT(INOUT), DIMENSION(:) :: fs_p_temp, fs_p_tvel
1577 REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: fs_count_landed_all, fs_fire_area, fs_fuel_spotting_risk ! Use memory bounds so positons match indices
1578 REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: fs_count_landed_hist, fs_spotting_lkhd, fs_frac_landed
1579 INTEGER, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: fs_landing_mask, fs_gen_inst
1580 REAL, INTENT(INOUT), DIMENSION(ifps:ifpe,jfps:jfpe) :: fs_fire_ROSdt
1582 !-------------------------------------------------------------------------------
1583 ! Parameters
1584 !-------------------------------------------------------------------------------
1586 !-------------------------------------------------------------------------------
1587 ! Local Variables - Arrays
1588 !-------------------------------------------------------------------------------
1590 ! Subroutine Control
1591 LOGICAL, DIMENSION(fs_array_maxsize) :: sparse_mask
1593 TYPE(p_properties), DIMENSION(fs_array_maxsize) :: fs_p_prop
1595 ! Physics from WRF
1596 REAL, DIMENSION(ims:ime, ks:ke, jms:jme):: w, z, p_hyd, dz8w, z_at_w ! z_at_w: height at interface of mdl level
1597 REAL, DIMENSION(ims:ime, ks:ke, jms:jme):: u, v, rho, th_phy ! *_phy: half levels
1598 REAL, DIMENSION(ims:ime, ks:ke, jms:jme):: qtot, th8w, p_phy, p8w ! *8w: full levels
1599 REAL, DIMENSION(ims:ime, jms:jme) :: msft, w1
1600 REAL, DIMENSION(ks:ke) :: znw
1601 REAL, ALLOCATABLE, DIMENSION(:) :: xout, yout, zout ! firebrand position after advection
1603 !-------------------------------------------------------------------------------
1604 ! Likelihood
1605 !-------------------------------------------------------------------------------
1607 REAL, ALLOCATABLE, DIMENSION(:,:) :: fuel_spotting_risk_fr, spotthreat_fr ! only need it for history intervals
1609 ! only need it in hist interval - do not hold memory
1610 LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: landing_mask
1611 INTEGER, ALLOCATABLE, DIMENSION(:) :: landing_cnt
1612 INTEGER, ALLOCATABLE, DIMENSION(:) :: np_landed_cnt, np_nlanded ! (MPI: only Master allocates these)
1613 REAL, ALLOCATABLE, DIMENSION(:) :: np_landed_risk, np_frac_landed, np_lkhd ! (MPI: only Master allocates these)
1614 REAL, ALLOCATABLE, DIMENSION(:) :: landed_risk, recv_spot_lkhd, recv_frac_landed
1615 REAL, ALLOCATABLE, DIMENSION(:) :: tmp1
1616 INTEGER :: ndep, ndep_total
1617 !-------------------------------------------------------------------------------
1619 REAL, DIMENSION(ims:ime, jms:jme) :: dt_sum_landed ! sum timestep's deposited fb in 2-D memory array
1620 LOGICAL, SAVE :: hist_flag = .FALSE. ! time of history output
1622 !-------------------------------------------------------------------------------
1623 ! Release
1624 !-------------------------------------------------------------------------------
1626 INTEGER, ALLOCATABLE, DIMENSION(:) :: np_Ngrdpts, loc_ij
1627 REAL, ALLOCATABLE, DIMENSION(:) :: np_MeanPot
1628 REAL, ALLOCATABLE, DIMENSION(:,:) :: RelPot, fcat_factor, fr_rate_dummy
1629 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: fcat_fr, n_rel_fr2d
1630 LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: rel_mask
1632 INTEGER :: Ngrdpts, Npts_total, relpts, rankpts
1633 REAL :: TotPot, MeanPot, PotThr, rel_ratio, cdf_coeff, LimPot
1635 ! Firebrand Release and Transport
1636 TYPE(p_properties), ALLOCATABLE, DIMENSION(:) :: release_prop
1637 REAL, ALLOCATABLE, DIMENSION(:) :: release_i, release_j, release_k
1638 INTEGER, ALLOCATABLE, DIMENSION(:) :: release_n
1640 !-------------------------------------------------------------------------------
1641 ! Local Variables - Scalars
1642 !-------------------------------------------------------------------------------
1644 ! Note:
1645 ! Variables for tile, memory, and domain bounds are in module scope & are assigned in firebrand_spotting_init
1647 INTEGER :: kk, pp, ii ! counters
1648 INTEGER :: kk1, kk2, k_end
1649 INTEGER :: active_br, prior_br, nresets, idmax, firepts, nbr_sum
1650 REAL :: rdx, rdy, dt, sr_xy ! grid length (1/dx) and wrf's time-step
1652 INTEGER :: firebrand_max_life_dt
1653 REAL :: firebrand_land_hgt, firebrand_gen_maxhgt
1654 LOGICAL :: firebrand_gen_levrand
1655 REAL :: firebrand_gen_prop_diam, firebrand_gen_prop_effd
1656 REAL :: firebrand_gen_prop_temp, firebrand_gen_prop_tvel
1657 INTEGER :: levrand_seed
1659 !-------------------------------------------------------------------------------
1660 ! MPI local Variables (some variables are in module scope)
1661 !-------------------------------------------------------------------------------
1663 INTEGER, SAVE :: MasterId = 0
1664 INTEGER, SAVE :: ntiles = 0
1665 INTEGER, SAVE :: myprocid = 0
1666 LOGICAL, SAVE :: IamMaster = .FALSE.
1668 INTEGER, ALLOCATABLE, DIMENSION(:) :: nbr_id, r_id, r_dt, r_src
1669 REAL, ALLOCATABLE, DIMENSION(:) :: r_x, r_y, r_z
1671 ! MPI receive p_properties
1672 REAL, ALLOCATABLE, DIMENSION(:) :: r_p_m,r_p_d,r_p_e,r_p_t,r_p_v
1673 TYPE(p_properties), ALLOCATABLE, DIMENSION(:) :: fb_mdetv
1676 !-------------------------------------------------------------------------------
1677 ! It was a dummy variable back when fuel moisture was not implemented in wrf-v4.0.1
1678 !-------------------------------------------------------------------------------
1680 REAL, DIMENSION(ifps:ifpe,jfps:jfpe) :: fmcg_2df
1682 !===============================================================================
1683 ! Begin
1684 !===============================================================================
1687 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
1688 IF (mpiprocs == 0) &
1689 CALL fs_mpi_init(grid) ! should be called in fspotting_init (?)
1691 ntiles = mpiprocs
1692 myprocid = my_id * 1000000
1693 IF ( wrf_dm_on_monitor() ) THEN
1694 IamMaster = .TRUE.
1695 MasterId = my_id
1696 ENDIF
1697 #else
1698 IamMaster = .TRUE. ! initialized with FALSE. It should be TRUE in serial compilation
1699 ntiles = 1 ! used to allocate the arrays for incomming mpi comms. ntiles=1 in serial builds
1700 myprocid = 1* 1000000
1701 #endif
1703 hist_flag = .FALSE.
1704 IF ( Is_alarm_tstep(grid%domain_clock, grid%alarms(HISTORY_ALARM)) ) hist_flag = .TRUE.
1705 !IF ( Is_alarm_tstep(grid%domain_clock, grid%alarms(AUXHIST23_ALARM)) ) hist_flag = .TRUE.
1706 IF (hist_flag) &
1707 CALL wrf_debug (wrfdbg, 'SPFire_driver: History output on next tstep')
1709 fs_last_gen_dt = fs_last_gen_dt + 1
1710 active_br = COUNT( fs_p_id > 0 )
1712 IF (grid%itimestep == 1) &
1713 fs_fire_ROSdt(ifps:ifpe,jfps:jfpe) =0.0_dp
1715 !===============================================================================
1716 ! Reset array after writting history file
1717 !===============================================================================
1719 IF (fs_count_reset) THEN ! if true
1720 !CALL wrf_debug (wrfdbg, 'SPFire_lkhd: resetting fs_count_reset!')
1721 fs_count_landed_hist(:,:) = ZERO_dp
1722 fs_fuel_spotting_risk(:,:) = ZERO_dp
1723 fs_count_reset = .FALSE.
1725 fs_gen_inst(is:ie,js:je) = 0
1726 ENDIF
1728 fs_landing_mask(ims:ime,jms:jme) = 0
1730 !===============================================================================
1731 ! Assign properties to derived type (used locally to shorten argument list)
1732 !===============================================================================
1734 fs_p_prop%p_mass = ZERO_dp
1735 fs_p_prop%p_diam = ZERO_dp
1736 fs_p_prop%p_effd = ZERO_dp
1737 fs_p_prop%p_temp = ZERO_dp
1738 fs_p_prop%p_tvel = ZERO_dp
1740 fs_p_prop(:active_br)%p_mass = fs_p_mass(:active_br)
1741 fs_p_prop(:active_br)%p_diam = fs_p_diam(:active_br)
1742 fs_p_prop(:active_br)%p_effd = fs_p_effd(:active_br)
1743 fs_p_prop(:active_br)%p_temp = fs_p_temp(:active_br)
1744 fs_p_prop(:active_br)%p_tvel = fs_p_tvel(:active_br)
1746 !===============================================================================
1747 ! Assign user defined values from namelist
1748 !===============================================================================
1750 firebrand_max_life_dt = cf%fs_firebrand_max_life_dt
1751 firebrand_land_hgt = cf%fs_firebrand_land_hgt
1752 firebrand_gen_maxhgt = cf%fs_firebrand_gen_maxhgt
1753 firebrand_gen_levrand = cf%fs_firebrand_gen_levrand
1755 firebrand_gen_prop_diam = cf%fs_firebrand_gen_prop_diam
1756 firebrand_gen_prop_effd = cf%fs_firebrand_gen_prop_effd
1757 firebrand_gen_prop_temp = cf%fs_firebrand_gen_prop_temp
1758 firebrand_gen_prop_tvel = cf%fs_firebrand_gen_prop_tvel
1760 firebrand_dens = cf%fs_firebrand_dens ! 513000.0_dp ! density of firebrand (g / m3)
1761 firebrand_dens_char = cf%fs_firebrand_dens_char ! 299000.0_dp ! density of char (g m-3) [gupta et al 2002; fuel]
1763 !===============================================================================
1764 ! Firebrand Generation
1765 !===============================================================================
1767 !-------------------------------------------------------------------------------
1768 ! Accumulate burn rate between generation intervals
1769 !-------------------------------------------------------------------------------
1771 firepts = 0
1772 ! fs_fire_ROSdt(ifps:ifpe, jfps:jfpe) = fs_fire_ROSdt(ifps:ifpe, jfps:jfpe) + grid%burnt_area_dt(ifps:ifpe,jfps:jfpe)
1773 fs_fire_ROSdt(ifps:ifpe, jfps:jfpe) = grid%burnt_area_dt(ifps:ifpe,jfps:jfpe)
1774 ! Accumulation prioritises firebrands from gridpts mostly burnt vs recently ignited
1776 !-------------------------------------------------------------------------------
1777 ! Updated fuel moisture from WRF-fire
1778 !-------------------------------------------------------------------------------
1780 fmcg_2df(ifps:ifpe,jfps:jfpe) = grid%fmc_g(ifps:ifpe,jfps:jfpe)
1782 !===============================================================================
1783 ! Firebrand Generation - Release Firebrands - Part 1
1784 !===============================================================================
1786 !-------------------------------------------------------------------------------
1787 ! When a new release is due, calculate fraction of release by fuel category
1788 !-------------------------------------------------------------------------------
1790 IF (fs_last_gen_dt >= fs_gen_dt) THEN
1792 firepts = COUNT(fs_fire_ROSdt(ifps:ifpe,jfps:jfpe) > ZERO_dp)
1793 Ngrdpts = 0 ! Number of active fire points of valid fuel categories
1795 LimPot = 0.000001_dp ! 1e-6
1796 MeanPot = ZERO_dp
1798 IF (firepts > 0) THEN
1800 !WRITE (msg,'(1(i6,1x))') firepts
1801 !CALL wrf_debug (wrfdbg, 'SPFire_rel firepts: '//msg)
1803 !-------------------------------------------------------------------------------
1804 ! Get fuel category from WRF-FIRE and calculate release potential
1805 !-------------------------------------------------------------------------------
1807 ALLOCATE(fcat_fr(ifps:ifpe,jfps:jfpe), fcat_factor(ifps:ifpe,jfps:jfpe), RelPot(ifps:ifpe,jfps:jfpe))!, rel_fcat(ifps:ifpe,jfps:jfpe))
1808 fcat_fr(:,:) = 0
1809 RelPot(:,:) = ZERO_dp
1810 fcat_factor(:,:) = ZERO_dp
1812 fcat_fr(ifps:ifpe,jfps:jfpe) = RESHAPE(get_fuel_cat(crosswalk=grid%fuel_crosswalk, &
1813 cat=grid%nfuel_cat(ifps:ifpe,jfps:jfpe)), &
1814 SHAPE(fcat_fr))
1816 fcat_factor(ifps:ifpe,jfps:jfpe)= firebrand_gen_factor(fcat_fr(ifps:ifpe,jfps:jfpe)) !RESHAPE, SHAPE(fcat_fr))
1818 RelPot(ifps:ifpe,jfps:jfpe) = RESHAPE(firebrand_gen_potential(&
1819 fuel_fgi=grid%fgip(ifps:ifpe,jfps:jfpe), &
1820 fuel_mcg=fmcg_2df(ifps:ifpe,jfps:jfpe), &
1821 !fuel_mcg=grid%fmcg2df(ifps:ifpe,jfps:jfpe), &
1822 factor=fcat_factor(ifps:ifpe,jfps:jfpe),&
1823 fire_rate=fs_fire_ROSdt(ifps:ifpe,jfps:jfpe)), &
1824 SHAPE(fcat_fr))
1826 ! Calculate Release Potential for each burning fuel cat
1827 !-------------------------------------------------------------------------------
1829 Ngrdpts= COUNT(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)! RelPot do not count cat with factor of 0
1831 ! In Large Fires
1832 ! Try to increase minimum LimPot to remove very low values of Release Potential
1833 !-------------------------------------------------------------------------------
1834 IF (Ngrdpts > fs_gen_lim) THEN
1835 LimPot = 10.0_dp * LimPot ! 1e-5
1836 Ngrdpts= COUNT(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)
1838 IF (Ngrdpts > fs_gen_lim) THEN
1839 LimPot = 10.0_dp * LimPot ! 1e-4
1840 Ngrdpts= COUNT(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)
1842 ENDIF
1843 ENDIF
1844 ! WRITE (msg,'(1(i8,1x), 1(f14.8,1x))') Ngrdpts, LimPot
1845 ! CALL wrf_debug (wrfdbg, 'SPFire_rel_lim Ngrdpts, LimP: '//msg)
1848 WRITE (msg,'(1(i6,1x))') Ngrdpts
1849 CALL wrf_debug (wrfdbg, 'SPFire_rel Ngrdpts: '//msg)
1851 ! Set a threshold for Potentials - use only highest values
1852 !-------------------------------------------------------------------------------
1853 IF (Ngrdpts > fs_gen_lim) THEN
1855 ! Find threshold value to select fs_gen_lim points
1856 !-------------------------------------------------------------------------------
1858 ALLOCATE(tmp1(Ngrdpts))
1859 tmp1(1:Ngrdpts) = PACK(RelPot(ifps:ifpe,jfps:jfpe), &
1860 mask=(RelPot(ifps:ifpe,jfps:jfpe) > LimPot))
1862 ! Find a limit value that yields an array size close to fs_gen_lim
1863 LimPot = order_val(tmp1, ORD=fs_gen_lim)
1864 !maxpot = MAXVAL(tmp1)
1865 DEALLOCATE(tmp1)
1867 ! WRITE (msg,'(3(i8,1x), 2(f14.8,1x))') Ngrdpts, COUNT(RelPot > LimPot), SIZE(tmp1), LimPot, maxpot
1868 ! CALL wrf_debug (wrfdbg, 'SPFire_rel_lim Ngrdpts, Cnt>Lim, LimP, max: '//msg)
1870 ! Find the median rank among remaining points
1871 !-------------------------------------------------------------------------------
1873 ALLOCATE(tmp1(COUNT(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)))
1874 tmp1 = PACK(RelPot(ifps:ifpe,jfps:jfpe), &
1875 mask=(RelPot(ifps:ifpe,jfps:jfpe) > LimPot))
1877 MeanPot = order_val(tmp1, ORD=INT(REAL(SIZE(tmp1))/2.0))
1878 !maxpot = MAXVAL(tmp1)
1879 DEALLOCATE(tmp1)
1881 ! WRITE (msg,'(4(i8,1x), 3(f14.8,1x))') Ngrdpts, &
1882 ! COUNT(RelPot > MeanPot), SIZE(tmp1), INT(REAL(SIZE(tmp1))/2.0), &
1883 ! MeanPot, MINVAL(tmp1), maxpot
1884 ! CALL wrf_debug (wrfdbg, 'SPFire_rel_lim Ngrdpts, Cnt>Mean, MPot, min, max: '//msg)
1886 Ngrdpts= COUNT(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)
1888 !-------------------------------------------------------------------------------
1890 ELSEIF (Ngrdpts > 0) THEN
1892 ! Find the median rank among points
1893 !-------------------------------------------------------------------------------
1894 MeanPot = order_val(PACK(RelPot(ifps:ifpe,jfps:jfpe), &
1895 mask=(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)), &
1896 ORD=INT(dp05*Ngrdpts))
1898 ! WRITE (msg,'(2(i8,1x), 2(f14.8,1x))') Ngrdpts, COUNT(RelPot > MeanPot), MeanPot, LimPot
1899 ! CALL wrf_debug (wrfdbg, 'SPFire_rel_lim Ngrdpts, Cnt>Mean, MPot, LimP: '//msg)
1901 ENDIF
1903 ! Prepare for MPI
1904 !-------------------------------------------------------------------------------
1906 IF (Ngrdpts == 0) firepts = 0
1908 !-------------------------------------------------------------------------------
1910 ENDIF ! (firepts > 0)
1911 !-------------------------------------------------------------------------------
1913 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
1914 IF (.NOT. (IamMaster)) THEN
1916 ! nbr > 0 indicates there's fire on tile and RelPot > 0, not the array size to receive in next mpi comm
1917 ! if nbr > 0, next array will be one (array: MeanPot)
1918 CALL fs_mpi_sendbuff1_int(sendto=MasterId, nbr=Ngrdpts)
1920 IF (Ngrdpts > 0) THEN
1921 CALL fs_mpi_sendbuff1_real(sendto=MasterId, nbr=MeanPot)
1922 ! WRITE (msg,'(1(i6,1x), (f14.6,1x))') Ngrdpts, MeanPot
1923 ! CALL wrf_debug (wrfdbg, 'SPFire_rel mpi_send Npts,Meanpot: '//msg)
1924 ENDIF
1925 ENDIF
1926 #endif
1928 !===============================================================================
1929 ! Firebrand Generation - Release Firebrands - Part 2
1930 !===============================================================================
1932 !-------------------------------------------------------------------------------
1933 ! MPI Master: Calculate Release Field (fire gridpt cnt by fuel category)
1934 !-------------------------------------------------------------------------------
1936 IF (IamMaster) THEN
1938 CALL wrf_debug (wrfdbg, 'SPFire_rel_master --------------------------------------------------- ')
1939 !-------------------------------------------------------------------------------
1940 ! Declare arrays to receive data
1942 Npts_total = 0
1943 PotThr = ZERO_dp
1944 TotPot = ZERO_dp
1945 rel_ratio = ZERO_dp
1946 cdf_coeff = ZERO_dp
1949 ALLOCATE(np_Ngrdpts(ntiles)) ! Number of fuel cat actively burning in each mpiproc
1950 ALLOCATE(np_MeanPot(ntiles))
1951 np_MeanPot(:) = ZERO_dp
1952 np_Ngrdpts(:) = 0
1954 np_Ngrdpts(1) = Ngrdpts ! Array index for Master is 1
1955 np_MeanPot(1) = MeanPot
1957 !-------------------------------------------------------------------------------
1958 ! Receive release potential data from procs
1960 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
1962 DO kk=2, mpiprocs ! kk=1 is Master, MPI-IDs begin at 0
1963 np_Ngrdpts(kk) = fs_mpi_recvbuff1_int(fromid=kk-1) ! Comm receive: Number of pts with RelPot>0
1965 IF (np_Ngrdpts(kk) > 0) THEN
1966 np_MeanPot(kk) = fs_mpi_recvbuff1_real(fromid=kk-1)
1967 ! WRITE(msg, '(3(i8,1x), 2(f14.6,1x))') my_id, kk-1, np_Ngrdpts(kk), np_Meanpot(kk), TotPot
1968 ! CALL wrf_debug (wrfdbg, 'SPFire_rel master proc Ngrdpts, MeanPot, TotPot: '//msg)
1969 ENDIF
1970 ENDDO
1971 #endif
1973 DO kk=1, ntiles ! at least 1 for serial compile
1974 TotPot = np_MeanPot(kk) * REAL(np_Ngrdpts(kk)) + TotPot
1975 ENDDO
1977 !-------------------------------------------------------------------------------
1978 ! Communicate release sources between Master and procs
1980 Npts_total = SUM( np_Ngrdpts )
1981 IF ( Npts_total > 0) THEN
1983 TotPot = TotPot/(REAL(Npts_total))
1985 !-------------------------------------------------------------------------------
1986 ! Calculate potential threshold for grid pts
1988 ! All elements fit in the array
1989 IF (Npts_total <= fs_gen_lim) THEN
1990 PotThr = ZERO_dp
1992 ! Estimate threshold with a linear regression
1993 ELSE ! IF (Npts_total > fs_gen_lim)
1994 rel_ratio = REAL(fs_gen_lim)/(REAL(Npts_total))
1995 cdf_coeff = (dp1 - rel_ratio)/rel_ratio
1996 PotThr = cdf_coeff * TotPot
1997 ENDIF
1999 ! WRITE(msg, '(2(i8,1x), 4(f12.6,1x))') fs_gen_lim, Npts_total, PotThr, rel_ratio, cdf_coeff, TotPot
2000 ! CALL wrf_debug (wrfdbg, 'SPFire_rel_master max, Npts, PotThr: '//msg)
2002 !-------------------------------------------------------------------------------
2003 ! Send dat to procs
2005 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2007 DO kk=2, mpiprocs ! kk=1 is Master
2008 IF (np_Ngrdpts(kk) > 0) THEN
2009 CALL fs_mpi_sendbuff1_real(sendto=kk-1, nbr=PotThr)
2010 ! WRITE(msg, '(1(i8,1x), 1(f12.6,1x))') kk-1, PotThr
2011 ! CALL wrf_debug (wrfdbg, 'SPFire_rel_master mpi_send PotThr: '//msg)
2012 ENDIF
2013 ENDDO
2014 #endif
2015 ENDIF ! npts_total > 0
2016 !-------------------------------------------------------------------------------
2018 ENDIF ! IamMaster
2020 !-------------------------------------------------------------------------------
2022 !===============================================================================
2023 ! Firebrand Generation - Release Firebrands - Part 3
2024 !===============================================================================
2026 IF (Ngrdpts > 0) THEN
2028 !-------------------------------------------------------------------------------
2029 ! Receive Potential Threshold from Master
2031 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2032 IF (.NOT.(IamMaster)) THEN
2033 PotThr = ZERO_dp
2034 PotThr = fs_mpi_recvbuff1_real(fromid=MasterId)
2035 ! WRITE(msg, '((i4,1x), (f14.6,1x))') my_id, PotThr
2036 ! CALL wrf_debug (wrfdbg, 'SPFire_rel receive potential threshold: '//msg)
2037 ENDIF
2038 #endif
2040 !-------------------------------------------------------------------------------
2041 ! Assign release from the various fuel cat to corresponding gridpts
2043 !CALL wrf_debug (wrfdbg, 'SPFire_rel2fr relpts>0 --------------------------------------------- ')
2045 ALLOCATE(n_rel_fr2d(ifps:ifpe, jfps:jfpe))
2046 n_rel_fr2d(ifps:ifpe, jfps:jfpe) = 0
2048 IF (PotThr > ZERO_dp) &
2049 n_rel_fr2d(ifps:ifpe,jfps:jfpe) = n_rel_fr2d(ifps:ifpe,jfps:jfpe) + &
2050 UNPACK(SPREAD(1, DIM=1, NCOPIES=Ngrdpts), &
2051 mask=(RelPot(ifps:ifpe,jfps:jfpe) >= PotThr), FIELD=0)
2053 IF (PotThr == ZERO_dp) &
2054 n_rel_fr2d(ifps:ifpe,jfps:jfpe) = n_rel_fr2d(ifps:ifpe,jfps:jfpe) + &
2055 UNPACK(SPREAD(1, DIM=1, NCOPIES=Ngrdpts), &
2056 mask=(fs_fire_ROSdt(ifps:ifpe,jfps:jfpe) > ZERO_dp), FIELD=0)
2058 !-------------------------------------------------------------------------------
2059 ! update relpts to represent value assigned to gridpts
2061 relpts = COUNT(n_rel_fr2d(ifps:ifpe,jfps:jfpe)>0) ! Limits to 1 release per refined gridpt, at various heights
2063 WRITE(msg, '(1(i8,1x))') relpts
2064 CALL wrf_debug (wrfdbg, 'SPFire_rel2fr relpts: '//msg)
2066 !-------------------------------------------------------------------------------
2067 ! Convert release locations to met grid
2069 ALLOCATE(release_i(relpts*fs_gen_levels), &
2070 release_j(relpts*fs_gen_levels), &
2071 release_k(relpts*fs_gen_levels), &
2072 release_n(relpts))
2074 release_i(:) = 0.0_dp
2075 release_j(:) = 0.0_dp
2076 release_k(:) = 0.0_dp
2077 release_n(:) = 0
2079 ! Does not allow more than one release from the same refined gridpt
2080 CALL releaseijk2atm(nij_2d=n_rel_fr2d(ifps:ifpe,jfps:jfpe), &
2081 fcat=fcat_fr(ifps:ifpe,jfps:jfpe), &
2082 maxhgt_usr = firebrand_gen_maxhgt, &
2083 sr_x=grid%sr_x, sr_y=grid%sr_y, &
2084 pi=release_i(1:relpts), &
2085 pj=release_j(1:relpts), &
2086 pk=release_k(1:relpts), & ! max height
2087 nij=release_n)
2090 ! Adjust position to corresponding refined grid center
2091 release_i(1:relpts) = release_i(1:relpts) + dp05/grid%sr_x
2092 release_j(1:relpts) = release_j(1:relpts) + dp05/grid%sr_y
2094 DO kk=1,relpts ! track release points between history outputs
2096 fs_gen_inst(INT(release_i(kk)),INT(release_j(kk))) = &
2097 fs_gen_inst(INT(release_i(kk)),INT(release_j(kk))) + 1
2099 ENDDO
2101 CALL wrf_debug (wrfdbg, 'SPFire_rel2fr prep_release_hgt ---------------------------------------------')
2102 ALLOCATE(release_prop(relpts*fs_gen_levels))
2104 ! Create derived type
2105 release_prop(1:relpts) = firebrand_property([firebrand_gen_prop_diam, &
2106 firebrand_gen_prop_effd, &
2107 firebrand_gen_prop_temp, &
2108 firebrand_gen_prop_tvel])
2110 levrand_seed = 0 ! means no random levels
2111 IF (firebrand_gen_levrand) levrand_seed = cf%fs_firebrand_gen_levrand_seed + 1 ! in case user sets it to 0
2112 ! include more releases to span over various heights
2113 CALL prep_release_hgt( &
2114 release_i = release_i, &
2115 release_j = release_j, &
2116 release_k = release_k, & ! max height
2117 release_n = release_n, &
2118 landhgt = firebrand_land_hgt, &
2119 release_prop=release_prop, &
2120 levrand_seed = levrand_seed * grid%itimestep) ! force seed to vary in time
2122 prior_br = active_br
2123 idmax = fs_gen_idmax
2125 CALL wrf_debug (wrfdbg, 'SPFire_driver_call_generate_firebrands ----------------------------------------- ')
2127 IF (active_br + COUNT( INT(release_i) > 0 ) > fs_array_maxsize) &
2128 CALL wrf_debug (wrfdbg, 'SPFire_driver_release: brand array is full, cannot release new brands')
2130 CALL generate_firebrands(&
2131 fs_p_id = fs_p_id, &
2132 fs_p_src = fs_p_src, &
2133 fs_p_dt = fs_p_dt, &
2134 fs_p_z = fs_p_z, &
2135 fs_p_x = fs_p_x, &
2136 fs_p_y = fs_p_y, &
2137 fs_p_prop = fs_p_prop, &
2138 fs_gen_idmax= fs_gen_idmax, &
2139 myprocid = myprocid, &
2140 release_prop= release_prop, &
2141 release_i = release_i, &
2142 release_j = release_j, &
2143 release_k = release_k, &
2144 active_br = active_br)
2146 !-------------------------------------------------------------------------------
2147 ! Track released firebrands
2148 !-------------------------------------------------------------------------------
2150 WRITE (msg,'(i8)') active_br
2151 CALL wrf_debug (wrfdbg, 'SPFire_driver: Active brands: '// msg )
2153 IF (cf%trackember) THEN
2155 DO kk=prior_br+1,active_br
2156 WRITE (msg,'(3(i8,1x),i10,1x,8(f12.6,1x))') &
2157 kk, fs_p_id(kk), fs_p_dt(kk),fs_p_src(kk), &
2158 fs_p_x(kk), fs_p_y(kk), fs_p_z(kk), &
2159 fs_p_mass(kk), fs_p_diam(kk), fs_p_effd(kk), &
2160 fs_p_temp(kk), fs_p_tvel(kk)
2161 CALL wrf_debug (0, 'SPFire p,id,dt,src,x,y,z,m,d,e,t,v: '//msg)
2162 ENDDO
2163 ENDIF
2164 !-------------------------------------------------------------------------------
2166 IF (active_br /= COUNT( fs_p_id > 0 ) ) CALL wrf_error_fatal('SPFire_driver: Active brands do not match!')
2168 ENDIF ! Ngrdpts > 0
2170 fs_fire_ROSdt(:, :) = ZERO_dp
2171 relpts = 0
2172 fs_last_gen_dt = 0 ! Reset release timestep counter
2173 IF( ALLOCATED(fcat_fr)) DEALLOCATE(fcat_fr)
2174 IF( ALLOCATED(fcat_factor)) DEALLOCATE(fcat_factor)
2175 IF( ALLOCATED(RelPot)) DEALLOCATE(RelPot)
2176 IF( ALLOCATED(np_Ngrdpts)) DEALLOCATE(np_Ngrdpts)
2177 IF( ALLOCATED(np_MeanPot)) DEALLOCATE(np_MeanPot)
2178 IF( ALLOCATED(n_rel_fr2d)) DEALLOCATE(n_rel_fr2d)
2179 IF( ALLOCATED(release_prop)) DEALLOCATE(release_prop)
2181 ENDIF ! IF (fs_last_gen_dt >= fs_gen_dt) THEN
2183 !===============================================================================
2184 ! End of Releases
2185 !===============================================================================
2187 !-------------------------------------------------------------------------------
2188 ! Nothing to transport: MPI Wait
2189 !-------------------------------------------------------------------------------
2191 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2192 IF ((active_br == 0)) THEN
2194 CALL fs_mpi_nothing2send(sendto=left_id)
2195 CALL fs_mpi_nothing2send(sendto=right_id)
2196 CALL fs_mpi_nothing2send(sendto=up_id)
2197 CALL fs_mpi_nothing2send(sendto=down_id)
2199 CALL fs_mpi_nothing2send(sendto=upleft_id)
2200 CALL fs_mpi_nothing2send(sendto=upright_id)
2201 CALL fs_mpi_nothing2send(sendto=downleft_id)
2202 CALL fs_mpi_nothing2send(sendto=downright_id)
2203 !CALL wrf_debug (wrfdbg, 'SPFire_driver: no fires or active brands in patch. MPI signal sent.')
2205 ENDIF
2206 #endif
2207 !===============================================================================
2208 !===============================================================================
2212 !===============================================================================
2213 ! Transport active firebrands
2214 !===============================================================================
2216 IF (active_br > 0) THEN
2218 !-------------------------------------------------------------------------------
2219 ! Fetch wrf variables for advection
2220 !-------------------------------------------------------------------------------
2222 ! Use arrays computed for the physics scheme on half (_phy) and full (8w) levels
2223 ! will skip calc at model top - particle will not live during journey (module_big_step_utilities_em.F)
2225 !-------------------------------------------------------------------------------
2226 ! Local met for fb physics: based on variables calculated in module_big_step_utilities_em.F
2227 !-------------------------------------------------------------------------------
2229 !HALO_EM_E_5:48: u_1,u_2,v_1,v_2,w_1,w_2,t_1,t_2,ph_1,ph_2,tke_1,tke_2,;4:mu_1,mu_2 moist
2230 !-------------------------------------------------------------------------------
2231 ! P, Th, Rho, hgt at full levels
2233 ! Direct assignments from existing Halos:
2235 ! Half Levels
2236 k_end = MIN( ke, kde-1 )
2237 ks = 1
2239 p_phy(ims:, ks:, jms:) = ZERO_dp
2240 th_phy(ims:, ks:, jms:) = ZERO_dp
2243 p_phy(ims:, ks:k_end, jms:) = grid%p(ims:ime,ks:k_end,jms:jme) + grid%pb(ims:ime,ks:k_end,jms:jme)
2244 th_phy(ims:, ks:, jms:) = grid%t_2(ims:ime,ks:k_end,jms:jme) + t0
2246 !hypsometric_opt 2 (default) computes pressure in the model by using an alternative method
2247 ! Temperature is moist theta by default (v4+, https://github.com/wrf-model/WRF/releases)
2248 IF ( ( grid%use_theta_m == 1 ) .AND. (P_Qv >= PARAM_FIRST_SCALAR) ) &
2249 th_phy = th_phy/(dp1 + Rv/Rd * grid%moist(ims:ime,ks:k_end,jms:jme,P_QV))
2251 ! Full Levels
2252 z_at_w(ims:,ks:,jms:)= (grid%phb(ims:ime,ks:ke,jms:jme) + grid%ph_2(ims:ime,ks:ke,jms:jme))/grav ! ASL
2254 u(ims:, ks:, jms:) = grid%u_2(ims:ime,ks:ke,jms:jme)
2255 v(ims:, ks:, jms:) = grid%v_2(ims:ime,ks:ke,jms:jme)
2256 w(ims:, ks:, jms:) = grid%w_2(ims:ime,ks:ke,jms:jme)
2258 msft(ims:, jms:) = grid%msftx(ims:ime, jms:jme)
2260 znw = grid%znw
2261 DT = grid%dt
2262 rdx = grid%rdx ! inverse grid length [m]
2263 rdy = grid%rdy
2265 !-------------------------------------------------------------------------------
2266 ! Calculated variables
2268 p_hyd(ims:, ks:, jms:) = ZERO_dp
2269 dz8w(ims:, ks:, jms:) = ZERO_dp
2270 th8w(ims:, ks:, jms:) = ZERO_dp
2271 p8w(ims:, ks:, jms:) = ZERO_dp
2272 qtot(ims:, ks:, jms:) = ZERO_dp
2273 rho(ims:, ks:, jms:) = ZERO_dp
2274 z(ims:, ks:, jms:) = ZERO_dp
2276 dz8w(ims:,ks:ke-1,jms:) = z_at_w(:,2:,:)-z_at_w(:,:ke-1,:) ! z_at_w: Height above the ground at interfaces
2277 z(ims:,ks:ke-1,jms:) = dp05 * (z_at_w(:,2:,:) + z_at_w(:,:ke-1,:) )
2279 qtot(ims:, ks:k_end, jms:) = SUM(grid%moist(ims:ime,ks:k_end,jms:jme,PARAM_FIRST_SCALAR:num_moist),DIM=4)
2280 p_hyd(ims:, ke, jms:) = grid%p_top
2281 DO kk = ke-1, ks, -1
2282 p_hyd(:,kk,:) = p_hyd(:,kk+1,:) - &
2283 (dp1 + qtot(:,kk,:)) * (grid%c1h(kk) * grid%muts(:,:) + grid%c2h(kk)) * grid%dnw(kk)
2284 ENDDO
2286 rho(ims:,ks:k_end,jms:) = (dp1/&
2287 (grid%al(ims:ime,ks:k_end,jms:jme)+grid%alb(ims:ime,ks:k_end,jms:jme)))&
2288 *(dp1+grid%moist(ims:ime,ks:k_end,jms:jme,P_QV))
2289 DO kk = 2, k_end
2290 th8w(ims:,kk,jms:) = grid%fnm(kk) * th_phy(:,kk,:) + th_phy(:,kk-1,:) * grid%fnp(kk)
2291 p8w(ims:,kk,jms:) = grid%fnm(kk) * p_phy(:,kk,:) + p_phy(:,kk-1,:) * grid%fnp(kk)
2292 ENDDO
2294 w1 = (z(:,1,:) - z(:,2,:))
2295 WHERE(NINT(w1*100.0) == 0 ) w1 = dp1
2296 w1 = (z_at_w(:,1,:) - z(:,2,:)) / w1
2298 th8w(:,1,:) = w1 * th_phy(:,1,:) + (dp1 - w1) * th_phy(:,2,:) ! interp at bottom only - ignore top
2299 p8w(:,1,:) = w1 * p_phy(:,1,:) + (dp1 - w1) * p_phy(:,2,:) ! interp at bottom only - ignore top
2301 DO kk = 1, ke
2302 z_at_w(:,kk,:)= z_at_w(:,kk,:) - z_at_w(:,1,:)
2303 ENDDO
2305 !-------------------------------------------------------------------------------
2306 !-------------------------------------------------------------------------------
2308 !-------------------------------------------------------------------------------
2309 ! Advect particle
2310 !-------------------------------------------------------------------------------
2312 CALL wrf_debug (wrfdbg, 'SPFire_driver begin transport ---------------------------------------------------- ')
2314 ALLOCATE(xout(active_br), &
2315 yout(active_br), &
2316 zout(active_br))
2318 CALL advect_xyz_m(grid=grid, &
2319 xp= fs_p_x(:active_br), &
2320 yp= fs_p_y(:active_br), &
2321 hgt= fs_p_z(:active_br), & ! height[m] is the fixed position reference between wrf time steps
2322 dtp= fs_p_dt(:active_br),&
2323 idp= fs_p_id(:active_br),&
2324 znw= znw, &
2325 mf = rdx*msft, &
2326 z_at_w=z_at_w, &
2327 u= u, &
2328 v= v, &
2329 w= w, &
2330 dt=dt, &
2331 ims= ims, &
2332 jms= jms, &
2333 kms= ks, &
2334 phyd= p_hyd, & ! air pressure (Pa)
2335 thet= th8w, & ! temperature (K)
2336 rho = rho, & ! density (kg/m3)
2337 xout= xout, &
2338 yout= yout, &
2339 zout= zout, &
2340 fs_p_prop=fs_p_prop(:active_br), &
2341 land_hgt = firebrand_land_hgt, &
2342 start_mom3d_dt = cf%fs_firebrand_gen_mom3d_dt, &
2343 msg = msg)
2346 IF (LEN(TRIM(msg)) > 1) &
2347 CALL wrf_debug (wrfdbg, 'SPFire_transp: '//msg)
2349 !-------------------------------------------------------------------------------
2350 ! Update array for output before resetting deposited indices to zero
2351 !-------------------------------------------------------------------------------
2353 fs_p_x(:active_br) = xout
2354 fs_p_y(:active_br) = yout
2355 fs_p_z(:active_br) = zout
2356 fs_p_dt(:active_br)= fs_p_dt(:active_br) +1
2358 fs_p_mass(:active_br) = fs_p_prop(:active_br)%p_mass
2359 fs_p_diam(:active_br) = fs_p_prop(:active_br)%p_diam
2360 fs_p_effd(:active_br) = fs_p_prop(:active_br)%p_effd
2361 fs_p_temp(:active_br) = fs_p_prop(:active_br)%p_temp
2362 fs_p_tvel(:active_br) = fs_p_prop(:active_br)%p_tvel
2364 CALL wrf_debug (wrfdbg, 'SPFire_driver end transport ---------------------------------------------------- ')
2366 !=============================================================
2367 ! Transport done :)
2368 !=============================================================
2370 !-------------------------------------------------------------------------------
2371 ! Track firebrands
2372 !-------------------------------------------------------------------------------
2374 WRITE (msg,'(2(i8,x))') active_br, COUNT( fs_p_id > 0 )
2375 CALL wrf_debug (wrfdbg, 'SPFire_driver: Active brands, id>0: '// msg )
2377 IF (cf%trackember) THEN
2378 DO kk=1,active_br
2379 WRITE (msg,'(3(i8,1x),i10,1x,8(f12.6,1x))') &
2380 kk, fs_p_id(kk), fs_p_dt(kk),fs_p_src(kk), &
2381 fs_p_x(kk), fs_p_y(kk), fs_p_z(kk), &
2382 fs_p_mass(kk), fs_p_diam(kk), fs_p_effd(kk), &
2383 fs_p_temp(kk), fs_p_tvel(kk)
2384 CALL wrf_debug (0, 'SPFire p,id,dt,src,x,y,z,m,d,e,t,v: '//msg)
2385 ENDDO
2386 ENDIF
2387 !-------------------------------------------------------------------------------
2389 !===============================================================================
2390 ! Remove particles
2391 !===============================================================================
2393 !-------------------------------------------------------------------------------
2394 ! Set NaN properties for particle removal
2395 !-------------------------------------------------------------------------------
2396 sparse_mask = .FALSE.
2398 WHERE (IEEE_IS_NAN(fs_p_mass) .OR. IEEE_IS_NAN(fs_p_diam) .OR. &
2399 IEEE_IS_NAN(fs_p_effd) .OR. IEEE_IS_NAN(fs_p_temp) .OR. &
2400 IEEE_IS_NAN(fs_p_tvel)) fs_p_effd = ZERO_dp ! ediam=zero will be removed in remove_br
2403 !===============================================================================
2404 ! Remove particles - from burnout or out of domain
2405 !===============================================================================
2406 nresets = 0
2407 CALL remove_br(& !active_br= active_br, &
2408 fs_p_id = fs_p_id, &
2409 fs_p_x = fs_p_x, &
2410 fs_p_y = fs_p_y, &
2411 fs_p_z = fs_p_z, &
2412 fs_p_dt = fs_p_dt, &
2413 fs_p_effd = fs_p_effd, &
2414 cnt = nresets, &
2415 max_life_dt= firebrand_max_life_dt, &
2416 land_hgt = firebrand_land_hgt)
2417 WRITE (msg,'(i12)') nresets
2418 CALL wrf_debug (wrfdbg, 'SPFire_driver remove br: '//msg)
2421 !===============================================================================
2422 ! Spotfires: Deposited Brands
2423 !===============================================================================
2425 dt_sum_landed(ims:, jms:) = ZERO_dp
2426 CALL deposit_br(& !active_br= active_br, &
2427 fs_p_id = fs_p_id, &
2428 fs_p_x = fs_p_x, &
2429 fs_p_y = fs_p_y, &
2430 fs_p_z = fs_p_z, &
2431 sum_fbrand = dt_sum_landed,&
2432 land_hgt = firebrand_land_hgt)
2434 ! The mask is based on floor(fs_p_x,y) so deposits should not exist at ie or je
2435 fs_count_landed_all(is:ie, js:je) = fs_count_landed_all(is:ie, js:je) + dt_sum_landed(is:ie, js:je)
2436 fs_count_landed_hist(is:ie, js:je) = fs_count_landed_hist(is:ie, js:je) + dt_sum_landed(is:ie, js:je)
2437 WRITE (msg,'(2(i12,x))') COUNT(dt_sum_landed(is:ie, js:je) > 0), INT(SUM(dt_sum_landed(is:ie, js:je)))
2438 CALL wrf_debug (wrfdbg, 'SPFire_driver deposit br: '//msg)
2441 !===============================================================================
2442 !===============================================================================
2444 !******************************************************************
2445 !******************************************************************
2446 !* *
2447 !* MPI Send firebrands *
2448 !* *
2449 !******************************************************************
2450 !******************************************************************
2453 !-------------------------------------------------------------------------------
2454 !-------------------------------------------------------------------------------
2456 ! Re-PACK before incoming MPI
2458 !===============================================================================
2459 ! MPI: Send Brands (out of this patch)
2460 !===============================================================================
2462 ! CALL wrf_debug (wrfdbg, 'SPFire_driver mpi transfer particles')
2464 ! Send / Receive arrays:
2465 ! 3 REAL arrays, 2 INTEGER arrays: fs_p_x, fs_p_y, fs_p_z, fs_p_id, fs_p_dt:
2467 ! |____|____|____|
2468 ! | | U | |
2469 ! |____|____|____|je
2470 ! | L | Me | R |
2471 ! |____|____|____|js
2472 ! | | D | |
2473 ! |____|____|____|
2474 ! is ie
2476 ! Diagonal???
2478 !-------------------------------------------------------------------------------
2479 ! Send any particle?
2480 !-------------------------------------------------------------------------------
2482 sparse_mask = .FALSE.
2483 sparse_mask = (fs_p_id > 0 .AND. &
2484 (FLOOR(fs_p_x) < is .OR. & ! LEFT
2485 FLOOR(fs_p_x) > ie .OR. & ! RIGHT
2486 FLOOR(fs_p_y) > je .OR. & ! UP
2487 FLOOR(fs_p_y) < js )) ! DONW
2490 WRITE (msg,'(i6)') COUNT(sparse_mask)
2491 CALL wrf_debug (wrfdbg, 'SPFire_driver mpi send away: '//msg)
2493 !-------------------------------------------------------------------------------
2494 ! interrupt before MPI
2495 ! IF (COUNT(sparse_mask) > 0) CALL wrf_error_fatal( "MPI called")
2497 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2498 CALL fs_mpi_send2neighbors(task_id=task_id,&
2499 mask=sparse_mask,&
2500 p_x = fs_p_x, &
2501 p_y = fs_p_y, &
2502 p_z = fs_p_z, &
2503 p_id = fs_p_id, &
2504 p_src = fs_p_src, &
2505 p_dt = fs_p_dt, &
2506 fs_p_m = fs_p_mass, &
2507 fs_p_d = fs_p_diam, &
2508 fs_p_e = fs_p_effd, &
2509 fs_p_t = fs_p_temp, &
2510 fs_p_v = fs_p_tvel)
2511 #endif
2512 !-------------------------------------------------------------------------------
2513 ! Reset
2514 !-------------------------------------------------------------------------------
2515 ! CALL wrf_debug (wrfdbg, 'SPFire_driver reset particles sent away')
2517 ! In serial compilation, firebrands outside the domain limits will be removed
2518 WHERE(sparse_mask) fs_p_id= 0
2520 !===============================================================================
2521 ! Last step for IF active_br > 0
2522 !===============================================================================
2524 ! WRITE (msg,'(2(i8,1x))') active_br, COUNT( fs_p_id > 0 )
2525 ! CALL wrf_debug (wrfdbg, 'SPFire_driver pack BEFORE: active_br >>> '// msg)
2527 !-------------------------------------------------------------------------------
2528 ! Repack all zeros: PACK where mask is TRUE
2529 !-------------------------------------------------------------------------------
2530 sparse_mask = .FALSE.
2531 sparse_mask = (fs_p_id>0)
2532 fs_p_x = PACK(fs_p_x, sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2533 fs_p_y = PACK(fs_p_y, sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2534 fs_p_z = PACK(fs_p_z, sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2535 fs_p_id = PACK(fs_p_id, sparse_mask, SPREAD(0,1,fs_array_maxsize) )
2536 fs_p_src = PACK(fs_p_src, sparse_mask, SPREAD(0,1,fs_array_maxsize) )
2537 fs_p_dt = PACK(fs_p_dt, sparse_mask, SPREAD(0,1,fs_array_maxsize) )
2538 fs_p_mass = PACK(fs_p_mass , sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2539 fs_p_diam = PACK(fs_p_diam , sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2540 fs_p_effd = PACK(fs_p_effd, sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2541 fs_p_temp = PACK(fs_p_temp , sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2542 fs_p_tvel = PACK(fs_p_tvel , sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2544 active_br = COUNT( fs_p_id > 0 )
2546 WRITE (msg,'(1(i8,1x))') active_br
2547 CALL wrf_debug (wrfdbg, 'SPFire_driver pack AFTER mpi_send2neighbors active_br >>> '// msg)
2550 ENDIF
2552 !===============================================================================
2553 ! Active_br work done ;)
2554 !===============================================================================
2555 IF (active_br /= COUNT( fs_p_id > 0 ) ) CALL wrf_error_fatal('SPFire_driver: Active brands do not match!')
2557 !===============================================================================
2558 !===============================================================================
2562 !===============================================================================
2563 ! Calculate Spotfire Likelihood at history output time
2564 !===============================================================================
2566 ! Registry
2567 ! fs_count_landed_all
2568 ! fs_count_landed_hist
2569 ! fs_frac_landed
2570 ! fs_spotting_lkhd
2571 ! fs_count_reset
2572 ! fs_fire_area
2574 ! store array of data for each proc: Master only
2575 ! np_nlanded [np] ! ALLOCATABLE
2576 ! np_landed_cnt [SUM(np_nlanded)] ! ALLOCATABLE
2577 ! np_landed_risk [SUM(np_nlanded)] ! ALLOCATABLE
2578 ! np_frac_landed [SUM(np_nlanded)] ! ALLOCATABLE
2579 ! np_lkhd [SUM(np_nlanded)] ! ALLOCATABLE
2580 !
2582 ! store array of local data
2583 ! landing_mask ! ALLOCATABLE
2584 ! landing_cnt ! ALLOCATABLE (mpi send)
2585 ! landed_risk ! ALLOCATABLE (mpi send)
2586 ! recv_spot_lkhd ! ALLOCATABLE (mpi recv)
2587 ! recv_frac_landed ! ALLOCATABLE (mpi recv)
2589 !===============================================================================
2590 ! Check history interval (use diag_flag defined in solve_em for microphysics)
2591 !===============================================================================
2593 IF ((hist_flag) .AND. (grid%itimestep > 1)) THEN
2595 fs_count_reset = .TRUE. ! true
2596 fs_frac_landed(ims:ime, jms:jme) = ZERO_dp
2597 fs_spotting_lkhd(ims:ime, jms:jme) = ZERO_dp
2599 ! Was there any firebrand deposit? *Note that fire may exist on a different tile
2600 !-------------------------------------------------------------------------------
2602 ALLOCATE(landing_mask(ims:ime, jms:jme)) ! Allocate memory but only mask this tile
2603 landing_mask(ims:ime, jms:jme) = .FALSE.
2604 landing_mask(is:ie, js:je) = (fs_count_landed_hist(is:ie, js:je) > 0)
2605 ndep = COUNT(landing_mask)
2607 ! convert fire_area from fire grid to atm grid
2608 !-------------------------------------------------------------------------------
2609 fs_fire_area(is:ie,js:je) = ZERO_dp
2610 fs_fire_area(is:ie, js:je) = fire2tile(fr_arr=grid%fire_area(ifps:ifpe,jfps:jfpe), &
2611 dsx=grid%sr_x,dsy=grid%sr_y)
2613 ! WRITE(msg, '(i8,1x,i12)') ndep, COUNT(fs_fire_area(is:ie,js:je)> ZERO_dp)
2614 ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd: deposit points, fire points: '//msg)
2616 IF (ndep > 0) THEN ! deposited anywhere
2618 landing_mask(ims:ime, jms:jme) = .FALSE.
2619 landing_mask(is:ie, js:je) = ((fs_count_landed_hist(is:ie, js:je) > 0) .AND. &
2620 (fs_fire_area(is:ie, js:je) == ZERO_dp) )
2621 ndep = COUNT(landing_mask)
2623 ! deposited on nonburning gridpt
2624 fs_landing_mask(is:ie, js:je) = UNPACK(SPREAD(1,DIM=1,NCOPIES=ndep), MASK=landing_mask(is:ie, js:je), FIELD=0)
2626 ! WRITE(msg, '(4(i8,1x))') COUNT(fs_count_landed_hist > ZERO_dp), ndep, COUNT(fs_landing_mask == 1), COUNT((fs_count_landed_hist > ZERO_dp).AND.(fs_fire_area > ZERO_dp))
2627 ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd deposits, non fire-point deposits, fire-point deposits: '//msg)
2629 ALLOCATE(tmp1(ndep))
2630 tmp1 = PACK(fs_count_landed_hist(is:ie, js:je), landing_mask(is:ie, js:je)) ! WHERE doest work on memory bounds
2631 fs_count_landed_hist(is:ie, js:je) = ZERO_dp
2632 fs_count_landed_hist(is:ie, js:je) = UNPACK(tmp1, MASK=landing_mask(is:ie, js:je), FIELD=ZERO_dp)
2634 ! WRITE(msg, '(2(i8,1x), 1(f14.6,1x))') COUNT(fs_count_landed_hist > ZERO_dp), ndep, MAXVAL(fs_count_landed_hist)
2635 ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd deposits, non fire-point deposits, fire-point deposits: '//msg)
2636 DEALLOCATE(tmp1)
2638 ! Were these firebrand deposited outside the fire area?
2639 !-------------------------------------------------------------------------------
2641 IF (ndep > 0) THEN
2643 ! Calculate fuel risk
2644 !-------------------------------------------------------------------------------
2646 ALLOCATE(fuel_spotting_risk_fr(ifps:ifpe,jfps:jfpe), &!fs_fuel_spotting_risk(ims:ime,jms:jme), &
2647 spotthreat_fr(ifps:ifpe,jfps:jfpe))
2648 fuel_spotting_risk_fr(:,:) = ZERO_dp
2649 fs_fuel_spotting_risk(:,:) = ZERO_dp
2650 spotthreat_fr(:,:) = ZERO_dp
2652 IF (.NOT. ALLOCATED(fcat_fr)) THEN
2653 ALLOCATE(fcat_fr(ifps:ifpe,jfps:jfpe))
2654 fcat_fr(:,:) = ZERO_dp
2655 fcat_fr(ifps:ifpe,jfps:jfpe) = RESHAPE(&
2656 get_fuel_cat(crosswalk=grid%fuel_crosswalk, &
2657 cat=grid%nfuel_cat(ifps:ifpe,jfps:jfpe)), &
2658 SHAPE(fuel_spotting_risk_fr))
2659 ENDIF
2660 !WRITE(msg, *) MAXVAL(fcat_fr), MINVAL(fcat_fr)
2661 ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd fcat max,min '//msg)
2663 spotthreat_fr(ifps:ifpe,jfps:jfpe)= spotting_threat_factor(fcat_fr(ifps:ifpe,jfps:jfpe)) !RESHAPE(SHAPE(fuel_spotting_risk_fr))
2664 ! get_fuel_cat(crosswalk=grid%fuel_crosswalk, &
2665 ! cat=grid%nfuel_cat(ifps:ifpe,jfps:jfpe))), &
2666 ! WRITE(msg, *) MAXVAL(spotthreat_fr), MINVAL(spotthreat_fr)
2667 ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd spotting_threat_factor max,min'//msg)
2669 fuel_spotting_risk_fr(ifps:ifpe,jfps:jfpe) = RESHAPE(&
2670 fuel_spotting_risk(&
2671 fuel_fgi=grid%fgip(ifps:ifpe,jfps:jfpe), &
2672 fuel_mcg=fmcg_2df(ifps:ifpe,jfps:jfpe), &
2673 !fuel_mcg=grid%fmcg2df(ifps:ifpe,jfps:jfpe), &
2674 factor=spotthreat_fr(ifps:ifpe,jfps:jfpe)), &
2675 SHAPE(fuel_spotting_risk_fr)) !, fuel_mcg=grid%fmcg2df)
2677 fs_fuel_spotting_risk(is:ie,js:je) = fire2tile(fr_arr=fuel_spotting_risk_fr(ifps:ifpe,jfps:jfpe), &
2678 dsx=grid%sr_x,dsy=grid%sr_y)
2679 DEALLOCATE(fcat_fr)
2681 !-------------------------------------------------------------------------------
2684 ! Prepare for MPI
2685 !-------------------------------------------------------------------------------
2687 ALLOCATE(landing_cnt(ndep), landed_risk(ndep), recv_spot_lkhd(ndep), recv_frac_landed(ndep))
2688 landing_cnt = PACK(fs_count_landed_hist(is:ie,js:je), landing_mask(is:ie,js:je))
2689 landed_risk = PACK(fs_fuel_spotting_risk(is:ie,js:je), landing_mask(is:ie,js:je))
2691 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2692 IF (.NOT. (IamMaster)) THEN
2694 CALL fs_mpi_sendbuffsize(sendto=MasterId, nbr=ndep) ! Comm send >0: BuffSz
2695 CALL fs_mpi_sendbuff_int(sendto=MasterId, bsz=ndep, buff=landing_cnt(1:ndep))
2696 CALL fs_mpi_sendbuff_real(sendto=MasterId, bsz=ndep, buff=landed_risk(1:ndep))
2698 ENDIF
2699 #endif
2700 ENDIF ! (fs_count_landed_hist > 0 .AND. fire_area == 0)
2701 ENDIF ! SUM(fs_count_landed_hist) > 0
2703 !-------------------------------------------------------------------------------
2705 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2706 IF ((ndep == 0) .AND. .NOT.(IamMaster)) &
2707 CALL fs_mpi_nothing2send(sendto=MasterId) ! Comm send 0: BuffSz
2708 #endif
2709 !-------------------------------------------------------------------------------
2710 ! Master: Calculate Likelihood Field
2711 !-------------------------------------------------------------------------------
2713 IF (IamMaster) THEN
2715 !-------------------------------------------------------------------------------
2716 ! Declare arrays to receive data
2718 ALLOCATE(np_nlanded(ntiles))
2719 np_nlanded(:) = 0
2720 np_nlanded(1) = ndep ! Array index for Master is 1
2722 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2723 DO kk=2, mpiprocs ! kk=1 is Master, MPI-IDs begin at 0
2724 np_nlanded(kk) = fs_mpi_recvbuffsize(fromid=kk-1) ! Comm receive: BuffSz
2725 ! WRITE(msg, '(2(i8,1x))') kk,np_nlanded(kk)
2726 ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd proc deposit: '//msg)
2727 ENDDO
2728 #endif
2729 !-------------------------------------------------------------------------------
2730 ! Communicate deposits between Master and procs
2732 ndep_total = SUM(np_nlanded)
2734 IF ( ndep_total > 0) THEN
2735 ALLOCATE(np_landed_cnt(ndep_total), np_landed_risk(ndep_total))
2736 IF (ndep > 0) THEN
2737 np_landed_cnt(1:ndep) = landing_cnt
2738 np_landed_risk(1:ndep)= landed_risk
2739 ELSE
2740 np_landed_cnt(1) = ZERO_dp ! I don't think I need this
2741 np_landed_risk(1)= ZERO_dp
2742 ENDIF
2744 !-------------------------------------------------------------------------------
2745 ! Receive deposit data from procs
2747 kk1 = ndep + 1 ! index start after master's
2749 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2750 DO kk=2, mpiprocs ! kk=1 is Master
2752 IF (np_nlanded(kk) > 0) THEN
2754 kk2 = kk1 + np_nlanded(kk) -1
2756 ! WRITE(msg, '(4(i8,1x))') kk, np_nlanded(kk), kk1, kk2
2757 ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd master kk, np_dep, k1:k2: '//msg)
2758 np_landed_cnt(kk1:kk2) = fs_mpi_recvfrompatch_int(bsz=np_nlanded(kk), fromid=kk-1)
2759 np_landed_risk(kk1:kk2)= fs_mpi_recvfrompatch_real(bsz=np_nlanded(kk), fromid=kk-1)
2760 kk1 = kk1 + np_nlanded(kk)
2761 ENDIF
2762 ENDDO
2763 #endif
2765 !-------------------------------------------------------------------------------
2766 ! Calculate ratio and lkhd
2768 ALLOCATE(np_frac_landed(ndep_total), np_lkhd(ndep_total))
2769 np_frac_landed(:) = ZERO_dp
2770 np_lkhd(:) = ZERO_dp
2772 np_frac_landed(1:ndep_total) = REAL(np_landed_cnt(1:ndep_total))/REAL(SUM(np_landed_cnt))
2773 np_lkhd(1:ndep_total) = np_frac_landed(1:ndep_total) * np_landed_risk(1:ndep_total)
2774 np_lkhd(1:ndep_total) = np_lkhd(1:ndep_total)/MAXVAL(np_lkhd)
2776 ! DO pp=1,SIZE(np_landed_cnt)
2777 ! WRITE(msg, '(3(i8,1x),4(f12.6,1x))') pp, np_landed_cnt(pp), SUM(np_landed_cnt), np_landed_risk(pp), np_frac_landed(pp), MAXVAL(np_lkhd), np_lkhd(pp)
2778 ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd master np_frac_landed, np_lkhd, : '//msg)
2779 ! ENDDO
2781 !-------------------------------------------------------------------------------
2782 ! Send dat to procs
2784 kk1 = ndep + 1 ! index start after master's
2786 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2787 DO kk=2, mpiprocs ! kk=1 is Master
2789 IF (np_nlanded(kk) > 0) THEN
2791 kk2 = kk1 + np_nlanded(kk) -1
2793 CALL fs_mpi_sendbuff_real(sendto=kk-1, bsz=np_nlanded(kk), buff=np_frac_landed(kk1:kk2))
2794 CALL fs_mpi_sendbuff_real(sendto=kk-1, bsz=np_nlanded(kk), buff=np_lkhd(kk1:kk2))
2796 kk1 = kk1 + np_nlanded(kk)
2797 ENDIF
2798 ENDDO
2799 #endif
2800 ENDIF ! SUM(np_nlanded) > 0
2801 !-------------------------------------------------------------------------------
2803 ENDIF ! IamMaster
2805 !-------------------------------------------------------------------------------
2807 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2808 IF (.NOT.(IamMaster) .AND. (ndep > 0)) THEN
2810 recv_frac_landed = fs_mpi_recvfrompatch_real(bsz=ndep, fromid=MasterId)
2811 fs_frac_landed(is:ie,js:je) = UNPACK(recv_frac_landed, MASK=landing_mask(is:ie,js:je), FIELD=ZERO_dp)
2813 recv_spot_lkhd = fs_mpi_recvfrompatch_real(bsz=ndep, fromid=MasterId)
2814 fs_spotting_lkhd(is:ie,js:je) = UNPACK(recv_spot_lkhd, MASK=landing_mask(is:ie,js:je), FIELD=ZERO_dp)
2816 ENDIF
2817 #endif
2819 IF ((IamMaster) .AND. (ndep > 0)) THEN
2821 fs_frac_landed(is:ie,js:je) = UNPACK(np_frac_landed(1:ndep), MASK=landing_mask(is:ie,js:je), FIELD=ZERO_dp)
2822 fs_spotting_lkhd(is:ie,js:je) = UNPACK(np_lkhd(1:ndep), MASK=landing_mask(is:ie,js:je), FIELD=ZERO_dp)
2824 ENDIF
2826 ENDIF ! hist_flag
2828 !===============================================================================
2829 !===============================================================================
2832 !******************************************************************
2833 !******************************************************************
2834 !* *
2835 !* MPI Receive firebrands *
2836 !* *
2837 !******************************************************************
2838 !******************************************************************
2841 !===============================================================================
2842 ! MPI: Receive Brands (into this patch)
2843 !===============================================================================
2847 !-------------------------------------------------------------------------------
2848 ! Anything to receive? This part is only needed when compiled in parallel
2849 !-------------------------------------------------------------------------------
2851 #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
2853 ALLOCATE(nbr_id(neighbors))
2854 nbr_id(:)=0
2855 nbr_id(:)=fs_mpi_checkreceive(task_list=task_id, np=neighbors)
2857 nbr_sum = SUM(nbr_id)
2858 IF (nbr_sum > 0) THEN
2860 WRITE (msg,'(16(i4,1x))') ([task_id(ii), nbr_id(ii)], ii=1,neighbors)
2861 CALL wrf_debug (wrfdbg, 'SPFire_driver mpi_check_receive (taskID, N): '//msg)
2863 ! Receiving positions, id, dt
2864 ALLOCATE(r_x(nbr_sum), r_y(nbr_sum), r_z(nbr_sum), r_id(nbr_sum), r_src(nbr_sum), r_dt(nbr_sum))
2865 r_id(:) = 0
2866 r_src(:) = 0
2867 r_dt(:) = 0
2868 r_x(:) = ZERO_dp
2869 r_y(:) = ZERO_dp
2870 r_z(:) = ZERO_dp
2872 ! Receiving p_propertieserties
2873 ALLOCATE(r_p_m(nbr_sum), r_p_d(nbr_sum), r_p_e(nbr_sum), r_p_t(nbr_sum), r_p_v(nbr_sum))
2874 ALLOCATE(fb_mdetv(nbr_sum))
2876 CALL fs_mpi_recv(np_id=nbr_id, task_id=task_id, &
2877 r_x=r_x, r_y=r_y, r_z=r_z, &
2878 r_id=r_id, r_src=r_src, r_dt=r_dt, &
2879 r_p_m=r_p_m, &
2880 r_p_d=r_p_d, &
2881 r_p_e=r_p_e, &
2882 r_p_t=r_p_t, &
2883 r_p_v=r_p_v)
2885 ! Assign values to array and release incoming particles
2886 !release_mpi_ijk = RESHAPE([r_x,r_y,r_z], [nbr_sum,3])
2887 fb_mdetv = [(p_properties(r_p_m(kk), r_p_d(kk), r_p_e(kk), r_p_t(kk), r_p_v(kk)), kk=1,nbr_sum)]
2889 ! DO kk=1, nbr_sum
2890 ! WRITE(msg,'(3(i6,1x),3(f12.6,1x))') r_id(kk), r_src(kk), r_dt(kk), r_x(kk), r_y(kk), r_z(kk)
2891 ! CALL wrf_debug (wrfdbg, 'SPFire_mpi recv AFTER >>> '// msg)
2892 ! ENDDO
2894 prior_br = active_br + 1 ! Array is packed so new particles will be assgned at this position
2895 CALL generate_firebrands(&
2896 fs_p_id = fs_p_id, &
2897 fs_p_src = fs_p_src, &
2898 fs_p_dt = fs_p_dt, &
2899 fs_p_z = fs_p_z, &
2900 fs_p_x = fs_p_x, &
2901 fs_p_y = fs_p_y, &
2902 fs_gen_idmax = fs_gen_idmax, &
2903 myprocid = myprocid, &
2904 active_br = active_br, &
2905 release_i = r_x, &
2906 release_j = r_y, &
2907 release_k = r_z, &
2908 release_dt = r_dt, &
2909 release_src = r_src, &
2910 release_prop= fb_mdetv, &
2911 fs_p_prop = fs_p_prop)
2913 WRITE (msg,'(2(i8,1x))') active_br, fs_gen_idmax
2914 CALL wrf_debug (wrfdbg, 'SPFire_driver mpi recv AFTER : ii, fs_gen_idmax >>> '// msg)
2916 fs_p_mass (prior_br:active_br) = r_p_m
2917 fs_p_diam (prior_br:active_br) = r_p_d
2918 fs_p_effd(prior_br:active_br) = r_p_e
2919 fs_p_temp (prior_br:active_br) = r_p_t
2920 fs_p_tvel (prior_br:active_br) = r_p_v
2922 ! DO kk=prior_br,prior_br+MIN(nbr_sum,5)
2923 ! WRITE(msg,'(2(i6,1x),6(f12.6,1x))') fs_p_id(kk), fs_p_dt(kk), fs_p_x(kk), fs_p_y(kk), fs_p_z(kk), fs_p_diam(kk), fs_p_temp(kk), fs_p_tvel(kk)
2924 ! CALL wrf_debug (wrfdbg, 'SPFire_mpi recv AFTER2 >>> '// msg)
2925 ! ENDDO
2927 DEALLOCATE(r_x, r_y, r_z, r_id, r_src, r_dt)
2928 !DEALLOCATE(release_mpi_ijk)
2929 DEALLOCATE(fb_mdetv)
2930 DEALLOCATE(r_p_m, r_p_d, r_p_e, r_p_t, r_p_v)
2932 ! ELSE
2933 ! CALL wrf_debug (wrfdbg, 'SPFire_driver mpi nothing to receive')
2934 ENDIF
2935 DEALLOCATE(nbr_id)
2936 !-------------------------------------------------------------------------------
2937 ! DONE :)
2938 !-------------------------------------------------------------------------------
2939 #endif
2940 !===============================================================================
2941 ! Check history interval
2942 !===============================================================================
2944 hist_flag = .FALSE.
2947 END SUBROUTINE firebrand_spotting_em_driver
2948 !=============================================================
2950 !=============================================================
2954 !=============================================================
2955 !=============================================================
2957 !******************************************************************
2958 !******************************************************************
2959 !* *
2960 !* Particle Transport & Physics *
2961 !* *
2962 !******************************************************************
2963 !******************************************************************
2965 !=============================================================
2966 PURE &
2967 SUBROUTINE advect_xyz_m(grid, xp, yp, hgt, dtp, idp, &
2968 u, v, w, dt, mf, z_at_w, znw, ims, jms, kms, &
2969 phyd, thet, rho, &
2970 xout, yout, zout, fs_p_prop, land_hgt, &
2971 start_mom3d_dt, msg)
2972 !=============================================================
2974 ! Transport a particle at a given xp, yp, zp position in meters
2975 ! given the corresponding velocities [m/s] at the enclosing grid points
2977 ! arguments
2978 ! xp, yp, hgt: starting particle position
2979 ! u, v, w: wind components
2980 ! dt: time interval in sec
2981 ! ims, jms : start index of wind components
2982 ! mf: constant to convert from grid position to meters: rdx*msft(x,y)
2983 ! z_at_w: wrf height at vertical interface
2984 ! znw: eta levels
2986 !-------------------------------------------------------------
2988 IMPLICIT NONE
2990 TYPE(domain), INTENT(IN) :: grid
2991 INTEGER, INTENT(IN) :: ims, jms, kms, start_mom3d_dt
2992 REAL, INTENT(IN) :: dt
2993 REAL, INTENT(IN) :: land_hgt
2994 REAL, INTENT(IN), DIMENSION(:) :: xp, yp, hgt ! particle position
2995 INTEGER, INTENT(IN), DIMENSION(:) :: dtp, idp ! particle life-time and ID
2996 REAL, INTENT(IN), DIMENSION(ims:, kms:, jms:) :: u, v, w, z_at_w
2997 REAL, INTENT(IN), DIMENSION(ims:, kms:, jms:) :: phyd, thet, rho
2998 REAL, INTENT(IN), DIMENSION(:) :: znw
2999 REAL, INTENT(IN), DIMENSION(ims:, jms:) :: mf ! map factor, msft * rdx
3000 REAL, INTENT(OUT), DIMENSION(:) :: xout, yout, zout!, aglh
3001 TYPE(p_properties),INTENT(INOUT),DIMENSION(:):: fs_p_prop
3003 REAL, DIMENSION(3) :: uvw2p
3004 INTEGER :: pp, aux
3005 INTEGER :: hsz, ihs, jhs, ihe, jhe ! tile +- halo array bounds
3006 REAL :: xp_m0, yp_m0, zp_m0, zp0
3007 REAL :: xp_m1, yp_m1, zp_m1, xp1, yp1, zp1
3008 REAL :: xp_m2, yp_m2, zp_m2, xp2, yp2, zp2
3009 REAL :: zp, wp, brz
3010 REAL :: loc_p, loc_t, loc_d ! local met properties
3011 CHARACTER (LEN=200), INTENT(OUT) :: msg
3013 WRITE(msg, '(a100)') ' '
3014 !-------------------------------------------------------------------------------
3016 hsz = 4 ! halo size - get this value from grid%. Halo:48 = 4 points. How about staggering?
3017 ihs = is - hsz ! MAX(, ids)
3018 jhs = js - hsz ! MAX(, jds)
3019 ihe = ie -1 + hsz ! MIN(, ide-1) ! check for ide -1 bc we use i, i+1 to interpolate
3020 jhe = je -1 + hsz ! MIN(, jde-1)
3022 ! WRITE (msg,'(4(i4,1x))') ihs,ihe, jhs, jhe
3023 ! CALL wrf_debug (wrfdbg, 'SPFire_advect halo bounds ihs, ihe, jhs, jhe: >>> '//msg)
3025 !start_mom3d_dt = 4 ! min lifetime before calculating fb physics (number of timesteps for this nested grid)
3027 DO pp=1, SIZE(xp)
3029 ! CALL wrf_debug (wrfdbg, 'SPFire advect_xyz_m ---------------------------------------------------- ')
3031 ! WRITE (msg,'(3(i8,1x),3(f12.6,1x))') pp, idp(pp), dtp(pp), xp(pp), yp(pp), hgt(pp)
3032 ! CALL wrf_debug (wrfdbg, 'SPFire BEFORE advect: pp,id,dt,xp,yp,zp: >>> '//msg)
3034 !-------------------------------------------------------------------------------
3035 ! Convert horizontal positions to meters
3036 !-------------------------------------------------------------------------------
3038 xp_m0 = xp(pp) / mf(FLOOR(xp(pp)),FLOOR(yp(pp)))
3039 yp_m0 = yp(pp) / mf(FLOOR(xp(pp)),FLOOR(yp(pp)))
3040 zp_m0 = hgt(pp)
3042 ! zp changes slightly between wrf time steps, so it need to be calculated again
3043 zp0 = hgt2k(xp=xp(pp), yp=yp(pp), hgt=zp_m0, z_at_w=z_at_w, znw=znw)
3044 ! (Need fractional position for interpolation)
3046 ! WRITE (msg,'(3(i8,1x),4(f12.6,1x))') pp, idp(pp), dtp(pp), xp(pp), yp(pp), hgt(pp), zp0
3047 ! CALL wrf_debug (wrfdbg, 'SPFire_advect BEFORE: pp,id,dt,xp,yp,zp,zk: >>> '//msg)
3049 ! If this is a particle deposited by an adjacent tile, do not transport it
3050 IF (zp_m0 <= land_hgt) THEN
3051 ! CALL wrf_debug (wrfdbg, 'SPFire_advect particle at deposit threshold - will not advect')
3053 xout(pp) = xp(pp)
3054 yout(pp) = yp(pp)
3055 zout(pp) = zp_m0
3056 CYCLE
3058 ENDIF
3060 !-------------------------------------------------------------------------------
3061 !-------------------------------------------------------------------------------
3064 ! 1ST PASS
3065 !=============================================================
3068 !-------------------------------------------------------------------------------
3069 ! #1 interpolate velocities to particle position & advect particle
3070 !-------------------------------------------------------------------------------
3072 uvw2p = uvw_3d_interp(x=xp(pp), y=yp(pp), z=zp0, u=u, v=v, w=w, ihs=ihs, jhs=jhs, ihe=ihe, jhe=jhe) ! args: grid fractional positions
3074 xp_m1 = uvw2p(1)*DT + xp_m0
3075 yp_m1 = uvw2p(2)*DT + yp_m0
3076 zp_m1 = uvw2p(3)*DT + zp_m0
3078 !-------------------------------------------------------------------------------
3079 ! Convert horizontal positions from meters
3080 !-------------------------------------------------------------------------------
3082 xp1 = xp_m1 * mf(FLOOR(xp(pp)),FLOOR(yp(pp))) ! 2 pass will account for grid pt changes in mf & u
3083 yp1 = yp_m1 * mf(FLOOR(xp(pp)),FLOOR(yp(pp)))
3085 IF((ABS(xp1-xp(pp)) > 2.0) .OR. (ABS(yp1-yp(pp)) > 2.0)) THEN
3086 WRITE (msg,'(2(i4,1x),7(f10.6,1x))') pp, dtp(pp), xp(pp), yp(pp), zp0, &
3087 uvw2p(1), uvw2p(2), uvw2p(3), dt !
3088 zout(pp) = ZERO_dp
3089 xout(pp) = ZERO_dp
3090 yout(pp) = ZERO_dp
3091 CYCLE
3092 ENDIF
3094 !-------------------------------------------------------------------------------
3095 ! Convert z position from meters to k-fractional
3096 !-------------------------------------------------------------------------------
3098 ! CALL wrf_debug (wrfdbg, 'SPFire advect_xyz_m 1st pass --------------------------------------------- ')
3100 ! WRITE (msg,'(3(i8,1x),3(f12.6,1x))') pp, idp(pp), dtp(pp), xp1, yp1, zp_m1
3101 ! CALL wrf_debug (wrfdbg, 'SPFire after 1-pass adv: pp,id,dt,xp,yp,zp: >>> '//msg)
3102 ! WRITE (msg,*) pp, z_at_w(FLOOR(xp1),1,FLOOR(yp1))
3103 ! CALL wrf_debug (wrfdbg, 'SPFire after 1-pass adv: zn: >>> '//msg)
3104 ! WRITE (msg,*) pp, z_at_w(FLOOR(xp1-dp05),1,FLOOR(yp1-dp05))
3105 ! CALL wrf_debug (wrfdbg, 'SPFire after 1-pass adv: zn: >>> '//msg)
3107 zp1 = hgt2k(xp=xp1, yp=yp1, hgt=zp_m1, z_at_w=z_at_w, znw=znw)
3108 wp = uvw2p(3)
3110 !-------------------------------------------------------------------------------
3111 ! Is the new position not on the same memory?
3112 !-------------------------------------------------------------------------------
3113 ! check for floor(x)+1 bc we use i and i+1 to interpolate
3114 IF (FLOOR(xp1-dp05) < ihs .OR. FLOOR(xp1-dp05) +1 > ihe .OR. &
3115 FLOOR(yp1-dp05) < jhs .OR. FLOOR(yp1-dp05) +1 > jhe) THEN
3116 ! Checking staggered positions relative to unstaggered variables (+- dp05)
3118 !-------------------------------------------------------------------------------
3119 ! Use only one pass
3120 !-------------------------------------------------------------------------------
3122 xout(pp) = xp1
3123 yout(pp) = yp1
3124 zout(pp) = zp_m1
3125 zp = zp1
3127 ELSE
3129 ! 2ND PASS
3130 !=============================================================
3132 !-------------------------------------------------------------------------------
3133 ! #2 interpolate velocities to particle position & advect particle
3134 !-------------------------------------------------------------------------------
3136 uvw2p = uvw_3d_interp(x=xp1, y=yp1, z=zp1, u=u, v=v, w=w, ihs=ihs, jhs=jhs, ihe=ihe, jhe=jhe) ! args: grid fractional positions
3138 ! WRITE (msg,'(3(f14.9,1x))') uvw2p(1), uvw2p(2), uvw2p(3)
3139 ! CALL wrf_debug (wrfdbg, 'SPFire_advect effective vel m/s: u,v,w 2nd pass>>> '//TRIM(msg))
3141 xp_m2 = uvw2p(1)*DT + xp_m0
3142 yp_m2 = uvw2p(2)*DT + yp_m0
3143 zp_m2 = uvw2p(3)*DT + zp_m0
3145 xp2 = xp_m2 * mf(FLOOR(xp1),FLOOR(yp1))
3146 yp2 = yp_m2 * mf(FLOOR(xp1),FLOOR(yp1))
3148 IF((ABS(xp2-xp1) > 2.0) .OR. (ABS(yp2-yp1) > 2.0)) THEN
3149 WRITE (msg,'(2(i4,1x),7(f10.6,1x))') pp, dtp(pp), xp1, yp1, zp1, &
3150 uvw2p(1), uvw2p(2), uvw2p(3), dt !
3151 zout(pp) = ZERO_dp
3152 xout(pp) = ZERO_dp
3153 yout(pp) = ZERO_dp
3154 CYCLE
3155 ENDIF
3158 !-------------------------------------------------------------------------------
3159 ! Is the new position not on the same memory?
3160 !-------------------------------------------------------------------------------
3162 ! check against tile +- halo points instead
3164 IF (FLOOR(xp2-dp05) < ihs .OR. FLOOR(xp2-dp05) +1 > ihe .OR. &
3165 FLOOR(yp2-dp05) < jhs .OR. FLOOR(yp2-dp05) +1 > jhe) THEN
3166 ! Checking staggered positions relative to unstaggered variables (+- dp05)
3168 !-------------------------------------------------------------------------------
3169 ! Use xp1, yp1 to find zp2 (z_at_w will be out of bounds)
3170 !-------------------------------------------------------------------------------
3172 zp2 = hgt2k(xp=xp1, yp=yp1, hgt=zp_m2, z_at_w=z_at_w, znw=znw)
3174 ELSE
3176 zp2 = hgt2k(xp=xp2, yp=yp2, hgt=zp_m2, z_at_w=z_at_w, znw=znw)
3178 ENDIF
3180 !=============================================================
3182 !-------------------------------------------------------------------------------
3183 ! Average both passes
3184 !-------------------------------------------------------------------------------
3186 xp2 = (xp1 + xp2) * dp05
3187 yp2 = (yp1 + yp2) * dp05
3188 zp_m2 = (zp_m1 + zp_m2) * dp05
3189 zp2 = (zp1 + zp2) * dp05
3191 xout(pp) = xp2
3192 yout(pp) = yp2
3193 zout(pp) = zp_m2
3194 zp = zp2
3195 wp = (uvw2p(3) + wp) * dp05
3197 ENDIF
3199 !=============================================================
3200 ! Transport Done :)
3201 !=============================================================
3202 ! WRITE (msg, '(a100)') ' '
3203 ! WRITE (msg,'(2(i6,1x),4(f12.6,1x))') pp, dtp(pp), xout(pp), yout(pp), zout(pp), zp
3204 ! CALL wrf_debug (wrfdbg, 'SPFire_advect AFTER : pp,id,dt,xp,yp,zp,kp: >>> '//msg)
3206 !=============================================================
3207 ! Firebrand Physics
3208 !=============================================================
3210 loc_p = ZERO_dp
3211 loc_d = ZERO_dp
3212 loc_t = ZERO_dp
3214 IF (zout(pp) <= ZERO_dp) THEN
3215 fs_p_prop(pp) = p_properties(ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP)
3216 zout(pp) = ZERO_dp
3217 CYCLE ! Skip physics
3218 ENDIF
3220 IF (IEEE_IS_NAN(zout(pp)) .OR. ABS(zout(pp)) > HUGE(1.0)) THEN
3221 fs_p_prop(pp) = p_properties(ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP)
3222 zout(pp) = ZERO_dp
3223 xout(pp) = ZERO_dp
3224 yout(pp) = ZERO_dp
3226 CYCLE ! Skip physics
3227 ENDIF
3229 IF (dtp(pp) < start_mom3d_dt) THEN ! Temporarily testing transport: extending particle lifetime
3230 ! WRITE (msg,*) dt
3231 ! CALL wrf_debug (wrfdbg, 'SPFire_physics cycle dt: '//msg)
3232 CYCLE
3233 ENDIF
3236 CALL get_local_met( xp = xout(pp), &
3237 yp = yout(pp), &
3238 zp = zp, &
3239 ihs = ihs, &
3240 jhs = jhs, &
3241 ihe = ihe, &
3242 jhe = jhe, &
3243 p = phyd, & ! air pressure (Pa)
3244 t = thet, & ! temperature (K)
3245 d = rho, & ! density (kg/m3)
3246 loc_p = loc_p, & ! local air pressure 3d (Pa)
3247 loc_d = loc_d, & ! local density (g/m3)
3248 loc_t = loc_t) ! local temperature (K)
3250 ! WRITE (msg,'(3(f16.8,1x))') grid%p_hyd_w(FLOOR(xout(pp)),FLOOR(zp),FLOOR(yout(pp))), &
3251 ! grid%rho(FLOOR(xout(pp)),FLOOR(zp),FLOOR(yout(pp))), &
3252 ! grid%t_phy(FLOOR(xout(pp)),FLOOR(zp),FLOOR(yout(pp)))
3253 ! CALL wrf_debug (wrfdbg, 'SPFire grid%p, %rho, %t_phy : '//msg)
3255 ! WRITE (msg,'(4(f16.8,1x))') loc_p, loc_d, loc_t, wp
3256 ! CALL wrf_debug (wrfdbg, 'SPFire br met p, d, t, w : '//msg)
3258 !-------------------------------------------------------------------------------
3259 ! Firebrand Burnout and Terminal Velocity
3260 !-------------------------------------------------------------------------------
3262 brz = zout(pp)
3263 CALL firebrand_physics(dt = dt, &
3264 hgt = brz, &
3265 loc_w = wp, &
3266 loc_p = loc_p, &
3267 loc_t = loc_t, &
3268 loc_d = loc_d, &
3269 fbprop = fs_p_prop(pp))
3271 IF (IEEE_IS_NAN(fs_p_prop(pp)%p_tvel)) THEN
3272 fs_p_prop(pp) = p_properties(ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP)
3273 zout(pp) = ZERO_dp
3275 ! WRITE (msg,'(3(i8,1x),4(f12.6,1x))') pp, idp(pp), dtp(pp), xout(pp), yout(pp), brz, fs_p_prop(pp)%p_tvel
3276 ! CALL wrf_debug (wrfdbg, 'SPFire_physics AFTER: pp,id,dt,xp,yp,zp,fbvt: >>> '//msg)
3278 ELSE
3280 zout(pp) = brz
3282 ENDIF
3283 ! WRITE (msg,'(3(i8,1x),4(f12.6,1x))') pp, idp(pp), dtp(pp), xout(pp), yout(pp), brz, fs_p_prop(pp)%p_tvel*dt
3284 ! CALL wrf_debug (wrfdbg, 'SPFire_physics AFTER: pp,id,dt,xp,yp,zp,fbvt: >>> '//msg)
3286 ENDDO
3288 END SUBROUTINE advect_xyz_m
3289 !=============================================================
3292 ! Cannot be elemental: array
3293 !=============================================================
3294 PURE FUNCTION uvw_3d_interp(x, y, z, u, v, w, ihs, jhs, ihe, jhe)
3295 !=============================================================
3297 ! Wrapper to interpolate U,V,W to a point within a 3D grid box
3299 ! Interpolate velocities at grid point interfaces to particle position
3300 ! input range needed: i:i+1 | i < x(i) <= i+1
3301 ! j:j+1 | j < y(j) <= j+1
3302 ! k:k+1 | k < z(k) <= k+1
3304 ! arguments
3305 ! x, y, z: fractional position to interpolate velocities (e.g.: z is a fraction of the k index)
3306 ! u, v, w: wind components
3307 ! is, js : start index of wind components
3309 !-------------------------------------------------------------
3310 ! Interpolate between planes: U: xy(k), xy(k+1) stag x, dim=1
3311 ! V: xy(k), xy(k+1) stag y, dim=2
3312 ! W: xz(j), xz(j+1) stag z, dim=2
3313 !
3314 ! particle position is not staggered:
3315 ! grid center is given by (FLOOR(x), FLOOR(y))
3316 !-------------------------------------------------------------
3318 IMPLICIT NONE
3320 INTEGER, INTENT(IN) :: ihs, jhs, ihe, jhe ! tile +- halo array bounds
3321 REAL, INTENT(IN) :: x, y, z ! fractional positions relative to integral array indices
3322 REAL, INTENT(IN), DIMENSION(ims:, 1:, jms:) :: u, v ,w
3323 REAL, DIMENSION(3) :: uvw_3d_interp
3324 REAL :: tru, trv, trw, x0, y0, z0
3325 INTEGER :: i, j, k, i0, j0, k0
3327 ! Unstaggerecd indices
3328 k = MAX(FLOOR(z), 1)
3329 j = FLOOR(y)
3330 i = FLOOR(x)
3332 ! Variable Staggering:
3333 ! particle position is equivalent to a staggered dimension
3334 ! the staggered domain begins 1/2 pnt before and ends 1/2 pnt after the regular domain
3335 ! e.g: interpolate u (stag x): x_stag = xp, y_unst = yp -0.5, z_unst = zp -0.5
3336 ! By shifting the particle position by 1/2 point,
3337 ! the unstaggered velocities will represent the grid edge instead of grid center
3339 x0 = x - dp05
3340 y0 = y - dp05
3341 z0 = z - dp05 ! z can be < 1.0 when below the first half level, u_3d_interp will extrapolate values
3343 i0 = FLOOR(x0) !MIN(MAX(ihs, ), ihe) ! shift indices by 0.5 because unstaggered variables given at horizontal grid center
3344 j0 = FLOOR(y0) !MIN(MAX(jhs, ), ihe)
3345 k0 = MAX(1, FLOOR(z0))
3347 !-----------------------------------------------------------------------------
3348 ! u,v interpolates between xy planes: xy(k) xy(k+1)
3349 !-----------------------------------------------------------------------------
3351 ! staggered x
3352 tru= u_3d_interp( x=x, y=y0, z=z0, &
3353 u_i0j0_bot=u(i, k0, j0 ), u_i0j1_bot=u(i, k0, j0+1), u_i1j0_bot=u(i+1, k0, j0 ), u_i1j1_bot=u(i+1, k0, j0+1), &
3354 u_i0j0_top=u(i, k0+1, j0 ), u_i0j1_top=u(i, k0+1, j0+1), u_i1j0_top=u(i+1, k0+1, j0 ), u_i1j1_top=u(i+1, k0+1, j0+1))
3356 ! staggered y
3357 trv= u_3d_interp( x=x0, y=y, z=z0, &
3358 u_i0j0_bot=v(i0, k0, j ), u_i0j1_bot=v(i0, k0, j+1), u_i1j0_bot=v(i0+1, k0, j ), u_i1j1_bot=v(i0+1, k0, j+1), &
3359 u_i0j0_top=v(i0, k0+1, j ), u_i0j1_top=v(i0, k0+1, j+1), u_i1j0_top=v(i0+1, k0+1, j ), u_i1j1_top=v(i0+1, k0+1, j+1))
3361 !-----------------------------------------------------------------------------
3362 ! w interpolates between xz planes, not xy (so y corresponds to z): xz(j) xz(j+1)
3363 !-----------------------------------------------------------------------------
3365 ! z and w are both staggered (z is derived from grid%z_at_w)
3366 trw= u_3d_interp( x=x0, y=z, z=y, &
3367 u_i0j0_bot=w(i0, k, j ), u_i0j1_bot=w(i0, k+1, j ), u_i1j0_bot=w(i0+1, k, j ), u_i1j1_bot=w(i0+1, k+1, j ), &
3368 u_i0j0_top=w(i0, k, j+1), u_i0j1_top=w(i0, k+1, j+1), u_i1j0_top=w(i0+1, k, j+1), u_i1j1_top=w(i0+1, k+1, j+1))
3370 !-----------------------------------------------------------------------------
3371 ! Return U, V, W at x, y, z positions
3372 !-----------------------------------------------------------------------------
3374 uvw_3d_interp = [tru, trv, trw]
3376 END FUNCTION uvw_3d_interp
3377 !=============================================================
3381 !=============================================================
3382 ELEMENTAL FUNCTION u_3d_interp(x, y, z, &
3383 u_i0j0_bot, u_i0j1_bot, u_i1j0_bot, u_i1j1_bot, &
3384 u_i0j0_top, u_i0j1_top, u_i1j0_top, u_i1j1_top )
3385 !=============================================================
3387 ! Interpolate velocities at grid point interfaces to particle position
3388 ! input range needed: u( i: i+1, j: j+1, k: k+1) |
3389 ! i < x <= i+1
3390 ! j < y <= j+1
3391 ! k < z <= k+1
3393 ! arguments
3394 ! x, y, z: particle position to interpolate velocities
3395 ! u, v, w: wind components
3397 !-----------------------------------------------------------------------------
3398 ! Bilinear interpolation :
3399 ! The product of the value at the desired point (x,y,z) and the area enclosed by the corners
3400 ! is equal to the sum of the
3401 ! products of the value at each corner and the partial area diagonally opposite the corner
3402 !-----------------------------------------------------------------------------
3404 !-----------------------------------------------------------------------------
3405 ! The 3-D interpolation is achieved with a linear interpolation between
3406 ! the xy-interpolated points from two adjacent planes
3407 !-----------------------------------------------------------------------------
3409 !-----------------------------------------------------------------------------
3410 ! Z must be offset by 0.5 to correctly interpolate between vertically unstaggered planes
3411 !-----------------------------------------------------------------------------
3413 REAL, INTENT(IN) :: x, y, z ! fractional position
3414 REAL, INTENT(IN) :: u_i0j0_bot, u_i0j1_bot, u_i1j0_bot, u_i1j1_bot ! u_lower
3415 REAL, INTENT(IN) :: u_i0j0_top, u_i0j1_top, u_i1j0_top, u_i1j1_top ! u_upper
3416 REAL :: u_3d_interp
3418 REAL :: dbot, dtop, u_lower, u_upper
3419 INTEGER :: k, j, i
3421 !-----------------------------------------------------------------------------
3422 ! Get weighted U, V, W at particle position based on surrounding xy and xz planes
3423 !-----------------------------------------------------------------------------
3425 ! index
3426 k = MAX(1, FLOOR(z)) ! To interpolate w, z corresponds y: ( x=x, y=z, z=y, ...
3428 ! z can be < 1.0 when interpolating below the first level on half levels
3429 ! e.g.: z=0.75, dbot = 0.75 -1 =-0.25, dtop = 2 - 0.75 =1.25
3430 ! u_3d_interp = u_upper*-0.25 + u_lower*1.25
3432 dbot = z - REAL(k)
3433 dtop = REAL(k+1) - z
3435 !-----------------------------------------------------------------------------
3436 ! j = FLOOR(y)
3437 ! i = FLOOR(x)
3439 u_lower = u_2d_interp( u_i0j0=u_i0j0_bot, u_i0j1=u_i0j1_bot, &
3440 u_i1j0=u_i1j0_bot, u_i1j1=u_i1j1_bot, xp=x, yp=y)
3441 u_upper = u_2d_interp( u_i0j0=u_i0j0_top, u_i0j1=u_i0j1_top, &
3442 u_i1j0=u_i1j0_top, u_i1j1=u_i1j1_top, xp=x, yp=y)
3444 !-----------------------------------------------------------------------------
3445 ! Apply weights and calculate U at particle position
3446 !-----------------------------------------------------------------------------
3448 u_3d_interp = u_upper * dbot + & !
3449 u_lower * dtop
3451 !-----------------------------------------------------------------------------
3453 END FUNCTION u_3d_interp
3454 !=============================================================
3456 !=============================================================
3460 !=============================================================
3461 ELEMENTAL FUNCTION u_2d_interp(u_i0j0, u_i0j1, u_i1j0, u_i1j1, xp, yp)
3462 !=============================================================
3464 ! Interpolate any variable at grid point interfaces to a position in between
3465 ! Minimum range: u(i:i+1, j:j+1) | i < xp <= i+1 & j < yp <= j+1
3467 ! arguments
3468 ! u : 2D variable at grid 4 corners: (i0,j0); (i0,j1); (i1,j0); (i1,j1)
3469 ! xp, yp : position within grid box to interpolate velocity
3470 ! is, js : start index of variable array boundary
3472 !-------------------------------------------------------------
3473 ! 2D Interpolation sketch
3474 ! a, b, c, d: grid point corners
3475 ! da, db, dc, dd: distance between the nearest corner and point projection
3476 !
3477 ! a + + + + + + d
3478 ! + dc | db +
3479 ! +____|____+
3480 ! + dd | da +
3481 ! + | +
3482 ! b + + + + + + c
3483 !-------------------------------------------------------------
3485 IMPLICIT NONE
3487 !INTEGER, INTENT(IN) :: is, js
3488 REAL, INTENT(IN) :: xp, yp
3489 REAL, INTENT(IN) :: u_i0j1, u_i0j0, u_i1j0, u_i1j1 ! DIMENSION(is:,js:) :: uu
3490 REAL :: u_2d_interp
3491 REAL :: d_a, d_b, d_c, d_d ! grid vertices
3492 REAL :: uu_a, uu_b, uu_c, uu_d ! advection at vertices
3493 INTEGER :: i, j
3495 i = FLOOR(xp)
3496 j = FLOOR(yp)
3498 uu_a = u_i0j1 ! uu(i ,j+1)
3499 uu_b = u_i0j0 ! uu(i ,j )
3500 uu_c = u_i1j0 ! uu(i+1,j )
3501 uu_d = u_i1j1 ! uu(i+1,j+1)
3503 ! WRITE (msg,'(4(f14.9,1x))') uu_a, uu_b, uu_c, uu_d
3504 ! CALL wrf_debug (wrfdbg, 'SPFire_uu_2d: uu _a,_b,_c,_d >>> '//TRIM(msg))
3506 d_a = ABS( ( REAL(i+1) - xp ) * ( REAL(j ) - yp ) )
3507 d_b = ABS( ( REAL(i+1) - xp ) * ( REAL(j+1) - yp ) )
3508 d_c = ABS( ( REAL(i ) - xp ) * ( REAL(j+1) - yp ) )
3509 d_d = ABS( ( REAL(i ) - xp ) * ( REAL(j ) - yp ) )
3511 ! WRITE (msg,'(4(f14.9,1x))') d_a, d_b, d_c, d_d
3512 ! CALL wrf_debug (wrfdbg, 'SPFire_uu_2d: d _a,_b,_c,_d >>> '//TRIM(msg))
3515 u_2d_interp = ( uu_a * d_a + &
3516 uu_b * d_b + &
3517 uu_c * d_c + &
3518 uu_d * d_d ) / &
3519 (d_a + d_b + d_c + d_d)
3521 ! WRITE (msg,*) uu_2d_interp
3522 ! CALL wrf_debug (wrfdbg, 'SPFire_uu_2d: >>> '//TRIM(msg))
3524 END FUNCTION u_2d_interp
3525 !=============================================================
3527 !=============================================================
3532 !=============================================================
3533 PURE SUBROUTINE remove_br( fs_p_id, fs_p_x, fs_p_y, fs_p_z, fs_p_dt, fs_p_effd, cnt, max_life_dt, land_hgt)
3534 !=============================================================
3535 IMPLICIT NONE
3536 INTEGER, INTENT(INOUT), DIMENSION(:) :: fs_p_id
3537 INTEGER, INTENT(IN), DIMENSION(:) :: fs_p_dt
3538 !INTEGER, INTENT(IN) :: active_br
3539 REAL, INTENT(IN), DIMENSION(:) :: fs_p_x, fs_p_y, fs_p_z, fs_p_effd
3540 LOGICAL, DIMENSION(fs_array_maxsize) :: sparse_mask
3541 INTEGER, INTENT(OUT) :: cnt
3543 INTEGER, INTENT(IN) :: max_life_dt
3544 REAL, INTENT(IN) :: land_hgt
3546 !-------------------------------------------------------------------------------
3548 sparse_mask = .FALSE.
3549 !-------------------------------------------------------------------------------
3550 ! Flag to remove brands beyond domain bounds
3551 !-------------------------------------------------------------------------------
3553 sparse_mask = ( fs_p_id > 0 .AND. &
3554 (FLOOR(fs_p_x-dp05) < ids .OR. FLOOR(fs_p_y-dp05) < jds .OR. &
3555 FLOOR(fs_p_x-dp05)+1 > ide .OR. FLOOR(fs_p_y-dp05)+1 > jde))
3556 ! interpolation needs i, j, i+1, j+1
3557 ! First we unstagger fs_p_x, fs_p_y to align with met variables (x=fs_p_x-0.5, y=fs_p_y-0.5)
3558 ! Then we use FLOOR(x), FLOOR(y), FLOOR(x)+1, FLOOR(y)+1 to interpolate
3560 !-------------------------------------------------------------------------------
3561 ! Flag to remove brands living longer than firebrand_max_life_dt
3562 !-------------------------------------------------------------------------------
3564 sparse_mask = ( sparse_mask .OR. &
3565 (fs_p_id > 0 .AND. fs_p_dt > max_life_dt))
3567 !-------------------------------------------------------------------------------
3568 ! Flag to remove brands with near zero mass
3569 !-------------------------------------------------------------------------------
3571 ! sparse_mask = (sparse_mask .OR. &
3572 ! (fs_p_id > 0 .AND. fs_p_mass <= br_min_mass))
3574 !-------------------------------------------------------------------------------
3575 ! Flag to remove brands with near zero diameter that are not at deposit height
3576 !-------------------------------------------------------------------------------
3578 sparse_mask = (sparse_mask .OR. &
3579 (fs_p_id > 0 .AND. ((fs_p_effd <= TINY(1.0) .AND. (fs_p_z > land_hgt)))))
3581 !-------------------------------------------------------------------------------
3582 ! Reset array
3583 !-------------------------------------------------------------------------------
3584 cnt = 0
3585 cnt = COUNT(sparse_mask)
3586 IF (cnt > 0) THEN
3588 WHERE(sparse_mask) fs_p_id = 0
3590 ! WRITE (msg,'(3(i8,1x))') COUNT(sparse_mask), COUNT( fs_p_id > 0 )
3591 ! CALL wrf_debug (wrfdbg, 'SPFire_remove removed, remaining: '//msg)
3592 !CALL wrf_error_fatal( "end of test")
3593 ENDIF
3595 END SUBROUTINE remove_br
3597 !=============================================================
3598 !=============================================================
3601 !=============================================================
3602 PURE SUBROUTINE deposit_br(fs_p_id, fs_p_x, fs_p_y, fs_p_z, sum_fbrand, land_hgt)
3603 !=============================================================
3605 IMPLICIT NONE
3606 INTEGER, INTENT(INOUT), DIMENSION(:) :: fs_p_id
3607 REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: sum_fbrand
3608 REAL, INTENT(IN), DIMENSION(:) :: fs_p_x, fs_p_y, fs_p_z
3609 REAL, INTENT(IN) :: land_hgt
3611 !INTEGER, INTENT(IN) :: active_br
3613 !INTEGER, ALLOCATABLE, DIMENSION(:) :: mask_idx
3614 INTEGER, ALLOCATABLE, DIMENSION(:) :: rx, ry, rid
3615 LOGICAL, DIMENSION(fs_array_maxsize) :: sparse_mask
3616 LOGICAL, DIMENSION(fs_array_maxsize) :: bounds_mask
3617 !LOGICAL, ALLOCATABLE, DIMENSION(:) :: bounds_mask
3618 INTEGER :: nresets
3619 INTEGER :: x, y, i0,k ! counters
3620 !CHARACTER(LEN=10) fmt
3622 !CALL wrf_debug (wrfdbg, 'SPFire_deposit: check')
3623 sum_fbrand(:,:) = ZERO_dp
3625 !-------------------------------------------------------------------------------
3626 ! Flag brands near ground (deps_lev is a declared parameter)
3627 !-------------------------------------------------------------------------------
3629 sparse_mask = .FALSE.
3630 bounds_mask = .FALSE.
3632 ! Only deposit particles on this tile
3633 bounds_mask = ((FLOOR(fs_p_x) >= is) .OR. &
3634 (FLOOR(fs_p_x) <= ie) .OR. &
3635 (FLOOR(fs_p_y) >= js) .OR. &
3636 (FLOOR(fs_p_y) <= je))
3638 sparse_mask = ( bounds_mask .AND. (fs_p_id > 0 .AND. fs_p_z <= land_hgt))
3640 !-------------------------------------------------------------------------------
3641 ! Deposit flagged brands
3642 !-------------------------------------------------------------------------------
3643 nresets = 0
3644 nresets = COUNT(sparse_mask)
3646 ! WRITE (msg,'(i6)') nresets
3647 ! CALL wrf_debug (wrfdbg, 'SPFire_deposit total: '//msg)
3649 IF (nresets > 0) THEN
3651 !-------------------------------------------------------------------------------
3652 ! convert position to integer indices
3653 !-------------------------------------------------------------------------------
3655 ALLOCATE(rx(nresets),ry(nresets))!, mask_idx(nresets), rid(nresets), bounds_mask(nresets))
3657 ! Get indices of particles when sparse_mask = True
3658 rx(:) = 0
3659 ry(:) = 0
3661 rx(1:nresets) = FLOOR(PACK(fs_p_x, sparse_mask))
3662 ry(1:nresets) = FLOOR(PACK(fs_p_y, sparse_mask))
3664 !-------------------------------------------------------------------------------
3665 ! Increment sum_fbrand
3666 !-------------------------------------------------------------------------------
3668 DO k=1,nresets
3669 sum_fbrand(rx(k),ry(k)) = sum_fbrand(rx(k),ry(k)) + 1.0_dp
3670 ENDDO
3672 !-------------------------------------------------------------------------------
3673 ! Reset array
3674 !-------------------------------------------------------------------------------
3676 WHERE(sparse_mask) fs_p_id= 0
3678 ! WRITE (msg,'(5(i8,1x))') COUNT(sparse_mask), COUNT( fs_p_id > 0 ), COUNT(rx >= ie), COUNT(ry >= je)
3679 ! CALL wrf_debug (wrfdbg, 'SPFire_deposit deposited, remaining, outofbounds: '//msg)
3681 ENDIF
3683 END SUBROUTINE deposit_br
3686 !=============================================================
3687 !=============================================================
3691 !------------------------------------------------------------------------------
3692 ! UTILITIES
3693 !------------------------------------------------------------------------------
3696 !=============================================================
3698 SUBROUTINE get_local_ijk(grid, &
3699 ims, ime, jms, jme, kms, kme, &
3700 ips, ipe, jps, jpe, kps, kpe, &
3701 ifps, ifpe, jfps, jfpe, &
3702 ifms, ifme, jfms, jfme, &
3703 ids, ide, jds, jde, kds, kde, &
3704 m_idim, m_jdim, m_kdim, &
3705 p_idim, p_jdim, p_kdim, &
3706 d_idim, d_jdim, d_kdim)
3707 !p2m_is, p2m_ie, p2m_js, p2m_je)
3709 USE module_domain, ONLY: get_ijk_from_grid, get_ijk_from_subgrid
3711 IMPLICIT NONE
3713 TYPE(DOMAIN), INTENT(IN) :: grid
3714 INTEGER, INTENT(OUT), OPTIONAL :: ims, ime, jms, jme, kms, kme
3715 INTEGER, INTENT(OUT), OPTIONAL :: ips, ipe, jps, jpe, kps, kpe
3716 INTEGER, INTENT(OUT), OPTIONAL :: ifps, ifpe, jfps, jfpe
3717 INTEGER, INTENT(OUT), OPTIONAL :: ifms, ifme, jfms, jfme
3718 INTEGER, INTENT(OUT), OPTIONAL :: ids, ide, jds, jde, kds, kde
3719 INTEGER, INTENT(OUT), OPTIONAL :: m_idim, m_jdim, m_kdim
3720 INTEGER, INTENT(OUT), OPTIONAL :: p_idim, p_jdim, p_kdim
3721 INTEGER, INTENT(OUT), OPTIONAL :: d_idim, d_jdim, d_kdim
3722 !INTEGER, INTENT(OUT), OPTIONAL :: p2m_is, p2m_ie, p2m_js, p2m_je
3725 INTEGER :: iims, iime, jjms, jjme, kkms, kkme
3726 INTEGER :: iips, iipe, jjps, jjpe, kkps, kkpe
3727 INTEGER :: iids, iide, jjds, jjde, kkds, kkde
3729 INTEGER :: iifps, iifpe, jjfps, jjfpe, kkfps, kkfpe
3730 INTEGER :: iifms, iifme, jjfms, jjfme, kkfms, kkfme
3731 INTEGER :: iifds, iifde, jjfds, jjfde, kkfds, kkfde
3733 IF ((.NOT. PRESENT(ims)) .AND. &
3734 (.NOT. PRESENT(jms)) .AND. &
3735 (.NOT. PRESENT(kms)) .AND. &
3736 (.NOT. PRESENT(ime)) .AND. &
3737 (.NOT. PRESENT(jme)) .AND. &
3738 (.NOT. PRESENT(kme)) .AND. &
3739 !
3740 (.NOT. PRESENT(ips)) .AND. &
3741 (.NOT. PRESENT(jps)) .AND. &
3742 (.NOT. PRESENT(kps)) .AND. &
3743 (.NOT. PRESENT(ipe)) .AND. &
3744 (.NOT. PRESENT(jpe)) .AND. &
3745 (.NOT. PRESENT(kpe)) .AND. &
3746 !
3747 (.NOT. PRESENT(ifps)) .AND. &
3748 (.NOT. PRESENT(jfps)) .AND. &
3749 (.NOT. PRESENT(ifpe)) .AND. &
3750 (.NOT. PRESENT(jfpe)) .AND. &
3751 !
3752 (.NOT. PRESENT(ifms)) .AND. &
3753 (.NOT. PRESENT(jfms)) .AND. &
3754 (.NOT. PRESENT(ifme)) .AND. &
3755 (.NOT. PRESENT(jfme)) .AND. &
3756 !
3757 (.NOT. PRESENT(ids)) .AND. &
3758 (.NOT. PRESENT(jds)) .AND. &
3759 (.NOT. PRESENT(kds)) .AND. &
3760 (.NOT. PRESENT(ide)) .AND. &
3761 (.NOT. PRESENT(jde)) .AND. &
3762 (.NOT. PRESENT(kde)) .AND. &
3763 !
3764 ! (.NOT. PRESENT(p2m_is)) .AND. &
3765 ! (.NOT. PRESENT(p2m_ie)) .AND. &
3766 ! (.NOT. PRESENT(p2m_js)) .AND. &
3767 ! (.NOT. PRESENT(p2m_je)) .AND. &
3768 !
3769 (.NOT. PRESENT(m_idim)) .AND. &
3770 (.NOT. PRESENT(m_jdim)) .AND. &
3771 (.NOT. PRESENT(m_kdim)) .AND. &
3772 (.NOT. PRESENT(p_idim)) .AND. &
3773 (.NOT. PRESENT(p_jdim)) .AND. &
3774 (.NOT. PRESENT(p_kdim)) .AND. &
3775 (.NOT. PRESENT(d_idim)) .AND. &
3776 (.NOT. PRESENT(d_jdim)) .AND. &
3777 (.NOT. PRESENT(d_kdim))) &
3778 !
3779 CALL wrf_debug ( 1, 'get_local_ijk function is NOT requesting a result' )
3781 CALL get_ijk_from_grid ( grid , &
3782 iids, iide, jjds, jjde, kkds, kkde, &
3783 iims, iime, jjms, jjme, kkms, kkme, &
3784 iips, iipe, jjps, jjpe, kkps, kkpe)
3786 IF (PRESENT(ifps) .OR. &
3787 PRESENT(jfps) .OR. &
3788 PRESENT(ifpe) .OR. &
3789 PRESENT(jfpe) .OR. &
3790 PRESENT(ifms) .OR. &
3791 PRESENT(jfms) .OR. &
3792 PRESENT(ifme) .OR. &
3793 PRESENT(jfme)) CALL get_ijk_from_subgrid(grid , &
3794 iifds, iifde, jjfds, jjfde, kkfds, kkfde, &
3795 iifms, iifme, jjfms, jjfme, kkfms, kkfme, &
3796 iifps, iifpe, jjfps, jjfpe, kkfps, kkfpe)
3799 IF (PRESENT(ims)) ims = iims
3800 IF (PRESENT(jms)) jms = jjms
3801 IF (PRESENT(kms)) kms = kkms
3802 IF (PRESENT(ime)) ime = iime
3803 IF (PRESENT(jme)) jme = jjme
3804 IF (PRESENT(kme)) kme = kkme
3806 IF (PRESENT(ips)) ips = iips
3807 IF (PRESENT(jps)) jps = jjps
3808 IF (PRESENT(kps)) kps = kkps
3809 IF (PRESENT(ipe)) ipe = iipe
3810 IF (PRESENT(jpe)) jpe = jjpe
3811 IF (PRESENT(kpe)) kpe = kkpe
3813 IF (PRESENT(ifps)) ifps = iifps
3814 IF (PRESENT(jfps)) jfps = jjfps
3815 IF (PRESENT(ifpe)) ifpe = iifpe
3816 IF (PRESENT(jfpe)) jfpe = jjfpe
3818 IF (PRESENT(ifms)) ifms = iifms
3819 IF (PRESENT(jfms)) jfms = jjfms
3820 IF (PRESENT(ifme)) ifme = iifme
3821 IF (PRESENT(jfme)) jfme = jjfme
3823 IF (PRESENT(ids)) ids = iids
3824 IF (PRESENT(jds)) jds = jjds
3825 IF (PRESENT(kds)) kds = kkds
3826 IF (PRESENT(ide)) ide = iide
3827 IF (PRESENT(jde)) jde = jjde
3828 IF (PRESENT(kde)) kde = kkde
3830 IF (PRESENT(m_idim)) m_idim = iime - iims + 1
3831 IF (PRESENT(m_jdim)) m_jdim = jjme - jjms + 1
3832 IF (PRESENT(m_kdim)) m_kdim = kkme - kkms + 1
3833 IF (PRESENT(p_idim)) p_idim = iipe - iips + 1
3834 IF (PRESENT(p_jdim)) p_jdim = jjpe - jjps + 1
3835 IF (PRESENT(p_kdim)) p_kdim = kkpe - kkps + 1
3836 IF (PRESENT(d_idim)) d_idim = iide - iids + 1
3837 IF (PRESENT(d_jdim)) d_jdim = jjde - jjds + 1
3838 IF (PRESENT(d_kdim)) d_kdim = kkde - kkds + 1
3840 ! IF (PRESENT(p2m_is)) p2m_is = iips - iims
3841 ! IF (PRESENT(p2m_ie)) p2m_ie = iipe - iims
3842 ! IF (PRESENT(p2m_js)) p2m_js = jjps - jjms
3843 ! IF (PRESENT(p2m_je)) p2m_je = jjpe - jjms
3845 ! p2m_is = ips - ims
3846 ! p2m_ie = ips - ims + p_idim - 1
3847 ! p2m_js = jps - jms
3848 ! p2m_je = jps - jms + p_jdim - 1
3850 END SUBROUTINE get_local_ijk
3852 !=============================================================
3853 !=============================================================
3858 END MODULE module_firebrand_spotting