anchor.diff:              ! stage forcing. (Not too surprising, because the boundary is quite close to the
               ! coast, in 'not very deep' water).
   
  -            !! Approach 1
  -            !! Transmissive momentum. Unstable with flather?
  -            !!
  -            !! FIXME: Currently this doesn't carefully treat 'corners' (e.g. i==1, j==1).
  -            !! Could be an issue depending on boundary setup
  -            !if(j == domain%nx(2)) stage_uh_vh_elev(2:3) = domain%U(i,j-1,2:3)
  -            !if(j == 1) stage_uh_vh_elev(2:3) = domain%U(i,j+1,2:3)
  -            !if(i == 1) stage_uh_vh_elev(2:3) = domain%U(i+1,j,2:3)
  -            !if(i == domain%nx(1)) stage_uh_vh_elev(2:3) = domain%U(i-1,j,2:3)
  -            !
  -            !
  -            ! Approach 3:
  -            ! Assume the incoming wave is like a plane wave, so that vh = sqrt(g * depth) * stage
  -            ! This is wrong so close to shore, and the wave is over-amplified. But in 'really deep water'
  -            ! where the coast was > 1 wave-length away, I suppose this approach might work OK?
  -            !if(j == domain%nx(2)) then
  -            !    stage_uh_vh_elev(2) = 0.0_dp
  -            !    ! Linear wave, assuming MSL = 0
  -            !    stage_uh_vh_elev(3) = -sqrt(9.81_dp * (-1.0_dp*local_elev)) * stage_uh_vh_elev(1)
  -            !else
  -            !    stage_uh_vh_elev(2:3) = 0.0_dp
  -            !end if
  -
  -
               select case(boundary_type)
   
               case('boundary_stage_transmissive_normal_momentum')
  
positive.code:         md%domains(1)%boundary_subroutine => boundary_stage_transmissive_normal_momentum
      case('boundary_stage_transmissive_momentum')
          md%domains(1)%boundary_subroutine => boundary_stage_transmissive_momentum
      !case('boundary_stage_transmissive_momentum_sponge')
      !    md%domains(1)%boundary_subroutine => boundary_stage_transmissive_momentum_sponge
      case('flather_with_vh_from_continuity')
          md%domains(1)%boundary_subroutine => flather_boundary
      case('flather_with_vh_equal_zero')
  
negative.code: 
          !
          ! The DEM needs to be 'fixed' in a few places where bridges remain. Google earth    
          ! suggests the bridges should not strongly impede the flow. So based on checks of the DEM, ....
          !
          ! pol1 should have a value about -6.0
          !
  