anchor.diff:  
       !
       ! Store a 1 byte integer, denoting the cells in the current domain that are 
  -    ! priority_domain cells.
  +    ! priority_domain cells. Also store regular integer grids with the priority domain index
  +    ! and image.
       !
       subroutine store_priority_domain_cells(nc_grid_output, priority_domain_index, &
           priority_domain_image, my_index, my_image)
  
positive.code:             (/ nc_grid_output%dim_time_id /), nc_grid_output%var_time_id), __LINE__ )
  
          ! Stage
          call check( nf90_def_var(iNcid, "stage", output_prec, &
              [nc_grid_output%dim_x_id, nc_grid_output%dim_y_id, nc_grid_output%dim_time_id], &
              nc_grid_output%var_stage_id), __LINE__ )
          ! Compress
          if(using_netcdf4) call check(nf90_def_var_deflate(iNcid, nc_grid_output%var_stage_id, 1, 1, 3), __LINE__)
          ! UH
          call check( nf90_def_var(iNcid, "uh", output_prec, &
              [nc_grid_output%dim_x_id, nc_grid_output%dim_y_id, nc_grid_output%dim_time_id], &
              nc_grid_output%var_uh_id), __LINE__ )
          ! Compress
          if(using_netcdf4) call check(nf90_def_var_deflate(iNcid, nc_grid_output%var_uh_id, 1, 1, 3), __LINE__)
  
          ! VH
          call check( nf90_def_var(iNcid, "vh", output_prec, &
              [nc_grid_output%dim_x_id, nc_grid_output%dim_y_id, nc_grid_output%dim_time_id], &
              nc_grid_output%var_vh_id), __LINE__ )
          ! Compress
          if(using_netcdf4) call check(nf90_def_var_deflate(iNcid, nc_grid_output%var_vh_id, 1, 1, 3), __LINE__)
  
          ! Elev
          call check( nf90_def_var(iNcid, "elev", output_prec, &
              [nc_grid_output%dim_x_id, nc_grid_output%dim_y_id, nc_grid_output%dim_time_id], &
              nc_grid_output%var_elev_id), __LINE__ )
          ! Compress
          if(using_netcdf4) call check(nf90_def_var_deflate(iNcid, nc_grid_output%var_elev_id, 1, 1, 3), __LINE__)
  
          ! priority_domain -- Use 1 if the current domain is the priority domain, and zero otherwise
          call check(nf90_def_var(iNcid, "is_priority_domain", output_byte, &
  
negative.code: 
          ! Flag for whether we store max stage, etc
          logical :: record_max_U
  
          ! Options to store stage/uh/vh/elev over time.
          ! .true. or .false. for STG, UH, VH, ELV
          logical :: time_var_store_flag(4) = [.true. , .true., .true., .true.]
          ! For example, if we didn't want to store ELV every time-step, we
          ! would do "nc_grid_output%time_var_store_flag(ELV) = .false."
  
          ! Allow only storing every n'th point in x/y space.
          ! For instance, on a 1-arc-minute model, a value of 4 would
  