anchor.diff:      !! NTHMP benchmark problem 9 -- Okushiri tsunami field test case
       !!
   
  -    use global_mod, only: ip, dp, minimum_allowed_depth, charlen, default_nonlinear_timestepping_method,&
  +    use global_mod, only: ip, dp, charlen, default_nonlinear_timestepping_method,&
           default_linear_timestepping_method
  -    use domain_mod, only: domain_type
  -    use multidomain_mod, only: multidomain_type, setup_multidomain, test_multidomain_mod
  -    use boundary_mod, only: flather_boundary, transmissive_boundary
  -    use local_routines
  -    use timer_mod
  +    use multidomain_mod, only: multidomain_type
  +    use boundary_mod, only: flather_boundary
  +    use timer_mod, only: timer_type
       use logging_mod, only: log_output_unit, send_log_output_to_file
       use stop_mod, only: generic_stop
  -    use iso_c_binding, only: C_DOUBLE !, C_INT, C_LONG
       use coarray_intrinsic_alternatives, only : swals_mpi_init, swals_mpi_finalize
  +    use iso_c_binding, only: C_DOUBLE
  +    use local_routines
   
       implicit none
   
  
positive.code:     real(dp), parameter :: mesh_refine = 0.4_dp ! Good for more efficient testing
      !real(dp), parameter :: mesh_refine = 1.0_dp ! Standard
  
      ! Optionally put a "very high res" domain around monai. Slows things down and the code requires care
      ! If this is true, then also use "mesh_refine = 1.0_dp"
      !logical, parameter:: very_high_res_monai = .true. 
      logical, parameter:: very_high_res_monai = .false.
  
      ! If using very high res domain, then reduce the entire model timestep when it becomes active
      real(dp), parameter :: very_high_res_timestep_reduction = 10.0_dp
      ! Assume(require) the very high res domain is dry before this time. If not set correctly, expect the wrong answer!
      real(dp), parameter :: very_high_res_static_before_time = 235.0_dp
  
      ! Approx timestep between outputs
      real(dp), parameter :: approximate_writeout_frequency = 7.50_dp !0.6_dp * (1.0_dp/0.4_dp ) * 1.0_dp/3.0_dp!7.50_dp
      real(dp), parameter :: final_time = 3600.0_dp * 1.0_dp
      !real(dp), parameter :: final_time = 30.0_dp * 1.0_dp
      !real(dp), parameter :: approximate_writeout_frequency = 0.2_dp
      !real(dp), parameter :: final_time = 10.0_dp
  
      ! Length/width
      real(dp), parameter, dimension(2):: global_lw = [3.9_dp, 4.6_dp]
      ! Lower-left corner coordinate
      real(dp), parameter, dimension(2):: global_ll = [137.6_dp, 39.6_dp]
      ! grid size (number of x/y cells). Note that dx is in lon/lat, so an increment of dx(1) is smaller than dx(2)
      integer(ip), parameter, dimension(2):: global_nx = nint(global_lw*[77, 100]*1.12_dp*mesh_refine)
  
      ! Refinement factor for domains
      integer(ip), parameter :: nest_ratio = 5_ip
      ! The global (i.e. outer-domain) time-step in the multidomain 
      real(dp) ::  global_dt = 0.60_dp * (1.0_dp/mesh_refine) * (1.0_dp / 3.0_dp)
  
      ! Useful misc variables
      integer(ip):: j, nd
  
negative.code:         ! Ensure stage >= elevation
          domain%U(:,:,STG) = max(domain%U(:,:,STG), domain%U(:,:,ELV) + 1.0e-07_dp)
  
          write(log_output_unit,*) 'Stage range is: ', minval(domain%U(:,:,STG)), maxval(domain%U(:,:,STG))
          write(log_output_unit,*) 'Elev range is: ', minval(domain%U(:,:,ELV)), maxval(domain%U(:,:,ELV))
  
  