anchor.diff:          ! Ensure stage >= elevation
           domain%U(:,:,STG) = max(domain%U(:,:,STG), domain%U(:,:,ELV) + 1.0e-07_dp)
   
  -        !write(log_output_unit,*) 'DEBUG: SETTING STAGE TO ZERO'
  -        !domain%U(:,:,STG) = max(0.0_dp, 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))
   
  
positive.code:     end do
      call md%make_initial_conditions_consistent()
      
      ! NOTE: 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'
  
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)
  
      !
  