anchor.diff:          print*, 'Elevation range: ', minval(domain%U(:,:,ELV)), &
               maxval(domain%U(:,:,ELV))
   
  -        ! A bit unclear how this boundary should be treated -- but works ok
  +        ! A bit unclear how this boundary should be treated -- a wall works ok
           wall = 20._dp
           domain%U((domain%nx(1)-1):domain%nx(1),:,ELV) = wall
   
           ! Friction 
           if(domain%timestepping_method /= 'linear') then
  -            domain%manning_squared = 0.025_dp * 0.025_dp
  +            domain%manning_squared = 0.025_dp**2
           end if
   
           ! Ensure stage >= elevation
  
positive.code:         end select
          md%domains(j)%boundary_function => boundary_function
      end do
     
      ! Print the gravity-wave CFL limit, to guide timestepping
      do j = 1, size(md%domains)
  
negative.code:     ! Build boundary conditions
      call setup_boundary_information(bc_file, bc_elev)
  
      do j = 1, size(md%domains) 
          select case(boundary_type)
          case('flather_with_vh_from_continuity')
  