anchor.diff:          call md%evolve_one_step(global_dt)
   
           ! Add discharge inflow
  -        md%domains(1)%U(channel_x_indices, md%domains(1)%nx(2)-1, 1) = &
  -            md%domains(1)%U(channel_x_indices, md%domains(1)%nx(2)-1, 1) + &
  +        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
  
positive.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,4) = y*slope
                  if(x > -channel_width/2.0_dp .and. x < channel_width/2.0_dp) then
                      domain%U(i,j,4) = domain%U(i,j,4) - channel_depth 
                  end if
              end do
          end do
          ! Wall boundaries along the sides and top
          domain%U(1,:,4) = wall_elevation 
          domain%U(domain%nx(1),:,4) = wall_elevation 
          domain%U(:,domain%nx(2),4) = wall_elevation 
          domain%U(:,1,4) = wall_elevation 
  
          ! Stage
          domain%U(:,:,1) = domain%U(:,:,4)
  
          ! Ensure stage >= elevation
          domain%U(:,:,1 ) = max(domain%U(:,:,1), domain%U(:,:,4))
  
          ! Friction
          domain%manning_squared = 0.03_dp**2
  
negative.code: 
      ! Mass conservation check
      theoretical_vol = Qin * md%domains(1)%time * channel_width
      model_vol = sum(md%domains(1)%U(:,:,STG) - md%domains(1)%U(:,:,ELV)) * md%domains(1)%dx(1) * md%domains(1)%dx(2)
  
      ! Analytical solution for uniform slope
      ! vd = discharge_per_unit_width
  