anchor.diff:          boundary_information%boundary_data(nr - skip + extra,2) = &
               boundary_information%boundary_data(nr - skip + extra - 1, 2)
   
  -        !print*, 'ZEROING STAGE '
  -        !boundary_information$boundary_data(:,2) = 0.0_dp
           call boundary_information%gauge4_ts_function%initialise(&
  -                boundary_information%boundary_data(:,1), boundary_information%boundary_data(:,2))
  +            boundary_information%boundary_data(:,1), boundary_information%boundary_data(:,2))
           boundary_information%t0 = boundary_information%boundary_data(1,1)
   
       end subroutine
  -    
  -    ! Make a function to evaluate the boundary at the domain
  +   
  +    ! 
  +    ! Evaluate the boundary condition -- this will be called inside the domain%boundary_subroutine
       !
       function boundary_function(domain, t, i, j) result(stage_uh_vh_elev)
           type(domain_type), intent(in):: domain
  
positive.code:         dx_refinement_factor=nest_ratio, &
          timestepping_refinement_factor=nest_ratio,&
          rounding_method='nearest')
      md%domains(2)%timestepping_method = default_nonlinear_timestepping_method !'cliffs' !'midpoint' !'rk2'
      !md%domains(2)%theta = 4.0_dp ! No big impact for this problem.
  
      print*, 2, ' lw: ', md%domains(2)%lw, ' ll: ', md%domains(2)%lower_left, ' dx: ', md%domains(2)%dx, &
          ' nx: ', md%domains(2)%nx
  
negative.code: 
      !
      ! Parameter controlling the northern boundary treatment, which substantially
      ! affects the tsunami in this case. Possible values are:
      !   'boundary_stage_radiation_momentum' -- Works pretty well without tuning for rk2 and related. 
      !   'flather_with_vh_from_continuity' -- about right? But depends on a scale for the boundary 
      !       VH term, which we can estimate from continuity, but needs tuning. OK for leapfrog_nonlinear.
  