anchor.diff:      use global_mod, only: ip, dp, minimum_allowed_depth, default_nonlinear_timestepping_method
       use domain_mod, only: domain_type
       use multidomain_mod, only: multidomain_type, setup_multidomain, test_multidomain_mod
  -    use boundary_mod, only: boundary_stage_transmissive_normal_momentum, flather_boundary, &
  -        boundary_stage_transmissive_momentum, boundary_stage_radiation_momentum
  +    use boundary_mod, only: flather_boundary, boundary_stage_radiation_momentum
       use local_routines
       use timer_mod
       use logging_mod, only: log_output_unit
  
positive.code:             stage_uh_vh_elev(2:3) = 0.0_dp
          else
  
              ! Much experimentation was conducted here. This problem is sensitive to
              ! imperfections in our semi-transmissive boundary conditions which allow a 
              ! stage forcing. (Not too surprising, because the boundary is quite close to the
              ! coast, in 'not very deep' water).
  
              select case(boundary_type)
  
              case('boundary_stage_transmissive_normal_momentum')
                  ! Do nothing, because we do not need uh/vh.
                  stage_uh_vh_elev(2:3) = 0.0_dp 
  
              case('boundary_stage_transmissive_momentum')
                  ! Do nothing, because we do not need uh/vh.
                  stage_uh_vh_elev(2:3) = 0.0_dp 
  
              case('boundary_stage_radiation_momentum') 
                  ! These will never be used for this boundary
                  stage_uh_vh_elev(2:3) = 0.0_dp
  
              case('flather_with_vh_equal_zero') 
                  ! This absorbs, but also distorts the stage at the boundary too much when
                  ! wave frequencies are lower
                  stage_uh_vh_elev(2:3) = 0.0_dp
  
              case('flather_with_vh_from_continuity')
  
negative.code:     call md%set_point_gauges_from_csv("point_gauges.csv", skip_header=1_ip)
     
      do j = 1, size(md%domains) 
          select case(boundary_type)
          case('flather_with_vh_from_continuity')
              md%domains(j)%boundary_subroutine => flather_boundary
          case('boundary_stage_radiation_momentum')
              md%domains(j)%boundary_subroutine => boundary_stage_radiation_momentum
  