anchor.diff:          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.
  +    md%domains(2)%timestepping_method = default_nonlinear_timestepping_method
  +    !md%domains(2)%use_eddy_viscosity = .true.
   
       print*, 2, ' lw: ', md%domains(2)%lw, ' ll: ', md%domains(2)%lower_left, ' dx: ', md%domains(2)%dx, &
           ' nx: ', md%domains(2)%nx
  
positive.code:     ! nd domains in this model
      nd = 2
      allocate(md%domains(nd))
      !md%load_balance_file = 'load_balance_partition.txt'
  
      !
      ! Setup basic metadata
  
negative.code:         boundary_information%boundary_data(nr - skip + extra,2) = &
              boundary_information%boundary_data(nr - skip + extra - 1, 2)
  
          call boundary_information%gauge4_ts_function%initialise(&
              boundary_information%boundary_data(:,1), boundary_information%boundary_data(:,2))
          boundary_information%t0 = boundary_information%boundary_data(1,1)
  
      end subroutine
     
      ! 
      ! 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
  