anchor.diff:      end do
       call md%make_initial_conditions_consistent()
       
  -    ! NOTE: For stability in 'null' regions, we set them to 'high land' that
  -    ! should be inactive. 
  +    ! For stability in 'null' regions, we set them to 'high land' that should be inactive. 
       call md%set_null_regions_to_dry()
      
       write(log_output_unit,*) 'End setup'
  
positive.code:         real(dp), allocatable:: x(:), y(:)
          type(multi_raster_type):: elevation_data, stage_data
          real(dp), allocatable :: random_uniform(:)
          ! Use this to check the effect of very small perturbations of the initial condition
          ! Eventually  a tiny perturbation will lead to differences in the results. This seems comparable
          ! to the differences that eventually emerge in DEFAULT openmp/coarray versions, which seem
          ! to begin with differences of this order in elevation/stage. That case is due to differences in
          ! domain partitioning, and can be eliminated by providing a load_balance_partition.txt file
          real(dp), parameter :: random_perturbation_scale = 0.0e-10_dp
  
          ! Stage
  
negative.code:     call generic_stop
  #endif
      
      !! nd domains in this model
      !if(very_high_res_monai) then
      !    ! In this case we put a 30cm x 30cm cell domain around the monai inundation peak (with mesh_refine=1.0_dp)
      !    !
      !    ! It leads to a peak of > 30.0m with mesh_refine=1.0_dp (vs 31.7 obs, although the latter varies depending on which
      !    ! dataset is used).
      !    !
      !    ! The model needs various modifications to do this stably, and it takes more than twice as long.
      !    nd = 7
      !else
          ! This case has a domain res ~ 1.5x1.5m around the monai inundation peak (with mesh_refine = 1.0_dp)
          !
          ! It leads to a peak of >28.0m with mesh_refine=1.0_dp (vs 31.7 obs, although the latter varies depending on which
          ! dataset is used). So the very_high_res_monai case does better.
          nd = 6
      !end if
  
      allocate(md%domains(nd))
  
      ! Need to provide a load_balance_partition.txt file to get exact reproducibility 
      ! with varying number of openmp threads & MPI ranks
      call get_command_argument(1, md%load_balance_file)
  
      !
  