anchor.diff:  module local_routines 
       use global_mod, only: dp, ip, wall_elevation
  -    use domain_mod, only: domain_type
  +    use domain_mod, only: domain_type, STG, UH, VH, ELV
       implicit none
       
   
  
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
              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))
  
  
          domain%manning_squared = 0.03_dp**2
  
  
negative.code: 
  
      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(:,:,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
  