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:         boundary_type = 'boundary_stage_radiation_momentum'
      end if
          
  
      ! Initial conditions
      do j = 1, size(md%domains)
          call set_initial_conditions(md%domains(j))
          ! Even for linear solver, allow the (g depth dstage/dx) term to have depth varying
          md%domains(j)%linear_solver_is_truely_linear = .false.
      end do
      call md%make_initial_conditions_consistent 
      ! NOTE: For stability in 'null' regions, we set them to 'high land' that
      ! should be inactive. 
      call md%set_null_regions_to_dry()
  
  
negative.code:     md%domains(1)%timestepping_refinement_factor = 1_ip
      md%domains(1)%dx_refinement_factor = 1.0_dp
      md%domains(1)%timestepping_method = default_nonlinear_timestepping_method
      !md%domains(1)%use_eddy_viscosity = .true.
  
      print*, 1, ' lw: ', md%domains(1)%lw, ' ll: ', md%domains(1)%lower_left, ' dx: ', md%domains(1)%dx, &
          ' nx: ', md%domains(1)%nx
  
      ! A detailed domain
      call md%domains(2)%match_geometry_to_parent(&
          parent_domain=md%domains(1), &
          lower_left=[28500.0_dp, 11500.0_dp], &
  