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:     end subroutine
  
      subroutine initialise(nc_grid_output, filename, output_precision, &
          xs, ys, &
          time_grids_to_store, nontemporal_grids_to_store, &
          attribute_names, attribute_values)
          !!
          !! Create nc file to hold gridded outputs
          !!
  
          class(nc_grid_output_type), intent(inout) :: nc_grid_output 
              !!type with nc_grid_output info
          character(len=charlen), intent(in) :: filename 
              !! Name of the output netcdf file
          integer(ip), intent(in) :: output_precision 
              !! Fortran real kind denoting precision of output (either C_FLOAT or C_DOUBLE)
          real(dp), intent(in) :: xs(:), ys(:) 
              !! array of grid x/y coordinates
          character(charlen), allocatable, intent(inout) :: time_grids_to_store(:) 
              !! Array with names of time-variables to store. Values should match those in nc_grid_output%time_var_names, or ''.
              !! To support some legacy behaviour we allow this to be unallocated on entry. Some older code specified
              !! variables to store by directly manipulating the nc_grid_output variable, and that approach is still supported
              !! if time_grids_to_store is unallocated (although it is not recommended -- in future code please just specify
              !! the variables to store in domain%time_grids_to_store).
          character(charlen), intent(inout) :: nontemporal_grids_to_store(:) 
              !! Array with names of variables without a time dimension. Note that irrespective of the value of this variable, 
              !! we also store other "priority-domain" variables. Values should match those in nc_grid_output%static_var_names, or ''.
          character(charlen), optional, intent(in):: attribute_names(:) 
              !! global attribute names for netcdf file
          character(charlen), optional, intent(in):: attribute_values(:)
              !! values for global attributes of netcdf file (same length as attribute_names)
  
  