anchor.diff:  
   
       theoretical_vol = Qin * domain%time * (domain%nx(1) - 2)!* domain%lw(1) * (global_nx(1) - 2) * 1.0_dp / global_nx(1)
  -    model_vol = sum(domain%U(:,:,1) - domain%U(:,:,4)) * domain%dx(1) * domain%dx(2)
  +    model_vol = sum(domain%U(:,:,STG) - domain%U(:,:,ELV)) * domain%dx(1) * domain%dx(2)
       print*, 'Cell area: ', domain%dx, product(domain%dx)
       ! Analytical solution
       ! vd = discharge_per_unit_width
  
positive.code: 
  program uniform_slope
      use global_mod, only: ip, dp, charlen, gravity
      use domain_mod, only: domain_type
      use file_io_mod, only: read_csv_into_array
      use local_routines
      implicit none
  
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
              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))
  
          domain%manning_squared = 0.03_dp**2
  
  