anchor.diff:  
       ! 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)
  +    call which(md%domains(1)%U(:,2,ELV) < (minval(md%domains(1)%U(:,2,ELV)) + 1.0e-06), channel_x_indices)
       print*, channel_x_indices
   
       ! Time-step at which we evolve the solution
  
positive.code:     !! Setup uniform channel problem (aligned to grid)
      !!
      use global_mod, only: dp, ip, wall_elevation
      use domain_mod, only: domain_type
      use which_mod, only: which
      implicit none
      
  
negative.code:         call md%evolve_one_step(global_dt)
  
          ! Add discharge inflow
          md%domains(1)%U(channel_x_indices, md%domains(1)%nx(2)-1, STG) = &
              md%domains(1)%U(channel_x_indices, md%domains(1)%nx(2)-1, STG) + &
              Qin*global_dt/product(md%domains(1)%dx)
  
      end do
  