anchor.diff:              [nc_grid_output%dim_x_id, nc_grid_output%dim_y_id], &
               nc_grid_output%var_priority_img), __LINE__ )
   
  -        ! Peak stage
  -        if(record_max_U) then
  -
  -            nc_grid_output%record_max_U = .TRUE.
  -
  -            call check( nf90_def_var(iNcid, "max_stage", output_prec, &
  -                [nc_grid_output%dim_x_id, nc_grid_output%dim_y_id], &
  -                nc_grid_output%var_max_stage_id), __LINE__ )
  -
  -            call check( nf90_def_var(iNcid, "elevation0", output_prec, &
  -                [nc_grid_output%dim_x_id, nc_grid_output%dim_y_id], &
  -                nc_grid_output%var_elev0_id), __LINE__ )
  -
  -            call check( nf90_def_var(iNcid, "manning_squared", output_prec, &
  -                [nc_grid_output%dim_x_id, nc_grid_output%dim_y_id], &
  -                nc_grid_output%var_manningsq_id), __LINE__ )
  -        else
  -
  -            nc_grid_output%record_max_U = .FALSE.
  -
  -        end if
  +
  +        ! Use nontemporal_grids_to_store to define which variables to store.
  +        !
  +        ! Check for typos. It can contain either an empty string '', or
  +        ! a value in nc_grid_output%static_var_names
  +        do i = 1, size(nontemporal_grids_to_store)
  +            if(.not. (nontemporal_grids_to_store(i) == '' .or. &
  +                  any(nontemporal_grids_to_store(i) == nc_grid_output%static_var_names))) then
  +                write(log_output_unit, *)&
  +                    "Unknown value in nontemporal_grids_to_store:", trim(time_grids_to_store(i))
  +                call generic_stop
  +            end if
  +        end do
  +        ! Set the value of static_var_store_flag so that the variables that we want to store are stored.
  +        do i = 1, size(nc_grid_output%static_var_names)
  +            nc_grid_output%static_var_store_flag(i) = &
  +                any(nc_grid_output%static_var_names(i) == nontemporal_grids_to_store)
  +        end do
  +
  +        ! Storage for the static variables
  +        do i = 1, size(nc_grid_output%static_var_names)
  +
  +            if(.not. nc_grid_output%static_var_store_flag(i)) cycle
  +
  +            call check( &
  +                nf90_def_var(iNcid, &
  +                    trim(nc_grid_output%static_var_names(i)), &
  +                    output_prec, &
  +                    [nc_grid_output%dim_x_id, nc_grid_output%dim_y_id], &
  +                    nc_grid_output%static_var_id(i)), &
  +                __LINE__ )
  +        end do
   
           !
           ! Add some attributes
  
positive.code:         call check(nf90_put_var(iNcid, nc_grid_output%var_time_id, time, &
              start=[nc_grid_output%num_output_steps]), __LINE__ )
  
          ! Save the stage
          if(nc_grid_output%time_var_store_flag(STG)) then
              call check(nf90_put_var(&
                  iNcid, &
                  nc_grid_output%var_stage_id, &
                  U(spatial_start(1):nxy(1):spatial_stride(1),spatial_start(2):nxy(2):spatial_stride(2),STG), &
                  start=[1,1,nc_grid_output%num_output_steps]),&
                  __LINE__)
          end if
  
          ! Save uh
          if(nc_grid_output%time_var_store_flag(UH)) then
              call check(nf90_put_var(&
                  iNcid, &
                  nc_grid_output%var_uh_id, &
                  U(spatial_start(1):nxy(1):spatial_stride(1),spatial_start(2):nxy(2):spatial_stride(2),UH), &
                  start=[1,1,nc_grid_output%num_output_steps]), &
                  __LINE__)
          end if
  
          ! Save vh
          if(nc_grid_output%time_var_store_flag(VH)) then
              call check(nf90_put_var(&
                  iNcid, &
                  nc_grid_output%var_vh_id, &
                  U(spatial_start(1):nxy(1):spatial_stride(1),spatial_start(2):nxy(2):spatial_stride(2),VH), &
                  start=[1,1,nc_grid_output%num_output_steps]), &
                  __LINE__)
          end if
  
          ! Save elev
          if(nc_grid_output%time_var_store_flag(ELV)) then
              call check(nf90_put_var(&
                  iNcid, &
                  nc_grid_output%var_elev_id, &
                  U(spatial_start(1):nxy(1):spatial_stride(1),spatial_start(2):nxy(2):spatial_stride(2),ELV), &
                  start=[1,1,nc_grid_output%num_output_steps]), &
                  __LINE__)
          end if
  
  #ifdef DEBUG_ARRAY
          if(present(debug_array)) then
  
negative.code:             !! If compiled with -DDEBUG_ARRAY, this array will also be written to the netcdf file.
              !! Provides a useful means of debugging or reporting non-standard quantities.
  
          integer:: iNcid, nxy(2), spatial_start(2), spatial_stride(2), i
          real(C_DOUBLE) :: current_cpu_time
  
          ! Gracefully compile in case we don't link with netcdf
  