anchor.diff:          ! there, leading to tiny numerical differences in the time.
           ! To avoid any issues, we explicitly set them to the same number
   
  -        !DEBUG
  -        !call check_for_overflow(md, 'start')
  +#ifdef TRACK_MULTIDOMAIN_STABILITY
  +        call check_multidomain_stability(md, 'start')
  +#endif
   
           ! Evolve every domain, sending the data to communicate straight after
           ! the evolve
  
positive.code:         !$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
  
negative.code: 
  
      ! Scan domains for unreasonably large numbers or NaN issues
      subroutine check_multidomain_stability(md, flag, domain_ind)
          type(multidomain_type), intent(inout) :: md
          character(len=*), intent(in) :: flag
          integer(ip), optional, intent(in):: domain_ind
  
          integer :: throw_error
          integer(ip) :: i1, i2, i3, j, d1, d2, ecw
          real(dp), parameter :: bignum = HUGE(1.0)/2.0 ! Deliberate single-precision number
  
          if(present(domain_ind)) then
  