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:     ! A more detailed Monai domain 
      call md%domains(5)%match_geometry_to_parent(&
          parent_domain=md%domains(4), &
          !lower_left  = [139.4208_dp, 42.0968_dp], &
          lower_left  = [139.4150_dp, 42.0968_dp], &
          upper_right = [139.4264_dp, 42.1020_dp], &
          dx_refinement_factor = nest_ratio, &
  
negative.code:     
      ! nd domains in this model
      if(very_high_res_monai) then
          ! In this case we put a 30cm x 30cm cell domain around the monai inundation peak. (with mesh_refine=1.0_dp)
          ! It leads to a peak of > 30.0m (vs 31.7 obs, although the latter varies depending on which
          ! dataset is used).
          !
  