anchor.diff:      !! Setup uniform channel problem (aligned to grid)
       !!
       use global_mod, only: dp, ip, wall_elevation
  -    use domain_mod, only: domain_type
  +    use domain_mod, only: domain_type, STG, UH, VH, ELV
       use which_mod, only: which
       implicit none
       
  
positive.code: 
      ! Find channel indices 
      ! This corresponds to the lowest points on a horizontal cross-section
      call which(md%domains(1)%U(:,2,4) < (minval(md%domains(1)%U(:,2,4)) + 1.0e-06), channel_x_indices)
      print*, channel_x_indices
  
      ! Time-step at which we evolve the solution
  
negative.code:             do j = 1, domain%nx(2)
                  x = (i-0.5_dp)*domain%dx(1) - cx
                  y = (j-0.5_dp)*domain%dx(2) - cy 
                  domain%U(i,j,ELV) = y*slope
                  if(x > -channel_width/2.0_dp .and. x < channel_width/2.0_dp) then
                      domain%U(i,j,ELV) = domain%U(i,j,ELV) - channel_depth 
                  end if
              end do
          end do
          ! Wall boundaries along the sides and top
          domain%U(1,:,ELV) = wall_elevation 
          domain%U(domain%nx(1),:,ELV) = wall_elevation 
          domain%U(:,domain%nx(2),ELV) = wall_elevation 
          domain%U(:,1,ELV) = wall_elevation 
  
          ! Stage
          domain%U(:,:,STG) = domain%U(:,:,ELV)
  
          ! Ensure stage >= elevation
          domain%U(:,:,STG) = max(domain%U(:,:,STG), domain%U(:,:,ELV))
  
          ! Friction
          domain%manning_squared = 0.03_dp**2
  