anchor.diff:          allocate(domain%energy_kinetic_on_rho_work(ny))
   
           !
  +        ! If we specified the time_grids_to_store, check for spelling errors
  +        !
  +        if(allocated(domain%time_grids_to_store)) then
  +            do i = 1, size(domain%time_grids_to_store)
  +                if(.not. ((domain%time_grids_to_store(i) == '') .or. &
  +                    any(domain%time_grids_to_store(i) == gridded_output_variables_time))) then
  +
  +                    write(log_output_unit, *) 'Unsupported variable in time_grids_to_store: ', &
  +                        trim(domain%time_grids_to_store(i))
  +                    call generic_stop
  +                end if
  +            end do
  +        end if 
  +        !
  +        ! Determine which flow statistics we store that are not time-series.
  +        !
  +        if(.not. allocated(domain%nontemporal_grids_to_store)) then
  +            ! Mimic 'old behaviour' before we allowed specification of domain%nontempral_grids_to_store.
  +            ! In new code please avoid setting domain%record_max_U -- instead set domain%nontemporal_grids_to_store
  +            if(domain%record_max_U) then
  +                ! Old default where we store maxima
  +                domain%nontemporal_grids_to_store = [character(len=charlen):: 'max_stage', 'elevation0', 'manning_squared'] 
  +            else
  +                ! Old default where we don't store maxima or elevation or manning
  +                domain%nontemporal_grids_to_store = [''] 
  +            end if
  +        end if
  +        ! Check for spelling errors in nontemporal_grids_to_store 
  +        do i = 1, size(domain%nontemporal_grids_to_store)
  +            if(.not. ((domain%nontemporal_grids_to_store(i) == '') .or. &
  +                any(domain%nontemporal_grids_to_store(i) == gridded_output_variables_max_U) .or. &
  +                any(domain%nontemporal_grids_to_store(i) == gridded_output_variables_other) )) then
  +                
  +                write(log_output_unit, *) 'Unsupported variable in nontemporal_grids_to_store: ', &
  +                    trim(domain%nontemporal_grids_to_store(i))
  +                call generic_stop
  +            end if 
  +        end do
  +        ! Define the names of nontemporal variables to store that require tracking flow statistics
  +        allocate(domain%max_U_variables(0))
  +        do i = 1, size(gridded_output_variables_max_U)
  +            if(any(domain%nontemporal_grids_to_store == gridded_output_variables_max_U(i))) &
  +                domain%max_U_variables = [domain%max_U_variables, gridded_output_variables_max_U(i)]
  +        end do
  +        if(size(domain%max_U_variables) > 0) then
  +            ! Ensure the first max_U variable is 'max_stage'
  +            if(domain%max_U_variables(1) /= 'max_stage') then
  +                write(log_output_unit, *) 'Error: "max_stage" MUST be added to domain%nontemporal_grids_to_store.'
  +                write(log_output_unit, *) 'While not logically essential, some older code assumes that if'
  +                write(log_output_unit, *) 'domain%max_U is allocated, then domain%max_U(:,:,1) has max-stage'
  +                call generic_stop
  +            end if
  +        end if
  +
  +        !
           ! Below we allocate the main arrays, using loops to promote openmp memory affinity (based on the first-touch principle).
           !
           allocate(domain%U(nx, ny, 1:nvar))
  
positive.code: 
      subroutine compute_depth_and_velocity(domain, min_depth)
          !! Compute domain%depth and domain%velocity, so they are consistent with domain%U.
          !! We also implement a mass conservation check in this routine, and zero velocities
          !! at sites with depth < minimum_allowed_depth.
  
  
negative.code:             !! 'nesting_boundary_flux_*')
  
          !
          ! Output file content -- bespoke binary format (that is no longer properly supported, consider removing)
          !
          character(len=charlen):: metadata_ascii_filename
              !! Name of ascii file where we output metadata
  