diff --git a/physics/fire_model_mod.F90 b/physics/fire_model_mod.F90 index e5b9ae5..4efa9d7 100644 --- a/physics/fire_model_mod.F90 +++ b/physics/fire_model_mod.F90 @@ -2,7 +2,7 @@ module fire_model_mod use fire_physics_mod, only: Calc_flame_length, Calc_fire_fluxes, Calc_smoke_emissions use level_set_mod, only: Calc_fuel_left, Update_ignition_times, Reinit_level_set, Prop_level_set, Extrapol_var_at_bdys, & - Stop_if_close_to_bdy, Copy_lfnout_to_lfn, Reinit_level_set_fast_dist + Stop_if_close_to_bdy, Copy_lfnout_to_lfn, Reinit_level_set_fast_dist, Check_isolated_negative_lfn use namelist_mod, only : namelist_t use ros_mod, only : ros_t use state_mod, only: state_fire_t, N_POINTS_IN_HALO @@ -29,6 +29,7 @@ subroutine Advance_fire_model (config_flags, grid) integer :: ij, ifds, ifde, jfds, jfde, ifts, ifte, jfts, jfte, ifms, ifme, jfms, jfme real :: tbound, time_start + integer, parameter :: CHECK_ISOLATED_NEG_LFN = 0 logical, parameter :: DEBUG_LOCAL = .false. @@ -130,6 +131,8 @@ subroutine Advance_fire_model (config_flags, grid) #ifdef DM_PARALLEL call Do_halo_exchange_with_corners (grid%lfn, ifms, ifme, jfms, jfme, grid%ifps, grid%ifpe, grid%jfps, grid%jfpe, N_POINTS_IN_HALO, grid%cart_comm) #endif + + if (CHECK_ISOLATED_NEG_LFN == 1) call Check_isolated_negative_lfn (grid) if (DEBUG_LOCAL) call Print_message ('calling Ignite_prescribed_fires...') !$OMP PARALLEL DO & diff --git a/physics/level_set_mod.F90 b/physics/level_set_mod.F90 index fb1077b..0f3d35e 100644 --- a/physics/level_set_mod.F90 +++ b/physics/level_set_mod.F90 @@ -20,7 +20,7 @@ module level_set_mod use ros_mod, only : ros_t #ifdef DM_PARALLEL - use mpi_mod, only : Do_halo_exchange + use mpi_mod, only : Do_halo_exchange, Sum_across_mpi_tasks, Max_across_mpi_tasks #endif implicit none @@ -28,7 +28,7 @@ module level_set_mod private public :: Calc_fuel_left, Update_ignition_times, Reinit_level_set, Prop_level_set, Extrapol_var_at_bdys, Stop_if_close_to_bdy, & - Copy_lfnout_to_lfn, Reinit_level_set_fast_dist + Copy_lfnout_to_lfn, Reinit_level_set_fast_dist, Check_isolated_negative_lfn integer, parameter :: BDY_ENO1 = 10, FAST_DIST_REINIT_FSM = 1 integer, parameter :: FAR = 0, TRIAL = 1, KNOWN = 2 @@ -291,6 +291,137 @@ pure subroutine Copy_lfnout_to_lfn (ifts, ifte, jfts, jfte, ifms, ifme, jfms, jf end subroutine Copy_lfnout_to_lfn + subroutine Check_isolated_negative_lfn (grid, threshold, min_size, max_size) + + implicit none + + type (state_fire_t), intent (in out) :: grid + real, intent (in), optional :: threshold + integer, intent (in), optional :: min_size, max_size + + character (len = 256) :: msg + integer :: cluster_i, cluster_j, cluster_size, flagged_size + integer :: head, tail, i, ic, inb, j, jc, jnb, max_queue, min_size_value, max_size_value + integer :: nx_loc, ny_loc + logical :: found_cluster, touches_boundary + logical, dimension(:, :), allocatable :: visited + integer, dimension(:), allocatable :: queue_i, queue_j + real :: detection_flag, detection_sum, global_location_i, global_location_j, global_size_value + real :: location_i, location_j, size_value, threshold_value + + + threshold_value = 0.0 + if (present (threshold)) threshold_value = threshold + + min_size_value = 1 + if (present (min_size)) min_size_value = min_size + + max_size_value = 6 + if (present (max_size)) max_size_value = max_size + + if (max_size_value < min_size_value) return + + nx_loc = grid%ifpe - grid%ifps + 1 + ny_loc = grid%jfpe - grid%jfps + 1 + if (nx_loc <= 0 .or. ny_loc <= 0) return + max_queue = nx_loc * ny_loc + + allocate (visited(grid%ifps:grid%ifpe, grid%jfps:grid%jfpe)) + allocate (queue_i(max_queue)) + allocate (queue_j(max_queue)) + + visited = .false. + found_cluster = .false. + flagged_size = 0 + cluster_i = grid%ifps + cluster_j = grid%jfps + + Cluster_search: do j = grid%jfps, grid%jfpe + do i = grid%ifps, grid%ifpe + if (visited(i, j)) cycle + visited(i, j) = .true. + if (grid%lfn(i, j) >= threshold_value) cycle + + head = 1 + tail = 1 + queue_i(1) = i + queue_j(1) = j + cluster_size = 0 + touches_boundary = .false. + + do while (head <= tail) + ic = queue_i(head) + jc = queue_j(head) + head = head + 1 + + cluster_size = cluster_size + 1 + if (ic == grid%ifps .or. ic == grid%ifpe .or. jc == grid%jfps .or. jc == grid%jfpe) touches_boundary = .true. + + do jnb = jc - 1, jc + 1 + do inb = ic - 1, ic + 1 + if (inb == ic .and. jnb == jc) cycle + if (inb < grid%ifps .or. inb > grid%ifpe) cycle + if (jnb < grid%jfps .or. jnb > grid%jfpe) cycle + if (visited(inb, jnb)) cycle + visited(inb, jnb) = .true. + if (grid%lfn(inb, jnb) >= threshold_value) cycle + + tail = tail + 1 + queue_i(tail) = inb + queue_j(tail) = jnb + end do + end do + end do + + if (cluster_size >= min_size_value .and. cluster_size <= max_size_value .and. .not. touches_boundary) then + found_cluster = .true. + flagged_size = cluster_size + cluster_i = i + cluster_j = j + exit Cluster_search + end if + end do + end do Cluster_search + + deallocate (visited) + deallocate (queue_i) + deallocate (queue_j) + + detection_flag = 0.0 + location_i = -1.0 + location_j = -1.0 + size_value = -1.0 + if (found_cluster) then + detection_flag = 1.0 + location_i = real (cluster_i) + location_j = real (cluster_j) + size_value = real (flagged_size) + end if + +#ifdef DM_PARALLEL + call Sum_across_mpi_tasks (detection_flag, grid%cart_comm, detection_sum) + call Max_across_mpi_tasks (location_i, grid%cart_comm, global_location_i) + call Max_across_mpi_tasks (location_j, grid%cart_comm, global_location_j) + call Max_across_mpi_tasks (size_value, grid%cart_comm, global_size_value) +#else + detection_sum = detection_flag + global_location_i = location_i + global_location_j = location_j + global_size_value = size_value +#endif + + if (detection_sum > 0.0) then + grid%datetime_now = grid%datetime_start + call grid%datetime_now%Add_seconds (grid%itimestep * grid%dt) + call grid%Save_state () + + write (msg, '(a, i6, a, i6, a, i6, a)') 'Isolated negative LFN cluster (size', int (global_size_value), ') near (i=', & + int (global_location_i), ', j=', int (global_location_j), '). Simulation stopped after saving state.' + call Stop_simulation (trim (msg)) + end if + + end subroutine Check_isolated_negative_lfn + subroutine Prop_level_set (ifds, ifde, jfds, jfde, ifms, ifme, jfms, jfme, & num_tiles, i_start, i_end, j_start, j_end, ts, dt, dx, dy, fire_upwinding, fire_viscosity, & fire_viscosity_bg, fire_viscosity_band, fire_viscosity_ngp, fire_lsm_band_ngp, & @@ -1667,34 +1798,25 @@ end function Select_godunov pure function Select_eno (diff_lx, diff_rx) result (return_value) - ! 1st order ENO scheme + ! 1st order ENO scheme: choose the smaller magnitude derivative + ! Both positive: pick smaller + ! Both negative: pick larger (less negative) + ! Mixed signs: zero (derivative crosses zero) implicit none - real, intent (in):: diff_lx, diff_rx + real, intent(in) :: diff_lx, diff_rx real :: return_value - real :: diff2x - - if (.not. diff_lx > 0.0 .and. .not. diff_rx > 0.0) then - diff2x = diff_rx - else if (.not. diff_lx < 0.0 .and. .not. diff_rx < 0.0) then - diff2x = diff_lx - else if (.not. diff_lx < 0.0 .and. .not. diff_rx > 0.0) then - if (.not. abs (diff_rx) < abs(diff_lx)) then - diff2x = diff_rx - else - diff2x = diff_lx - end if + if (diff_lx * diff_rx > 0.0) then + return_value = sign(1.0, diff_lx) * min(abs(diff_lx), abs(diff_rx)) else - diff2x = 0.0 + return_value = 0.0 end if - return_value = diff2x - end function Select_eno - + pure function Select_2nd (dx, lfn_i, lfn_im1, lfn_ip1) result (return_value) ! 2nd-order advection scheme in the x,y-direction (DME) diff --git a/physics/ros_wrffire_mod.F90 b/physics/ros_wrffire_mod.F90 index 788e89f..b8cc25f 100644 --- a/physics/ros_wrffire_mod.F90 +++ b/physics/ros_wrffire_mod.F90 @@ -88,15 +88,7 @@ pure function Calc_ros_wrffire (this, ifms, ifme, jfms, jfme, i, j, nvx, nvy, uf ros_wind = ros_wind * cor_wind ros_slope = ros_slope * cor_slope - ! Limit the ros - excess = ros_base + ros_wind + ros_slope - ROS_MAX - if (excess > 0.0) then - ! take excess out of wind and slope in proportion - ros_wind = ros_wind - excess * ros_wind / (ros_wind + ros_slope) - ros_slope = ros_slope - excess * ros_slope/ (ros_wind + ros_slope) - end if - - return_value = ros_base + ros_wind + SLOPE_FACTOR * ros_slope + return_value = min (ros_base + ros_wind + SLOPE_FACTOR * ros_slope, ROS_MAX) if (FIRE_GROWS_ONLY) return_value = max (return_value, 0.0) end function Calc_ros_wrffire