anchor.diff:          allocate(random_uniform(domain%nx(1)))
           
           ! Set stage and elevation row-by-row.
  -        ! This saves memory compared to doing it all at once.
           do j = 1, domain%nx(2)
               y = domain%y(j)
   
  -            ! Set stage
  -            call stage_data%get_xy(x,y, domain%U(:,j,STG), domain%nx(1), bilinear=1_ip)
  -            ! Clip 'NA' regions (since the stage raster does no cover the entire domain)
  -            do i = 1, domain%nx(1)
  -                if(domain%U(i,j,STG) < (-HUGE(1.0_dp)*0.99_dp) ) domain%U(i,j,STG) = 0.0_dp 
  -            end do
  -            ! Set elevation
  +            ! Set stage - and clip 'NA' regions (since the stage raster does no cover the entire domain)
  +            call stage_data%get_xy(x,y, domain%U(:,j,STG), domain%nx(1), bilinear=1_ip, na_below_limit=-1.0e-20_dp)
  +            where(domain%U(:,j,STG) <= -1.0e-20_dp ) domain%U(:,j,STG) = 0.0_dp 
  +
  +            ! Set elevation -- no need to clip NA.
               call elevation_data%get_xy(x,y, domain%U(:,j,ELV), domain%nx(1), bilinear=1_ip)
   
  -            !! Test only
  -            !domain%U(:,j,STG) = 1.0e+02_dp*(x + y) + 1.0e-03_dp
  -            !domain%U(:,j,ELV) = 1.0e+02_dp*(x + y)
  +            ! Here we can experiment with the random perturbation
               call random_number(random_uniform)
               domain%U(:,j,STG) = domain%U(:,j,STG) + random_perturbation_scale * (random_uniform-0.5_dp)
               call random_number(random_uniform)
  
positive.code:     md%domains(3)%timestepping_method = 'rk2'
  
      ! The monai domain 
      ! (Peak stage seems to be affected by the westward extent
      ! of this, possibly mesh design needs more tuning)
      call md%domains(4)%match_geometry_to_parent(&
          parent_domain=md%domains(3), &
          lower_left  = [139.38_dp, 42.0828_dp], &
  
negative.code:             dx_refinement_factor = nest_ratio, &
              timestepping_refinement_factor = 6_ip,&
              rounding_method='nearest')
          ! Older versions of the code needed euler for stability. Does it still?
          md%domains(7)%timestepping_method = 'euler'  
          ! Should be dry before this time -- note also we change the time-step 
          ! in the evolve loop on the basis of this value
  