phys/module_firebrand_spotting.F

   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 !-------------------------------------------------------------------------------
   6
   7 !=============================================================
   8 !=============================================================
   9
  10 MODULE module_firebrand_spotting
  11
  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
  19
  20     IMPLICIT NONE
  21
  22     PRIVATE
  23     PUBLIC firebrand_spotting_em_init, firebrand_spotting_em_driver, get_local_ijk
  24
  25     ! THESE VARIABLES ARE IN MODULE SCOPE ! Careful with reassignments - don't reassign
  26     ! SAVE attribute is default
  27
  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)
  37
  38     !-------------------------------------------------------------------------------
  39     ! Fixed parameters ***MODULE SCOPE***
  40     !-------------------------------------------------------------------------------
  41
  42     INTEGER, PARAMETER :: dp = KIND(0.d0)      ! double precision
  43
  44     INTEGER, PARAMETER :: wrfdbg = 0
  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
  49
  50     INTEGER, PARAMETER :: sp_fuel_src_typ_n = 1
  51
  52     REAL,    PARAMETER :: grav = 9.80616_dp    ! gravity (m/s2)
  53     REAL,    PARAMETER :: rdry = 287.04_dp     ! dry air (J/Kg-K)
  54     REAL,    PARAMETER :: p2jm = 100.0_dp      ! mb to j/m3
  55     REAL,    PARAMETER :: rcd  = 0.45_dp       ! drag constant
  56     REAL,    PARAMETER :: rd   = 287.15_dp    ! gas constant dry air (J / kg K)
  57     REAL,    PARAMETER :: rv   = 461.6_dp
  58     REAL,    PARAMETER :: t0   = 300.0_dp
  59     REAL,    PARAMETER :: cp   = 7.0_dp*rd/2.0_dp
  60     REAL,    PARAMETER :: rcp  = rd/cp
  61     REAL,    PARAMETER :: p1000mb= 100000.0_dp
  62     REAL,    PARAMETER :: r1o3   = 1.0_dp/3.0_dp  ! ratio 1 over 3
  63     REAL,    PARAMETER :: sboltz = 5.67E-5_dp     ! stefan-boltzmann (g / s3-K4)
  64     REAL,    PARAMETER :: emiss  = 0.9_dp         ! emissivity
  65     REAL,    PARAMETER :: cpw    = 1466.0_dp      ! specific heat wood (J / kg-K)
  66     REAL,    PARAMETER :: cpc    = 712.0_dp       ! specific heat char (J / kg-K)
  67     REAL,    PARAMETER :: beta0  = 4.8E-7_dp      ! burning rate constant (m2 / s)
  68     REAL,    PARAMETER :: s_coeff= 110.4_dp       ! Sutherland's law coefficient (K)
  69     REAL,    PARAMETER :: b_coeff= 1.458E-3_dp    ! Sutherland's law coefficient [g / m-s-K-(1/2)]
  70     REAL,    PARAMETER :: shmt   = 0.7_dp         ! schmidt number
  71     REAL,    PARAMETER :: thcon  = 27.0_dp        ! thermal conductivity air
  72
  73     !USE module_state_description, ONLY: p_qv
  74     !INTEGER, PARAMETER :: p_qv = 1
  75     REAL,    PARAMETER :: pr     = 0.7_dp         ! Prandtl
  76
  77     REAL,    PARAMETER :: NEGLIGIBLE = 10*EPSILON(1.0)
  78     REAL,    PARAMETER :: ZERO_dp = 0.0_dp ! this is a real type variable, not a double precision type
  79     REAL,    PARAMETER :: dp05 = 0.5_dp
  80     REAL,    PARAMETER :: dp1 = 1.0_dp
  81
  82
  83     !-------------------------------------------------------------------------------
  84     ! ***MODULE SCOPE***
  85     !-------------------------------------------------------------------------------
  86
  87     ! Mass threshold to consider burnout (g)
  88     REAL, PARAMETER :: br_min_mass = EPSILON(dp1) ! max precision for real(4)
  89
  90     ! Diameter threshold to consider burnout (mm)
  91     REAL, PARAMETER :: br_min_diam = 0.0000001_dp
  92
  93     !-------------------------------------------------------------------------------
  94     ! Generic variables for multiple use within module ***MODULE SCOPE***
  95     !-------------------------------------------------------------------------------
  96
  97     CHARACTER (LEN=200), SAVE     :: msg
  98     CHARACTER (LEN=256), SAVE     :: fmt
  99     CHARACTER (LEN=200), DIMENSION(10) :: amsg
 100     INTEGER, SAVE :: imsg ! loop counters
 101
 102     !-------------------------------------------------------------------------------
 103     ! variables from namelist ***MODULE SCOPE***
 104     !-------------------------------------------------------------------------------
 105
 106     ! These will be available to all subroutines and will be set in init routine
 107     INTEGER, SAVE :: fs_array_maxsize   ! maximum number of particles carried in simulation
 108     INTEGER, SAVE :: fs_gen_levels      ! one per grid pt (may change if releasing at multiple levels)
 109     INTEGER, SAVE :: fs_gen_lim
 110     INTEGER, SAVE :: fs_gen_dt
 111
 112     REAL,    SAVE :: firebrand_dens      ! 513000.0_dp     ! density of firebrand (g / m3)
 113     REAL,    SAVE :: firebrand_dens_char ! 299000.0_dp     ! density of char (g m-3) [gupta et al 2002; fuel]
 114
 115
 116     !-------------------------------------------------------------------------------
 117     ! Fixed indices, ranks ***MODULE SCOPE***
 118     !-------------------------------------------------------------------------------
 119
 120     LOGICAL              :: this_is_ideal = .FALSE.
 121
 122     TYPE p_properties  ! fb_prop
 123         REAL :: p_mass ! fbmass
 124         REAL :: p_diam ! fbdiam
 125         REAL :: p_effd ! fbediam
 126         REAL :: p_temp ! fbtemp
 127         REAL :: p_tvel ! fbvelo
 128     END TYPE p_properties
 129
 130     !-------------------------------------------------------------------------------
 131     ! grid and cf are not here because dimensions are given at runtime (derived types)
 132     ! Also, grid values change with timestep, so they need to be passed as dummy arguments
 133     !-------------------------------------------------------------------------------
 134
 135     !-------------------------------------------------------------------------------
 136     ! Variable bounds - Initialized in init, used in dummy arguments in driver
 137     ! ***MODULE SCOPE***
 138     !-------------------------------------------------------------------------------
 139     INTEGER, SAVE :: ids, jds, ide, jde, kde      ! domain bounds
 140     INTEGER, SAVE :: ims, jms, ime, jme, kms, kme ! memory bounds
 141     INTEGER, SAVE :: is, ie, js, je, ks, ke       ! patch start/end
 142     INTEGER, SAVE :: ifps, jfps, ifpe, jfpe       ! refined fire grid bounds
 143
 144 CONTAINS
 145
 146 !******************************************************************
 147 !******************************************************************
 148 !*                                                                *
 149 !*                       Module Routines                          *
 150 !*                                                                *
 151 !******************************************************************
 152 !******************************************************************
 153
 154
 155
 156 !=============================================================
 157     SUBROUTINE firebrand_spotting_em_init ( &
 158 !=============================================================
 159         grid,      &
 160         cf,        &
 161         fs_p_id,     &
 162         fs_p_dt,     &
 163         fs_p_x,      &
 164         fs_p_y,      &
 165         fs_p_z,      &
 166         fs_last_gen_dt,&
 167         fs_gen_idmax,  &
 168         fs_count_reset,&
 169         fs_p_mass, &
 170         fs_p_diam, &
 171         fs_p_effd,&
 172         fs_p_temp, &
 173         fs_p_tvel, &
 174         fs_count_landed_all,&
 175         fs_count_landed_hist,&
 176         fs_landing_mask,&
 177         fs_spotting_lkhd,&
 178         fs_frac_landed)
 179
 180 !=============================================================
 181 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
 182         USE module_firebrand_spotting_mpi,            ONLY: fs_mpi_init
 183 #endif
 184         !-------------------------------------------------------------------------------
 185         ! Initialize all I/O and variables declared in the registry
 186         !-------------------------------------------------------------------------------
 187
 188         IMPLICIT NONE
 189
 190         !-------------------------------------------------------------------------------
 191         ! Arguments
 192         !-------------------------------------------------------------------------------
 193
 194         TYPE(domain),               INTENT(IN) :: grid !
 195         TYPE(grid_config_rec_type), INTENT(IN) :: cf   !  run-time configuration (namelist) for domain
 196
 197         INTEGER,  INTENT(INOUT) :: fs_last_gen_dt, fs_gen_idmax
 198         LOGICAL,  INTENT(INOUT) :: fs_count_reset
 199
 200         ! Firebrand Spotting Particle Properties (fs_p_*)
 201         INTEGER,  INTENT(INOUT), DIMENSION(:)    :: fs_p_id, fs_p_dt
 202         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)
 203         REAL,     INTENT(INOUT), DIMENSION(:)    :: fs_p_mass, fs_p_diam, fs_p_effd
 204         REAL,     INTENT(INOUT), DIMENSION(:)    :: fs_p_temp, fs_p_tvel
 205         REAL,     INTENT(INOUT), DIMENSION(ims:,jms:)  :: fs_count_landed_all, fs_count_landed_hist, fs_spotting_lkhd, fs_frac_landed
 206         INTEGER,  INTENT(INOUT), DIMENSION(ims:,jms:) :: fs_landing_mask
 207
 208         !-------------------------------------------------------------------------------
 209
 210         !-------------------------------------------------------------------------------
 211         ! Module variables from namelist/default
 212
 213         fs_array_maxsize = cf%fs_array_maxsize
 214         fs_gen_levels    = cf%fs_firebrand_gen_levels
 215         fs_gen_lim       = cf%fs_firebrand_gen_lim
 216         fs_gen_dt        = cf%fs_firebrand_gen_dt
 217
 218         fmt = '(A,1x,I6)'
 219         WRITE (amsg(1),*)   'SPFire_init: fspotting_em_init'
 220         WRITE (amsg(2),fmt) 'SPFire_init:  firebrand limit =', fs_gen_lim
 221         WRITE (amsg(3),fmt) 'SPFire_init:               dt =', fs_gen_dt
 222         WRITE (amsg(4),fmt) 'SPFire_init: fs_array_maxsize =', fs_array_maxsize
 223         WRITE (amsg(5),fmt) 'SPFire_init:    fs_gen_levels =', fs_gen_levels
 224         DO imsg=1,5
 225             CALL wrf_debug (wrfdbg, TRIM(amsg(imsg)) )
 226         ENDDO
 227
 228         !-------------------------------------------------------------------------------
 229         ! Set bounds to be used as dummy arguments in driver
 230         ! ***variables declared in MODULE SCOPE***
 231         !-------------------------------------------------------------------------------
 232
 233         CALL get_local_ijk(grid, &
 234                            ifps=ifps, jfps=jfps, ifpe=ifpe, jfpe=jfpe, &
 235                            ids=ids, jds=jds, ide=ide, jde=jde, kde=kde, &
 236                            ims=ims, jms=jms, ime=ime, jme=jme, kms=kms, kme=kme, &
 237                            ips=is,  jps=js,  ipe=ie,  jpe=je,  kps=ks,  kpe=ke)
 238
 239         WRITE (msg,'(6(i6,1x))') is, ie, js, je, ks, ke
 240         CALL wrf_debug (wrfdbg, 'SPFire_init tile bounds: '//msg)
 241
 242         WRITE (msg,'(6(i6,1x))') ims, ime, jms, jme, kms, kme
 243         CALL wrf_debug (wrfdbg, 'SPFire_init memory bounds: '//msg)
 244
 245         WRITE (msg,'(4(i6,1x))') ifps, ifpe, jfps, jfpe
 246         CALL wrf_debug (wrfdbg, 'SPFire_init fire refined bounds: '//msg)
 247
 248
 249         !-------------------------------------------------------------------------------
 250         ! Initialize registry arrays
 251         !-------------------------------------------------------------------------------
 252         fs_count_reset=.FALSE.
 253
 254         fs_last_gen_dt = 0
 255         fs_gen_idmax   = 0
 256         fs_p_id    = 0
 257         fs_p_dt    = 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
 266
 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
 272
 273         !fs_fire_ROSdt(:,:) = 0.0_dp
 274         !fs_gen_inst(:,:) = 0
 275
 276         !-------------------------------------------------------------------------------
 277
 278         IF ( grid%this_is_an_ideal_run ) THEN
 279
 280             this_is_ideal = .TRUE.
 281             CALL wrf_debug (wrfdbg, 'SPFire_init: Ideal Run detected' )
 282
 283         ELSE
 284
 285             this_is_ideal = .FALSE.
 286             CALL wrf_debug (wrfdbg, 'SPFire_init: Not an Ideal Run' )
 287
 288         ENDIF
 289
 290
 291 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
 292         CALL fs_mpi_init(grid)
 293 #endif
 294
 295     END SUBROUTINE firebrand_spotting_em_init
 296 !=============================================================
 297 !=============================================================
 298
 299
 300 !=============================================================
 301     PURE &
 302     FUNCTION order_val(arr, ord)
 303 !=============================================================
 304
 305         IMPLICIT NONE
 306
 307         REAL,    INTENT(IN), DIMENSION(:) :: arr
 308         INTEGER, INTENT(IN) :: ord
 309         REAL    :: order_val
 310         INTEGER :: diff, besti, iord, i, error
 311
 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         !-----------------------------------------------------------
 324
 325         error = 2
 326
 327         diff = SIZE(arr) -1
 328         besti = 1
 329
 330         DO i=1,SIZE(arr)
 331
 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
 337
 338         ! WRITE(msg, '(2(i6,1x))') iord, diff
 339         ! CALL wrf_debug (wrfdbg, 'SPFire_rel_order_val: '//msg)
 340
 341         order_val = arr(besti)
 342
 343 END FUNCTION order_val
 344 !=============================================================
 345 !=============================================================
 346
 347
 348
 349 !=============================================================
 350     PURE &
 351     FUNCTION fire2tile(fr_arr, dsx,dsy) RESULT(new_arr)
 352 !=============================================================
 353
 354         IMPLICIT NONE
 355
 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
 361
 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         !-----------------------------------------------------------
 366
 367         fshp = SHAPE(fr_arr)
 368
 369         ALLOCATE(new_arr(1:fshp(1)/dsx, 1:fshp(2)/dsy))
 370         new_arr(:,:) = ZERO_dp
 371
 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])
 378
 379     END FUNCTION fire2tile
 380 !=============================================================
 381 !=============================================================
 382
 383
 384
 385 !=============================================================
 386     ELEMENTAL FUNCTION fuel_spotting_risk(fuel_fgi, factor, fuel_mcg)!, fuel_mcg)
 387 !=============================================================
 388
 389         ! Returns the refined grid array with the fuel_spotting_risk for spotfires
 390         !-------------------------------------------------------------------------------
 391
 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
 414
 415         IMPLICIT NONE
 416
 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)
 420
 421         fuel_spotting_risk = factor * (dp1 - MIN(dp1, fuel_mcg/fuel_fgi)) ! *(fuel_fgi-fuel_mcg)
 422
 423         END FUNCTION fuel_spotting_risk
 424 !=============================================================
 425 !=============================================================
 426
 427
 428
 429 !=============================================================
 430     ELEMENTAL FUNCTION spotting_threat_factor(fcat)
 431 !=============================================================
 432         IMPLICIT NONE
 433
 434         INTEGER, INTENT(IN)    :: fcat
 435         REAL,    DIMENSION(nfuelcats+1) :: factor
 436         REAL                            :: spotting_threat_factor
 437
 438         ! ****************               read this from external input              ****************
 439         ! **************** call it in firebrand_spotting_driver once, when grid%itimestep =1  ****************
 440         ! ****************              and load a registry array                   ****************
 441
 442         ! User defined coefficient to increase spotfire likelihood;
 443         ! threat = 0 yields zero likelihood
 444
 445         factor(1)     =  ZERO_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
 449
 450         spotting_threat_factor = factor(fcat)
 451
 452         !CALL read_table_static
 453
 454     END FUNCTION spotting_threat_factor
 455
 456 !=============================================================
 457 !=============================================================
 458
 459
 460 ! !=============================================================
 461 !     SUBROUTINE read_table_static ! Not debuggeed - from module_mp_thompson
 462 ! !=============================================================
 463
 464 !         USE module_domain
 465 !         USE module_dm
 466 !         IMPLICIT NONE
 467
 468 !         LOGICAL, EXTERNAL:: wrf_dm_on_monitor
 469 !         INTEGER:: iunit_fs1, i
 470 !         LOGICAL:: opened
 471 !         CHARACTER (LEN=64) :: filename
 472
 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
 483
 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')
 489
 490 !       IF ( wrf_dm_on_monitor() ) THEN
 491
 492 !           WRITE(msg, '(i2,1x,a20)') iunit_fs1, filename
 493 !           CALL wrf_debug(wrfdbg,'SPFire_read_table_static: opening table on unit ' // msg)
 494
 495 !           OPEN(iunit_fs1,FILE=tablename, FORM='UNFORMATTED',STATUS='OLD',ERR=9009)
 496
 497 !       ENDIF
 498
 499 ! #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes(A, size(A)*R4SIZE)
 500
 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
 505
 506 !       RETURN
 507
 508 !  9009 CONTINUE
 509 !       CALL wrf_error_fatal('SPFire_read_table_static: error openinig static table')
 510 !       RETURN
 511
 512 !  9010 CONTINUE
 513 !       CALL wrf_error_fatal('SPFire_read_table_static: error reading static table')
 514
 515
 516 !       END SUBROUTINE read_table_static
 517
 518 ! !=============================================================
 519 ! !=============================================================
 520
 521
 522 !=============================================================
 523     ELEMENTAL FUNCTION get_fuel_cat(crosswalk, cat)
 524 !=============================================================
 525
 526         ! Use with Scott & Burgen fuel model (WRF CO-FSP)
 527
 528         IMPLICIT NONE
 529
 530         INTEGER :: get_fuel_cat
 531         LOGICAL, INTENT(IN):: crosswalk
 532         REAL,    INTENT(IN):: cat
 533
 534         INTEGER, DIMENSION(1:nf_sb) :: ksb ! Anderson82 + S&B2005 fuel categories array
 535         INTEGER :: i, icat
 536         !-------------------------------------------------------------------------------
 537
 538         icat = INT(cat)
 539
 540         ksb = no_fuel_cat ! 14 ! Initialize all values to no-fuel
 541
 542         ! Anderson 1982
 543         ksb(1:13) = [(i, i=1, 13)]
 544
 545         IF ( crosswalk ) THEN ! Scott & Burgan crosswalks to Anderson
 546
 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
 561
 562         ELSE ! full Scott and Burgan (2005)
 563
 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)
 570
 571         ENDIF
 572
 573         get_fuel_cat = ksb(icat)
 574
 575     END FUNCTION get_fuel_cat
 576 !=============================================================
 577 !=============================================================
 578
 579
 580 !=============================================================
 581     ELEMENTAL FUNCTION firebrand_gen_factor(fcat)
 582 !=============================================================
 583         IMPLICIT NONE
 584
 585         REAL                :: firebrand_gen_factor
 586         INTEGER, INTENT(IN) :: fcat
 587         REAL,    DIMENSION(nfuelcats+1) :: factor
 588
 589
 590         ! ****************               read this from external input              ****************
 591         ! **************** call it in firebrand_spotting_driver once, when grid%itimestep =1  ****************
 592         ! ****************              and load a registry array                   ****************
 593
 594         ! User defined coefficient to set release priority
 595         ! factor = 0 yields zero releases from gridpoint
 596
 597         factor(1:nfuelcats+1)=   1.0_dp
 598         firebrand_gen_factor = factor(fcat)
 599
 600     END FUNCTION firebrand_gen_factor
 601
 602 !=============================================================
 603 !=============================================================
 604
 605
 606
 607 ! !=============================================================
 608 !     ELEMENTAL FUNCTION firebrand_gen_maxhgt(fcat)
 609 ! !=============================================================
 610 !         IMPLICIT NONE
 611
 612 !         REAL                :: firebrand_gen_maxhgt
 613 !         INTEGER, INTENT(IN) :: fcat
 614 !         REAL,    DIMENSION(nfuelcats+1) :: maxhgt
 615
 616
 617 !         ! ****************               read this from external input              ****************
 618 !         ! **************** call it in firebrand_spotting_driver once, when grid%itimestep =1  ****************
 619 !         ! ****************              and load a registry array                   ****************
 620
 621 !         ! User defined max height to release firebrand
 622
 623 !         maxhgt(1:nfuelcats+1)=  50.0_dp
 624 !         firebrand_gen_maxhgt = maxhgt(fcat)
 625
 626 !     END FUNCTION firebrand_gen_maxhgt
 627
 628 !=============================================================
 629 !=============================================================
 630
 631
 632
 633 !=============================================================
 634     ELEMENTAL FUNCTION firebrand_gen_potential(fuel_fgi, factor, fire_rate, fuel_mcg)!, fuel_mcg)
 635 !=============================================================
 636
 637         ! Returns the refined grid array with the potential for firebrand generation
 638         !-------------------------------------------------------------------------------
 639
 640         IMPLICIT NONE
 641
 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)
 645
 646         firebrand_gen_potential = factor * fire_rate * fuel_fgi !*(dp1 - fuel_mcg/(dp1+fuel_mcg)) ! fuelloadm; dry load from fr_fire_phys
 647
 648         END FUNCTION firebrand_gen_potential
 649 !=============================================================
 650 !=============================================================
 651
 652
 653
 654 !=============================================================
 655 !    PURE & ! can't call wrf_fatal in pure routines
 656     SUBROUTINE generate_firebrands( fs_p_id, fs_p_dt, fs_p_z, fs_p_x, fs_p_y,   &
 657                                release_i, release_j, release_k, &
 658                                active_br, fs_gen_idmax, &
 659                                release_prop, fs_p_prop)
 660 !=============================================================
 661
 662         IMPLICIT NONE
 663
 664         TYPE(p_properties), INTENT(IN), DIMENSION(:):: release_prop
 665         TYPE(p_properties), INTENT(OUT),DIMENSION(:):: fs_p_prop
 666
 667         INTEGER, INTENT(INOUT), DIMENSION(:)   :: fs_p_id, fs_p_dt
 668         REAL,    INTENT(INOUT), DIMENSION(:)   :: fs_p_z, fs_p_x, fs_p_y
 669         REAL,    INTENT(INOUT), DIMENSION(:)   :: release_i, release_j, release_k
 670         INTEGER, INTENT(INOUT)  :: fs_gen_idmax
 671         INTEGER, INTENT(OUT)    :: active_br
 672
 673         LOGICAL :: release_true
 674         INTEGER :: br, ii
 675         REAL    :: rx, ry, rz
 676         INTEGER :: new_release
 677         LOGICAL, DIMENSION(fs_array_maxsize) :: flag_true
 678
 679         !-------------------------------------------------------------------------------
 680         release_true = .FALSE.  ! scalar: true when brands are ready to be released
 681         flag_true = .FALSE.     ! logical array: flag to find active elements in array
 682
 683         active_br = 0           ! number of active brands
 684         new_release = 0         ! number of brands to release
 685
 686         ! when all brand IDs of array elements are zero: nothing active
 687         active_br   = COUNT( fs_p_id > 0 )
 688         new_release = COUNT( INT(release_i) > 0 )
 689
 690         IF (new_release > 0 ) release_true = .TRUE. ! to-do: include dependency on release_dt
 691
 692         IF ( .NOT. release_true) THEN
 693             RETURN !CALL wrf_debug (wrfdbg,  'SPFire_generate_firebrands: Found nothing to release')
 694         ENDIF
 695
 696         !-------------------------------------------------------------------------------
 697         ! New brands fit in array? Release particles: set lat/lon/lev of new brands
 698         !-------------------------------------------------------------------------------
 699
 700         IF (active_br + new_release <= fs_array_maxsize) THEN
 701
 702             !-------------------------------------------------------------------------------
 703             ! Add new grid positions to array: create a subroutine for this: brand_release
 704             !-------------------------------------------------------------------------------
 705
 706             ! Flag is True where positions are available to be filled
 707             WHERE(fs_p_id == 0) flag_true = .TRUE. ! id_indx holds brand ID (unique integer)
 708
 709             ii = active_br + 1
 710             DO br=1, new_release
 711
 712                 ! Find an ii position available to store a new release
 713                 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)
 714                     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?')
 715                 ENDIF
 716
 717                 ! Sanity check
 718                 IF (INT(release_i(br)) == 0 .OR. &
 719                     INT(release_j(br)) == 0 ) &
 720                     CALL wrf_error_fatal('SPFire_generate_firebrands: release_ijk is zero! Positions cannot be zero.')
 721
 722                 IF (fs_gen_idmax + 10 >= HUGE(1)) fs_gen_idmax = 0
 723
 724                 ! Assign values and convert height to z index
 725                 fs_p_x(ii) = release_i(br)
 726                 fs_p_y(ii) = release_j(br)
 727                 fs_p_z(ii) = release_k(br)
 728                 fs_p_id(ii) = fs_gen_idmax + 1              ! Assign new IDmax to particle
 729                 fs_p_dt(ii) = 0 ! firebrand has not been advected yet so dt is zero
 730                 fs_p_prop(ii)%p_mass  = release_prop(br)%p_mass
 731                 fs_p_prop(ii)%p_diam  = release_prop(br)%p_diam
 732                 fs_p_prop(ii)%p_effd = release_prop(br)%p_effd
 733                 fs_p_prop(ii)%p_temp  = release_prop(br)%p_temp
 734                 fs_p_prop(ii)%p_tvel  = release_prop(br)%p_tvel
 735
 736                 fs_gen_idmax = fs_p_id(ii)
 737
 738                 flag_true(ii) = .FALSE.
 739                 ii = ii + 1
 740
 741             ENDDO
 742
 743             active_br = active_br + new_release
 744             release_i = ZERO_dp
 745             release_j = ZERO_dp
 746             release_k = ZERO_dp
 747
 748 !             WRITE (msg,'(2(i6,1x))') new_release, fs_gen_idmax
 749 !             CALL wrf_debug (wrfdbg, 'SPFire_release:  Released new brands fs_gen_idmax: '//msg)
 750
 751 !         ELSE
 752
 753 !             WRITE (msg,*) active_br + new_release, fs_array_maxsize
 754 !             CALL wrf_debug (wrfdbg, 'SPFire_release: brand array is full, cannot release new brands'// msg)
 755
 756         ENDIF
 757
 758
 759     END SUBROUTINE generate_firebrands
 760
 761 !=============================================================
 762 !=============================================================
 763
 764
 765
 766
 767 !=============================================================
 768     PURE &
 769     FUNCTION hgt2k(xp, yp, hgt, z_at_w, znw)
 770 !=============================================================
 771
 772     !-------------------------------------------------------------------------------
 773     ! Converts height in meters to vertical level
 774     !-------------------------------------------------------------------------------
 775
 776         IMPLICIT NONE
 777         !INTEGER, INTENT(IN) :: ims, jms
 778         REAL,    INTENT(IN) :: hgt, xp, yp
 779         REAL,    INTENT(IN), DIMENSION(ims:,:,jms:) :: z_at_w
 780         REAL,    INTENT(IN), DIMENSION(:)         :: znw
 781
 782         REAL    :: zdz, z_lower_at_p, z_upper_at_p, x0,y0
 783         REAL    :: hgt2k
 784         INTEGER :: k, i, j
 785
 786         ! An (xp,yp)=(1.5, 1.5) places the particle at the gridpoint center
 787         ! but z_at_w is horizontally unstaggered,
 788         ! so z_at_w(i,j) is at the horizontal grid center when (i,j)=(1,1)
 789         !
 790         ! Hence, we adjust the particle position to fetch the horizontally adjacent heights for the linear interp.
 791         ! We subtract the deltax between
 792         ! position at grid center and position at grid edge from the particle position xp, yp, such that
 793         ! z_at_w(i,1,j) for (i,j)=(1,1) corresponds to the height at (xp,yp)=(1.5,1.5):
 794         ! (i = xp - 0.5, j = yp - 0.5) = (i=1, j=1), for (xp=1.5, yp=1.5)
 795
 796         x0 = xp - dp05
 797         y0 = yp - dp05
 798
 799         i = FLOOR(x0) ! shift indices by 0.5 because z_at_w is given at horizontal grid center
 800         j = FLOOR(y0)
 801
 802         ! Minloc where hgt - z_k is positive: level right below brand release
 803         k = MINLOC(hgt - z_at_w(i,:,j), dim=1, &
 804                    mask=( hgt - z_at_w(i,:,j) >= 0.0_dp ))
 805
 806         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),&
 807                                                  u_i1j0=z_at_w(i+1,k,j), u_i1j1=z_at_w(i+1,k,j+1))
 808
 809        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),&
 810                                                  u_i1j0=z_at_w(i+1,k+1,j), u_i1j1=z_at_w(i+1,k+1,j+1))
 811
 812         zdz = (hgt - z_lower_at_p) / (z_upper_at_p-z_lower_at_p)
 813         hgt2k    = k + zdz
 814
 815     END FUNCTION hgt2k
 816 !=============================================================
 817 !=============================================================
 818
 819
 820
 821
 822 !=============================================================
 823     PURE &
 824     SUBROUTINE releaseijk2atm (nij_2d, sr_x, sr_y, fcat, maxhgt_usr, pi, pj, pk, nij)
 825 !=============================================================
 826
 827         IMPLICIT NONE
 828
 829         INTEGER, INTENT(IN),  DIMENSION(:,:) :: nij_2d ! do i need ifps here? No
 830         INTEGER, INTENT(IN),  DIMENSION(:,:) :: fcat   ! do i need ifps here? No
 831         INTEGER, INTENT(IN)   :: sr_x, sr_y
 832         REAL,    INTENT(IN)   :: maxhgt_usr
 833
 834         REAL,    INTENT(INOUT), DIMENSION(:)   :: pj, pi, pk
 835         INTEGER, INTENT(INOUT), DIMENSION(:)   :: nij
 836
 837         INTEGER, ALLOCATABLE, DIMENSION(:)   :: fpi,fpj
 838         INTEGER :: i, j, cnt
 839
 840         !-------------------------------------------------------------------------------
 841         ! Returns a 1-D array with the release positions (i, j, k[m])
 842         !   and number of particles to release (nij)
 843         !   converted from fire refine grid to wrf grid
 844         !   (i,j k are type float and represent the original refined grid position)
 845         !-------------------------------------------------------------------------------
 846
 847         cnt = COUNT(nij_2d > 0)
 848
 849         pi(:) = ZERO_dp
 850         pj(:) = ZERO_dp
 851         pk(:) = ZERO_dp
 852         nij(:) = 0
 853
 854         ALLOCATE(fpi(cnt), fpj(cnt))
 855         fpi(:) = 0
 856         fpj(:) = 0
 857
 858
 859         ! fpi = RESHAPE([( [(i, i=ifps, ifpe)], j=jfps, jfpe)], SHAPE(fr_arr))
 860         ! fpj = RESHAPE([( [(j, i=ifps, ifpe)], j=jfps, jfpe)], SHAPE(fr_arr))
 861
 862         fpi = PACK(RESHAPE([( [(i, i=ifps, ifpe)], j=jfps, jfpe)], SHAPE(nij_2d)), mask=(nij_2d > 0))
 863         fpj = PACK(RESHAPE([( [(j, i=ifps, ifpe)], j=jfps, jfpe)], SHAPE(nij_2d)), mask=(nij_2d > 0))
 864         nij = PACK(nij_2d, mask=(nij_2d > 0))
 865
 866         !-------------------------------------------------------------------------------
 867         ! convert from refined fire grid to atm grid
 868         !-------------------------------------------------------------------------------
 869         pi = fire2atm_ij(fpi, sr_x )
 870         pj = fire2atm_ij(fpj, sr_y )
 871
 872         ! Release height and firebrand properties can depend on fuel category (in future developments)
 873         ! Because fcat is on the refined grid, hgt and fprop must be specified here
 874
 875         ! pk = firebrand_gen_maxhgt(PACK(fcat, mask=(nij_2d > 0))) ! Use this to specify hgts for fuel types
 876         pk = maxhgt_usr
 877
 878
 879     END SUBROUTINE  releaseijk2atm
 880 !=============================================================
 881 !=============================================================
 882
 883
 884
 885 !=============================================================
 886     ELEMENTAL FUNCTION fire2atm_ij (fp_ij, sr)
 887 !=============================================================
 888 ! Convert refined subgrid to default atm grid
 889 ! Use function separately for i and j:
 890 !     pi = fire2atm_ij(fpi, grid%sr_x )
 891 !     pj = fire2atm_ij(fpj, grid%sr_y )
 892
 893         IMPLICIT NONE
 894
 895         REAL  :: fire2atm_ij
 896         INTEGER, INTENT(IN) :: fp_ij
 897         INTEGER, INTENT(IN) :: sr
 898
 899         fire2atm_ij = ( dp1+ REAL(fp_ij - dp1)/REAL(sr) )
 900
 901 ! fpij = 298, sr = 3: 297/3+1 = 100
 902 ! fpij = 299, sr = 3: 298/3+1 = 100.333
 903 ! fpij = 300, sr = 3: 299/3+1 = 100.666
 904 ! fpij = 301, sr = 3: 300/3+1 = 101
 905
 906 ! fpij = 297, sr = 4: 296/4+1 = 75.00
 907 ! fpij = 298, sr = 4: 297/4+1 = 75.25
 908 ! fpij = 299, sr = 4: 298/4+1 = 75.50
 909 ! fpij = 300, sr = 4: 299/4+1 = 75.75
 910 ! fpij = 301, sr = 4: 299/4+1 = 76.00
 911
 912     END FUNCTION fire2atm_ij
 913 !=============================================================
 914 !=============================================================
 915
 916
 917
 918 !=============================================================
 919     ELEMENTAL FUNCTION atm2fire_ij (pij, sr)
 920 !=============================================================
 921 ! Convert default atm grid to lower left index of refined grid
 922 ! (atm_i, atm_j) = (f_i : f_i + srx , f_j : f_j + sry)
 923
 924         IMPLICIT NONE
 925
 926         INTEGER  :: atm2fire_ij
 927         INTEGER, INTENT(IN) :: pij
 928         INTEGER, INTENT(IN) :: sr
 929
 930         atm2fire_ij = (pij - 1) * sr + 1
 931
 932 ! pij = 100, sr = 3: 99*3+1 = 298
 933 ! pij = 101, sr = 3: 100*3+1 = 301
 934
 935     END FUNCTION atm2fire_ij
 936 !=============================================================
 937 !=============================================================
 938
 939
 940
 941 ! !=============================================================
 942 !     SUBROUTINE subgrid_fire_property(fire_ij, fuel_frac, grid)
 943 ! !=============================================================
 944
 945 ! ! Not used - may be delted from final code version
 946
 947 !         IMPLICIT NONE
 948
 949 !         TYPE(domain), INTENT(IN) :: grid ! input data **Note: Intent IN**
 950 !         REAL,    INTENT(OUT), DIMENSION(:) :: fuel_frac
 951 !         INTEGER, INTENT(IN),  DIMENSION(:,:) :: fire_ij
 952
 953 !         INTEGER, ALLOCATABLE, DIMENSION(:) :: fpi, fpj
 954 !         INTEGER :: k, cnt, sri, srj, srx, sry
 955
 956 !         srx = grid%sr_x
 957 !         sry = grid%sr_y
 958
 959 !         cnt = SIZE(fire_ij(:,1))
 960 !         fpi = atm2fire_ij( fire_ij(:,1), srx)
 961 !         fpj = atm2fire_ij( fire_ij(:,2), sry)
 962 !         !ALLOCATE( flame(cnt), fuel(cnt))
 963
 964 !         ! Get flame and fuel_frac maxvals from subgrids
 965 !         DO k=1,SIZE(fire_ij(:,1))
 966
 967 !             sri = fpi(k)
 968 !             srj = fpj(k)
 969
 970 !             fuel_frac(k)  = MAXVAL(grid%fuel_frac( sri: sri+srx-1, srj : srj+sry-1))
 971
 972 !             !WRITE (msg,*) MAXVAL(grid%zsf(sri:sri+dsx,srj:srj+dsy)) always zero??
 973
 974 !         ENDDO
 975
 976 !     END SUBROUTINE subgrid_fire_property
 977 ! !=============================================================
 978 ! !=============================================================
 979
 980
 981
 982 ! !=============================================================
 983 !     PURE FUNCTION uniq_ij (fpi, fpj)
 984 ! !=============================================================
 985 ! ! remove duplicated grid points (after converting from refined fire grid to grid)
 986
 987 ! ! Not used but quite handy - do not delete it.
 988
 989 !         IMPLICIT NONE
 990
 991 !         INTEGER, INTENT(IN), DIMENSION(:) :: fpi, fpj
 992 !         INTEGER, ALLOCATABLE, DIMENSION(:,:) :: uniq_ij
 993 !         LOGICAL, ALLOCATABLE, DIMENSION(:) :: mask
 994 !         INTEGER :: cnt, i
 995
 996 !         cnt = SIZE(fpi)
 997
 998 !         ALLOCATE(mask(cnt))
 999 !         mask = .TRUE.
1000
1001 !         DO i = cnt, 2, -1
1002 !             mask(i) = .NOT.( ANY(fpi(:i-1) == fpi(i) .AND. &
1003 !                                  fpj(:i-1) == fpj(i)))
1004 !         END DO
1005
1006 !         ALLOCATE( uniq_ij(COUNT(mask),2) )
1007 !         uniq_ij(:,1) = REAL(PACK(fpi, mask))
1008 !         uniq_ij(:,2) = REAL(PACK(fpj, mask))
1009
1010 !     END FUNCTION uniq_ij
1011
1012 ! !=============================================================
1013 ! !=============================================================
1014
1015
1016
1017 !=============================================================
1018     SUBROUTINE prep_release_hgt(release_i, release_j, release_k, release_n, release_prop, levrand_seed)
1019 !=============================================================
1020
1021         !----------------------------------------------------------
1022         ! Set release locations for new particles
1023         ! and distribute particles over fs_gen_levels heights
1024         !----------------------------------------------------------
1025
1026         IMPLICIT NONE
1027
1028         TYPE(p_properties), INTENT(INOUT),   DIMENSION(:) :: release_prop
1029         REAL,          INTENT(INOUT), DIMENSION(:) :: release_i, release_j
1030         REAL,          INTENT(INOUT), DIMENSION(:) :: release_k
1031         INTEGER,       INTENT(IN),    DIMENSION(:) :: release_n
1032         INTEGER,       INTENT(IN) :: levrand_seed
1033
1034         INTEGER :: ii, kk, cnt                   ! counters
1035         INTEGER :: nseeds
1036         REAL,    ALLOCATABLE, DIMENSION(:) :: rand_arr, frachgt
1037         INTEGER, ALLOCATABLE, DIMENSION(:) :: seeds
1038
1039         !-------------------------------------------------------------------------------
1040         ! Calculate number of releases and heights from this fire
1041         !-------------------------------------------------------------------------------
1042
1043         ALLOCATE(frachgt(SIZE(release_n)*fs_gen_levels))
1044         ALLOCATE(rand_arr(SIZE(release_n)*fs_gen_levels))
1045
1046         IF (fs_gen_levels == 1) THEN
1047             frachgt(:) = dp1
1048         ELSE
1049             frachgt(:) = [( [(REAL(ii-1)* (dp1/REAL(fs_gen_levels-1)), ii=fs_gen_levels, 1, -1)], &
1050                                        kk=1,SIZE(release_n))]
1051         ENDIF
1052
1053         IF (levrand_seed > 0) THEN
1054
1055             nseeds = SIZE(rand_arr)
1056             CALL random_seed(size=nseeds)
1057             ALLOCATE(seeds(nseeds))
1058             seeds = [(( release_i(ii) * release_j(ii) * levrand_seed * kk, &
1059                         kk=1, fs_gen_levels -1), &
1060                         ii=1, SIZE(release_n))] ! force seed to vary across space
1061             CALL random_seed(put = seeds)
1062             DEALLOCATE(seeds)
1063             CALL random_number(rand_arr)
1064
1065         ENDIF
1066
1067
1068         ii = SIZE(release_n)+1           ! array index for release at various height (append to end of the array)
1069         DO cnt=1,SIZE(release_n)         ! loop over release gridpoints and set a new release at the same ij for 1 < hgt < release_k
1070
1071             DO kk=1, fs_gen_levels -1 ! release at multiple levels
1072
1073                 !-------------------------------------------------------------------------------
1074                 ! Fill release array
1075                 !-------------------------------------------------------------------------------
1076
1077                 release_i(ii) = release_i(cnt)    ! ii increments inside this loop
1078                 release_j(ii) = release_j(cnt)
1079                 release_k(ii) = (release_k(cnt)-dp1)*frachgt(ii)+dp1 ! keep a minimum height of 1-m
1080
1081                 ! if we want to perturb the max_release_hgt in releases_k(cnt), we can change its value now
1082                 IF (levrand_seed > 0) &
1083                     release_k(ii) = (release_k(cnt)-dp1) * rand_arr(ii) + dp1 ! keep a minimum height of 1-m
1084
1085                 !release_prop(cnt) = firebrand_property(pprop(cnt))
1086                 release_prop(ii)  = release_prop(cnt)
1087
1088                 ! IF (cnt < 3) THEN
1089                 !     WRITE (msg,'(3(i4,1x),5(f8.2,1x))') cnt, kk, ii, &
1090                 !         release_i(ii), release_j(ii), release_k(ii), release_k(cnt), rand_arr(ii)
1091                 !     CALL wrf_debug (wrfdbg, 'SPFire_prep_release c,k,i, ijk: '//msg)
1092                 ! ENDIF
1093
1094                 ii = ii + 1 ! increment release array index
1095
1096             ENDDO
1097
1098         ENDDO
1099
1100     END SUBROUTINE prep_release_hgt
1101 !=============================================================
1102 !=============================================================
1103
1104
1105
1106 !=============================================================
1107     FUNCTION firebrand_property(arr) RESULT(prop)
1108 !=============================================================
1109
1110         IMPLICIT NONE
1111
1112         TYPE (p_properties) :: prop
1113         REAL, INTENT(IN),  DIMENSION(:) :: arr
1114         REAL, PARAMETER :: pi    = 4.0_dp * ATAN (1.0_dp)
1115
1116         prop%p_diam  = arr(1)   ! 10.0_dp    ! mm
1117         prop%p_effd  = arr(2)   ! 10.0_dp    ! mm
1118         prop%p_temp  = arr(3)   ! 900.0_dp   ! K
1119         ! termvel should start with 0m/s because firebrand is released at top of trajectory
1120         prop%p_tvel  = arr(4)   ! 0.0_dp     ! m/s
1121         prop%p_mass  =  (firebrand_dens * pi * (prop%p_effd/1000.0_dp)**2)/6.0_dp  ! g
1122
1123     END FUNCTION firebrand_property
1124
1125 !=============================================================
1126 !=============================================================
1127
1128
1129
1130 !=============================================================
1131     FUNCTION idx_packed_1d(mask) RESULT(mask_idx)
1132 !=============================================================
1133
1134         ! Return an array of indices where mask is true
1135         IMPLICIT NONE
1136         LOGICAL,  INTENT(IN),  DIMENSION(:) :: mask
1137         !INTEGER,  INTENT(IN) :: active_br
1138         INTEGER, ALLOCATABLE, DIMENSION(:) :: mask_idx
1139         INTEGER :: nresets, ii
1140
1141         nresets = COUNT(mask)
1142         ALLOCATE(mask_idx(nresets))
1143         mask_idx = PACK([(ii, ii=1,SIZE(mask))], mask)   ! Get flag indices
1144
1145     END FUNCTION idx_packed_1d
1146 !=============================================================
1147 !=============================================================
1148
1149
1150
1151
1152 !=============================================================
1153     PURE &
1154     SUBROUTINE firebrand_physics(dt, hgt, loc_p, loc_t, loc_d, loc_w, fbprop)
1155 !=============================================================
1156
1157         IMPLICIT NONE
1158
1159         TYPE(p_properties), INTENT(INOUT) :: fbprop
1160         REAL, INTENT(IN)    :: loc_p, loc_t, loc_d, loc_w  ! local pressure, temp, density, w
1161         REAL, INTENT(IN)    :: dt
1162         REAL, INTENT(INOUT) :: hgt
1163         !REAL, INTENT(INOUT) :: p_mass, p_diam, p_effd, p_temp, p_tvel
1164
1165         !-------------------------------------------------------------------------------
1166         ! firebrand burnout
1167         !-------------------------------------------------------------------------------
1168
1169         ! WRITE (msg,'(5(f12.8,1x))') fbprop%p_mass, fbprop%p_diam, &
1170         !                             fbprop%p_effd, fbprop%p_temp, fbprop%p_tvel
1171         ! CALL wrf_debug (wrfdbg, 'SPFire br physics1 m, d, e, t, v: '//msg)
1172
1173         CALL burnout(p_mass  = fbprop%p_mass,   & ! firebrand mass (g)
1174                      p_diam  = fbprop%p_diam,   & ! firebrand diameter (mm)
1175                      p_effd  = fbprop%p_effd,   & ! firebrand effective diameter (mm)
1176                      p_temp  = fbprop%p_temp,   & ! firebrand temperature (K)
1177                      p_tvel  = fbprop%p_tvel,   & ! firebrand terminal velocity (m/s)
1178                      aird    = loc_d,    & ! g/m3
1179                      pres    = loc_p,    & ! Pa
1180                      temp    = loc_t,    & ! K
1181                      loc_w   = loc_w,    & ! m/s
1182                      dt      = dt)
1183
1184 !         WRITE (msg,'(4(f12.8,1x))') fbprop%p_mass, fbprop%p_diam, &
1185 !                                     fbprop%p_effd, fbprop%p_temp
1186 !         CALL wrf_debug (wrfdbg, 'SPFire br physics2 m, d, e, t: '//msg)
1187
1188         !-------------------------------------------------------------------------------
1189         ! firebrand terminal velocity
1190         !-------------------------------------------------------------------------------
1191
1192         CALL termvel(p_diam  = fbprop%p_diam, &
1193                      p_effd  = fbprop%p_effd,&
1194                      p_temp  = fbprop%p_temp, &
1195                      p_tvel  = fbprop%p_tvel, &
1196                      aird    = loc_d,  & ! g/m3
1197                      pres    = loc_p,  & ! Pa
1198                      temp    = loc_t,  & ! K
1199                      dt      = dt,     & ! sec
1200                      hgt     = hgt)      ! m
1201
1202         ! WRITE (msg,'(6(f12.8,1x))') fbprop%p_mass, fbprop%p_diam, &
1203         !                             fbprop%p_effd, fbprop%p_temp, fbprop%p_tvel, hgt
1204         ! CALL wrf_debug (wrfdbg, 'SPFire br physics2 m,d,e,t,v,h : '//msg)
1205
1206         ! **********
1207         ! NOTE:
1208         ! Should we calculate in double precision when diam is in meters?
1209         ! double: real(8), 15 precision digits
1210         ! single: real(4), 7  precision digits
1211
1212     END SUBROUTINE firebrand_physics
1213
1214 !=============================================================
1215 !=============================================================
1216
1217
1218
1219 !=============================================================
1220      ELEMENTAL SUBROUTINE burnout(pres, aird, temp, loc_w, dt, p_mass, p_diam, p_effd, p_temp, p_tvel)
1221 !=============================================================
1222
1223 !-------------------------------------------------------------------------------
1224 ! This subroutine was developed in RAL-NCAR/UCAR, 2019-2020 by:
1225 ! *** T. W. Juliano ***
1226 ! based on
1227 ! Tse & Fernandez-Pello 1998 (DOI: 10.1016/ S0379-7112(97)00050-7)
1228 ! Bhutia, Jenkins & Sun 2010 (DOI:10.3894/JAMES.2010.2.4)
1229 !-------------------------------------------------------------------------------
1230
1231         IMPLICIT NONE
1232
1233         REAL,         INTENT(IN)    :: pres, aird, temp, loc_w, dt
1234         REAL,         INTENT(INOUT) :: p_mass    ! firebrand mass (g)
1235         REAL,         INTENT(INOUT) :: p_diam    ! firebrand diameter (mm)
1236         REAL,         INTENT(INOUT) :: p_effd   ! firebrand effective diameter (mm)
1237         REAL,         INTENT(INOUT) :: p_temp    ! firebrand temperature (K)
1238         REAL,         INTENT(IN)    :: p_tvel    ! new firebrand terminal velocity (m/s)
1239
1240         !-------------------------------------------------------------------------------
1241         ! Constants are at module top
1242         !-------------------------------------------------------------------------------
1243
1244         REAL, PARAMETER :: pi    = 4.0_dp * ATAN (1.0_dp)
1245
1246         REAL :: aird2, wind, gama
1247         REAL :: pdia, pedia
1248         REAL :: p_effd2, pdia_new4, reyn, beta, deff, dmvc, dmvc2
1249         REAL :: parta, partb, dtemp, nuss, hbar, fbvol, qcon, qrad, partc
1250         REAL :: pratio, cpmix
1251
1252         p_mass = MAX(p_mass, EPSILON(dp1))
1253         !-------------------------------------------------------------------------------
1254
1255         ! air density surrounding firebrand
1256         aird2 = 1000.0_dp * pres / ( rd * p_temp )
1257
1258         ! Assume vel diff between firebrand and wind is due only to
1259         ! vertical; that is, assume particle moves with horiz wind
1260         wind = ABS(loc_w - p_tvel)
1261
1262         ! particle diameter (mm) to (m)
1263         pdia = p_diam / 1000.0_dp
1264
1265         ! particle effective diameter (mm) to (m)
1266         pedia = p_effd / 1000.0_dp
1267
1268         ! local atmospheric dynamic viscosity according to Sutherland's law
1269         dmvc = (b_coeff * temp**1.5_dp) / (temp + s_coeff)
1270
1271         ! dynamic viscosity surrounding firebrand according to Sutherland's law
1272         dmvc2 = (b_coeff * p_temp**1.5_dp) / (p_temp + s_coeff)
1273
1274         ! kinematic viscosity average for atmosphere and firebrand
1275         gama = ((dmvc / aird) + (dmvc2 / aird2))/2.0_dp
1276
1277         ! reynolds number
1278         reyn = wind * pdia / gama
1279
1280         ! nusselt number
1281         nuss = 2.0_dp + 0.6_dp * reyn**(dp05) * pr**(r1o3)
1282
1283         ! convection heat transfer coefficient
1284         hbar = nuss * thcon / pdia
1285
1286         ! burning rate
1287         beta = beta0 + beta0 * (0.276_dp * reyn**(dp05) * shmt**(r1o3))
1288
1289         ! WRITE (msg,'(9(f12.8,1x),f12.3,1x,f12.8)') aird2, wind, pdia, &
1290         !                              pedia, dmvc, dmvc2, &
1291         !                               gama, reyn, nuss, hbar, beta
1292         ! CALL wrf_debug (wrfdbg, 'SPFire burnout1: '//msg)
1293
1294         ! change in firebrand mass diameter (m)
1295         parta = pdia**4
1296         partb = SQRT(3.0_dp) * (beta**2) * (dt**2)
1297         pdia_new4 = parta - partb
1298         p_diam = pdia_new4**(0.25_dp)
1299
1300         ! change in firebrand effective mass diameter (m)
1301         p_effd2 = (pedia**2) - beta * dt
1302         p_effd = SQRT(p_effd2)
1303
1304         ! new firebrand mass (g)
1305         p_mass = (firebrand_dens * pi * p_effd2 * p_effd)/6.0_dp
1306
1307         ! new firebrand volume (m3)
1308         fbvol = p_mass / firebrand_dens
1309
1310         ! new firebrand temp (K)
1311         qcon = hbar * (p_temp - temp)
1312         qrad = sboltz * emiss * (p_temp**4 - temp**4)
1313         pratio = p_effd / p_diam
1314         cpmix = (pratio * cpw) + ((1.0_dp - pratio) * cpc)
1315         partc = 6.0_dp / (p_mass * cpmix * p_diam)
1316         dtemp = dt * fbvol * partc * (qcon + qrad)
1317         p_temp = p_temp - dtemp
1318
1319         ! WRITE (msg,'(8(f14.12,1x))') parta, partb, pdia_new4, p_diam, p_effd2, p_effd, p_mass, fbvol
1320         ! CALL wrf_debug (wrfdbg, 'SPFire burnout2: '//msg)
1321         ! 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
1322         ! CALL wrf_debug (wrfdbg, 'SPFire burnout3: '//msg)
1323
1324         ! firebrand diameter from (m) to (mm)
1325         p_diam = 1000.0_dp * p_diam
1326         p_effd = 1000.0_dp * p_effd
1327
1328     END SUBROUTINE burnout
1329
1330 !=============================================================
1331 !=============================================================
1332
1333
1334
1335 !=============================================================
1336     ELEMENTAL  &
1337     SUBROUTINE termvel(p_diam, p_effd, p_tvel, p_temp, hgt, dt, pres, aird, temp)
1338 !=============================================================
1339
1340 !-------------------------------------------------------------------------------
1341 ! This subroutine was developed in RAL-NCAR/UCAR, 2019-2020 by:
1342 ! *** T. W. Juliano ***
1343 ! based on Bhutia, Jenkins & Sun 2010 (DOI:10.3894/JAMES.2010.2.4)
1344 !-------------------------------------------------------------------------------
1345
1346         ! Constants are at module top
1347
1348         IMPLICIT NONE
1349
1350         REAL, INTENT(IN) :: pres, aird, temp, dt
1351         REAL, INTENT(IN) :: p_diam    ! firebrand diameter (mm)
1352         REAL, INTENT(IN) :: p_effd    ! firebrand effective diameter (mm)
1353         REAL, INTENT(IN) :: p_temp    ! firebrand temperature (K)
1354         REAL, INTENT(INOUT):: hgt     ! particle position at height agl (m)
1355         REAL, INTENT(INOUT):: p_tvel  ! new firebrand terminal velocity (m/s)
1356
1357         REAL    :: pbot, ptop, aird2, pdia, parta, partb
1358         REAL    :: drop, vt, pratio ! fbtype,
1359
1360         !-------------------------------------------------------------------------------
1361
1362         ! air density around particle
1363         aird2 = 1000.0_dp * pres / (rd * p_temp)
1364
1365         ! particle diameter (mm) to (m)
1366         pdia = p_diam /1000.0_dp
1367
1368         ! terminal velocity
1369         pratio = p_effd / p_diam
1370         parta = ( (pratio*firebrand_dens) + ( (dp1-pratio) * firebrand_dens_char) ) *pdia*grav
1371         partb = 3.0_dp * ((aird + aird2) / 2.0_dp)*rcd
1372
1373         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
1374
1375         ! WRITE (msg,'(4(f12.6,1x))') parta, partb, pratio, SQRT(pratio)
1376         ! CALL wrf_debug (wrfdbg, 'SPFire termvel : '//msg)
1377
1378         vt = SQRT(pratio)
1379
1380         ! change in vertical position due to settling
1381         drop = vt * ABS(dt)
1382         pbot = MAX(0.0_dp, hgt-drop)
1383         hgt  = pbot
1384
1385         p_tvel=vt
1386
1387     END SUBROUTINE termvel
1388
1389 !=============================================================
1390 !=============================================================
1391
1392
1393 !=============================================================
1394     PURE  &
1395     SUBROUTINE get_local_met(xp, yp, zp, loc_p, loc_d, loc_t, p, t, d, ihs, jhs, ihe, jhe)
1396 !=============================================================
1397
1398         !-------------------------------------------------------------------------------
1399         !
1400         ! Calculate met from variables updated in HALO after the call to microphysics:
1401         !
1402         ! 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
1403
1404
1405
1406         IMPLICIT NONE
1407
1408         !TYPE(domain), INTENT(IN) :: grid
1409         INTEGER, INTENT(IN)  :: ihs, jhs, ihe, jhe  ! tile +- halo array bounds
1410         REAL,    INTENT(IN)  :: xp, yp, zp          ! particle position
1411         REAL,    INTENT(OUT) :: loc_p, loc_t, loc_d ! local met properties
1412         REAL,    INTENT(IN), DIMENSION(ims:,kms:,jms:):: p,t,d
1413
1414         INTEGER :: i, j, k , kk
1415         REAL :: tmp1, tmp2, zph ! corresponding zp at half levels
1416
1417         k = FLOOR(zp)
1418         j = MIN(MAX(FLOOR(yp-dp05), jhs), jhe)
1419         i = MIN(MAX(FLOOR(xp-dp05), ihs), ihe)
1420
1421         ! physics variables with (i,j) = (1,1) corresponds to (x,y) = (1.5,1.5)
1422         ! T(i=1): x=1.5
1423         ! T(i=2): x=2.5
1424         ! To interpolate T(i=1:2), x=[1.5:2.5[
1425         ! See comment on function hgt2k
1426
1427         ! boundaries ims, ..., ke are declared in module scope
1428
1429         ! p_hydrostatic is on full levels: no vertical (de)staggering needed
1430         loc_p = u_3d_interp( x=xp-dp05, y=yp-dp05, z=zp, &
1431             u_i0j0_bot = p(i,  k, j),   u_i0j1_bot = p(i,  k, j+1),&
1432             u_i1j0_bot = p(i+1,k, j),   u_i1j1_bot = p(i+1,k, j+1),&
1433             u_i0j0_top = p(i,  k+1, j), u_i0j1_top = p(i,  k+1,j+1),&
1434             u_i1j0_top = p(i+1,k+1, j), u_i1j1_top = p(i+1,k+1,j+1))
1435
1436         ! theta moist is on full levels: no vertical (de)staggering needed: th8w
1437         loc_t = u_3d_interp( x=xp-dp05, y=yp-dp05, z=zp, &
1438             u_i0j0_bot = t(i,  k, j),   u_i0j1_bot = t(i,  k, j+1),&
1439             u_i1j0_bot = t(i+1,k, j),   u_i1j1_bot = t(i+1,k, j+1),&
1440             u_i0j0_top = t(i,  k+1, j), u_i0j1_top = t(i,  k+1,j+1),&
1441             u_i1j0_top = t(i+1,k+1, j), u_i1j1_top = t(i+1,k+1,j+1))
1442
1443         ! Variables on half levels: provide the k index corresponding to particle bottom & top levels
1444         ! Is a linear interpolation ok?
1445
1446         ! when zp = 1.0, particle is at the surface on the stag grid,
1447         !                and at k=0.5 on half levels (below the first level)
1448         ! when zp = 1.5, particle is in the middle of the stag grid, and at k=1.0 on half level
1449         !                This is the location where the unstaggered variable value is valid,
1450         !                so the distance between the particle and this point must be 0.
1451
1452         ! Physics variables are all unstaggered in x and y.
1453
1454         zph = zp - dp05
1455         k = MAX(1, FLOOR(zph))
1456
1457         loc_d = u_3d_interp( x=xp-dp05, y=yp-dp05, z=zph, &
1458             u_i0j0_bot = d(i,  k, j),   u_i0j1_bot = d(i,  k, j+1),&
1459             u_i1j0_bot = d(i+1,k, j),   u_i1j1_bot = d(i+1,k, j+1),&
1460             u_i0j0_top = d(i,  k+1, j), u_i0j1_top = d(i,  k+1,j+1),&
1461             u_i1j0_top = d(i+1,k+1, j), u_i1j1_top = d(i+1,k+1,j+1))
1462
1463         loc_d = loc_d * 1000.0_dp  ! Convert from kg/m3 to g/m3
1464
1465
1466     END SUBROUTINE get_local_met
1467
1468 !=============================================================
1469
1470 !=============================================================
1471
1472
1473
1474 !=============================================================
1475     SUBROUTINE firebrand_spotting_em_driver( &
1476 !=============================================================
1477         cf,                    &
1478         grid,                  &
1479         fs_p_id,               &
1480         fs_p_dt,               &
1481         fs_p_x,                &
1482         fs_p_y,                &
1483         fs_p_z,                &
1484         fs_gen_inst,           &
1485         fs_p_mass,             &
1486         fs_p_diam,             &
1487         fs_p_effd,             &
1488         fs_p_temp,             &
1489         fs_p_tvel,             &
1490         fs_last_gen_dt,        &
1491         fs_gen_idmax,          &
1492         fs_fire_ROSdt,         &
1493         fs_fire_area,          &
1494         fs_count_landed_all,   &
1495         fs_count_landed_hist,  &
1496         fs_landing_mask,       &
1497         fs_spotting_lkhd,      &
1498         fs_frac_landed,        &
1499         fs_fuel_spotting_risk, &
1500         fs_count_reset)
1501
1502
1503         !-------------------------------------------------------------------------------
1504         ! Modules
1505         !-------------------------------------------------------------------------------
1506
1507         USE module_domain,            ONLY: domain_get_time_since_sim_start, &
1508                                             domain_clock_get, is_alarm_tstep
1509         USE module_domain_type,       ONLY: HISTORY_ALARM, restart_alarm, AUXHIST23_ALARM
1510         USE module_state_description, ONLY: p_qv, num_moist, param_first_scalar
1511         USE module_utility,           ONLY: WRFU_Alarm !WRFU_Clock,
1512 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
1513         USE module_firebrand_spotting_mpi,  ONLY: fs_mpi_init, &
1514                                             fs_mpi_sendbuff1_int,      &
1515                                             fs_mpi_sendbuff1_real,     &
1516                                             fs_mpi_recvbuff1_int,      &
1517                                             fs_mpi_recvbuff1_real,     &
1518                                             fs_mpi_nothing2send,       &
1519                                             fs_mpi_send2neighbors,     &
1520                                             fs_mpi_sendbuffsize,       &
1521                                             fs_mpi_sendbuff_int,       &
1522                                             fs_mpi_sendbuff_real,      &
1523                                             fs_mpi_recvbuffsize,       &
1524                                             fs_mpi_recvfrompatch_int,  &
1525                                             fs_mpi_recvfrompatch_real, &
1526                                             fs_mpi_sendbuff_real,      &
1527                                             fs_mpi_checkreceive,       &
1528                                             fs_mpi_recv
1529         USE module_firebrand_spotting_mpi,  ONLY: neighbors, my_id, task_id, mpiprocs,  &
1530                                             left_id, right_id, up_id, down_id, &
1531                                             upleft_id, upright_id, downleft_id, downright_id
1532 #endif
1533
1534
1535         IMPLICIT NONE
1536 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
1537         LOGICAL, EXTERNAL :: wrf_dm_on_monitor
1538 #endif
1539         !-------------------------------------------------------------------------------
1540         ! WRF variables
1541         !-------------------------------------------------------------------------------
1542
1543         TYPE(domain),               INTENT(IN)    :: grid ! input data **Note: Intent IN**
1544         TYPE (GRID_config_rec_type),INTENT(IN)    :: cf   ! type available by module host association
1545
1546         !-------------------------------------------------------------------------------
1547         ! Dummy Arguments
1548         !-------------------------------------------------------------------------------
1549
1550         LOGICAL, INTENT(INOUT)  :: fs_count_reset
1551         INTEGER, INTENT(INOUT)  :: fs_last_gen_dt, fs_gen_idmax
1552
1553         ! Firebrand Spotting Particle Properties (fs_p_*)
1554         INTEGER, INTENT(INOUT), DIMENSION(:)    :: fs_p_id, fs_p_dt
1555         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)
1556         REAL,    INTENT(INOUT), DIMENSION(:)    :: fs_p_mass, fs_p_diam, fs_p_effd
1557         REAL,    INTENT(INOUT), DIMENSION(:)    :: fs_p_temp, fs_p_tvel
1558         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
1559         REAL,    INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: fs_count_landed_hist, fs_spotting_lkhd, fs_frac_landed
1560         INTEGER, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: fs_landing_mask, fs_gen_inst
1561         REAL,    INTENT(INOUT), DIMENSION(ifps:ifpe,jfps:jfpe)  :: fs_fire_ROSdt
1562
1563         !-------------------------------------------------------------------------------
1564         ! Parameters
1565         !-------------------------------------------------------------------------------
1566
1567         !-------------------------------------------------------------------------------
1568         ! Local Variables - Arrays
1569         !-------------------------------------------------------------------------------
1570
1571         ! Subroutine Control
1572         LOGICAL, DIMENSION(fs_array_maxsize) :: sparse_mask
1573
1574         TYPE(p_properties), DIMENSION(fs_array_maxsize)   :: fs_p_prop
1575
1576         ! Physics from WRF
1577         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
1578         REAL, DIMENSION(ims:ime, ks:ke, jms:jme):: u, v, rho, th_phy         ! *_phy: half levels
1579         REAL, DIMENSION(ims:ime, ks:ke, jms:jme):: qtot, th8w, p_phy, p8w    ! *8w: full levels
1580         REAL, DIMENSION(ims:ime, jms:jme)       :: msft, w1
1581         REAL, DIMENSION(ks:ke)                  :: znw
1582         REAL, ALLOCATABLE, DIMENSION(:)         :: xout, yout, zout ! firebrand position after advection
1583
1584         !-------------------------------------------------------------------------------
1585         ! Likelihood
1586         !-------------------------------------------------------------------------------
1587
1588         REAL,    ALLOCATABLE, DIMENSION(:,:) :: fuel_spotting_risk_fr, spotthreat_fr ! only need it for history intervals
1589
1590         ! only need it in hist interval - do not hold memory
1591         LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: landing_mask
1592         INTEGER, ALLOCATABLE, DIMENSION(:)   :: landing_cnt
1593         INTEGER, ALLOCATABLE, DIMENSION(:)   :: np_landed_cnt, np_nlanded  ! (MPI: only Master allocates these)
1594         REAL,    ALLOCATABLE, DIMENSION(:)   :: np_landed_risk, np_frac_landed, np_lkhd ! (MPI: only Master allocates these)
1595         REAL,    ALLOCATABLE, DIMENSION(:)   :: landed_risk, recv_spot_lkhd, recv_frac_landed
1596         REAL,    ALLOCATABLE, DIMENSION(:)   :: tmp1
1597         INTEGER :: ndep, ndep_total
1598         !-------------------------------------------------------------------------------
1599
1600         REAL, DIMENSION(ims:ime, jms:jme) :: dt_sum_landed ! sum timestep's deposited fb in 2-D memory array
1601         LOGICAL, SAVE :: hist_flag = .FALSE.  ! time of history output
1602
1603         !-------------------------------------------------------------------------------
1604         ! Release
1605         !-------------------------------------------------------------------------------
1606
1607         INTEGER, ALLOCATABLE, DIMENSION(:)    :: np_Ngrdpts, loc_ij
1608         REAL,    ALLOCATABLE, DIMENSION(:)    :: np_MeanPot
1609         REAL,    ALLOCATABLE, DIMENSION(:,:)  :: RelPot, fcat_factor, fr_rate_dummy
1610         INTEGER, ALLOCATABLE, DIMENSION(:,:)  :: fcat_fr, n_rel_fr2d
1611         LOGICAL, ALLOCATABLE, DIMENSION(:,:)  :: rel_mask
1612
1613         INTEGER :: Ngrdpts, Npts_total, relpts, rankpts
1614         REAL    :: TotPot, MeanPot, PotThr, rel_ratio, cdf_coeff, LimPot
1615
1616         ! Firebrand Release and Transport
1617         TYPE(p_properties), ALLOCATABLE, DIMENSION(:)   :: release_prop
1618         REAL,          ALLOCATABLE, DIMENSION(:)   :: release_i, release_j, release_k
1619         INTEGER,       ALLOCATABLE, DIMENSION(:)   :: release_n
1620
1621         !-------------------------------------------------------------------------------
1622         ! Local Variables - Scalars
1623         !-------------------------------------------------------------------------------
1624
1625         ! Note:
1626         ! Variables for tile, memory, and domain bounds are in module scope & are assigned in firebrand_spotting_init
1627
1628         INTEGER :: kk, pp, ii ! counters
1629         INTEGER :: kk1, kk2, k_end
1630         INTEGER :: active_br, prior_br, nresets, idmax, firepts, nbr_sum
1631         REAL    :: rdx, rdy, dt, sr_xy ! grid length (1/dx) and wrf's time-step
1632
1633         INTEGER :: firebrand_max_life_dt
1634         REAL    :: firebrand_land_hgt, firebrand_gen_maxhgt
1635         LOGICAL :: firebrand_gen_levrand
1636         REAL    :: firebrand_gen_prop_diam, firebrand_gen_prop_effd
1637         REAL    :: firebrand_gen_prop_temp, firebrand_gen_prop_tvel
1638         INTEGER :: levrand_seed
1639
1640         !-------------------------------------------------------------------------------
1641         ! MPI local Variables (some variables are in module scope)
1642         !-------------------------------------------------------------------------------
1643
1644         INTEGER, SAVE :: MasterId = 0
1645         INTEGER, SAVE :: ntiles = 0
1646         LOGICAL, SAVE :: IamMaster = .FALSE.
1647
1648         INTEGER, ALLOCATABLE, DIMENSION(:) :: nbr_id, r_id, r_dt
1649         REAL,    ALLOCATABLE, DIMENSION(:) :: r_x, r_y, r_z
1650
1651         ! MPI receive p_properties
1652         REAL, ALLOCATABLE, DIMENSION(:) :: r_p_m,r_p_d,r_p_e,r_p_t,r_p_v
1653         TYPE(p_properties), ALLOCATABLE, DIMENSION(:) :: fb_mdetv
1654
1655
1656         !-------------------------------------------------------------------------------
1657         ! It was a dummy variable back when fuel moisture was not implemented in wrf-v4.0.1
1658         !-------------------------------------------------------------------------------
1659
1660         REAL,  DIMENSION(ifps:ifpe,jfps:jfpe)  :: fmcg_2df
1661
1662         !===============================================================================
1663         ! Begin
1664         !===============================================================================
1665
1666
1667 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
1668         IF (mpiprocs == 0) &
1669             CALL fs_mpi_init(grid) ! should be called in fspotting_init (?)
1670
1671         ntiles = mpiprocs
1672         IF ( wrf_dm_on_monitor() ) THEN
1673             IamMaster = .TRUE.
1674             MasterId = my_id
1675         ENDIF
1676 #else
1677         IamMaster = .TRUE.  ! initialized with FALSE. It should be TRUE in serial compilation
1678         ntiles = 1 ! used to allocate the arrays for incomming mpi comms. ntiles=1 in serial builds
1679 #endif
1680
1681         hist_flag = .FALSE.
1682         IF ( Is_alarm_tstep(grid%domain_clock, grid%alarms(HISTORY_ALARM)) ) hist_flag = .TRUE.
1683         !IF ( Is_alarm_tstep(grid%domain_clock, grid%alarms(AUXHIST23_ALARM)) ) hist_flag = .TRUE.
1684         IF (hist_flag) &
1685             CALL wrf_debug (wrfdbg, 'SPFire_driver: History output on next tstep')
1686
1687         fs_last_gen_dt = fs_last_gen_dt + 1
1688         active_br = COUNT( fs_p_id > 0 )
1689
1690         IF (grid%itimestep == 1) &
1691             fs_fire_ROSdt(ifps:ifpe,jfps:jfpe) =0.0_dp
1692
1693         !===============================================================================
1694         ! Reset array after writting history file
1695         !===============================================================================
1696
1697         IF (fs_count_reset) THEN ! if true
1698             !CALL wrf_debug (wrfdbg, 'SPFire_lkhd: resetting fs_count_reset!')
1699             fs_count_landed_hist(:,:)  = ZERO_dp
1700             fs_fuel_spotting_risk(:,:) = ZERO_dp
1701             fs_count_reset      = .FALSE.
1702
1703             fs_gen_inst(is:ie,js:je) = 0
1704         ENDIF
1705
1706         fs_landing_mask(ims:ime,jms:jme) = 0
1707
1708         !===============================================================================
1709         ! Assign properties to derived type (used locally to shorten argument list)
1710         !===============================================================================
1711
1712         fs_p_prop%p_mass  = ZERO_dp
1713         fs_p_prop%p_diam  = ZERO_dp
1714         fs_p_prop%p_effd  = ZERO_dp
1715         fs_p_prop%p_temp  = ZERO_dp
1716         fs_p_prop%p_tvel  = ZERO_dp
1717
1718         fs_p_prop(:active_br)%p_mass  = fs_p_mass(:active_br)
1719         fs_p_prop(:active_br)%p_diam  = fs_p_diam(:active_br)
1720         fs_p_prop(:active_br)%p_effd  = fs_p_effd(:active_br)
1721         fs_p_prop(:active_br)%p_temp  = fs_p_temp(:active_br)
1722         fs_p_prop(:active_br)%p_tvel  = fs_p_tvel(:active_br)
1723
1724         !===============================================================================
1725         ! Assign user defined values from namelist
1726         !===============================================================================
1727
1728         firebrand_max_life_dt = cf%fs_firebrand_max_life_dt
1729         firebrand_land_hgt = cf%fs_firebrand_land_hgt
1730         firebrand_gen_maxhgt  = cf%fs_firebrand_gen_maxhgt
1731         firebrand_gen_levrand = cf%fs_firebrand_gen_levrand
1732
1733         firebrand_gen_prop_diam = cf%fs_firebrand_gen_prop_diam
1734         firebrand_gen_prop_effd = cf%fs_firebrand_gen_prop_effd
1735         firebrand_gen_prop_temp = cf%fs_firebrand_gen_prop_temp
1736         firebrand_gen_prop_tvel = cf%fs_firebrand_gen_prop_tvel
1737
1738         firebrand_dens      = cf%fs_firebrand_dens      ! 513000.0_dp     ! density of firebrand (g / m3)
1739         firebrand_dens_char = cf%fs_firebrand_dens_char ! 299000.0_dp     ! density of char (g m-3) [gupta et al 2002; fuel]
1740
1741         !===============================================================================
1742         ! Firebrand Generation
1743         !===============================================================================
1744
1745         !-------------------------------------------------------------------------------
1746         ! Accumulate burn rate between generation intervals
1747         !-------------------------------------------------------------------------------
1748
1749         firepts = 0
1750         fs_fire_ROSdt(ifps:ifpe, jfps:jfpe) = fs_fire_ROSdt(ifps:ifpe, jfps:jfpe) + grid%burnt_area_dt(ifps:ifpe,jfps:jfpe)
1751
1752         !-------------------------------------------------------------------------------
1753         ! Updated fuel moisture from WRF-fire
1754         !-------------------------------------------------------------------------------
1755
1756         fmcg_2df(ifps:ifpe,jfps:jfpe) = grid%fmc_g(ifps:ifpe,jfps:jfpe)
1757
1758         !===============================================================================
1759         ! Firebrand Generation - Release Firebrands - Part 1
1760         !===============================================================================
1761
1762         !-------------------------------------------------------------------------------
1763         ! When a new release is due, calculate fraction of release by fuel category
1764         !-------------------------------------------------------------------------------
1765
1766         IF (fs_last_gen_dt >= fs_gen_dt) THEN
1767
1768             firepts = COUNT(fs_fire_ROSdt(ifps:ifpe,jfps:jfpe) > ZERO_dp)
1769             Ngrdpts = 0 ! Number of active fire points of valid fuel categories
1770
1771             LimPot = 0.000001_dp ! 1e-6
1772             MeanPot = ZERO_dp
1773
1774             IF (firepts > 0) THEN
1775
1776                 !WRITE (msg,'(1(i6,1x))') firepts
1777                 !CALL wrf_debug (wrfdbg, 'SPFire_rel firepts: '//msg)
1778
1779                 !-------------------------------------------------------------------------------
1780                 ! Get fuel category from WRF-FIRE and calculate release potential
1781                 !-------------------------------------------------------------------------------
1782
1783                 ALLOCATE(fcat_fr(ifps:ifpe,jfps:jfpe), fcat_factor(ifps:ifpe,jfps:jfpe), RelPot(ifps:ifpe,jfps:jfpe))!, rel_fcat(ifps:ifpe,jfps:jfpe))
1784                 fcat_fr(:,:)  = 0
1785                 RelPot(:,:) = ZERO_dp
1786                 fcat_factor(:,:) = ZERO_dp
1787
1788                 fcat_fr(ifps:ifpe,jfps:jfpe) = RESHAPE(get_fuel_cat(crosswalk=grid%fuel_crosswalk, &
1789                                                                     cat=grid%nfuel_cat(ifps:ifpe,jfps:jfpe)), &
1790                                                        SHAPE(fcat_fr))
1791
1792                 fcat_factor(ifps:ifpe,jfps:jfpe)= firebrand_gen_factor(fcat_fr(ifps:ifpe,jfps:jfpe)) !RESHAPE, SHAPE(fcat_fr))
1793
1794                 RelPot(ifps:ifpe,jfps:jfpe) = RESHAPE(firebrand_gen_potential(&
1795                                                              fuel_fgi=grid%fgip(ifps:ifpe,jfps:jfpe), &
1796                                                              fuel_mcg=fmcg_2df(ifps:ifpe,jfps:jfpe), &
1797                                                              !fuel_mcg=grid%fmcg2df(ifps:ifpe,jfps:jfpe), &
1798                                                              factor=fcat_factor(ifps:ifpe,jfps:jfpe),&
1799                                                              fire_rate=fs_fire_ROSdt(ifps:ifpe,jfps:jfpe)), &
1800                                                          SHAPE(fcat_fr))
1801
1802                 ! Calculate Release Potential for each burning fuel cat
1803                 !-------------------------------------------------------------------------------
1804
1805                 Ngrdpts= COUNT(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)! RelPot do not count cat with factor of 0
1806
1807                 ! In Large Fires
1808                 ! Try to increase minimum LimPot to remove very low values of Release Potential
1809                 !-------------------------------------------------------------------------------
1810                 IF (Ngrdpts > fs_gen_lim) THEN
1811                     LimPot = 10.0_dp * LimPot ! 1e-5
1812                     Ngrdpts= COUNT(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)
1813
1814                     IF (Ngrdpts > fs_gen_lim) THEN
1815                         LimPot = 10.0_dp * LimPot ! 1e-4
1816                         Ngrdpts= COUNT(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)
1817
1818                     ENDIF
1819                 ENDIF
1820                 ! WRITE (msg,'(1(i8,1x), 1(f14.8,1x))') Ngrdpts, LimPot
1821                 ! CALL wrf_debug (wrfdbg, 'SPFire_rel_lim Ngrdpts, LimP: '//msg)
1822
1823
1824                 WRITE (msg,'(1(i6,1x))') Ngrdpts
1825                 CALL wrf_debug (wrfdbg, 'SPFire_rel Ngrdpts: '//msg)
1826
1827                 ! Set a threshold for Potentials - use only highest values
1828                 !-------------------------------------------------------------------------------
1829                 IF (Ngrdpts > fs_gen_lim) THEN
1830
1831                     ! Find threshold value to select fs_gen_lim points
1832                     !-------------------------------------------------------------------------------
1833
1834                     ALLOCATE(tmp1(Ngrdpts))
1835                     tmp1(1:Ngrdpts) = PACK(RelPot(ifps:ifpe,jfps:jfpe), &
1836                                 mask=(RelPot(ifps:ifpe,jfps:jfpe) > LimPot))
1837
1838                     ! Find a limit value that yields an array size close to fs_gen_lim
1839                     LimPot = order_val(tmp1, ORD=fs_gen_lim)
1840                     !maxpot = MAXVAL(tmp1)
1841                     DEALLOCATE(tmp1)
1842
1843                     ! WRITE (msg,'(3(i8,1x), 2(f14.8,1x))') Ngrdpts, COUNT(RelPot > LimPot), SIZE(tmp1), LimPot, maxpot
1844                     ! CALL wrf_debug (wrfdbg, 'SPFire_rel_lim Ngrdpts, Cnt>Lim, LimP, max: '//msg)
1845
1846                     ! Find the median rank among remaining points
1847                     !-------------------------------------------------------------------------------
1848
1849                     ALLOCATE(tmp1(COUNT(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)))
1850                     tmp1 = PACK(RelPot(ifps:ifpe,jfps:jfpe), &
1851                                 mask=(RelPot(ifps:ifpe,jfps:jfpe) > LimPot))
1852
1853                     MeanPot = order_val(tmp1, ORD=INT(REAL(SIZE(tmp1))/2.0))
1854                     !maxpot = MAXVAL(tmp1)
1855                     DEALLOCATE(tmp1)
1856
1857                     ! WRITE (msg,'(4(i8,1x), 3(f14.8,1x))') Ngrdpts, &
1858                     !     COUNT(RelPot > MeanPot), SIZE(tmp1), INT(REAL(SIZE(tmp1))/2.0), &
1859                     !     MeanPot, MINVAL(tmp1), maxpot
1860                     ! CALL wrf_debug (wrfdbg, 'SPFire_rel_lim Ngrdpts, Cnt>Mean, MPot, min, max: '//msg)
1861
1862                     Ngrdpts= COUNT(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)
1863
1864                 !-------------------------------------------------------------------------------
1865
1866                 ELSEIF (Ngrdpts > 0) THEN
1867
1868                     ! Find the median rank among points
1869                     !-------------------------------------------------------------------------------
1870                     MeanPot = order_val(PACK(RelPot(ifps:ifpe,jfps:jfpe), &
1871                                       mask=(RelPot(ifps:ifpe,jfps:jfpe) > LimPot)), &
1872                                       ORD=INT(dp05*Ngrdpts))
1873
1874                     ! WRITE (msg,'(2(i8,1x), 2(f14.8,1x))') Ngrdpts, COUNT(RelPot > MeanPot), MeanPot, LimPot
1875                     ! CALL wrf_debug (wrfdbg, 'SPFire_rel_lim Ngrdpts, Cnt>Mean, MPot, LimP: '//msg)
1876
1877                 ENDIF
1878
1879                 ! Prepare for MPI
1880                 !-------------------------------------------------------------------------------
1881
1882                 IF (Ngrdpts == 0) firepts = 0
1883
1884                 !-------------------------------------------------------------------------------
1885
1886             ENDIF ! (firepts > 0)
1887             !-------------------------------------------------------------------------------
1888
1889 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
1890             IF (.NOT. (IamMaster)) THEN
1891
1892                 ! nbr > 0 indicates there's fire on tile and RelPot > 0, not the array size to receive in next mpi comm
1893                 ! if nbr > 0, next array will be one (array: MeanPot)
1894                 CALL fs_mpi_sendbuff1_int(sendto=MasterId, nbr=Ngrdpts)
1895
1896                 IF (Ngrdpts > 0) THEN
1897                     CALL fs_mpi_sendbuff1_real(sendto=MasterId, nbr=MeanPot)
1898                     ! WRITE (msg,'(1(i6,1x), (f14.6,1x))') Ngrdpts, MeanPot
1899                     ! CALL wrf_debug (wrfdbg, 'SPFire_rel mpi_send Npts,Meanpot: '//msg)
1900                 ENDIF
1901             ENDIF
1902 #endif
1903
1904             !===============================================================================
1905             ! Firebrand Generation - Release Firebrands - Part 2
1906             !===============================================================================
1907
1908             !-------------------------------------------------------------------------------
1909             ! MPI Master: Calculate Release Field (fire gridpt cnt by fuel category)
1910             !-------------------------------------------------------------------------------
1911
1912             IF (IamMaster) THEN
1913
1914                 CALL wrf_debug (wrfdbg, 'SPFire_rel_master --------------------------------------------------- ')
1915                 !-------------------------------------------------------------------------------
1916                 ! Declare arrays to receive data
1917
1918                 Npts_total = 0
1919                 PotThr = ZERO_dp
1920                 TotPot = ZERO_dp
1921                 rel_ratio = ZERO_dp
1922                 cdf_coeff = ZERO_dp
1923
1924
1925                 ALLOCATE(np_Ngrdpts(ntiles)) ! Number of fuel cat actively burning in each mpiproc
1926                 ALLOCATE(np_MeanPot(ntiles))
1927                 np_MeanPot(:) = ZERO_dp
1928                 np_Ngrdpts(:) = 0
1929
1930                 np_Ngrdpts(1) = Ngrdpts ! Array index for Master is 1
1931                 np_MeanPot(1) = MeanPot
1932
1933                 !-------------------------------------------------------------------------------
1934                 ! Receive release potential data from procs
1935
1936 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
1937
1938                 DO kk=2, mpiprocs ! kk=1 is Master, MPI-IDs begin at 0
1939                     np_Ngrdpts(kk) = fs_mpi_recvbuff1_int(fromid=kk-1)  ! Comm receive: Number of pts with RelPot>0
1940
1941                     IF (np_Ngrdpts(kk) > 0) THEN
1942                         np_MeanPot(kk) = fs_mpi_recvbuff1_real(fromid=kk-1)
1943                         ! WRITE(msg, '(3(i8,1x), 2(f14.6,1x))') my_id, kk-1, np_Ngrdpts(kk), np_Meanpot(kk), TotPot
1944                         ! CALL wrf_debug (wrfdbg, 'SPFire_rel master proc Ngrdpts, MeanPot, TotPot: '//msg)
1945                     ENDIF
1946                 ENDDO
1947 #endif
1948
1949                 DO kk=1, ntiles ! at least 1 for serial compile
1950                    TotPot = np_MeanPot(kk) * REAL(np_Ngrdpts(kk)) + TotPot
1951                 ENDDO
1952
1953                 !-------------------------------------------------------------------------------
1954                 ! Communicate release sources between Master and procs
1955
1956                 Npts_total = SUM( np_Ngrdpts )
1957                 IF ( Npts_total > 0) THEN
1958
1959                     TotPot = TotPot/(REAL(Npts_total))
1960
1961                     !-------------------------------------------------------------------------------
1962                     ! Calculate potential threshold for grid pts
1963
1964                     ! All elements fit in the array
1965                     IF (Npts_total <= fs_gen_lim) THEN
1966                         PotThr = ZERO_dp
1967
1968                     ! Estimate threshold with a linear regression
1969                     ELSE ! IF (Npts_total > fs_gen_lim)
1970                         rel_ratio = REAL(fs_gen_lim)/(REAL(Npts_total))
1971                         cdf_coeff = (dp1 - rel_ratio)/rel_ratio
1972                         PotThr = cdf_coeff * TotPot
1973                     ENDIF
1974
1975                     ! WRITE(msg, '(2(i8,1x), 4(f12.6,1x))') fs_gen_lim, Npts_total, PotThr, rel_ratio, cdf_coeff, TotPot
1976                     ! CALL wrf_debug (wrfdbg, 'SPFire_rel_master max, Npts, PotThr: '//msg)
1977
1978                     !-------------------------------------------------------------------------------
1979                     ! Send dat to procs
1980
1981 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
1982
1983                     DO kk=2, mpiprocs ! kk=1 is Master
1984                         IF (np_Ngrdpts(kk) > 0) THEN
1985                             CALL fs_mpi_sendbuff1_real(sendto=kk-1, nbr=PotThr)
1986                             ! WRITE(msg, '(1(i8,1x), 1(f12.6,1x))') kk-1, PotThr
1987                             ! CALL wrf_debug (wrfdbg, 'SPFire_rel_master mpi_send PotThr: '//msg)
1988                         ENDIF
1989                     ENDDO
1990 #endif
1991                 ENDIF ! npts_total > 0
1992                 !-------------------------------------------------------------------------------
1993
1994             ENDIF ! IamMaster
1995
1996             !-------------------------------------------------------------------------------
1997
1998             !===============================================================================
1999             ! Firebrand Generation - Release Firebrands - Part 3
2000             !===============================================================================
2001
2002             IF (Ngrdpts > 0) THEN
2003
2004                 !-------------------------------------------------------------------------------
2005                 ! Receive Potential Threshold from Master
2006
2007 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
2008                 IF (.NOT.(IamMaster)) THEN
2009                     PotThr = ZERO_dp
2010                     PotThr = fs_mpi_recvbuff1_real(fromid=MasterId)
2011                     ! WRITE(msg, '((i4,1x), (f14.6,1x))')  my_id, PotThr
2012                     ! CALL wrf_debug (wrfdbg, 'SPFire_rel receive potential threshold: '//msg)
2013                 ENDIF
2014 #endif
2015
2016                 !-------------------------------------------------------------------------------
2017                 ! Assign release from the various fuel cat to corresponding gridpts
2018
2019                 !CALL wrf_debug (wrfdbg, 'SPFire_rel2fr relpts>0 --------------------------------------------- ')
2020
2021                 ALLOCATE(n_rel_fr2d(ifps:ifpe, jfps:jfpe))
2022                 n_rel_fr2d(ifps:ifpe, jfps:jfpe) = 0
2023
2024                 IF (PotThr > ZERO_dp) &
2025                     n_rel_fr2d(ifps:ifpe,jfps:jfpe) = n_rel_fr2d(ifps:ifpe,jfps:jfpe)  + &
2026                                                       UNPACK(SPREAD(1, DIM=1, NCOPIES=Ngrdpts), &
2027                                                                     mask=(RelPot(ifps:ifpe,jfps:jfpe) >= PotThr), FIELD=0)
2028
2029                 IF (PotThr == ZERO_dp) &
2030                     n_rel_fr2d(ifps:ifpe,jfps:jfpe) = n_rel_fr2d(ifps:ifpe,jfps:jfpe)  + &
2031                                                       UNPACK(SPREAD(1, DIM=1, NCOPIES=Ngrdpts), &
2032                                                                     mask=(fs_fire_ROSdt(ifps:ifpe,jfps:jfpe) > ZERO_dp), FIELD=0)
2033
2034                 !-------------------------------------------------------------------------------
2035                 ! update relpts to represent value assigned to gridpts
2036
2037                 relpts = COUNT(n_rel_fr2d(ifps:ifpe,jfps:jfpe)>0) ! Limits to 1 release per refined gridpt, at various heights
2038
2039                 WRITE(msg, '(1(i8,1x))') relpts
2040                 CALL wrf_debug (wrfdbg, 'SPFire_rel2fr relpts: '//msg)
2041
2042                 !-------------------------------------------------------------------------------
2043                 ! Convert release locations to met grid
2044
2045                 ALLOCATE(release_i(relpts*fs_gen_levels), &
2046                          release_j(relpts*fs_gen_levels), &
2047                          release_k(relpts*fs_gen_levels), &
2048                          release_n(relpts))
2049
2050                 release_i(:)  = 0.0_dp
2051                 release_j(:)  = 0.0_dp
2052                 release_k(:)  = 0.0_dp
2053                 release_n(:)  = 0
2054
2055                 ! Does not allow more than one release from the same refined gridpt
2056                 CALL releaseijk2atm(nij_2d=n_rel_fr2d(ifps:ifpe,jfps:jfpe), &
2057                                     fcat=fcat_fr(ifps:ifpe,jfps:jfpe), &
2058                                     maxhgt_usr = firebrand_gen_maxhgt, &
2059                                     sr_x=grid%sr_x, sr_y=grid%sr_y, &
2060                                     pi=release_i(1:relpts), &
2061                                     pj=release_j(1:relpts), &
2062                                     pk=release_k(1:relpts), & ! max height
2063                                     nij=release_n)
2064
2065
2066                 ! Adjust position to corresponding refined grid center
2067                 release_i(1:relpts) = release_i(1:relpts) + dp05/grid%sr_x
2068                 release_j(1:relpts) = release_j(1:relpts) + dp05/grid%sr_y
2069
2070                 DO kk=1,relpts
2071
2072                     fs_gen_inst(INT(release_i(kk)),INT(release_j(kk))) = &
2073                         fs_gen_inst(INT(release_i(kk)),INT(release_j(kk))) + 1
2074
2075                 ENDDO
2076
2077                 CALL wrf_debug (wrfdbg, 'SPFire_rel2fr prep_release_hgt ---------------------------------------------')
2078                 ALLOCATE(release_prop(relpts*fs_gen_levels))
2079
2080                 ! Create derived type
2081                 release_prop(1:relpts) = firebrand_property([firebrand_gen_prop_diam, &
2082                                                              firebrand_gen_prop_effd, &
2083                                                              firebrand_gen_prop_temp, &
2084                                                              firebrand_gen_prop_tvel])
2085
2086                 levrand_seed = 0 ! means no random levels
2087                 IF (firebrand_gen_levrand) levrand_seed = cf%fs_firebrand_gen_levrand_seed + 1 ! in case user sets it to 0
2088                 ! include more releases to span over various heights
2089                 CALL prep_release_hgt( &
2090                                   release_i = release_i, &
2091                                   release_j = release_j, &
2092                                   release_k = release_k, & ! max height
2093                                   release_n = release_n, &
2094                                   release_prop=release_prop, &
2095                                   levrand_seed = levrand_seed * grid%itimestep) ! force seed to vary in time
2096
2097                 prior_br = active_br
2098                 idmax = fs_gen_idmax
2099
2100                 CALL wrf_debug (wrfdbg, 'SPFire_driver_call_generate_firebrands ----------------------------------------- ')
2101
2102                 IF (active_br + COUNT( INT(release_i) > 0 )  > fs_array_maxsize) &
2103                     CALL wrf_debug (wrfdbg, 'SPFire_driver_release: brand array is full, cannot release new brands')
2104
2105                 CALL generate_firebrands(&
2106                     fs_p_id = fs_p_id, &
2107                     fs_p_dt = fs_p_dt, &
2108                     fs_p_z  = fs_p_z,  &
2109                     fs_p_x  = fs_p_x,  &
2110                     fs_p_y  = fs_p_y,  &
2111                     fs_p_prop   = fs_p_prop,    &
2112                     fs_gen_idmax= fs_gen_idmax, &
2113                     release_prop= release_prop, &
2114                     release_i = release_i, &
2115                     release_j = release_j, &
2116                     release_k = release_k, &
2117                     active_br = active_br)
2118
2119                 IF (active_br /= COUNT( fs_p_id > 0 ) ) CALL wrf_error_fatal('SPFire_driver: Active brands do not match!')
2120
2121             ENDIF ! Ngrdpts > 0
2122
2123             fs_fire_ROSdt(:, :) = ZERO_dp
2124             relpts = 0
2125             fs_last_gen_dt = 0 ! Reset release timestep counter
2126             IF( ALLOCATED(fcat_fr))     DEALLOCATE(fcat_fr)
2127             IF( ALLOCATED(fcat_factor)) DEALLOCATE(fcat_factor)
2128             IF( ALLOCATED(RelPot))      DEALLOCATE(RelPot)
2129             IF( ALLOCATED(np_Ngrdpts))  DEALLOCATE(np_Ngrdpts)
2130             IF( ALLOCATED(np_MeanPot))  DEALLOCATE(np_MeanPot)
2131             IF( ALLOCATED(n_rel_fr2d))  DEALLOCATE(n_rel_fr2d)
2132             IF( ALLOCATED(release_prop)) DEALLOCATE(release_prop)
2133
2134         ENDIF ! IF (fs_last_gen_dt >= fs_gen_dt) THEN
2135
2136         !===============================================================================
2137         ! End of Releases
2138         !===============================================================================
2139
2140         !-------------------------------------------------------------------------------
2141         ! Nothing to transport: MPI Wait
2142         !-------------------------------------------------------------------------------
2143
2144 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
2145         IF ((active_br == 0))  THEN
2146
2147             CALL fs_mpi_nothing2send(sendto=left_id)
2148             CALL fs_mpi_nothing2send(sendto=right_id)
2149             CALL fs_mpi_nothing2send(sendto=up_id)
2150             CALL fs_mpi_nothing2send(sendto=down_id)
2151
2152             CALL fs_mpi_nothing2send(sendto=upleft_id)
2153             CALL fs_mpi_nothing2send(sendto=upright_id)
2154             CALL fs_mpi_nothing2send(sendto=downleft_id)
2155             CALL fs_mpi_nothing2send(sendto=downright_id)
2156             !CALL wrf_debug (wrfdbg, 'SPFire_driver: no fires or active brands in patch. MPI signal sent.')
2157
2158         ENDIF
2159 #endif
2160         !===============================================================================
2161         !===============================================================================
2162
2163
2164
2165         !===============================================================================
2166         ! Transport active firebrands
2167         !===============================================================================
2168
2169         IF (active_br > 0) THEN
2170
2171             WRITE (msg,'(i8)') active_br
2172             CALL wrf_debug (wrfdbg, 'SPFire_driver: Active brands: '// msg )
2173
2174             !-------------------------------------------------------------------------------
2175             ! Fetch wrf variables for advection
2176             !-------------------------------------------------------------------------------
2177
2178             ! Use arrays computed for the physics scheme on half (_phy) and full (8w) levels
2179             ! will skip calc at model top - particle will not live during journey (module_big_step_utilities_em.F)
2180
2181             !-------------------------------------------------------------------------------
2182             ! Local met for fb physics: based on variables calculated in module_big_step_utilities_em.F
2183             !-------------------------------------------------------------------------------
2184
2185             !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
2186             !-------------------------------------------------------------------------------
2187             ! P, Th, Rho, hgt at full levels
2188
2189             ! Direct assignments from existing Halos:
2190
2191             ! Half Levels
2192             k_end = MIN( ke, kde-1 )
2193             ks = 1
2194
2195             p_phy(ims:, ks:, jms:) = ZERO_dp
2196             th_phy(ims:, ks:, jms:) = ZERO_dp
2197
2198
2199             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)
2200             th_phy(ims:, ks:, jms:) = grid%t_2(ims:ime,ks:k_end,jms:jme) + t0
2201
2202             !hypsometric_opt 2 (default) computes pressure in the model by using an alternative method
2203             ! Temperature is moist theta by default (v4+, https://github.com/wrf-model/WRF/releases)
2204             IF ( ( grid%use_theta_m == 1 ) .AND. (P_Qv >= PARAM_FIRST_SCALAR) ) &
2205                 th_phy = th_phy/(dp1 + Rv/Rd * grid%moist(ims:ime,ks:k_end,jms:jme,P_QV))
2206
2207             ! Full Levels
2208             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
2209
2210             u(ims:, ks:, jms:) = grid%u_2(ims:ime,ks:ke,jms:jme)
2211             v(ims:, ks:, jms:) = grid%v_2(ims:ime,ks:ke,jms:jme)
2212             w(ims:, ks:, jms:) = grid%w_2(ims:ime,ks:ke,jms:jme)
2213
2214             msft(ims:,   jms:) = grid%msftx(ims:ime,    jms:jme)
2215
2216             znw = grid%znw
2217             DT  = grid%dt
2218             rdx = grid%rdx ! inverse grid length [m]
2219             rdy = grid%rdy
2220
2221             !-------------------------------------------------------------------------------
2222             ! Calculated variables
2223
2224             p_hyd(ims:, ks:, jms:) = ZERO_dp
2225             dz8w(ims:, ks:, jms:) = ZERO_dp
2226             th8w(ims:, ks:, jms:) = ZERO_dp
2227             p8w(ims:, ks:, jms:) = ZERO_dp
2228             qtot(ims:, ks:, jms:) = ZERO_dp
2229             rho(ims:, ks:, jms:) = ZERO_dp
2230             z(ims:, ks:, jms:) = ZERO_dp
2231
2232             dz8w(ims:,ks:ke-1,jms:) = z_at_w(:,2:,:)-z_at_w(:,:ke-1,:) ! z_at_w: Height above the ground at interfaces
2233             z(ims:,ks:ke-1,jms:) = dp05 * (z_at_w(:,2:,:) + z_at_w(:,:ke-1,:) )
2234
2235             qtot(ims:, ks:k_end, jms:) = SUM(grid%moist(ims:ime,ks:k_end,jms:jme,PARAM_FIRST_SCALAR:num_moist),DIM=4)
2236             p_hyd(ims:, ke, jms:) = grid%p_top
2237             DO kk = ke-1, ks, -1
2238                 p_hyd(:,kk,:) = p_hyd(:,kk+1,:) - &
2239                     (dp1 + qtot(:,kk,:)) * (grid%c1h(kk) * grid%muts(:,:) + grid%c2h(kk)) * grid%dnw(kk)
2240             ENDDO
2241
2242             rho(ims:,ks:k_end,jms:) = (dp1/&
2243                                         (grid%al(ims:ime,ks:k_end,jms:jme)+grid%alb(ims:ime,ks:k_end,jms:jme)))&
2244                                         *(dp1+grid%moist(ims:ime,ks:k_end,jms:jme,P_QV))
2245             DO kk = 2, k_end
2246                 th8w(ims:,kk,jms:) = grid%fnm(kk) * th_phy(:,kk,:) +  th_phy(:,kk-1,:) * grid%fnp(kk)
2247                 p8w(ims:,kk,jms:)  = grid%fnm(kk) *  p_phy(:,kk,:) +   p_phy(:,kk-1,:) * grid%fnp(kk)
2248             ENDDO
2249
2250             w1 = (z(:,1,:) - z(:,2,:))
2251             WHERE(NINT(w1*100.0) == 0 ) w1 = dp1
2252             w1 = (z_at_w(:,1,:) - z(:,2,:)) / w1
2253
2254             th8w(:,1,:) = w1 * th_phy(:,1,:) + (dp1 - w1) * th_phy(:,2,:) ! interp at bottom only - ignore top
2255             p8w(:,1,:)  = w1 *  p_phy(:,1,:) + (dp1 - w1) *  p_phy(:,2,:) ! interp at bottom only - ignore top
2256
2257             DO kk = 1, ke
2258                 z_at_w(:,kk,:)= z_at_w(:,kk,:) - z_at_w(:,1,:)
2259             ENDDO
2260
2261             !-------------------------------------------------------------------------------
2262             !-------------------------------------------------------------------------------
2263
2264             !-------------------------------------------------------------------------------
2265             ! Advect particle
2266             !-------------------------------------------------------------------------------
2267
2268             CALL wrf_debug (wrfdbg, 'SPFire_driver begin transport ---------------------------------------------------- ')
2269
2270             ALLOCATE(xout(active_br), &
2271                      yout(active_br), &
2272                      zout(active_br))
2273
2274             CALL advect_xyz_m(grid=grid,           &
2275                             xp=  fs_p_x(:active_br), &
2276                             yp=  fs_p_y(:active_br), &
2277                             hgt= fs_p_z(:active_br), & ! height[m] is the fixed position reference between wrf time steps
2278                             dtp= fs_p_dt(:active_br),&
2279                             idp= fs_p_id(:active_br),&
2280                             znw= znw,      &
2281                             mf = rdx*msft, &
2282                             z_at_w=z_at_w, &
2283                             u= u,     &
2284                             v= v,     &
2285                             w= w,     &
2286                             dt=dt,    &
2287                             ims= ims,  &
2288                             jms= jms,  &
2289                             kms= ks,  &
2290                             phyd= p_hyd, & ! air pressure (Pa)
2291                             thet= th8w, &  ! temperature (K)
2292                             rho = rho,  &  ! density (kg/m3)
2293                             xout= xout, &
2294                             yout= yout, &
2295                             zout= zout, &
2296                             fs_p_prop=fs_p_prop(:active_br), &
2297                             land_hgt = firebrand_land_hgt,   &
2298                             start_mom3d_dt = cf%fs_firebrand_gen_mom3d_dt, &
2299                             msg = msg)
2300
2301
2302             IF (LEN(TRIM(msg)) > 1) &
2303                 CALL wrf_debug (wrfdbg, 'SPFire_transp: '//msg)
2304
2305             !-------------------------------------------------------------------------------
2306             ! Update array for output before resetting deposited indices to zero
2307             !-------------------------------------------------------------------------------
2308
2309             fs_p_x(:active_br) = xout
2310             fs_p_y(:active_br) = yout
2311             fs_p_z(:active_br) = zout
2312             fs_p_dt(:active_br)= fs_p_dt(:active_br) +1
2313
2314             fs_p_mass(:active_br) = fs_p_prop(:active_br)%p_mass
2315             fs_p_diam(:active_br) = fs_p_prop(:active_br)%p_diam
2316             fs_p_effd(:active_br) = fs_p_prop(:active_br)%p_effd
2317             fs_p_temp(:active_br) = fs_p_prop(:active_br)%p_temp
2318             fs_p_tvel(:active_br) = fs_p_prop(:active_br)%p_tvel
2319
2320             CALL wrf_debug (wrfdbg, 'SPFire_driver end transport   ---------------------------------------------------- ')
2321
2322             !=============================================================
2323             ! Transport done :)
2324             !=============================================================
2325
2326
2327             !===============================================================================
2328             ! Remove particles
2329             !===============================================================================
2330
2331             !-------------------------------------------------------------------------------
2332             ! Set NaN properties for particle removal
2333             !-------------------------------------------------------------------------------
2334             sparse_mask = .FALSE.
2335
2336             WHERE (IEEE_IS_NAN(fs_p_mass) .OR. IEEE_IS_NAN(fs_p_diam) .OR. &
2337                    IEEE_IS_NAN(fs_p_effd) .OR. IEEE_IS_NAN(fs_p_temp) .OR. &
2338                    IEEE_IS_NAN(fs_p_tvel)) fs_p_effd = ZERO_dp ! ediam=zero will be removed in remove_br
2339
2340
2341             !===============================================================================
2342             ! Remove particles - from burnout or out of domain
2343             !===============================================================================
2344
2345             CALL remove_br(& !active_br= active_br, &
2346                            fs_p_id    = fs_p_id,     &
2347                            fs_p_x     = fs_p_x,      &
2348                            fs_p_y     = fs_p_y,      &
2349                            fs_p_z     = fs_p_z,      &
2350                            fs_p_dt    = fs_p_dt,     &
2351                            fs_p_effd  = fs_p_effd,   &
2352                            cnt        = nresets,     &
2353                            max_life_dt= firebrand_max_life_dt, &
2354                            land_hgt   = firebrand_land_hgt)
2355             WRITE (msg,'(i12)') nresets
2356             CALL wrf_debug (wrfdbg, 'SPFire_driver remove br: '//msg)
2357
2358
2359             !===============================================================================
2360             ! Spotfires: Deposited Brands
2361             !===============================================================================
2362
2363             dt_sum_landed(ims:, jms:) = ZERO_dp
2364             CALL deposit_br(& !active_br= active_br, &
2365                            fs_p_id    = fs_p_id,      &
2366                            fs_p_x     = fs_p_x,       &
2367                            fs_p_y     = fs_p_y,       &
2368                            fs_p_z     = fs_p_z,       &
2369                            sum_fbrand = dt_sum_landed,&
2370                            land_hgt   = firebrand_land_hgt)
2371
2372             ! The mask is based on floor(fs_p_x,y) so deposits should not exist at ie or je
2373             fs_count_landed_all(is:ie, js:je)  = fs_count_landed_all(is:ie, js:je)  + dt_sum_landed(is:ie, js:je)
2374             fs_count_landed_hist(is:ie, js:je) = fs_count_landed_hist(is:ie, js:je) + dt_sum_landed(is:ie, js:je)
2375             WRITE (msg,'(i12)') COUNT(dt_sum_landed(is:ie, js:je) > 0)
2376             CALL wrf_debug (wrfdbg, 'SPFire_driver deposit br: '//msg)
2377
2378
2379             !===============================================================================
2380             !===============================================================================
2381
2382             !******************************************************************
2383             !******************************************************************
2384             !*                                                                *
2385             !*                  MPI Send firebrands                   *
2386             !*                                                                *
2387             !******************************************************************
2388             !******************************************************************
2389
2390
2391             !-------------------------------------------------------------------------------
2392             !-------------------------------------------------------------------------------
2393
2394             ! Re-PACK before incoming MPI
2395
2396             !===============================================================================
2397             ! MPI: Send Brands (out of this patch)
2398             !===============================================================================
2399
2400             ! CALL wrf_debug (wrfdbg, 'SPFire_driver mpi transfer particles')
2401
2402             ! Send / Receive arrays:
2403             !     3 REAL arrays, 2 INTEGER arrays: fs_p_x, fs_p_y, fs_p_z, fs_p_id, fs_p_dt:
2404
2405             !   |____|____|____|
2406             !   |    | U  |    |
2407             !   |____|____|____|je
2408             !   |  L | Me | R  |
2409             !   |____|____|____|js
2410             !   |    | D  |    |
2411             !   |____|____|____|
2412             !        is   ie
2413
2414             ! Diagonal???
2415
2416             !-------------------------------------------------------------------------------
2417             ! Send any particle?
2418             !-------------------------------------------------------------------------------
2419
2420             sparse_mask = .FALSE.
2421             sparse_mask = (fs_p_id > 0 .AND. &
2422                           (FLOOR(fs_p_x) < is .OR. & ! LEFT
2423                            FLOOR(fs_p_x) > ie .OR. & ! RIGHT
2424                            FLOOR(fs_p_y) > je .OR. & ! UP
2425                            FLOOR(fs_p_y) < js ))     ! DONW
2426
2427
2428             WRITE (msg,'(i6)') COUNT(sparse_mask)
2429             CALL wrf_debug (wrfdbg, 'SPFire_driver mpi send away: '//msg)
2430
2431             !-------------------------------------------------------------------------------
2432             ! interrupt before MPI
2433             ! IF (COUNT(sparse_mask) > 0) CALL wrf_error_fatal( "MPI called")
2434
2435 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
2436             CALL fs_mpi_send2neighbors(task_id=task_id,&
2437                     mask=sparse_mask,&
2438                     p_x  = fs_p_x,   &
2439                     p_y  = fs_p_y,   &
2440                     p_z  = fs_p_z,   &
2441                     p_id = fs_p_id,  &
2442                     p_dt = fs_p_dt,  &
2443                     fs_p_m = fs_p_mass, &
2444                     fs_p_d = fs_p_diam, &
2445                     fs_p_e = fs_p_effd, &
2446                     fs_p_t = fs_p_temp, &
2447                     fs_p_v = fs_p_tvel)
2448 #endif
2449             !-------------------------------------------------------------------------------
2450             ! Reset
2451             !-------------------------------------------------------------------------------
2452             ! CALL wrf_debug (wrfdbg, 'SPFire_driver reset particles sent away')
2453
2454             ! In serial compilation, firebrands outside the domain limits will be removed
2455             WHERE(sparse_mask) fs_p_id= 0
2456
2457         !===============================================================================
2458         ! Last step for IF active_br > 0
2459         !===============================================================================
2460
2461             ! WRITE (msg,'(2(i8,1x))') active_br, COUNT( fs_p_id > 0 )
2462             ! CALL wrf_debug (wrfdbg, 'SPFire_driver pack BEFORE: active_br >>> '// msg)
2463
2464             !-------------------------------------------------------------------------------
2465             ! Repack all zeros: PACK where mask is TRUE
2466             !-------------------------------------------------------------------------------
2467             sparse_mask = .FALSE.
2468             sparse_mask = (fs_p_id>0)
2469             fs_p_x  = PACK(fs_p_x,  sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2470             fs_p_y  = PACK(fs_p_y,  sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2471             fs_p_z  = PACK(fs_p_z,  sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2472             fs_p_id = PACK(fs_p_id, sparse_mask, SPREAD(0,1,fs_array_maxsize) )
2473             fs_p_dt = PACK(fs_p_dt, sparse_mask, SPREAD(0,1,fs_array_maxsize) )
2474             fs_p_mass  = PACK(fs_p_mass , sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2475             fs_p_diam  = PACK(fs_p_diam , sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2476             fs_p_effd = PACK(fs_p_effd, sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2477             fs_p_temp  = PACK(fs_p_temp , sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2478             fs_p_tvel  = PACK(fs_p_tvel , sparse_mask, SPREAD(ZERO_dp,1,fs_array_maxsize) )
2479
2480             active_br = COUNT( fs_p_id > 0 )
2481
2482             WRITE (msg,'(1(i8,1x))') active_br
2483             CALL wrf_debug (wrfdbg, 'SPFire_driver pack AFTER mpi_send2neighbors active_br >>> '// msg)
2484
2485
2486         ENDIF
2487
2488         !===============================================================================
2489         ! Active_br work done ;)
2490         !===============================================================================
2491         IF (active_br /= COUNT( fs_p_id > 0 ) ) CALL wrf_error_fatal('SPFire_driver: Active brands do not match!')
2492
2493         !===============================================================================
2494         !===============================================================================
2495
2496
2497
2498         !===============================================================================
2499         ! Calculate Spotfire Likelihood at history output time
2500         !===============================================================================
2501
2502         ! Registry
2503         ! fs_count_landed_all
2504         ! fs_count_landed_hist
2505         ! fs_frac_landed
2506         ! fs_spotting_lkhd
2507         ! fs_count_reset
2508         ! fs_fire_area
2509
2510         ! store array of data for each proc: Master only
2511         ! np_nlanded      [np]              ! ALLOCATABLE
2512         ! np_landed_cnt   [SUM(np_nlanded)] ! ALLOCATABLE
2513         ! np_landed_risk  [SUM(np_nlanded)] ! ALLOCATABLE
2514         ! np_frac_landed  [SUM(np_nlanded)] ! ALLOCATABLE
2515         ! np_lkhd         [SUM(np_nlanded)] ! ALLOCATABLE
2516         !
2517
2518         ! store array of local data
2519         ! landing_mask    ! ALLOCATABLE
2520         ! landing_cnt     ! ALLOCATABLE (mpi send)
2521         ! landed_risk     ! ALLOCATABLE (mpi send)
2522         ! recv_spot_lkhd  ! ALLOCATABLE (mpi recv)
2523         ! recv_frac_landed ! ALLOCATABLE (mpi recv)
2524
2525         !===============================================================================
2526         ! Check history interval (use diag_flag defined in solve_em for microphysics)
2527         !===============================================================================
2528
2529         IF ((hist_flag) .AND. (grid%itimestep > 1)) THEN
2530
2531             fs_count_reset = .TRUE. ! true
2532             fs_frac_landed(ims:ime, jms:jme) = ZERO_dp
2533             fs_spotting_lkhd(ims:ime, jms:jme) = ZERO_dp
2534
2535             ! Was there any firebrand deposit? *Note that fire may exist on a different tile
2536             !-------------------------------------------------------------------------------
2537
2538             ALLOCATE(landing_mask(ims:ime, jms:jme)) ! Allocate memory but only mask this tile
2539             landing_mask(ims:ime, jms:jme) = .FALSE.
2540             landing_mask(is:ie, js:je) = (fs_count_landed_hist(is:ie, js:je) > 0)
2541             ndep = COUNT(landing_mask)
2542
2543             ! convert fire_area from fire grid to atm grid
2544             !-------------------------------------------------------------------------------
2545             fs_fire_area(is:ie,js:je)  = ZERO_dp
2546                 fs_fire_area(is:ie, js:je) = fire2tile(fr_arr=grid%fire_area(ifps:ifpe,jfps:jfpe), &
2547                                                    dsx=grid%sr_x,dsy=grid%sr_y)
2548
2549             ! WRITE(msg, '(i8,1x,i12)') ndep, COUNT(fs_fire_area(is:ie,js:je)> ZERO_dp)
2550             ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd: Spotfire Likelihood, deposit points: '//msg)
2551
2552             IF (ndep > 0) THEN
2553
2554                 landing_mask(ims:ime, jms:jme) = .FALSE.
2555                 landing_mask(is:ie, js:je) = ((fs_count_landed_hist(is:ie, js:je) > 0) .AND. &
2556                                           (fs_fire_area(is:ie, js:je) == ZERO_dp) )
2557                 ndep = COUNT(landing_mask)
2558                 fs_landing_mask(is:ie, js:je) = UNPACK(SPREAD(1,DIM=1,NCOPIES=ndep), MASK=landing_mask(is:ie, js:je), FIELD=0)
2559
2560                 ! 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))
2561                 ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd deposits, non fire-point deposits, fire-point deposits: '//msg)
2562
2563                 ALLOCATE(tmp1(ndep))
2564                 tmp1 = PACK(fs_count_landed_hist(is:ie, js:je), landing_mask(is:ie, js:je)) ! WHERE doest work on memory bounds
2565                 fs_count_landed_hist(is:ie, js:je) = ZERO_dp
2566                 fs_count_landed_hist(is:ie, js:je) = UNPACK(tmp1, MASK=landing_mask(is:ie, js:je), FIELD=ZERO_dp)
2567
2568                 ! WRITE(msg, '(2(i8,1x), 1(f14.6,1x))') COUNT(fs_count_landed_hist > ZERO_dp), ndep, MAXVAL(fs_count_landed_hist)
2569                 ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd deposits, non fire-point deposits, fire-point deposits: '//msg)
2570                 DEALLOCATE(tmp1)
2571
2572                 ! Were these firebrand deposited outside the fire area?
2573                 !-------------------------------------------------------------------------------
2574
2575                 IF (ndep > 0) THEN
2576
2577                     ! Calculate fuel risk
2578                     !-------------------------------------------------------------------------------
2579
2580                     ALLOCATE(fuel_spotting_risk_fr(ifps:ifpe,jfps:jfpe), &!fs_fuel_spotting_risk(ims:ime,jms:jme), &
2581                              spotthreat_fr(ifps:ifpe,jfps:jfpe))
2582                     fuel_spotting_risk_fr(:,:) = ZERO_dp
2583                     fs_fuel_spotting_risk(:,:) = ZERO_dp
2584                     spotthreat_fr(:,:) = ZERO_dp
2585
2586                     IF (.NOT. ALLOCATED(fcat_fr)) THEN
2587                         ALLOCATE(fcat_fr(ifps:ifpe,jfps:jfpe))
2588                         fcat_fr(:,:) = ZERO_dp
2589                         fcat_fr(ifps:ifpe,jfps:jfpe) = RESHAPE(&
2590                                                        get_fuel_cat(crosswalk=grid%fuel_crosswalk, &
2591                                                                     cat=grid%nfuel_cat(ifps:ifpe,jfps:jfpe)), &
2592                                                        SHAPE(fuel_spotting_risk_fr))
2593                     ENDIF
2594                     !WRITE(msg, *) MAXVAL(fcat_fr), MINVAL(fcat_fr)
2595                     ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd fcat max,min       '//msg)
2596
2597                     spotthreat_fr(ifps:ifpe,jfps:jfpe)= spotting_threat_factor(fcat_fr(ifps:ifpe,jfps:jfpe)) !RESHAPE(SHAPE(fuel_spotting_risk_fr))
2598                                                              ! get_fuel_cat(crosswalk=grid%fuel_crosswalk, &
2599                                                              !     cat=grid%nfuel_cat(ifps:ifpe,jfps:jfpe))), &
2600                     ! WRITE(msg, *) MAXVAL(spotthreat_fr), MINVAL(spotthreat_fr)
2601                     ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd spotting_threat_factor max,min'//msg)
2602
2603                     fuel_spotting_risk_fr(ifps:ifpe,jfps:jfpe) = RESHAPE(&
2604                                                         fuel_spotting_risk(&
2605                                                              fuel_fgi=grid%fgip(ifps:ifpe,jfps:jfpe), &
2606                                                              fuel_mcg=fmcg_2df(ifps:ifpe,jfps:jfpe), &
2607                                                              !fuel_mcg=grid%fmcg2df(ifps:ifpe,jfps:jfpe), &
2608                                                              factor=spotthreat_fr(ifps:ifpe,jfps:jfpe)), &
2609                                                              SHAPE(fuel_spotting_risk_fr)) !, fuel_mcg=grid%fmcg2df)
2610
2611                     fs_fuel_spotting_risk(is:ie,js:je) = fire2tile(fr_arr=fuel_spotting_risk_fr(ifps:ifpe,jfps:jfpe), &
2612                                                          dsx=grid%sr_x,dsy=grid%sr_y)
2613                     DEALLOCATE(fcat_fr)
2614
2615                     !-------------------------------------------------------------------------------
2616
2617
2618                     ! Prepare for MPI
2619                     !-------------------------------------------------------------------------------
2620
2621                     ALLOCATE(landing_cnt(ndep), landed_risk(ndep), recv_spot_lkhd(ndep), recv_frac_landed(ndep))
2622                     landing_cnt = PACK(fs_count_landed_hist(is:ie,js:je), landing_mask(is:ie,js:je))
2623                     landed_risk = PACK(fs_fuel_spotting_risk(is:ie,js:je), landing_mask(is:ie,js:je))
2624
2625 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
2626                     IF (.NOT. (IamMaster)) THEN
2627
2628                         CALL fs_mpi_sendbuffsize(sendto=MasterId, nbr=ndep) ! Comm send >0: BuffSz
2629                         CALL fs_mpi_sendbuff_int(sendto=MasterId, bsz=ndep, buff=landing_cnt(1:ndep))
2630                         CALL fs_mpi_sendbuff_real(sendto=MasterId, bsz=ndep, buff=landed_risk(1:ndep))
2631
2632                     ENDIF
2633 #endif
2634                 ENDIF ! (fs_count_landed_hist > 0 .AND. fire_area == 0)
2635             ENDIF ! SUM(fs_count_landed_hist) > 0
2636
2637             !-------------------------------------------------------------------------------
2638
2639 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
2640             IF ((ndep == 0) .AND. .NOT.(IamMaster)) &
2641                 CALL fs_mpi_nothing2send(sendto=MasterId) ! Comm send 0: BuffSz
2642 #endif
2643             !-------------------------------------------------------------------------------
2644             ! Master: Calculate Likelihood Field
2645             !-------------------------------------------------------------------------------
2646
2647             IF (IamMaster) THEN
2648
2649                 !-------------------------------------------------------------------------------
2650                 ! Declare arrays to receive data
2651
2652                 ALLOCATE(np_nlanded(ntiles))
2653                 np_nlanded(:) = 0
2654                 np_nlanded(1) = ndep ! Array index for Master is 1
2655
2656 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
2657                 DO kk=2, mpiprocs ! kk=1 is Master, MPI-IDs begin at 0
2658                     np_nlanded(kk) = fs_mpi_recvbuffsize(fromid=kk-1)  ! Comm receive: BuffSz
2659                     ! WRITE(msg, '(2(i8,1x))') kk,np_nlanded(kk)
2660                     ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd proc deposit: '//msg)
2661                 ENDDO
2662 #endif
2663                 !-------------------------------------------------------------------------------
2664                 ! Communicate deposits between Master and procs
2665
2666                 ndep_total = SUM(np_nlanded)
2667
2668                 IF ( ndep_total > 0) THEN
2669                     ALLOCATE(np_landed_cnt(ndep_total), np_landed_risk(ndep_total))
2670                     IF (ndep > 0) THEN
2671                         np_landed_cnt(1:ndep) = landing_cnt
2672                         np_landed_risk(1:ndep)= landed_risk
2673                     ELSE
2674                         np_landed_cnt(1) = ZERO_dp ! I don't think I need this
2675                         np_landed_risk(1)= ZERO_dp
2676                     ENDIF
2677
2678                     !-------------------------------------------------------------------------------
2679                     ! Receive deposit data from procs
2680
2681                     kk1 = ndep + 1 ! index start after master's
2682
2683 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
2684                     DO kk=2, mpiprocs ! kk=1 is Master
2685
2686                         IF (np_nlanded(kk) > 0) THEN
2687
2688                             kk2 = kk1 + np_nlanded(kk) -1
2689
2690                             ! WRITE(msg, '(4(i8,1x))') kk, np_nlanded(kk), kk1, kk2
2691                             ! CALL wrf_debug (wrfdbg, 'SPFire_lkhd master kk, np_dep, k1:k2: '//msg)
2692                             np_landed_cnt(kk1:kk2) = fs_mpi_recvfrompatch_int(bsz=np_nlanded(kk), fromid=kk-1)
2693                             np_landed_risk(kk1:kk2)= fs_mpi_recvfrompatch_real(bsz=np_nlanded(kk), fromid=kk-1)
2694                             kk1 = kk1 + np_nlanded(kk)
2695                         ENDIF
2696                     ENDDO
2697 #endif
2698
2699                     !-------------------------------------------------------------------------------
2700                     ! Calculate ratio and lkhd
2701
2702                     ALLOCATE(np_frac_landed(ndep_total), np_lkhd(ndep_total))
2703                     np_frac_landed(:) = ZERO_dp
2704                     np_lkhd(:) = ZERO_dp
2705
2706                     np_frac_landed(1:ndep_total) = REAL(np_landed_cnt(1:ndep_total))/REAL(SUM(np_landed_cnt))
2707                     np_lkhd(1:ndep_total) = np_frac_landed(1:ndep_total) * np_landed_risk(1:ndep_total)
2708                     np_lkhd(1:ndep_total) = np_lkhd(1:ndep_total)/MAXVAL(np_lkhd)
2709
2710                     ! DO pp=1,SIZE(np_landed_cnt)
2711                     !     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)
2712                     !     CALL wrf_debug (wrfdbg, 'SPFire_lkhd master np_frac_landed, np_lkhd, : '//msg)
2713                     ! ENDDO
2714
2715                     !-------------------------------------------------------------------------------
2716                     ! Send dat to procs
2717
2718                     kk1 = ndep + 1 ! index start after master's
2719
2720 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
2721                     DO kk=2, mpiprocs ! kk=1 is Master
2722
2723                         IF (np_nlanded(kk) > 0) THEN
2724
2725                             kk2 = kk1 + np_nlanded(kk) -1
2726
2727                             CALL fs_mpi_sendbuff_real(sendto=kk-1, bsz=np_nlanded(kk), buff=np_frac_landed(kk1:kk2))
2728                             CALL fs_mpi_sendbuff_real(sendto=kk-1, bsz=np_nlanded(kk), buff=np_lkhd(kk1:kk2))
2729
2730                             kk1 = kk1 + np_nlanded(kk)
2731                         ENDIF
2732                     ENDDO
2733 #endif
2734                 ENDIF ! SUM(np_nlanded) > 0
2735                 !-------------------------------------------------------------------------------
2736
2737             ENDIF ! IamMaster
2738
2739             !-------------------------------------------------------------------------------
2740
2741 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
2742             IF (.NOT.(IamMaster) .AND. (ndep > 0)) THEN
2743
2744                 recv_frac_landed = fs_mpi_recvfrompatch_real(bsz=ndep, fromid=MasterId)
2745                 fs_frac_landed(is:ie,js:je) = UNPACK(recv_frac_landed, MASK=landing_mask(is:ie,js:je), FIELD=ZERO_dp)
2746
2747                 recv_spot_lkhd = fs_mpi_recvfrompatch_real(bsz=ndep, fromid=MasterId)
2748                 fs_spotting_lkhd(is:ie,js:je) = UNPACK(recv_spot_lkhd, MASK=landing_mask(is:ie,js:je), FIELD=ZERO_dp)
2749
2750             ENDIF
2751 #endif
2752
2753             IF ((IamMaster) .AND. (ndep > 0)) THEN
2754
2755                 fs_frac_landed(is:ie,js:je) = UNPACK(np_frac_landed(1:ndep), MASK=landing_mask(is:ie,js:je), FIELD=ZERO_dp)
2756                 fs_spotting_lkhd(is:ie,js:je) = UNPACK(np_lkhd(1:ndep), MASK=landing_mask(is:ie,js:je), FIELD=ZERO_dp)
2757
2758             ENDIF
2759
2760         ENDIF ! hist_flag
2761
2762         !===============================================================================
2763         !===============================================================================
2764
2765
2766         !******************************************************************
2767         !******************************************************************
2768         !*                                                                *
2769         !*                  MPI Receive firebrands                   *
2770         !*                                                                *
2771         !******************************************************************
2772         !******************************************************************
2773
2774
2775         !===============================================================================
2776         ! MPI: Receive Brands (into this patch)
2777         !===============================================================================
2778
2779
2780
2781         !-------------------------------------------------------------------------------
2782         ! Anything to receive? This part is only needed when compiled in parallel
2783         !-------------------------------------------------------------------------------
2784
2785 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
2786
2787         ALLOCATE(nbr_id(neighbors))
2788         nbr_id(:)=fs_mpi_checkreceive(task_list=task_id, np=neighbors)
2789
2790         nbr_sum = SUM(nbr_id)
2791         IF (nbr_sum > 0) THEN
2792
2793             WRITE (msg,'(16(i4,1x))') ([task_id(ii), nbr_id(ii)], ii=1,neighbors)
2794             CALL wrf_debug (wrfdbg,  'SPFire_driver mpi_check_receive: '//msg)
2795
2796             ! Receiving positions, id, dt
2797             ALLOCATE(r_x(nbr_sum), r_y(nbr_sum), r_z(nbr_sum), r_id(nbr_sum), r_dt(nbr_sum))
2798             !ALLOCATE(release_mpi_ijk(nbr_sum,3))
2799
2800             ! Receiving p_propertieserties
2801             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))
2802             ALLOCATE(fb_mdetv(nbr_sum))
2803
2804             CALL fs_mpi_recv(np_id=nbr_id, task_id=task_id, &
2805                                    r_x=r_x, r_y=r_y, r_z=r_z, &
2806                                    r_id=r_id, r_dt=r_dt, &
2807                                    r_p_m=r_p_m, &
2808                                    r_p_d=r_p_d, &
2809                                    r_p_e=r_p_e, &
2810                                    r_p_t=r_p_t, &
2811                                    r_p_v=r_p_v)
2812
2813             ! Assign values to array and release incoming particles
2814             !release_mpi_ijk = RESHAPE([r_x,r_y,r_z], [nbr_sum,3])
2815             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)]
2816
2817             ! DO kk=1,MIN(nbr_sum,5)
2818             !     WRITE(msg,'(6(f12.6,1x))') (release_mpi_ijk(kk,ii), ii=1,3), fb_mdetv(kk)%p_diam, fb_mdetv(kk)%p_temp, fb_mdetv(kk)%p_tvel
2819             !     CALL wrf_debug (wrfdbg, 'SPFire_mpi recv AFTER  >>> '// msg)
2820             ! ENDDO
2821
2822             prior_br = active_br + 1 ! Array is packed so new particles will be assgned at this position
2823             CALL generate_firebrands(&
2824                         fs_p_id = fs_p_id, &
2825                         fs_p_dt = fs_p_dt, &
2826                         fs_p_z  = fs_p_z,  &
2827                         fs_p_x  = fs_p_x,  &
2828                         fs_p_y  = fs_p_y,  &
2829                         fs_gen_idmax    = fs_gen_idmax,    &
2830                         active_br   = active_br,   &
2831                         !release_ijk = release_mpi_ijk, &
2832                         release_i = r_x, &
2833                         release_j = r_y, &
2834                         release_k = r_z, &
2835                         release_prop= fb_mdetv, &
2836                         fs_p_prop   = fs_p_prop)
2837
2838             WRITE (msg,'(2(i8,1x))') active_br, fs_gen_idmax
2839             CALL wrf_debug (wrfdbg, 'SPFire_driver mpi recv AFTER : ii, fs_gen_idmax >>> '// msg)
2840
2841             fs_p_mass (prior_br:active_br) = r_p_m
2842             fs_p_diam (prior_br:active_br) = r_p_d
2843             fs_p_effd(prior_br:active_br) = r_p_e
2844             fs_p_temp (prior_br:active_br) = r_p_t
2845             fs_p_tvel (prior_br:active_br) = r_p_v
2846
2847             ! DO kk=prior_br,prior_br+MIN(nbr_sum,5)
2848             !     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)
2849             !     CALL wrf_debug (wrfdbg, 'SPFire_mpi recv AFTER2  >>> '// msg)
2850             ! ENDDO
2851
2852             DEALLOCATE(r_x, r_y, r_z, r_id, r_dt)
2853             !DEALLOCATE(release_mpi_ijk)
2854             DEALLOCATE(fb_mdetv)
2855             DEALLOCATE(r_p_m, r_p_d, r_p_e, r_p_t, r_p_v)
2856
2857         ! ELSE
2858         !     CALL wrf_debug (wrfdbg,  'SPFire_driver mpi nothing to receive')
2859         ENDIF
2860         DEALLOCATE(nbr_id)
2861         !-------------------------------------------------------------------------------
2862         ! DONE :)
2863         !-------------------------------------------------------------------------------
2864 #endif
2865         !===============================================================================
2866         ! Check history interval
2867         !===============================================================================
2868
2869         hist_flag = .FALSE.
2870
2871
2872     END SUBROUTINE firebrand_spotting_em_driver
2873 !=============================================================
2874
2875 !=============================================================
2876
2877
2878
2879 !=============================================================
2880 !=============================================================
2881
2882 !******************************************************************
2883 !******************************************************************
2884 !*                                                                *
2885 !*                Particle Transport & Physics                    *
2886 !*                                                                *
2887 !******************************************************************
2888 !******************************************************************
2889
2890 !=============================================================
2891     PURE &
2892     SUBROUTINE advect_xyz_m(grid, xp, yp, hgt, dtp, idp, &
2893                             u, v, w, dt, mf, z_at_w, znw, ims, jms, kms, &
2894                             phyd, thet, rho, &
2895                             xout, yout, zout, fs_p_prop, land_hgt, &
2896                             start_mom3d_dt, msg)
2897 !=============================================================
2898
2899     ! Transport a particle at a given xp, yp, zp position in meters
2900     ! given the corresponding velocities [m/s] at the enclosing grid points
2901
2902     ! arguments
2903     ! xp, yp, hgt: starting particle position
2904     ! u, v, w: wind components
2905     ! dt: time interval in sec
2906     ! ims, jms : start index of wind components
2907     ! mf: constant to convert from grid position to meters: rdx*msft(x,y)
2908     ! z_at_w: wrf height at vertical interface
2909     ! znw: eta levels
2910
2911     !-------------------------------------------------------------
2912
2913         IMPLICIT NONE
2914
2915         TYPE(domain), INTENT(IN) :: grid
2916         INTEGER, INTENT(IN)   :: ims, jms, kms, start_mom3d_dt
2917         REAL,    INTENT(IN)   :: dt
2918         REAL, INTENT(IN)      :: land_hgt
2919         REAL,    INTENT(IN),  DIMENSION(:) :: xp, yp, hgt ! particle position
2920         INTEGER, INTENT(IN),  DIMENSION(:) :: dtp, idp    ! particle life-time and ID
2921         REAL,    INTENT(IN),  DIMENSION(ims:, kms:, jms:) :: u, v, w, z_at_w
2922         REAL,    INTENT(IN),  DIMENSION(ims:, kms:, jms:) :: phyd, thet, rho
2923         REAL,    INTENT(IN),  DIMENSION(:) :: znw
2924         REAL,    INTENT(IN),  DIMENSION(ims:, jms:) :: mf ! map factor, msft * rdx
2925         REAL,    INTENT(OUT), DIMENSION(:) :: xout, yout, zout!, aglh
2926         TYPE(p_properties),INTENT(INOUT),DIMENSION(:):: fs_p_prop
2927
2928         REAL,    DIMENSION(3) :: uvw2p
2929         INTEGER :: pp, aux
2930         INTEGER :: hsz, ihs, jhs, ihe, jhe  ! tile +- halo array bounds
2931         REAL    :: xp_m0, yp_m0, zp_m0, zp0
2932         REAL    :: xp_m1, yp_m1, zp_m1, xp1, yp1, zp1
2933         REAL    :: xp_m2, yp_m2, zp_m2, xp2, yp2, zp2
2934         REAL    :: zp, wp, brz
2935         REAL    :: loc_p, loc_t, loc_d ! local met properties
2936         CHARACTER (LEN=100), INTENT(OUT) :: msg
2937
2938         WRITE(msg, '(a100)') ' '
2939     !-------------------------------------------------------------------------------
2940
2941         hsz = 4 ! halo size - get this value from grid%. Halo:48 = 4 points. How about staggering?
2942         ihs = is - hsz    ! MAX(, ids)
2943         jhs = js - hsz    ! MAX(, jds)
2944         ihe = ie -1 + hsz ! MIN(, ide-1) ! check for ide -1 bc we use i, i+1 to interpolate
2945         jhe = je -1 + hsz ! MIN(, jde-1)
2946
2947         ! WRITE (msg,'(4(i4,1x))') ihs,ihe, jhs, jhe
2948         ! CALL wrf_debug (wrfdbg, 'SPFire_advect halo bounds ihs, ihe, jhs, jhe: >>> '//msg)
2949
2950         !start_mom3d_dt = 4 ! min lifetime before calculating fb physics (number of timesteps for this nested grid)
2951
2952         DO pp=1, SIZE(xp)
2953
2954             ! CALL wrf_debug (wrfdbg, 'SPFire advect_xyz_m ---------------------------------------------------- ')
2955
2956             ! WRITE (msg,'(3(i8,1x),3(f12.6,1x))') pp, idp(pp), dtp(pp), xp(pp), yp(pp), hgt(pp)
2957             ! CALL wrf_debug (wrfdbg, 'SPFire BEFORE advect:    pp,id,dt,xp,yp,zp:    >>> '//msg)
2958
2959             !-------------------------------------------------------------------------------
2960             ! Convert horizontal positions to meters
2961             !-------------------------------------------------------------------------------
2962
2963             xp_m0 = xp(pp) / mf(FLOOR(xp(pp)),FLOOR(yp(pp)))
2964             yp_m0 = yp(pp) / mf(FLOOR(xp(pp)),FLOOR(yp(pp)))
2965             zp_m0 = hgt(pp)
2966
2967             ! zp changes slightly between wrf time steps, so it need to be calculated again
2968             zp0 = hgt2k(xp=xp(pp), yp=yp(pp), hgt=zp_m0, z_at_w=z_at_w, znw=znw)
2969             ! (Need fractional position for interpolation)
2970
2971             ! WRITE (msg,'(3(i8,1x),4(f12.6,1x))') pp, idp(pp), dtp(pp), xp(pp), yp(pp), hgt(pp), zp0
2972             ! CALL wrf_debug (wrfdbg, 'SPFire_advect BEFORE:    pp,id,dt,xp,yp,zp,zk:      >>> '//msg)
2973
2974             ! If this is a particle deposited by an adjacent tile, do not transport it
2975             IF (zp_m0 <= land_hgt) THEN
2976                 ! CALL wrf_debug (wrfdbg, 'SPFire_advect particle at deposit threshold - will not advect')
2977
2978                 xout(pp) = xp(pp)
2979                 yout(pp) = yp(pp)
2980                 zout(pp) = zp_m0
2981                 CYCLE
2982
2983             ENDIF
2984
2985             !-------------------------------------------------------------------------------
2986             !-------------------------------------------------------------------------------
2987
2988
2989         ! 1ST PASS
2990         !=============================================================
2991
2992
2993             !-------------------------------------------------------------------------------
2994             ! #1 interpolate velocities to particle position & advect particle
2995             !-------------------------------------------------------------------------------
2996
2997             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
2998
2999             xp_m1 = uvw2p(1)*DT + xp_m0
3000             yp_m1 = uvw2p(2)*DT + yp_m0
3001             zp_m1 = uvw2p(3)*DT + zp_m0
3002
3003             !-------------------------------------------------------------------------------
3004             ! Convert horizontal positions from meters
3005             !-------------------------------------------------------------------------------
3006
3007             xp1 = xp_m1 * mf(FLOOR(xp(pp)),FLOOR(yp(pp))) ! 2 pass will account for grid pt changes in mf & u
3008             yp1 = yp_m1 * mf(FLOOR(xp(pp)),FLOOR(yp(pp)))
3009
3010             IF((ABS(xp1-xp(pp)) > 2.0) .OR. (ABS(yp1-yp(pp)) > 2.0)) THEN
3011                 WRITE (msg,'(2(i4,1x),7(f10.6,1x))') pp, dtp(pp), xp(pp), yp(pp), zp0, &
3012                                                      uvw2p(1), uvw2p(2), uvw2p(3), dt !
3013                 zout(pp) = ZERO_dp
3014                 xout(pp) = ZERO_dp
3015                 yout(pp) = ZERO_dp
3016                 CYCLE
3017             ENDIF
3018
3019             !-------------------------------------------------------------------------------
3020             ! Convert z position from meters to k-fractional
3021             !-------------------------------------------------------------------------------
3022
3023             ! CALL wrf_debug (wrfdbg, 'SPFire advect_xyz_m 1st pass --------------------------------------------- ')
3024
3025             ! WRITE (msg,'(3(i8,1x),3(f12.6,1x))') pp, idp(pp), dtp(pp), xp1, yp1, zp_m1
3026             ! CALL wrf_debug (wrfdbg, 'SPFire after 1-pass adv: pp,id,dt,xp,yp,zp:     >>> '//msg)
3027             ! WRITE (msg,*) pp, z_at_w(FLOOR(xp1),1,FLOOR(yp1))
3028             ! CALL wrf_debug (wrfdbg, 'SPFire after 1-pass adv: zn:                    >>> '//msg)
3029             ! WRITE (msg,*) pp, z_at_w(FLOOR(xp1-dp05),1,FLOOR(yp1-dp05))
3030             ! CALL wrf_debug (wrfdbg, 'SPFire after 1-pass adv: zn:                    >>> '//msg)
3031
3032             zp1 = hgt2k(xp=xp1, yp=yp1, hgt=zp_m1, z_at_w=z_at_w, znw=znw)
3033             wp = uvw2p(3)
3034
3035             !-------------------------------------------------------------------------------
3036             ! Is the new position not on the same memory?
3037             !-------------------------------------------------------------------------------
3038             ! check for floor(x)+1 bc we use i and i+1 to interpolate
3039             IF (FLOOR(xp1-dp05) < ihs .OR. FLOOR(xp1-dp05) +1 > ihe .OR. &
3040                 FLOOR(yp1-dp05) < jhs .OR. FLOOR(yp1-dp05) +1 > jhe) THEN
3041                 ! Checking staggered positions relative to unstaggered variables (+- dp05)
3042
3043                 !-------------------------------------------------------------------------------
3044                 ! Use only one pass
3045                 !-------------------------------------------------------------------------------
3046
3047                 xout(pp) = xp1
3048                 yout(pp) = yp1
3049                 zout(pp) = zp_m1
3050                 zp = zp1
3051
3052             ELSE
3053
3054             ! 2ND PASS
3055             !=============================================================
3056
3057                 !-------------------------------------------------------------------------------
3058                 ! #2 interpolate velocities to particle position & advect particle
3059                 !-------------------------------------------------------------------------------
3060
3061                 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
3062
3063                 ! WRITE (msg,'(3(f14.9,1x))') uvw2p(1), uvw2p(2), uvw2p(3)
3064                 ! CALL wrf_debug (wrfdbg, 'SPFire_advect effective vel m/s: u,v,w 2nd pass>>> '//TRIM(msg))
3065
3066                 xp_m2 = uvw2p(1)*DT + xp_m0
3067                 yp_m2 = uvw2p(2)*DT + yp_m0
3068                 zp_m2 = uvw2p(3)*DT + zp_m0
3069
3070                 xp2 = xp_m2 * mf(FLOOR(xp1),FLOOR(yp1))
3071                 yp2 = yp_m2 * mf(FLOOR(xp1),FLOOR(yp1))
3072
3073                 IF((ABS(xp2-xp1) > 2.0) .OR. (ABS(yp2-yp1) > 2.0)) THEN
3074                     WRITE (msg,'(2(i4,1x),7(f10.6,1x))') pp, dtp(pp), xp1, yp1, zp1, &
3075                                                          uvw2p(1), uvw2p(2), uvw2p(3), dt !
3076                     zout(pp) = ZERO_dp
3077                     xout(pp) = ZERO_dp
3078                     yout(pp) = ZERO_dp
3079                 CYCLE
3080             ENDIF
3081
3082
3083                 !-------------------------------------------------------------------------------
3084                 ! Is the new position not on the same memory?
3085                 !-------------------------------------------------------------------------------
3086
3087                 ! check against tile +- halo points instead
3088
3089                 IF (FLOOR(xp2-dp05) < ihs .OR. FLOOR(xp2-dp05) +1 > ihe .OR. &
3090                     FLOOR(yp2-dp05) < jhs .OR. FLOOR(yp2-dp05) +1 > jhe) THEN
3091                     ! Checking staggered positions relative to unstaggered variables (+- dp05)
3092
3093                     !-------------------------------------------------------------------------------
3094                     ! Use xp1, yp1 to find zp2 (z_at_w will be out of bounds)
3095                     !-------------------------------------------------------------------------------
3096
3097                     zp2 = hgt2k(xp=xp1, yp=yp1, hgt=zp_m2, z_at_w=z_at_w, znw=znw)
3098
3099                 ELSE
3100
3101                     zp2 = hgt2k(xp=xp2, yp=yp2, hgt=zp_m2, z_at_w=z_at_w, znw=znw)
3102
3103                 ENDIF
3104
3105             !=============================================================
3106
3107                 !-------------------------------------------------------------------------------
3108                 ! Average both passes
3109                 !-------------------------------------------------------------------------------
3110
3111                 xp2 = (xp1 + xp2) * dp05
3112                 yp2 = (yp1 + yp2) * dp05
3113                 zp_m2 = (zp_m1 + zp_m2) * dp05
3114                 zp2 = (zp1 + zp2) * dp05
3115
3116                 xout(pp) = xp2
3117                 yout(pp) = yp2
3118                 zout(pp) = zp_m2
3119                 zp = zp2
3120                 wp = (uvw2p(3) + wp) * dp05
3121
3122             ENDIF
3123
3124             !=============================================================
3125             ! Transport Done :)
3126             !=============================================================
3127             ! WRITE (msg, '(a100)') ' '
3128             ! WRITE (msg,'(2(i6,1x),4(f12.6,1x))') pp, dtp(pp), xout(pp), yout(pp), zout(pp), zp
3129             ! CALL wrf_debug (wrfdbg, 'SPFire_advect AFTER :    pp,id,dt,xp,yp,zp,kp:      >>> '//msg)
3130
3131             !=============================================================
3132             ! Firebrand Physics
3133             !=============================================================
3134
3135             loc_p = ZERO_dp
3136             loc_d = ZERO_dp
3137             loc_t = ZERO_dp
3138
3139             IF (zout(pp) <= ZERO_dp) THEN
3140                 fs_p_prop(pp) = p_properties(ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP)
3141                 zout(pp) = ZERO_dp
3142                 CYCLE ! Skip physics
3143             ENDIF
3144
3145             IF (IEEE_IS_NAN(zout(pp)) .OR. ABS(zout(pp)) > HUGE(1.0)) THEN
3146                 fs_p_prop(pp) = p_properties(ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP)
3147                 zout(pp) = ZERO_dp
3148                 xout(pp) = ZERO_dp
3149                 yout(pp) = ZERO_dp
3150
3151                 CYCLE ! Skip physics
3152             ENDIF
3153
3154             IF (dtp(pp) < start_mom3d_dt) THEN ! Temporarily testing transport: extending particle lifetime
3155                 ! WRITE (msg,*) dt
3156                 ! CALL wrf_debug (wrfdbg, 'SPFire_physics cycle dt: '//msg)
3157                 CYCLE
3158             ENDIF
3159
3160
3161             CALL get_local_met( xp  = xout(pp),   &
3162                                 yp  = yout(pp),   &
3163                                 zp  = zp,         &
3164                                 ihs = ihs,  &
3165                                 jhs = jhs,  &
3166                                 ihe = ihe,  &
3167                                 jhe = jhe,  &
3168                                 p = phyd, & ! air pressure (Pa)
3169                                 t = thet, & ! temperature (K)
3170                                 d = rho,  & ! density (kg/m3)
3171                                 loc_p = loc_p,    & ! local air pressure 3d (Pa)
3172                                 loc_d = loc_d,    & ! local density (g/m3)
3173                                 loc_t = loc_t)      ! local temperature (K)
3174
3175             ! WRITE (msg,'(3(f16.8,1x))') grid%p_hyd_w(FLOOR(xout(pp)),FLOOR(zp),FLOOR(yout(pp))), &
3176             !                             grid%rho(FLOOR(xout(pp)),FLOOR(zp),FLOOR(yout(pp))), &
3177             !                             grid%t_phy(FLOOR(xout(pp)),FLOOR(zp),FLOOR(yout(pp)))
3178             ! CALL wrf_debug (wrfdbg, 'SPFire grid%p, %rho, %t_phy      : '//msg)
3179
3180             ! WRITE (msg,'(4(f16.8,1x))') loc_p, loc_d, loc_t, wp
3181             ! CALL wrf_debug (wrfdbg, 'SPFire br met p, d, t, w : '//msg)
3182
3183             !-------------------------------------------------------------------------------
3184             ! Firebrand Burnout and Terminal Velocity
3185             !-------------------------------------------------------------------------------
3186
3187             brz = zout(pp)
3188             CALL firebrand_physics(dt = dt,         &
3189                                      hgt = brz,  &
3190                                      loc_w = wp,      &
3191                                      loc_p = loc_p,   &
3192                                      loc_t = loc_t,   &
3193                                      loc_d = loc_d,   &
3194                                      fbprop  = fs_p_prop(pp))
3195
3196             IF (IEEE_IS_NAN(fs_p_prop(pp)%p_tvel)) THEN
3197                 fs_p_prop(pp) = p_properties(ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP, ZERO_DP)
3198                 zout(pp) = ZERO_dp
3199
3200                 ! 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
3201                 ! CALL wrf_debug (wrfdbg, 'SPFire_physics AFTER:   pp,id,dt,xp,yp,zp,fbvt:    >>> '//msg)
3202
3203             ELSE
3204
3205                 zout(pp) = brz
3206
3207             ENDIF
3208             ! 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
3209             ! CALL wrf_debug (wrfdbg, 'SPFire_physics AFTER:   pp,id,dt,xp,yp,zp,fbvt:    >>> '//msg)
3210
3211         ENDDO
3212
3213     END SUBROUTINE advect_xyz_m
3214 !=============================================================
3215
3216
3217 ! Cannot be elemental: array
3218 !=============================================================
3219 PURE FUNCTION uvw_3d_interp(x, y, z, u, v, w, ihs, jhs, ihe, jhe)
3220 !=============================================================
3221
3222     ! Wrapper to interpolate U,V,W to a point within a 3D grid box
3223
3224     ! Interpolate velocities at grid point interfaces to particle position
3225     ! input range needed: i:i+1 | i < x(i) <= i+1
3226     !                     j:j+1 | j < y(j) <= j+1
3227     !                     k:k+1 | k < z(k) <= k+1
3228
3229     ! arguments
3230     ! x, y, z: fractional position to interpolate velocities (e.g.: z is a fraction of the k index)
3231     ! u, v, w: wind components
3232     ! is, js : start index of wind components
3233
3234     !-------------------------------------------------------------
3235     ! Interpolate between planes: U: xy(k), xy(k+1) stag x, dim=1
3236     !                             V: xy(k), xy(k+1) stag y, dim=2
3237     !                             W: xz(j), xz(j+1) stag z, dim=2
3238     !
3239     ! particle position is not staggered:
3240     !          grid center is given by (FLOOR(x), FLOOR(y))
3241     !-------------------------------------------------------------
3242
3243     IMPLICIT NONE
3244
3245     INTEGER, INTENT(IN)   :: ihs, jhs, ihe, jhe  ! tile +- halo array bounds
3246     REAL,    INTENT(IN)   :: x, y, z ! fractional positions relative to integral array indices
3247     REAL,    INTENT(IN), DIMENSION(ims:, 1:, jms:) :: u, v ,w
3248     REAL,    DIMENSION(3) :: uvw_3d_interp
3249     REAL     :: tru, trv, trw, x0, y0, z0
3250     INTEGER  :: i, j, k, i0, j0, k0
3251
3252     ! Unstaggerecd indices
3253     k = MAX(FLOOR(z), 1)
3254     j = FLOOR(y)
3255     i = FLOOR(x)
3256
3257     ! Variable Staggering:
3258     ! particle position is equivalent to a staggered dimension
3259     ! the staggered domain begins 1/2 pnt before and ends 1/2 pnt after the regular domain
3260     ! e.g: interpolate u (stag x): x_stag = xp, y_unst = yp -0.5, z_unst = zp -0.5
3261     ! By shifting the particle position by 1/2 point,
3262     ! the unstaggered velocities will represent the grid edge instead of grid center
3263
3264     x0 = x - dp05
3265     y0 = y - dp05
3266     z0 = z - dp05 ! z can be < 1.0 when below the first half level, u_3d_interp will extrapolate values
3267
3268     i0 = FLOOR(x0) !MIN(MAX(ihs, ), ihe) ! shift indices by 0.5 because unstaggered variables given at horizontal grid center
3269     j0 = FLOOR(y0) !MIN(MAX(jhs, ), ihe)
3270     k0 = MAX(1, FLOOR(z0))
3271
3272     !-----------------------------------------------------------------------------
3273     ! u,v interpolates between xy planes: xy(k) xy(k+1)
3274     !-----------------------------------------------------------------------------
3275
3276     ! staggered x
3277     tru= u_3d_interp( x=x, y=y0, z=z0, &
3278         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), &
3279         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))
3280
3281     ! staggered y
3282     trv= u_3d_interp( x=x0, y=y, z=z0, &
3283         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), &
3284         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))
3285
3286     !-----------------------------------------------------------------------------
3287     ! w interpolates between xz planes, not xy (so y corresponds to z): xz(j) xz(j+1)
3288     !-----------------------------------------------------------------------------
3289
3290     ! z and w are both staggered (z is derived from grid%z_at_w)
3291     trw= u_3d_interp( x=x0, y=z, z=y, &
3292         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  ), &
3293         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))
3294
3295     !-----------------------------------------------------------------------------
3296     ! Return U, V, W at x, y, z positions
3297     !-----------------------------------------------------------------------------
3298
3299     uvw_3d_interp = [tru, trv, trw]
3300
3301 END FUNCTION uvw_3d_interp
3302 !=============================================================
3303
3304
3305
3306 !=============================================================
3307 ELEMENTAL FUNCTION u_3d_interp(x, y, z, &
3308     u_i0j0_bot, u_i0j1_bot, u_i1j0_bot, u_i1j1_bot, &
3309     u_i0j0_top, u_i0j1_top, u_i1j0_top, u_i1j1_top )
3310 !=============================================================
3311
3312     ! Interpolate velocities at grid point interfaces to particle position
3313     ! input range needed: u( i: i+1, j: j+1, k: k+1) |
3314     !                       i < x <= i+1
3315     !                       j < y <= j+1
3316     !                       k < z <= k+1
3317
3318     ! arguments
3319     ! x, y, z: particle position to interpolate velocities
3320     ! u, v, w: wind components
3321
3322     !-----------------------------------------------------------------------------
3323     ! Bilinear interpolation :
3324     ! The product of the value at the desired point (x,y,z) and the area enclosed by the corners
3325     ! is equal to the sum of the
3326     ! products of the value at each corner and the partial area diagonally opposite the corner
3327     !-----------------------------------------------------------------------------
3328
3329     !-----------------------------------------------------------------------------
3330     ! The 3-D interpolation is achieved with a linear interpolation between
3331     ! the xy-interpolated points from two adjacent planes
3332     !-----------------------------------------------------------------------------
3333
3334     !-----------------------------------------------------------------------------
3335     ! Z must be offset by 0.5 to correctly interpolate between vertically unstaggered planes
3336     !-----------------------------------------------------------------------------
3337
3338     REAL, INTENT(IN)   :: x, y, z ! fractional position
3339     REAL, INTENT(IN) :: u_i0j0_bot, u_i0j1_bot, u_i1j0_bot, u_i1j1_bot ! u_lower
3340     REAL, INTENT(IN) :: u_i0j0_top, u_i0j1_top, u_i1j0_top, u_i1j1_top ! u_upper
3341     REAL :: u_3d_interp
3342
3343     REAL :: dbot, dtop, u_lower, u_upper
3344     INTEGER :: k, j, i
3345
3346     !-----------------------------------------------------------------------------
3347     ! Get weighted U, V, W at particle position based on surrounding xy and xz planes
3348     !-----------------------------------------------------------------------------
3349
3350     ! index
3351     k = MAX(1, FLOOR(z)) ! To interpolate w, z corresponds y:  ( x=x, y=z, z=y, ...
3352
3353     ! z can be < 1.0 when interpolating below the first level on half levels
3354     ! e.g.: z=0.75, dbot = 0.75 -1 =-0.25, dtop = 2 - 0.75 =1.25
3355     !       u_3d_interp = u_upper*-0.25 + u_lower*1.25
3356
3357     dbot = z - REAL(k)
3358     dtop = REAL(k+1) - z
3359
3360     !-----------------------------------------------------------------------------
3361     ! j = FLOOR(y)
3362     ! i = FLOOR(x)
3363
3364     u_lower = u_2d_interp( u_i0j0=u_i0j0_bot, u_i0j1=u_i0j1_bot, &
3365                            u_i1j0=u_i1j0_bot, u_i1j1=u_i1j1_bot, xp=x, yp=y)
3366     u_upper = u_2d_interp( u_i0j0=u_i0j0_top, u_i0j1=u_i0j1_top, &
3367                            u_i1j0=u_i1j0_top, u_i1j1=u_i1j1_top, xp=x, yp=y)
3368
3369     !-----------------------------------------------------------------------------
3370     ! Apply weights and calculate U at particle position
3371     !-----------------------------------------------------------------------------
3372
3373     u_3d_interp = u_upper * dbot + & !
3374                   u_lower * dtop
3375
3376     !-----------------------------------------------------------------------------
3377
3378 END FUNCTION u_3d_interp
3379 !=============================================================
3380
3381 !=============================================================
3382
3383
3384
3385 !=============================================================
3386  ELEMENTAL FUNCTION u_2d_interp(u_i0j0, u_i0j1, u_i1j0, u_i1j1, xp, yp)
3387 !=============================================================
3388
3389     ! Interpolate any variable at grid point interfaces to a position in between
3390     ! Minimum range: u(i:i+1, j:j+1) | i < xp <= i+1  &  j < yp <= j+1
3391
3392     ! arguments
3393     ! u      : 2D variable at grid 4 corners: (i0,j0); (i0,j1); (i1,j0); (i1,j1)
3394     ! xp, yp : position within grid box to interpolate velocity
3395     ! is, js : start index of variable array boundary
3396
3397     !-------------------------------------------------------------
3398     ! 2D Interpolation sketch
3399     ! a, b, c, d: grid point corners
3400     ! da, db, dc, dd: distance between the nearest corner and point projection
3401     !
3402     ! a + + + + + + d
3403     !   + dc | db +
3404     !   +____|____+
3405     !   + dd | da +
3406     !   +    |    +
3407     ! b + + + + + + c
3408     !-------------------------------------------------------------
3409
3410     IMPLICIT NONE
3411
3412     !INTEGER, INTENT(IN) :: is, js
3413     REAL,    INTENT(IN) :: xp, yp
3414     REAL,    INTENT(IN) :: u_i0j1, u_i0j0, u_i1j0, u_i1j1 ! DIMENSION(is:,js:) :: uu
3415     REAL    :: u_2d_interp
3416     REAL    :: d_a, d_b, d_c, d_d      ! grid vertices
3417     REAL    :: uu_a, uu_b, uu_c, uu_d  ! advection at vertices
3418     INTEGER :: i, j
3419
3420     i = FLOOR(xp)
3421     j = FLOOR(yp)
3422
3423     uu_a = u_i0j1 ! uu(i  ,j+1)
3424     uu_b = u_i0j0 ! uu(i  ,j  )
3425     uu_c = u_i1j0 ! uu(i+1,j  )
3426     uu_d = u_i1j1 ! uu(i+1,j+1)
3427
3428     ! WRITE (msg,'(4(f14.9,1x))') uu_a, uu_b, uu_c, uu_d
3429     ! CALL wrf_debug (wrfdbg, 'SPFire_uu_2d: uu _a,_b,_c,_d                   >>> '//TRIM(msg))
3430
3431     d_a = ABS( ( REAL(i+1) - xp ) * ( REAL(j  ) - yp ) )
3432     d_b = ABS( ( REAL(i+1) - xp ) * ( REAL(j+1) - yp ) )
3433     d_c = ABS( ( REAL(i  ) - xp ) * ( REAL(j+1) - yp ) )
3434     d_d = ABS( ( REAL(i  ) - xp ) * ( REAL(j  ) - yp ) )
3435
3436     ! WRITE (msg,'(4(f14.9,1x))') d_a, d_b, d_c, d_d
3437     ! CALL wrf_debug (wrfdbg, 'SPFire_uu_2d: d _a,_b,_c,_d                    >>> '//TRIM(msg))
3438
3439
3440    u_2d_interp = ( uu_a * d_a + &
3441                uu_b * d_b + &
3442                uu_c * d_c + &
3443                uu_d * d_d ) / &
3444                (d_a + d_b + d_c + d_d)
3445
3446     ! WRITE (msg,*) uu_2d_interp
3447     ! CALL wrf_debug (wrfdbg, 'SPFire_uu_2d:                              >>> '//TRIM(msg))
3448
3449 END FUNCTION u_2d_interp
3450 !=============================================================
3451
3452 !=============================================================
3453
3454
3455
3456
3457 !=============================================================
3458 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)
3459 !=============================================================
3460     IMPLICIT NONE
3461     INTEGER, INTENT(INOUT), DIMENSION(:) :: fs_p_id
3462     INTEGER, INTENT(IN),    DIMENSION(:) :: fs_p_dt
3463     !INTEGER, INTENT(IN) :: active_br
3464     REAL,    INTENT(IN), DIMENSION(:) :: fs_p_x, fs_p_y, fs_p_z, fs_p_effd
3465     LOGICAL, DIMENSION(fs_array_maxsize) :: sparse_mask
3466     INTEGER, INTENT(OUT) :: cnt
3467
3468     INTEGER, INTENT(IN) :: max_life_dt
3469     REAL, INTENT(IN)    :: land_hgt
3470
3471     !-------------------------------------------------------------------------------
3472
3473     sparse_mask = .FALSE.
3474     !-------------------------------------------------------------------------------
3475     ! Flag to remove brands beyond domain bounds
3476     !-------------------------------------------------------------------------------
3477
3478     sparse_mask = ( fs_p_id > 0 .AND. &
3479                     (FLOOR(fs_p_x-dp05)   < ids .OR. FLOOR(fs_p_y-dp05)   < jds .OR. &
3480                      FLOOR(fs_p_x-dp05)+1 > ide .OR. FLOOR(fs_p_y-dp05)+1 > jde))
3481     ! interpolation needs i, j, i+1, j+1
3482     ! 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)
3483     ! Then we use FLOOR(x), FLOOR(y), FLOOR(x)+1, FLOOR(y)+1 to interpolate
3484
3485     !-------------------------------------------------------------------------------
3486     ! Flag to remove brands living longer than firebrand_max_life_dt
3487     !-------------------------------------------------------------------------------
3488
3489     sparse_mask = ( sparse_mask .OR. &
3490                    (fs_p_id > 0 .AND. fs_p_dt > max_life_dt))
3491
3492     !-------------------------------------------------------------------------------
3493     ! Flag to remove brands with near zero mass
3494     !-------------------------------------------------------------------------------
3495
3496     ! sparse_mask = (sparse_mask .OR. &
3497     !                (fs_p_id > 0 .AND. fs_p_mass <= br_min_mass))
3498
3499     !-------------------------------------------------------------------------------
3500     ! Flag to remove brands with near zero diameter that are not at deposit height
3501     !-------------------------------------------------------------------------------
3502
3503     sparse_mask = (sparse_mask .OR. &
3504                    (fs_p_id > 0 .AND. ((fs_p_effd <= TINY(1.0) .AND. (fs_p_z > land_hgt)))))
3505
3506     !-------------------------------------------------------------------------------
3507     ! Reset array
3508     !-------------------------------------------------------------------------------
3509     cnt = 0
3510     cnt = COUNT(sparse_mask)
3511     IF (cnt > 0) THEN
3512
3513         WHERE(sparse_mask) fs_p_id = 0
3514
3515         ! WRITE (msg,'(3(i8,1x))') COUNT(sparse_mask), COUNT( fs_p_id > 0 )
3516         ! CALL wrf_debug (wrfdbg, 'SPFire_remove removed, remaining: '//msg)
3517         !CALL wrf_error_fatal( "end of test")
3518     ENDIF
3519
3520 END SUBROUTINE remove_br
3521
3522 !=============================================================
3523 !=============================================================
3524
3525
3526 !=============================================================
3527 PURE SUBROUTINE deposit_br(fs_p_id, fs_p_x, fs_p_y, fs_p_z,  sum_fbrand, land_hgt)
3528 !=============================================================
3529
3530     IMPLICIT NONE
3531     INTEGER, INTENT(INOUT), DIMENSION(:) :: fs_p_id
3532     REAL,    INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: sum_fbrand
3533     REAL,    INTENT(IN),    DIMENSION(:) :: fs_p_x, fs_p_y, fs_p_z
3534     REAL,    INTENT(IN)    :: land_hgt
3535
3536     !INTEGER, INTENT(IN) :: active_br
3537
3538     !INTEGER, ALLOCATABLE, DIMENSION(:) :: mask_idx
3539     INTEGER, ALLOCATABLE, DIMENSION(:) :: rx, ry, rid
3540     LOGICAL, DIMENSION(fs_array_maxsize)    :: sparse_mask
3541     LOGICAL, DIMENSION(fs_array_maxsize)    :: bounds_mask
3542     !LOGICAL, ALLOCATABLE, DIMENSION(:) :: bounds_mask
3543     INTEGER :: nresets
3544     INTEGER :: x, y, i0,k ! counters
3545     !CHARACTER(LEN=10) fmt
3546
3547     !CALL wrf_debug (wrfdbg, 'SPFire_deposit: check')
3548     sum_fbrand(:,:) = ZERO_dp
3549
3550     !-------------------------------------------------------------------------------
3551     ! Flag brands near ground (deps_lev is a declared parameter)
3552     !-------------------------------------------------------------------------------
3553
3554     sparse_mask = .FALSE.
3555     bounds_mask = .FALSE.
3556
3557     ! Only deposit particles on this tile
3558     bounds_mask = ((FLOOR(fs_p_x) >= is) .OR. &
3559                    (FLOOR(fs_p_x) <= ie) .OR. &
3560                    (FLOOR(fs_p_y) >= js) .OR. &
3561                    (FLOOR(fs_p_y) <= je))
3562
3563     sparse_mask = ( bounds_mask .AND. (fs_p_id > 0 .AND. fs_p_z <= land_hgt))
3564
3565     !-------------------------------------------------------------------------------
3566     ! Deposit flagged brands
3567     !-------------------------------------------------------------------------------
3568
3569     nresets = COUNT(sparse_mask)
3570
3571     ! WRITE (msg,'(i6)') nresets
3572     ! CALL wrf_debug (wrfdbg, 'SPFire_deposit total: '//msg)
3573
3574     IF (nresets > 0) THEN
3575
3576         !-------------------------------------------------------------------------------
3577         ! convert position to integer indices
3578         !-------------------------------------------------------------------------------
3579
3580         ALLOCATE(rx(nresets),ry(nresets))!, mask_idx(nresets), rid(nresets), bounds_mask(nresets))
3581
3582         ! Get indices of particles when sparse_mask = True
3583         !mask_idx = idx_packed_1d(mask=sparse_mask)
3584         rx = 0
3585         ry = 0
3586
3587         ! rx(1:nresets) = MIN(ie, NINT(PACK(fs_p_x, sparse_mask)) ) ! Using Nint to destager grid
3588         ! ry(1:nresets) = MIN(je, NINT(PACK(fs_p_y, sparse_mask)) )
3589         rx(1:nresets) = FLOOR(PACK(fs_p_x, sparse_mask))
3590         ry(1:nresets) = FLOOR(PACK(fs_p_y, sparse_mask))
3591
3592
3593         ! Use indices to fetch x, y positions and id of depositing particles
3594         ! rid = 0
3595         ! rx(1:nresets) = FLOOR(fs_p_x(mask_idx))
3596         ! ry(1:nresets) = FLOOR(fs_p_y(mask_idx))
3597         ! rid(1:nresets) = fs_p_id(mask_idx)
3598
3599         ! !-------------------------------------------------------------------------------
3600         ! ! adjust depositing position to inside this tile
3601         ! !  (incrementing sum_fbrand at memory buffer leads to memory error)
3602         ! !-------------------------------------------------------------------------------
3603
3604         ! bounds_mask = .FALSE.
3605         ! bounds_mask = (rid > 0    .AND. &
3606         !                ((rx < is) .OR. (rx > ie) .OR. &
3607         !                 (ry < js) .OR. (ry > je)))
3608
3609         ! IF (COUNT(bounds_mask) > 0) THEN
3610
3611         !     WRITE (msg,'(1(i6,1x))') COUNT(bounds_mask)
3612         !     CALL wrf_debug (wrfdbg, 'SPFire_deposit adjust out of tile bounds'//msg)
3613
3614         !     WHERE (rid > 0 .AND. rx < is) rx = is
3615         !     WHERE (rid > 0 .AND. rx > ie) rx = ie
3616         !     WHERE (rid > 0 .AND. ry < js) ry = js
3617         !     WHERE (rid > 0 .AND. ry > je) ry = je
3618         ! ENDIF
3619
3620         !-------------------------------------------------------------------------------
3621         ! Increment sum_fbrand
3622         !-------------------------------------------------------------------------------
3623
3624         DO k=1,nresets
3625             sum_fbrand(rx(k),ry(k)) = sum_fbrand(rx(k),ry(k)) + 1.0_dp
3626         ENDDO
3627
3628         !-------------------------------------------------------------------------------
3629         ! Reset array
3630         !-------------------------------------------------------------------------------
3631
3632         WHERE(sparse_mask) fs_p_id= 0
3633
3634         ! WRITE (msg,'(5(i8,1x))') COUNT(sparse_mask), COUNT( fs_p_id > 0 ), COUNT(rx >= ie), COUNT(ry >= je)
3635         ! CALL wrf_debug (wrfdbg, 'SPFire_deposit deposited, remaining, outofbounds: '//msg)
3636
3637     ENDIF
3638
3639 END SUBROUTINE deposit_br
3640
3641
3642 !=============================================================
3643 !=============================================================
3644
3645
3646
3647  !------------------------------------------------------------------------------
3648  ! UTILITIES
3649  !------------------------------------------------------------------------------
3650
3651
3652  !=============================================================
3653
3654     SUBROUTINE get_local_ijk(grid, &
3655         ims, ime, jms, jme, kms, kme, &
3656         ips, ipe, jps, jpe, kps, kpe, &
3657         ifps, ifpe, jfps, jfpe,       &
3658         ifms, ifme, jfms, jfme,       &
3659         ids, ide, jds, jde, kds, kde, &
3660         m_idim,   m_jdim,   m_kdim,   &
3661         p_idim,   p_jdim,   p_kdim,   &
3662         d_idim,   d_jdim,   d_kdim)
3663         !p2m_is, p2m_ie, p2m_js, p2m_je)
3664
3665         USE module_domain, ONLY: get_ijk_from_grid, get_ijk_from_subgrid
3666
3667         IMPLICIT NONE
3668
3669         TYPE(DOMAIN), INTENT(IN) :: grid
3670         INTEGER,      INTENT(OUT), OPTIONAL :: ims, ime, jms, jme, kms, kme
3671         INTEGER,      INTENT(OUT), OPTIONAL :: ips, ipe, jps, jpe, kps, kpe
3672         INTEGER,      INTENT(OUT), OPTIONAL :: ifps, ifpe, jfps, jfpe
3673         INTEGER,      INTENT(OUT), OPTIONAL :: ifms, ifme, jfms, jfme
3674         INTEGER,      INTENT(OUT), OPTIONAL :: ids, ide, jds, jde, kds, kde
3675         INTEGER,      INTENT(OUT), OPTIONAL :: m_idim,   m_jdim,   m_kdim
3676         INTEGER,      INTENT(OUT), OPTIONAL :: p_idim,   p_jdim,   p_kdim
3677         INTEGER,      INTENT(OUT), OPTIONAL :: d_idim,   d_jdim,   d_kdim
3678         !INTEGER,      INTENT(OUT), OPTIONAL :: p2m_is, p2m_ie, p2m_js, p2m_je
3679
3680
3681         INTEGER :: iims, iime, jjms, jjme, kkms, kkme
3682         INTEGER :: iips, iipe, jjps, jjpe, kkps, kkpe
3683         INTEGER :: iids, iide, jjds, jjde, kkds, kkde
3684
3685         INTEGER :: iifps, iifpe, jjfps, jjfpe, kkfps, kkfpe
3686         INTEGER :: iifms, iifme, jjfms, jjfme, kkfms, kkfme
3687         INTEGER :: iifds, iifde, jjfds, jjfde, kkfds, kkfde
3688
3689         IF ((.NOT. PRESENT(ims)) .AND. &
3690             (.NOT. PRESENT(jms)) .AND. &
3691             (.NOT. PRESENT(kms)) .AND. &
3692             (.NOT. PRESENT(ime)) .AND. &
3693             (.NOT. PRESENT(jme)) .AND. &
3694             (.NOT. PRESENT(kme)) .AND. &
3695             !
3696             (.NOT. PRESENT(ips)) .AND. &
3697             (.NOT. PRESENT(jps)) .AND. &
3698             (.NOT. PRESENT(kps)) .AND. &
3699             (.NOT. PRESENT(ipe)) .AND. &
3700             (.NOT. PRESENT(jpe)) .AND. &
3701             (.NOT. PRESENT(kpe)) .AND. &
3702             !
3703             (.NOT. PRESENT(ifps)) .AND. &
3704             (.NOT. PRESENT(jfps)) .AND. &
3705             (.NOT. PRESENT(ifpe)) .AND. &
3706             (.NOT. PRESENT(jfpe)) .AND. &
3707             !
3708             (.NOT. PRESENT(ifms)) .AND. &
3709             (.NOT. PRESENT(jfms)) .AND. &
3710             (.NOT. PRESENT(ifme)) .AND. &
3711             (.NOT. PRESENT(jfme)) .AND. &
3712             !
3713             (.NOT. PRESENT(ids)) .AND. &
3714             (.NOT. PRESENT(jds)) .AND. &
3715             (.NOT. PRESENT(kds)) .AND. &
3716             (.NOT. PRESENT(ide)) .AND. &
3717             (.NOT. PRESENT(jde)) .AND. &
3718             (.NOT. PRESENT(kde)) .AND. &
3719             !
3720             ! (.NOT. PRESENT(p2m_is)) .AND. &
3721             ! (.NOT. PRESENT(p2m_ie)) .AND. &
3722             ! (.NOT. PRESENT(p2m_js)) .AND. &
3723             ! (.NOT. PRESENT(p2m_je)) .AND. &
3724             !
3725             (.NOT. PRESENT(m_idim)) .AND. &
3726             (.NOT. PRESENT(m_jdim)) .AND. &
3727             (.NOT. PRESENT(m_kdim)) .AND. &
3728             (.NOT. PRESENT(p_idim)) .AND. &
3729             (.NOT. PRESENT(p_jdim)) .AND. &
3730             (.NOT. PRESENT(p_kdim)) .AND. &
3731             (.NOT. PRESENT(d_idim)) .AND. &
3732             (.NOT. PRESENT(d_jdim)) .AND. &
3733             (.NOT. PRESENT(d_kdim))) &
3734             !
3735             CALL wrf_debug ( 1, 'get_local_ijk function is NOT requesting a result' )
3736
3737         CALL get_ijk_from_grid (  grid ,        &
3738             iids, iide, jjds, jjde, kkds, kkde, &
3739             iims, iime, jjms, jjme, kkms, kkme, &
3740             iips, iipe, jjps, jjpe, kkps, kkpe)
3741
3742         IF (PRESENT(ifps) .OR. &
3743             PRESENT(jfps) .OR. &
3744             PRESENT(ifpe) .OR. &
3745             PRESENT(jfpe) .OR. &
3746             PRESENT(ifms) .OR. &
3747             PRESENT(jfms) .OR. &
3748             PRESENT(ifme) .OR. &
3749             PRESENT(jfme)) CALL get_ijk_from_subgrid(grid , &
3750                                     iifds, iifde, jjfds, jjfde, kkfds, kkfde, &
3751                                     iifms, iifme, jjfms, jjfme, kkfms, kkfme, &
3752                                     iifps, iifpe, jjfps, jjfpe, kkfps, kkfpe)
3753
3754
3755         IF (PRESENT(ims)) ims = iims
3756         IF (PRESENT(jms)) jms = jjms
3757         IF (PRESENT(kms)) kms = kkms
3758         IF (PRESENT(ime)) ime = iime
3759         IF (PRESENT(jme)) jme = jjme
3760         IF (PRESENT(kme)) kme = kkme
3761
3762         IF (PRESENT(ips)) ips = iips
3763         IF (PRESENT(jps)) jps = jjps
3764         IF (PRESENT(kps)) kps = kkps
3765         IF (PRESENT(ipe)) ipe = iipe
3766         IF (PRESENT(jpe)) jpe = jjpe
3767         IF (PRESENT(kpe)) kpe = kkpe
3768
3769         IF (PRESENT(ifps)) ifps = iifps
3770         IF (PRESENT(jfps)) jfps = jjfps
3771         IF (PRESENT(ifpe)) ifpe = iifpe
3772         IF (PRESENT(jfpe)) jfpe = jjfpe
3773
3774         IF (PRESENT(ifms)) ifms = iifms
3775         IF (PRESENT(jfms)) jfms = jjfms
3776         IF (PRESENT(ifme)) ifme = iifme
3777         IF (PRESENT(jfme)) jfme = jjfme
3778
3779         IF (PRESENT(ids)) ids = iids
3780         IF (PRESENT(jds)) jds = jjds
3781         IF (PRESENT(kds)) kds = kkds
3782         IF (PRESENT(ide)) ide = iide
3783         IF (PRESENT(jde)) jde = jjde
3784         IF (PRESENT(kde)) kde = kkde
3785
3786         IF (PRESENT(m_idim))  m_idim = iime - iims  + 1
3787         IF (PRESENT(m_jdim))  m_jdim = jjme - jjms  + 1
3788         IF (PRESENT(m_kdim))  m_kdim = kkme - kkms  + 1
3789         IF (PRESENT(p_idim))  p_idim = iipe - iips  + 1
3790         IF (PRESENT(p_jdim))  p_jdim = jjpe - jjps  + 1
3791         IF (PRESENT(p_kdim))  p_kdim = kkpe - kkps  + 1
3792         IF (PRESENT(d_idim))  d_idim = iide - iids  + 1
3793         IF (PRESENT(d_jdim))  d_jdim = jjde - jjds  + 1
3794         IF (PRESENT(d_kdim))  d_kdim = kkde - kkds  + 1
3795
3796         ! IF (PRESENT(p2m_is)) p2m_is = iips - iims
3797         ! IF (PRESENT(p2m_ie)) p2m_ie = iipe - iims
3798         ! IF (PRESENT(p2m_js)) p2m_js = jjps - jjms
3799         ! IF (PRESENT(p2m_je)) p2m_je = jjpe - jjms
3800
3801         ! p2m_is = ips - ims
3802         ! p2m_ie = ips - ims + p_idim - 1
3803         ! p2m_js = jps - jms
3804         ! p2m_je = jps - jms + p_jdim - 1
3805
3806     END SUBROUTINE get_local_ijk
3807
3808 !=============================================================
3809 !=============================================================
3810
3811
3812
3813
3814 END MODULE module_firebrand_spotting
3815