anchor.diff:          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')
  
positive.code:         ! Ensure stage >= elevation
          domain%U(:,:,STG) = max(domain%U(:,:,STG), domain%U(:,:,ELV) + 1.0e-07_dp)
  
          !! Gauges -- instead we set them from the csv file
          !gauge_xy(1:3, 1) = [2.724e+04_dp, 1.846e+04_dp, 1.0_dp]
          !gauge_xy(1:3, 2) = [3.085e+04_dp, 1.512e+04_dp, 2.0_dp]
          !gauge_xy(1:3, 3) = [3.200e+04_dp, 1.347e+04_dp, 3.0_dp]
          !gauge_xy(1:3, 4) = [3.005e+04_dp, 1.610e+04_dp, 4.0_dp]
          !! ADCP location -- reported
          !gauge_xy(1:3, 5) = [2.9250e+04_dp, 1.4660e+04_dp, 5.0_dp]
          !! ADCP location -- this is nearby, and the velocities seem to agree 
          !! better (this place has the highest velocities in the entrance)
          !gauge_xy(1:3, 6) = [2.9325e+04_dp, 1.4525e+04_dp, 6.0_dp]
          !call domain%setup_point_gauges(xy_coords = gauge_xy(1:2,:), gauge_ids=gauge_xy(3,:))
  
      end subroutine
  
  end module 
  
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
          !
  