MOM_debugging.F90

Go to the documentation of this file.

module mom_debugging

!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!   This module contains subroutines that perform various error checking and   !
! debugging functions for MOM6.  This routine is similar to it counterpart in  !
! the SIS2 code, except for the use of the ocean_grid_type and by keeping them !
! separate we retain the ability to set up MOM6 and SIS2 debugging separately. !
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!

use mom_checksums, only : hchksum, bchksum, qchksum, uvchksum
use mom_checksums, only : is_nan, chksum, mom_checksums_init
use mom_coms, only : pe_here, root_pe, num_pes, sum_across_pes
use mom_coms, only : min_across_pes, max_across_pes, reproducing_sum
use mom_domains, only : pass_vector, pass_var, pe_here
use mom_domains, only : bgrid_ne, agrid, to_all, scalar_pair
use mom_error_handler, only : mom_error, fatal, warning, is_root_pe
use mom_file_parser, only : log_version, param_file_type, get_param
use mom_grid, only : ocean_grid_type
use mom_hor_index, only : hor_index_type

implicit none ; private

public :: check_redundant_c, check_redundant_b, check_redundant_t, check_redundant
public :: vec_chksum, vec_chksum_c, vec_chksum_b, vec_chksum_a
public :: mom_debugging_init, totalstuff, totaltands
public :: check_column_integral, check_column_integrals

! These interfaces come from MOM_checksums.
public :: hchksum, bchksum, qchksum, is_nan, chksum
public :: uvchksum

interface check_redundant
  module procedure check_redundant_vc3d, check_redundant_vc2d
end interface check_redundant
interface check_redundant_c
  module procedure check_redundant_vc3d, check_redundant_vc2d
end interface check_redundant_c
interface check_redundant_b
  module procedure check_redundant_vb3d, check_redundant_vb2d
  module procedure check_redundant_sb3d, check_redundant_sb2d
end interface check_redundant_b
interface check_redundant_t
  module procedure check_redundant_st3d, check_redundant_st2d
  module procedure check_redundant_vt3d, check_redundant_vt2d
end interface check_redundant_t

interface vec_chksum
  module procedure chksum_vec_c3d, chksum_vec_c2d
end interface vec_chksum
interface vec_chksum_c
  module procedure chksum_vec_c3d, chksum_vec_c2d
end interface vec_chksum_c
interface vec_chksum_b
  module procedure chksum_vec_b3d, chksum_vec_b2d
end interface vec_chksum_b
interface vec_chksum_a
  module procedure chksum_vec_a3d, chksum_vec_a2d
end interface vec_chksum_a

integer :: max_redundant_prints = 100
integer :: redundant_prints(3) = 0
logical :: debug = .false.
logical :: debug_chksums = .true.
logical :: debug_redundant = .true.

contains

! =====================================================================

!> MOM_debugging_init initializes the MOM_debugging module, and sets
!! the parameterts that control which checks are active for MOM6.
subroutine mom_debugging_init(param_file)
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_debugging" ! This module's name.

  call log_version(param_file, mdl, version)
  call get_param(param_file, mdl, "DEBUG", debug, &
                 "If true, write out verbose debugging data.", default=.false.)
  call get_param(param_file, mdl, "DEBUG_CHKSUMS", debug_chksums, &
                 "If true, checksums are performed on arrays in the \n"//&
                 "various vec_chksum routines.", default=debug)
  call get_param(param_file, mdl, "DEBUG_REDUNDANT", debug_redundant, &
                 "If true, debug redundant data points during calls to \n"//&
                 "the various vec_chksum routines.", default=debug)

  call mom_checksums_init(param_file)

end subroutine mom_debugging_init

subroutine check_redundant_vc3d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                                direction)
  character(len=*),                    intent(in)    :: mesg
  type(ocean_grid_type),               intent(inout) :: G    !< The ocean's grid structure
  real, dimension(G%IsdB:,G%jsd:,:),   intent(in)    :: u_comp
  real, dimension(G%isd:,G%JsdB:,:),   intent(in)    :: v_comp
  integer,                   optional, intent(in)    :: is, ie, js, je
  integer,                   optional, intent(in)    :: direction
! Arguments: u_comp - The u-component of the vector being checked.
!  (in)      v_comp - The v-component of the vector being checked.
!  (in)      mesg - A message indicating what is being checked.
!  (in)      G - The ocean's grid structure.
!  (in/opt)  is, ie, js, je - the i- and j- range of indices to check.
!  (in/opt)  direction - the direction flag to be passed to pass_vector.

  character(len=24) :: mesg_k
  integer :: k

  do k=1,size(u_comp,3)
    if (k < 10) then ; write(mesg_k,'(" Layer",i2," ")') k
    elseif (k < 100) then ; write(mesg_k,'(" Layer",i3," ")') k
    elseif (k < 1000) then ; write(mesg_k,'(" Layer",i4," ")') k
    else ; write(mesg_k,'(" Layer",i9," ")') k ; endif

    call check_redundant_vc2d(trim(mesg)//trim(mesg_k), u_comp(:,:,k), &
             v_comp(:,:,k), g, is, ie, js, je, direction)
  enddo
end subroutine  check_redundant_vc3d

subroutine check_redundant_vc2d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                                direction)
  character(len=*),                intent(in)    :: mesg
  type(ocean_grid_type),           intent(inout) :: G    !< The ocean's grid structure
  real, dimension(G%IsdB:,G%jsd:), intent(in)    :: u_comp
  real, dimension(G%isd:,G%JsdB:), intent(in)    :: v_comp
  integer,               optional, intent(in)    :: is, ie, js, je
  integer,               optional, intent(in)    :: direction
! Arguments: u_comp - The u-component of the vector being checked.
!  (in)      v_comp - The v-component of the vector being checked.
!  (in)      mesg - A message indicating what is being checked.
!  (in)      G - The ocean's grid structure.
!  (in/opt)  is, ie, js, je - the i- and j- range of indices to check.
!  (in/opt)  direction - the direction flag to be passed to pass_vector.

  real :: u_nonsym(g%isd:g%ied,g%jsd:g%jed)
  real :: v_nonsym(g%isd:g%ied,g%jsd:g%jed)
  real :: u_resym(g%isdb:g%iedb,g%jsd:g%jed)
  real :: v_resym(g%isd:g%ied,g%jsdb:g%jedb)
  character(len=128) :: mesg2

  integer :: i, j, is_ch, ie_ch, js_ch, je_ch
  integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (.not.(present(is) .or. present(ie) .or. present(js) .or. present(je))) then
    ! This only works with symmetric memory, so otherwise return.
    if ((isd == isdb) .and. (jsd == jsdb)) return
  endif

  do i=isd,ied ; do j=jsd,jed
    u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j)
  enddo ; enddo

  if (.not.associated(g%Domain_aux)) call mom_error(fatal," check_redundant"//&
    " called with a non-associated auxiliary domain the grid type.")
  call pass_vector(u_nonsym, v_nonsym, g%Domain_aux, direction)

  do i=isdb,iedb ; do j=jsd,jed ; u_resym(i,j) = u_comp(i,j) ; enddo ; enddo
  do i=isd,ied ; do j=jsdb,jedb ; v_resym(i,j) = v_comp(i,j) ; enddo ; enddo
  do i=isd,ied ; do j=jsd,jed
    u_resym(i,j) = u_nonsym(i,j) ; v_resym(i,j) = v_nonsym(i,j)
  enddo ; enddo
  call pass_vector(u_resym, v_resym, g%Domain, direction)

  is_ch = isq ; ie_ch = ieq ; js_ch = jsq ; je_ch = jeq
  if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie
  if (present(js)) js_ch = js ; if (present(js)) je_ch = je

  do i=is_ch,ie_ch ; do j=js_ch+1,je_ch
    if (u_resym(i,j) /= u_comp(i,j) .and. &
        redundant_prints(3) < max_redundant_prints) then
      write(mesg2,'(" redundant u-components",2(1pe12.4)," differ by ", &
                    & 1pe12.4," at i,j = ",2i4," on pe ",i4)') &
           u_comp(i,j), u_resym(i,j),u_comp(i,j)-u_resym(i,j),i,j,pe_here()
      write(0,'(A130)') trim(mesg)//trim(mesg2)
      redundant_prints(3) = redundant_prints(3) + 1
    endif
  enddo ; enddo
  do i=is_ch+1,ie_ch ; do j=js_ch,je_ch
    if (v_resym(i,j) /= v_comp(i,j) .and. &
        redundant_prints(3) < max_redundant_prints) then
      write(mesg2,'(" redundant v-comps",2(1pe12.4)," differ by ", &
                    & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)" on pe ",i4)') &
           v_comp(i,j), v_resym(i,j),v_comp(i,j)-v_resym(i,j),i,j, &
           g%geoLonBu(i,j), g%geoLatBu(i,j), pe_here()
      write(0,'(A155)') trim(mesg)//trim(mesg2)
      redundant_prints(3) = redundant_prints(3) + 1
    endif
  enddo ; enddo

end subroutine  check_redundant_vc2d

subroutine check_redundant_sb3d(mesg, array, G, is, ie, js, je)
  character(len=*),                     intent(in)    :: mesg
  type(ocean_grid_type),                intent(inout) :: G    !< The ocean's grid structure
  real, dimension(G%IsdB:,G%JsdB:,:),   intent(in)    :: array
  integer,                    optional, intent(in)    :: is, ie, js, je
! Arguments: array - The array being checked.
!  (in)      mesg - A message indicating what is being checked.
!  (in)      G - The ocean's grid structure.
!  (in/opt)  is, ie, js, je - the i- and j- range of indices to check.

  character(len=24) :: mesg_k
  integer :: k

  do k=1,size(array,3)
    if (k < 10) then ; write(mesg_k,'(" Layer",i2," ")') k
    elseif (k < 100) then ; write(mesg_k,'(" Layer",i3," ")') k
    elseif (k < 1000) then ; write(mesg_k,'(" Layer",i4," ")') k
    else ; write(mesg_k,'(" Layer",i9," ")') k ; endif

    call check_redundant_sb2d(trim(mesg)//trim(mesg_k), array(:,:,k), &
                              g, is, ie, js, je)
  enddo
end subroutine  check_redundant_sb3d


subroutine check_redundant_sb2d(mesg, array, G, is, ie, js, je)
  character(len=*),                 intent(in)    :: mesg
  type(ocean_grid_type),            intent(inout) :: G    !< The ocean's grid structure
  real, dimension(G%IsdB:,G%JsdB:), intent(in)    :: array
  integer,                optional, intent(in)    :: is, ie, js, je
! Arguments: array - The array being checked.
!  (in)      mesg - A message indicating what is being checked.
!  (in)      G - The ocean's grid structure.
!  (in/opt)  is, ie, js, je - the i- and j- range of indices to check.

  real :: a_nonsym(g%isd:g%ied,g%jsd:g%jed)
  real :: a_resym(g%isdb:g%iedb,g%jsdb:g%jedb)
  character(len=128) :: mesg2

  integer :: i, j, is_ch, ie_ch, js_ch, je_ch
  integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (.not.(present(is) .or. present(ie) .or. present(js) .or. present(je))) then
    ! This only works with symmetric memory, so otherwise return.
    if ((isd == isdb) .and. (jsd == jsdb)) return
  endif

  do i=isd,ied ; do j=jsd,jed
    a_nonsym(i,j) = array(i,j)
  enddo ; enddo

  if (.not.associated(g%Domain_aux)) call mom_error(fatal," check_redundant"//&
    " called with a non-associated auxiliary domain the grid type.")
  call pass_vector(a_nonsym, a_nonsym, g%Domain_aux, &
                   direction=to_all+scalar_pair, stagger=bgrid_ne)

  do i=isdb,iedb ; do j=jsdb,jedb ; a_resym(i,j) = array(i,j) ; enddo ; enddo
  do i=isd,ied ; do j=jsd,jed
    a_resym(i,j) = a_nonsym(i,j)
  enddo ; enddo
  call pass_vector(a_resym, a_resym, g%Domain, direction=to_all+scalar_pair, &
                   stagger=bgrid_ne)

  is_ch = isq ; ie_ch = ieq ; js_ch = jsq ; je_ch = jeq
  if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie
  if (present(js)) js_ch = js ; if (present(js)) je_ch = je

  do i=is_ch,ie_ch ; do j=js_ch,je_ch
    if (a_resym(i,j) /= array(i,j) .and. &
        redundant_prints(2) < max_redundant_prints) then
      write(mesg2,'(" Redundant points",2(1pe12.4)," differ by ", &
                    & 1pe12.4," at i,j = ",2i4," on pe ",i4)') &
           array(i,j), a_resym(i,j),array(i,j)-a_resym(i,j),i,j,pe_here()
      write(0,'(A130)') trim(mesg)//trim(mesg2)
      redundant_prints(2) = redundant_prints(2) + 1
    endif
  enddo ; enddo

end subroutine  check_redundant_sb2d


subroutine check_redundant_vb3d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                                direction)
  character(len=*),                    intent(in)    :: mesg
  type(ocean_grid_type),               intent(inout) :: G    !< The ocean's grid structure
  real, dimension(G%IsdB:,G%JsdB:,:),  intent(in)    :: u_comp
  real, dimension(G%IsdB:,G%JsdB:,:),  intent(in)    :: v_comp
  integer,                   optional, intent(in)    :: is, ie, js, je
  integer,                   optional, intent(in)    :: direction
! Arguments: u_comp - The u-component of the vector being checked.
!  (in)      v_comp - The v-component of the vector being checked.
!  (in)      mesg - A message indicating what is being checked.
!  (in)      G - The ocean's grid structure.
!  (in/opt)  is, ie, js, je - the i- and j- range of indices to check.
!  (in/opt)  direction - the direction flag to be passed to pass_vector.

  character(len=24) :: mesg_k
  integer :: k

  do k=1,size(u_comp,3)
    if (k < 10) then ; write(mesg_k,'(" Layer",i2," ")') k
    elseif (k < 100) then ; write(mesg_k,'(" Layer",i3," ")') k
    elseif (k < 1000) then ; write(mesg_k,'(" Layer",i4," ")') k
    else ; write(mesg_k,'(" Layer",i9," ")') k ; endif

    call check_redundant_vb2d(trim(mesg)//trim(mesg_k), u_comp(:,:,k), &
             v_comp(:,:,k), g, is, ie, js, je, direction)
  enddo
end subroutine  check_redundant_vb3d

subroutine check_redundant_vb2d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                                direction)
  character(len=*),                intent(in)    :: mesg
  type(ocean_grid_type),            intent(inout) :: G    !< The ocean's grid structure
  real, dimension(G%IsdB:,G%JsdB:), intent(in)   :: u_comp
  real, dimension(G%IsdB:,G%JsdB:), intent(in)   :: v_comp
  integer,               optional, intent(in)    :: is, ie, js, je
  integer,               optional, intent(in)    :: direction
! Arguments: u_comp - The u-component of the vector being checked.
!  (in)      v_comp - The v-component of the vector being checked.
!  (in)      mesg - A message indicating what is being checked.
!  (in)      G - The ocean's grid structure.
!  (in/opt)  is, ie, js, je - the i- and j- range of indices to check.
!  (in/opt)  direction - the direction flag to be passed to pass_vector.

  real :: u_nonsym(g%isd:g%ied,g%jsd:g%jed)
  real :: v_nonsym(g%isd:g%ied,g%jsd:g%jed)
  real :: u_resym(g%isdb:g%iedb,g%jsdb:g%jedb)
  real :: v_resym(g%isdb:g%iedb,g%jsdb:g%jedb)
  character(len=128) :: mesg2

  integer :: i, j, is_ch, ie_ch, js_ch, je_ch
  integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (.not.(present(is) .or. present(ie) .or. present(js) .or. present(je))) then
    ! This only works with symmetric memory, so otherwise return.
    if ((isd == isdb) .and. (jsd == jsdb)) return
  endif

  do i=isd,ied ; do j=jsd,jed
    u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j)
  enddo ; enddo

  if (.not.associated(g%Domain_aux)) call mom_error(fatal," check_redundant"//&
    " called with a non-associated auxiliary domain the grid type.")
  call pass_vector(u_nonsym, v_nonsym, g%Domain_aux, direction, stagger=bgrid_ne)

  do i=isdb,iedb ; do j=jsdb,jedb
    u_resym(i,j) = u_comp(i,j) ; v_resym(i,j) = v_comp(i,j)
  enddo ; enddo
  do i=isd,ied ; do j=jsd,jed
    u_resym(i,j) = u_nonsym(i,j) ; v_resym(i,j) = v_nonsym(i,j)
  enddo ; enddo
  call pass_vector(u_resym, v_resym, g%Domain, direction, stagger=bgrid_ne)

  is_ch = isq ; ie_ch = ieq ; js_ch = jsq ; je_ch = jeq
  if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie
  if (present(js)) js_ch = js ; if (present(js)) je_ch = je

  do i=is_ch,ie_ch ; do j=js_ch,je_ch
    if (u_resym(i,j) /= u_comp(i,j) .and. &
        redundant_prints(2) < max_redundant_prints) then
      write(mesg2,'(" redundant u-components",2(1pe12.4)," differ by ", &
                    & 1pe12.4," at i,j = ",2i4," on pe ",i4)') &
           u_comp(i,j), u_resym(i,j),u_comp(i,j)-u_resym(i,j),i,j,pe_here()
      write(0,'(A130)') trim(mesg)//trim(mesg2)
      redundant_prints(2) = redundant_prints(2) + 1
    endif
  enddo ; enddo
  do i=is_ch,ie_ch ; do j=js_ch,je_ch
    if (v_resym(i,j) /= v_comp(i,j) .and. &
        redundant_prints(2) < max_redundant_prints) then
      write(mesg2,'(" redundant v-comps",2(1pe12.4)," differ by ", &
                    & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)" on pe ",i4)') &
           v_comp(i,j), v_resym(i,j),v_comp(i,j)-v_resym(i,j),i,j, &
           g%geoLonBu(i,j), g%geoLatBu(i,j), pe_here()
      write(0,'(A155)') trim(mesg)//trim(mesg2)
      redundant_prints(2) = redundant_prints(2) + 1
    endif
  enddo ; enddo

end subroutine  check_redundant_vb2d

subroutine check_redundant_st3d(mesg, array, G, is, ie, js, je)
  character(len=*),                     intent(in)    :: mesg
  type(ocean_grid_type),                intent(inout) :: G    !< The ocean's grid structure
  real, dimension(G%isd:,G%jsd:,:),     intent(in)    :: array
  integer,                    optional, intent(in)    :: is, ie, js, je
! Arguments: array - The array being checked.
!  (in)      mesg - A message indicating what is being checked.
!  (in)      G - The ocean's grid structure.
!  (in/opt)  is, ie, js, je - the i- and j- range of indices to check.

  character(len=24) :: mesg_k
  integer :: k

  do k=1,size(array,3)
    if (k < 10) then ; write(mesg_k,'(" Layer",i2," ")') k
    elseif (k < 100) then ; write(mesg_k,'(" Layer",i3," ")') k
    elseif (k < 1000) then ; write(mesg_k,'(" Layer",i4," ")') k
    else ; write(mesg_k,'(" Layer",i9," ")') k ; endif

    call check_redundant_st2d(trim(mesg)//trim(mesg_k), array(:,:,k), &
                              g, is, ie, js, je)
  enddo
end subroutine  check_redundant_st3d


subroutine check_redundant_st2d(mesg, array, G, is, ie, js, je)
  character(len=*),                 intent(in)    :: mesg
  type(ocean_grid_type),            intent(inout) :: G    !< The ocean's grid structure
  real, dimension(G%isd:,G%jsd:),   intent(in)    :: array
  integer,                optional, intent(in)    :: is, ie, js, je
! Arguments: array - The array being checked.
!  (in)      mesg - A message indicating what is being checked.
!  (in)      G - The ocean's grid structure.
!  (in/opt)  is, ie, js, je - the i- and j- range of indices to check.

  real :: a_nonsym(g%isd:g%ied,g%jsd:g%jed)
  character(len=128) :: mesg2

  integer :: i, j, is_ch, ie_ch, js_ch, je_ch
  integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  is_ch = g%isc ; ie_ch = g%iec ; js_ch = g%jsc ; je_ch = g%jec
  if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie
  if (present(js)) js_ch = js ; if (present(js)) je_ch = je

  ! This only works on points outside of the standard computational domain.
  if ((is_ch == g%isc) .and. (ie_ch == g%iec) .and. &
      (js_ch == g%jsc) .and. (je_ch == g%jec)) return

  do i=isd,ied ; do j=jsd,jed
    a_nonsym(i,j) = array(i,j)
  enddo ; enddo

  call pass_var(a_nonsym, g%Domain)

  do i=is_ch,ie_ch ; do j=js_ch,je_ch
    if (a_nonsym(i,j) /= array(i,j) .and. &
        redundant_prints(1) < max_redundant_prints) then
      write(mesg2,'(" Redundant points",2(1pe12.4)," differ by ", &
                    & 1pe12.4," at i,j = ",2i4," on pe ",i4)') &
           array(i,j), a_nonsym(i,j),array(i,j)-a_nonsym(i,j),i,j,pe_here()
      write(0,'(A130)') trim(mesg)//trim(mesg2)
     redundant_prints(1) = redundant_prints(1) + 1
   endif
  enddo ; enddo

end subroutine  check_redundant_st2d


subroutine check_redundant_vt3d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                               direction)
  character(len=*),                    intent(in)    :: mesg
  type(ocean_grid_type),               intent(inout) :: G    !< The ocean's grid structure
  real, dimension(G%isd:,G%jsd:,:),    intent(in)    :: u_comp
  real, dimension(G%isd:,G%jsd:,:),    intent(in)    :: v_comp
  integer,                   optional, intent(in)    :: is, ie, js, je
  integer,                   optional, intent(in)    :: direction
! Arguments: u_comp - The u-component of the vector being checked.
!  (in)      v_comp - The v-component of the vector being checked.
!  (in)      mesg - A message indicating what is being checked.
!  (in)      G - The ocean's grid structure.
!  (in/opt)  is, ie, js, je - the i- and j- range of indices to check.
!  (in/opt)  direction - the direction flag to be passed to pass_vector.

  character(len=24) :: mesg_k
  integer :: k

  do k=1,size(u_comp,3)
    if (k < 10) then ; write(mesg_k,'(" Layer",i2," ")') k
    elseif (k < 100) then ; write(mesg_k,'(" Layer",i3," ")') k
    elseif (k < 1000) then ; write(mesg_k,'(" Layer",i4," ")') k
    else ; write(mesg_k,'(" Layer",i9," ")') k ; endif

    call check_redundant_vt2d(trim(mesg)//trim(mesg_k), u_comp(:,:,k), &
             v_comp(:,:,k), g, is, ie, js, je, direction)
  enddo
end subroutine  check_redundant_vt3d

subroutine check_redundant_vt2d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                               direction)
  character(len=*),                intent(in)    :: mesg
  type(ocean_grid_type),           intent(inout) :: G    !< The ocean's grid structure
  real, dimension(G%isd:,G%jsd:),  intent(in)    :: u_comp
  real, dimension(G%isd:,G%jsd:),  intent(in)    :: v_comp
  integer,               optional, intent(in)    :: is, ie, js, je
  integer,               optional, intent(in)    :: direction
! Arguments: u_comp - The u-component of the vector being checked.
!  (in)      v_comp - The v-component of the vector being checked.
!  (in)      mesg - A message indicating what is being checked.
!  (in)      G - The ocean's grid structure.
!  (in/opt)  is, ie, js, je - the i- and j- range of indices to check.
!  (in/opt)  direction - the direction flag to be passed to pass_vector.

  real :: u_nonsym(g%isd:g%ied,g%jsd:g%jed)
  real :: v_nonsym(g%isd:g%ied,g%jsd:g%jed)
  character(len=128) :: mesg2

  integer :: i, j, is_ch, ie_ch, js_ch, je_ch
  integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  is_ch = g%isc ; ie_ch = g%iec ; js_ch = g%jsc ; je_ch = g%jec
  if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie
  if (present(js)) js_ch = js ; if (present(js)) je_ch = je

  ! This only works on points outside of the standard computational domain.
  if ((is_ch == g%isc) .and. (ie_ch == g%iec) .and. &
      (js_ch == g%jsc) .and. (je_ch == g%jec)) return

  do i=isd,ied ; do j=jsd,jed
    u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j)
  enddo ; enddo

  call pass_vector(u_nonsym, v_nonsym, g%Domain, direction, stagger=agrid)

  do i=is_ch,ie_ch ; do j=js_ch+1,je_ch
    if (u_nonsym(i,j) /= u_comp(i,j) .and. &
        redundant_prints(1) < max_redundant_prints) then
      write(mesg2,'(" redundant u-components",2(1pe12.4)," differ by ", &
                    & 1pe12.4," at i,j = ",2i4," on pe ",i4)') &
           u_comp(i,j), u_nonsym(i,j),u_comp(i,j)-u_nonsym(i,j),i,j,pe_here()
      write(0,'(A130)') trim(mesg)//trim(mesg2)
      redundant_prints(1) = redundant_prints(1) + 1
    endif
  enddo ; enddo
  do i=is_ch+1,ie_ch ; do j=js_ch,je_ch
    if (v_nonsym(i,j) /= v_comp(i,j) .and. &
        redundant_prints(1) < max_redundant_prints) then
      write(mesg2,'(" redundant v-comps",2(1pe12.4)," differ by ", &
                    & 1pe12.4," at i,j = ",2i4," x,y = ",2(1pe12.4)" on pe ",i4)') &
           v_comp(i,j), v_nonsym(i,j),v_comp(i,j)-v_nonsym(i,j),i,j, &
           g%geoLonBu(i,j), g%geoLatBu(i,j), pe_here()
      write(0,'(A155)') trim(mesg)//trim(mesg2)
      redundant_prints(1) = redundant_prints(1) + 1
    endif
  enddo ; enddo

end subroutine  check_redundant_vt2d

! =====================================================================

! This function does a checksum and redundant point check on a 3d C-grid vector.
subroutine chksum_vec_c3d(mesg, u_comp, v_comp, G, halos, scalars)
  character(len=*),                  intent(in)    :: mesg   !< An identifying message
  type(ocean_grid_type),             intent(inout) :: G      !< The ocean's grid structure
  real, dimension(G%IsdB:,G%jsd:,:), intent(in)    :: u_comp !< The u-component of the vector
  real, dimension(G%isd:,G%JsdB:,:), intent(in)    :: v_comp !< The v-component of the vector
  integer,                 optional, intent(in)    :: halos  !< The width of halos to check (default 0)
  logical,                 optional, intent(in)    :: scalars !< If true this is a pair of
                                                              !! scalars that are being checked.

  logical :: are_scalars
  are_scalars = .false. ; if (present(scalars)) are_scalars = scalars

  if (debug_chksums) then
    call uvchksum(mesg, u_comp, v_comp, g%HI, halos)
  endif
  if (debug_redundant) then
    if (are_scalars) then
      call check_redundant_c(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair)
    else
      call check_redundant_c(mesg, u_comp, v_comp, g)
    endif
  endif

end subroutine chksum_vec_c3d

! This function does a checksum and redundant point check on a 2d C-grid vector.
subroutine chksum_vec_c2d(mesg, u_comp, v_comp, G, halos, scalars)
  character(len=*),                intent(in)    :: mesg   !< An identifying message
  type(ocean_grid_type),           intent(inout) :: G      !< The ocean's grid structure
  real, dimension(G%IsdB:,G%jsd:), intent(in)    :: u_comp !< The u-component of the vector
  real, dimension(G%isd:,G%JsdB:), intent(in)    :: v_comp !< The v-component of the vector
  integer,               optional, intent(in)    :: halos  !< The width of halos to check (default 0)
  logical,               optional, intent(in)    :: scalars !< If true this is a pair of
                                                            !! scalars that are being checked.

  logical :: are_scalars
  are_scalars = .false. ; if (present(scalars)) are_scalars = scalars

  if (debug_chksums) then
    call uvchksum(mesg, u_comp, v_comp, g%HI, halos)
  endif
  if (debug_redundant) then
    if (are_scalars) then
      call check_redundant_c(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair)
    else
      call check_redundant_c(mesg, u_comp, v_comp, g)
    endif
  endif

end subroutine chksum_vec_c2d

! This function does a checksum and redundant point check on a 3d B-grid vector.
subroutine chksum_vec_b3d(mesg, u_comp, v_comp, G, halos, scalars)
  character(len=*),                   intent(in)    :: mesg   !< An identifying message
  type(ocean_grid_type),              intent(inout) :: G      !< The ocean's grid structure
  real, dimension(G%IsdB:,G%JsdB:,:), intent(in)    :: u_comp !< The u-component of the vector
  real, dimension(G%IsdB:,G%JsdB:,:), intent(in)    :: v_comp !< The v-component of the vector
  integer,                  optional, intent(in)    :: halos  !< The width of halos to check (default 0)
  logical,                  optional, intent(in)    :: scalars !< If true this is a pair of
                                                               !! scalars that are being checked.

  logical :: are_scalars
  are_scalars = .false. ; if (present(scalars)) are_scalars = scalars

  if (debug_chksums) then
    call bchksum(u_comp, mesg//"(u)", g%HI, halos)
    call bchksum(v_comp, mesg//"(v)", g%HI, halos)
  endif
  if (debug_redundant) then
    if (are_scalars) then
      call check_redundant_b(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair)
    else
      call check_redundant_b(mesg, u_comp, v_comp, g)
    endif
  endif

end subroutine chksum_vec_b3d

! This function does a checksum and redundant point check on a 2d B-grid vector.
subroutine chksum_vec_b2d(mesg, u_comp, v_comp, G, halos, scalars, symmetric)
  character(len=*),                 intent(in)    :: mesg   !< An identifying message
  type(ocean_grid_type),            intent(inout) :: G      !< The ocean's grid structure
  real, dimension(G%IsdB:,G%JsdB:), intent(in)    :: u_comp !< The u-component of the vector
  real, dimension(G%IsdB:,G%JsdB:), intent(in)    :: v_comp !< The v-component of the vector
  integer,                optional, intent(in)    :: halos  !< The width of halos to check (default 0)
  logical,                optional, intent(in)    :: scalars !< If true this is a pair of
                                                             !! scalars that are being checked.
  logical,                optional, intent(in)    :: symmetric !< If true, do the checksums on the
                                                               !! full symmetric computational domain.

  logical :: are_scalars
  are_scalars = .false. ; if (present(scalars)) are_scalars = scalars

  if (debug_chksums) then
    call bchksum(u_comp, mesg//"(u)", g%HI, halos, symmetric=symmetric)
    call bchksum(v_comp, mesg//"(v)", g%HI, halos, symmetric=symmetric)
  endif
  if (debug_redundant) then
    if (are_scalars) then
      call check_redundant_b(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair)
    else
      call check_redundant_b(mesg, u_comp, v_comp, g)
    endif
  endif

end subroutine chksum_vec_b2d

! This function does a checksum and redundant point check on a 3d C-grid vector.
subroutine chksum_vec_a3d(mesg, u_comp, v_comp, G, halos, scalars)
  character(len=*),                 intent(in)    :: mesg   !< An identifying message
  type(ocean_grid_type),            intent(inout) :: G      !< The ocean's grid structure
  real, dimension(G%isd:,G%jsd:,:), intent(in)    :: u_comp !< The u-component of the vector
  real, dimension(G%isd:,G%jsd:,:), intent(in)    :: v_comp !< The v-component of the vector
  integer,                optional, intent(in)    :: halos  !< The width of halos to check (default 0)
  logical,                optional, intent(in)    :: scalars !< If true this is a pair of
                                                             !! scalars that are being checked.

  logical :: are_scalars
  are_scalars = .false. ; if (present(scalars)) are_scalars = scalars

  if (debug_chksums) then
    call hchksum(u_comp, mesg//"(u)", g%HI, halos)
    call hchksum(v_comp, mesg//"(v)", g%HI, halos)
  endif
  if (debug_redundant) then
    if (are_scalars) then
      call check_redundant_t(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair)
    else
      call check_redundant_t(mesg, u_comp, v_comp, g)
    endif
  endif

end subroutine chksum_vec_a3d


! This function does a checksum and redundant point check on a 2d C-grid vector.
subroutine chksum_vec_a2d(mesg, u_comp, v_comp, G, halos, scalars)
  character(len=*),               intent(in)    :: mesg   !< An identifying message
  type(ocean_grid_type),          intent(inout) :: G      !< The ocean's grid structure
  real, dimension(G%isd:,G%jsd:), intent(in)    :: u_comp !< The u-component of the vector
  real, dimension(G%isd:,G%jsd:), intent(in)    :: v_comp !< The v-component of the vector
  integer,              optional, intent(in)    :: halos  !< The width of halos to check (default 0)
  logical,              optional, intent(in)    :: scalars !< If true this is a pair of
                                                           !! scalars that are being checked.

  logical :: are_scalars
  are_scalars = .false. ; if (present(scalars)) are_scalars = scalars

  if (debug_chksums) then
    call hchksum(u_comp, mesg//"(u)", g%HI, halos)
    call hchksum(v_comp, mesg//"(v)", g%HI, halos)
  endif
  if (debug_redundant) then
    if (are_scalars) then
      call check_redundant_t(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair)
    else
      call check_redundant_t(mesg, u_comp, v_comp, g)
    endif
  endif

end subroutine chksum_vec_a2d


! =====================================================================

!> This function returns the sum over computational domain of all
!! processors of hThick*stuff, where stuff is a 3-d array at tracer points.
function totalstuff(HI, hThick, areaT, stuff)
  type(hor_index_type),               intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: hThick !< The array of thicknesses to use as weights
  real, dimension(HI%isd:,HI%jsd:),   intent(in) :: areaT  !< The array of cell areas in m2
  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: stuff  !< The array of stuff to be summed
  real                                         :: totalStuff
  integer :: i, j, k, nz

  nz = size(hthick,3)
  totalstuff = 0.
  do k = 1, nz ; do j = hi%jsc, hi%jec ; do i = hi%isc, hi%iec
    totalstuff = totalstuff + hthick(i,j,k) * stuff(i,j,k) * areat(i,j)
  enddo ; enddo ; enddo
  call sum_across_pes(totalstuff)

end function totalstuff

! =====================================================================

!> This subroutine display the total thickness, temperature and salinity
!! as well as the change since the last call.
!! NOTE: This subroutine uses "save" data which is not thread safe and is purely
!! for extreme debugging without a proper debugger.
subroutine totaltands(HI, hThick, areaT, temperature, salinity, mesg)
  type(hor_index_type),               intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: hThick !< The array of thicknesses to use as weights
  real, dimension(HI%isd:,HI%jsd:),   intent(in) :: areaT  !< The array of cell areas in m2
  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: temperature !< The temperature field to sum
  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: salinity    !< The salinity field to sum
  character(len=*),                   intent(in) :: mesg        !< An identifying message

  ! NOTE: This subroutine uses "save" data which is not thread safe and is purely for
  ! extreme debugging without a proper debugger.
  real, save :: totalH = 0., totalt = 0., totals = 0.

  logical, save :: firstCall = .true.
  real :: thisH, thisT, thisS, delH, delT, delS
  integer :: i, j, k, nz

  nz = size(hthick,3)
  thish = 0.
  do k = 1, nz ; do j = hi%jsc, hi%jec ; do i = hi%isc, hi%iec
    thish = thish + hthick(i,j,k) * areat(i,j)
  enddo ; enddo ; enddo
  call sum_across_pes(thish)
  thist = totalstuff(hi, hthick, areat, temperature)
  thiss = totalstuff(hi, hthick, areat, salinity)

  if (is_root_pe()) then
    if (firstcall) then
      totalh = thish ; totalt = thist ; totals = thiss
      write(0,*) 'Totals H,T,S:',thish,thist,thiss,' ',mesg
      firstcall = .false.
    else
      delh = thish - totalh
      delt = thist - totalt
      dels = thiss - totals
      totalh = thish ; totalt = thist ; totals = thiss
      write(0,*) 'Tot/del H,T,S:',thish,thist,thiss,delh,delt,dels,' ',mesg
    endif
  endif

end subroutine totaltands

!> Returns false if the column integral of a given quantity is within roundoff
logical function check_column_integral(nk, field, known_answer)
  integer,             intent(in) :: nk           !< Number of levels in column
  real, dimension(nk), intent(in) :: field        !< Field to be summed
  real, optional,      intent(in) :: known_answer !< If present is the expected sum,
                                                  !! If missing, assumed zero
  ! Local variables
  real    :: u_sum, error, expected
  integer :: k

  u_sum = field(1)
  error = 0.

  ! Reintegrate and sum roundoff errors
  do k=2,nk
    u_sum = u_sum + field(k)
    error = error + epsilon(u_sum)*max(abs(u_sum),abs(field(k)))
  enddo

  ! Assign expected answer to either the optional input or 0
  if (present(known_answer)) then
    expected = known_answer
  else
    expected = 0.
  endif

  ! Compare the column integrals against calculated roundoff error
  if (abs(u_sum-expected) > error) then
    check_column_integral = .true.
  else
    check_column_integral = .false.
  endif

end function check_column_integral

!> Returns false if the column integrals of two given quantities are within roundoff of each other
logical function check_column_integrals(nk_1, field_1, nk_2, field_2, missing_value)
  integer,               intent(in) :: nk_1           !< Number of levels in field 1
  integer,               intent(in) :: nk_2           !< Number of levels in field 2
  real, dimension(nk_1), intent(in) :: field_1        !< First field to be summed
  real, dimension(nk_2), intent(in) :: field_2        !< Second field to be summed
  real, optional,        intent(in) :: missing_value  !< If column contains missing values,
                                                      !! mask them from the sum


  ! Local variables
  real    :: u1_sum, error1, u2_sum, error2, misval
  integer :: k

  ! Assign missing value
  if (present(missing_value)) then
    misval = missing_value
  else
    misval = 0.
  endif

  u1_sum = field_1(1)
  error1 = 0.

  ! Reintegrate and sum roundoff errors
  do k=2,nk_1
    if (field_1(k)/=misval) then
      u1_sum = u1_sum + field_1(k)
      error1 = error1 + epsilon(u1_sum)*max(abs(u1_sum),abs(field_1(k)))
    endif
  enddo

  u2_sum = field_2(1)
  error2 = 0.

  ! Reintegrate and sum roundoff errors
  do k=2,nk_2
    if (field_2(k)/=misval) then
      u2_sum = u2_sum + field_2(k)
      error2 = error2 + epsilon(u2_sum)*max(abs(u2_sum),abs(field_2(k)))
    endif
  enddo

  ! Compare the column integrals against calculated roundoff error
  if (abs(u1_sum-u2_sum) > (error1+error2)) then
    check_column_integrals = .true.
  else
    check_column_integrals = .false.
  endif

end function check_column_integrals

end module mom_debugging