anchor.diff:  
       md%domains(1)%lower_left = global_ll
       md%domains(1)%nx = global_nx
  -    md%domains(1)%lw = global_lw * ( ( 1.0_dp * md%domains(1)%nx ) / (1.0_dp * global_nx) )
  +    md%domains(1)%lw = global_lw
       md%domains(1)%dx = md%domains(1)%lw/md%domains(1)%nx
       md%domains(1)%timestepping_method = default_nonlinear_timestepping_method
  -    md%domains(1)%static_before_time = 1920.0_dp
  +    md%domains(1)%static_before_time = 1920.0_dp ! No need to evolve prior to this time
   
       print*, 1, ' lw: ', md%domains(1)%lw, ' ll: ', md%domains(1)%lower_left, &
           ' dx: ', md%domains(1)%dx, ' nx: ', md%domains(1)%nx
  
positive.code:     ! Build boundary conditions
      call setup_boundary_information(bc_file, bc_elev)
  
      ! Setup hazard points
      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')
  
negative.code:     type :: boundary_information_type
          character(charlen):: bc_file
          real(dp), allocatable:: boundary_data(:,:)
          type(linear_interpolator_type):: gauge_ts_function
          real(dp):: boundary_elev
          real(dp):: t0 = 0.0_dp
      end type
  