anchor.diff:              TIMER_STOP('nesting_flux_correction')
   
               !DEBUG
  -            !call md%check_for_overflow('step-after-flux_correct')
  +            !call check_for_overflow(md, 'step-after-flux_correct')
           end if
   
       end subroutine
  
positive.code: 
      ! Scan domains for unreasonably large numbers or NaN issues
      subroutine check_for_overflow(md, flag, domain_ind)
          class(multidomain_type), intent(inout) :: md
          character(len=*), intent(in) :: flag
          integer(ip), optional, intent(in):: domain_ind
  
          logical :: throw_error
          integer(ip) :: i1, i2, i3, j, d1, d2
          real(dp), parameter :: bignum = HUGE(1.0)/2.0 ! Deliberate single-precision number
  
  
negative.code:             d2 = size(md%domains)
          end if
  
          throw_error = 0
          !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(md, d1, d2) REDUCTION(+:throw_error)
          throw_error = 0
          do j = d1, d2
              !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
              do i1 = STG, ELV
                  do i2 = 1, md%domains(j)%nx(2)
                      do i3 = 1, md%domains(j)%nx(1)
                          if((md%domains(j)%U(i3, i2, i1) >  bignum) .or. &
                             (md%domains(j)%U(i3, i2, i1) < -bignum) .or. &
                             (md%domains(j)%U(i3, i2, i1) /= md%domains(j)%U(i3, i2, i1))) then
                             throw_error = throw_error + 1
                          end if
                      end do
                  end do
              end do
              !$OMP END DO
          end do
          !$OMP END PARALLEL
  
          if(throw_error > 0) then
              do j = d1, d2
                  do i1 = STG, ELV
                      do i2 = 1, md%domains(j)%nx(2)
                          do i3 = 1, md%domains(j)%nx(1)
                              if((md%domains(j)%U(i3, i2, i1) >  bignum) .or. &
                                 (md%domains(j)%U(i3, i2, i1) < -bignum) .or. &
                                 (md%domains(j)%U(i3, i2, i1) /= md%domains(j)%U(i3, i2, i1))) then
                                 write(log_output_unit,*) flag
                                 write(log_output_unit,*) 'error-count: ', throw_error
                                 write(log_output_unit,*) 'md%domains(', j,')%U(', i3, i2, i1, ')'
                                 write(log_output_unit,*) md%domains(j)%U(i3, i2, i1)
                                 write(log_output_unit,*) 'x: ', md%domains(j)%x(i3), '; y: ', md%domains(j)%y(i2)
                                 write(log_output_unit,*) 'time: ', md%domains(j)%time
                                 write(log_output_unit,*) 'stg:elv, ', md%domains(j)%U(i3, i2, STG:ELV)
                                 if(md%domains(j)%nesting%my_index > 0) then
                                     write(log_output_unit,*) 'priority domain index'
                                     write(log_output_unit,*) md%domains(j)%nesting%priority_domain_index(i3, i2)
                                     write(log_output_unit,*) 'priority domain image'
                                     write(log_output_unit,*) md%domains(j)%nesting%priority_domain_image(i3, i2)
                                 end if
                              end if
                          end do
                      end do
                  end do
              end do
  
              flush(log_output_unit)
              call md%finalise_and_print_timers()
              call generic_stop
          end if
  
      end subroutine
  
      ! Communicate max-stage array in halo regions.
      !
  