MOM_debugging.F90

!> Provides checksumming functions for debugging
!!
!! 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.
module mom_debugging

! This file is part of MOM6. See LICENSE.md for the license.

use mom_checksums, only : hchksum, bchksum, qchksum, uvchksum, hchksum_pair
use mom_checksums, only : is_nan, chksum, mom_checksums_init
use mom_coms, only : pe_here, root_pe, num_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, uvchksum, hchksum_pair

!> Check for consistency between the duplicated points of a C-grid vector
interface check_redundant
  module procedure check_redundant_vc3d, check_redundant_vc2d
end interface check_redundant
!> Check for consistency between the duplicated points of a C-grid vector
interface check_redundant_c
  module procedure check_redundant_vc3d, check_redundant_vc2d
end interface check_redundant_c
!> Check for consistency between the duplicated points of a B-grid vector or scalar
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
!> Check for consistency between the duplicated points of an A-grid vector or scalar
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

!> Do checksums on the components of a C-grid vector
interface vec_chksum
  module procedure chksum_vec_c3d, chksum_vec_c2d
end interface vec_chksum
!> Do checksums on the components of a C-grid vector
interface vec_chksum_c
  module procedure chksum_vec_c3d, chksum_vec_c2d
end interface vec_chksum_c
!> Do checksums on the components of a B-grid vector
interface vec_chksum_b
  module procedure chksum_vec_b3d, chksum_vec_b2d
end interface vec_chksum_b
!> Do checksums on the components of an A-grid vector
interface vec_chksum_a
  module procedure chksum_vec_a3d, chksum_vec_a2d
end interface vec_chksum_a

! Note: these parameters are module data but ONLY used when debugging and
!       so can violate the thread-safe requirement of no module/global data.
integer :: max_redundant_prints = 100 !< Maximum number of times to write redundant messages
integer :: redundant_prints(3) = 0 !< Counters for controlling redundant printing
logical :: debug = .false. !< Write out verbose debugging data
logical :: debug_chksums = .true. !< Perform checksums on arrays
logical :: debug_redundant = .true. !< Check redundant values on PE boundaries

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, debugging=.true.)
  call get_param(param_file, mdl, "DEBUG", debug, &
                 "If true, write out verbose debugging data.", &
                 default=.false., debuggingparam=.true.)
  call get_param(param_file, mdl, "DEBUG_CHKSUMS", debug_chksums, &
                 "If true, checksums are performed on arrays in the "//&
                 "various vec_chksum routines.", default=debug, &
                 debuggingparam=.true.)
  call get_param(param_file, mdl, "DEBUG_REDUNDANT", debug_redundant, &
                 "If true, debug redundant data points during calls to "//&
                 "the various vec_chksum routines.", default=debug, &
                 debuggingparam=.true.)

  call mom_checksums_init(param_file)

end subroutine mom_debugging_init

!> Check for consistency between the duplicated points of a 3-D C-grid vector
subroutine check_redundant_vc3d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                                direction)
  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
                                                               !! to be checked for consistency
  real, dimension(G%isd:,G%JsdB:,:),   intent(in)    :: v_comp !< The u-component of the vector
                                                               !! to be checked for consistency
  integer,                   optional, intent(in)    :: is     !< The starting i-index to check
  integer,                   optional, intent(in)    :: ie     !< The ending i-index to check
  integer,                   optional, intent(in)    :: js     !< The starting j-index to check
  integer,                   optional, intent(in)    :: je     !< The ending j-index to check
  integer,                   optional, intent(in)    :: direction !< the direction flag to be
                                                               !! passed to pass_vector
  ! Local variables
  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

!> Check for consistency between the duplicated points of a 2-D C-grid vector
subroutine check_redundant_vc2d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                                direction)
  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
                                                           !! to be checked for consistency
  real, dimension(G%isd:,G%JsdB:), intent(in)    :: v_comp !< The u-component of the vector
                                                           !! to be checked for consistency
  integer,               optional, intent(in)    :: is     !< The starting i-index to check
  integer,               optional, intent(in)    :: ie     !< The ending i-index to check
  integer,               optional, intent(in)    :: js     !< The starting j-index to check
  integer,               optional, intent(in)    :: je     !< The ending j-index to check
  integer,               optional, intent(in)    :: direction !< the direction flag to be
                                                           !! passed to pass_vector
  ! Local variables
  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

!> Check for consistency between the duplicated points of a 3-D scalar at corner points
subroutine check_redundant_sb3d(mesg, array, G, is, ie, js, je)
  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)    :: array !< The array to be checked for consistency
  integer,                    optional, intent(in)    :: is    !< The starting i-index to check
  integer,                    optional, intent(in)    :: ie    !< The ending i-index to check
  integer,                    optional, intent(in)    :: js    !< The starting j-index to check
  integer,                    optional, intent(in)    :: je    !< The ending j-index to check

  ! Local variables
  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

!> Check for consistency between the duplicated points of a 2-D scalar at corner points
subroutine check_redundant_sb2d(mesg, array, G, is, ie, js, je)
  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)    :: array !< The array to be checked for consistency
  integer,                optional, intent(in)    :: is    !< The starting i-index to check
  integer,                optional, intent(in)    :: ie    !< The ending i-index to check
  integer,                optional, intent(in)    :: js    !< The starting j-index to check
  integer,                optional, intent(in)    :: je    !< The ending j-index to check
  ! Local variables
  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

!> Check for consistency between the duplicated points of a 3-D B-grid vector
subroutine check_redundant_vb3d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                                direction)
  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
                                                           !! to be checked for consistency
  real, dimension(G%IsdB:,G%JsdB:,:),  intent(in)    :: v_comp !< The v-component of the vector
                                                           !! to be checked for consistency
  integer,                   optional, intent(in)    :: is     !< The starting i-index to check
  integer,                   optional, intent(in)    :: ie     !< The ending i-index to check
  integer,                   optional, intent(in)    :: js     !< The starting j-index to check
  integer,                   optional, intent(in)    :: je     !< The ending j-index to check
  integer,                   optional, intent(in)    :: direction !< the direction flag to be
                                                               !! passed to pass_vector
  ! Local variables
  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

!> Check for consistency between the duplicated points of a 2-D B-grid vector
subroutine check_redundant_vb2d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                                direction)
  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
                                                            !! to be checked for consistency
  real, dimension(G%IsdB:,G%JsdB:), intent(in)    :: v_comp !< The v-component of the vector
                                                            !! to be checked for consistency
  integer,                optional, intent(in)    :: is     !< The starting i-index to check
  integer,                optional, intent(in)    :: ie     !< The ending i-index to check
  integer,                optional, intent(in)    :: js     !< The starting j-index to check
  integer,                optional, intent(in)    :: je     !< The ending j-index to check
  integer,                optional, intent(in)    :: direction !< the direction flag to be
                                                            !! passed to pass_vector
  ! Local variables
  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

!> Check for consistency between the duplicated points of a 3-D scalar at tracer points
subroutine check_redundant_st3d(mesg, array, G, is, ie, js, je)
  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)    :: array !< The array to be checked for consistency
  integer,                    optional, intent(in)    :: is    !< The starting i-index to check
  integer,                    optional, intent(in)    :: ie    !< The ending i-index to check
  integer,                    optional, intent(in)    :: js    !< The starting j-index to check
  integer,                    optional, intent(in)    :: je    !< The ending j-index to check
  ! Local variables
  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


!> Check for consistency between the duplicated points of a 2-D scalar at tracer points
subroutine check_redundant_st2d(mesg, array, G, is, ie, js, je)
  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)    :: array !< The array to be checked for consistency
  integer,                optional, intent(in)    :: is    !< The starting i-index to check
  integer,                optional, intent(in)    :: ie    !< The ending i-index to check
  integer,                optional, intent(in)    :: js    !< The starting j-index to check
  integer,                optional, intent(in)    :: je    !< The ending j-index to check
  ! Local variables
  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

!> Check for consistency between the duplicated points of a 3-D A-grid vector
subroutine check_redundant_vt3d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                               direction)
  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
                                                           !! to be checked for consistency
  real, dimension(G%isd:,G%jsd:,:),    intent(in)    :: v_comp !< The v-component of the vector
                                                           !! to be checked for consistency
  integer,                   optional, intent(in)    :: is     !< The starting i-index to check
  integer,                   optional, intent(in)    :: ie     !< The ending i-index to check
  integer,                   optional, intent(in)    :: js     !< The starting j-index to check
  integer,                   optional, intent(in)    :: je     !< The ending j-index to check
  integer,                   optional, intent(in)    :: direction !< the direction flag to be
                                                           !! passed to pass_vector
  ! Local variables
  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

!> Check for consistency between the duplicated points of a 2-D A-grid vector
subroutine check_redundant_vt2d(mesg, u_comp, v_comp, G, is, ie, js, je, &
                               direction)
  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
                                                           !! to be checked for consistency
  real, dimension(G%isd:,G%jsd:),  intent(in)    :: v_comp !< The v-component of the vector
                                                           !! to be checked for consistency
  integer,               optional, intent(in)    :: is     !< The starting i-index to check
  integer,               optional, intent(in)    :: ie     !< The ending i-index to check
  integer,               optional, intent(in)    :: js     !< The starting j-index to check
  integer,               optional, intent(in)    :: je     !< The ending j-index to check
  integer,               optional, intent(in)    :: direction !< the direction flag to be
                                                           !! passed to pass_vector
  ! Local variables
  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

!> Do 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.
  ! Local variables
  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

!> Do 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.
  ! Local variables
  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

!> Do 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.
  ! Local variables
  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

! Do 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.
  ! Local variables
  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

!> Do 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.
  ! Local variables
  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

!> Do 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.
  ! Local variables
  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 [m2]
  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: stuff  !< The array of stuff to be summed
  real                                         :: totalstuff !< the globally integrated amoutn of stuff
  ! Local variables
  real, dimension(HI%isc:HI%iec, HI%jsc:HI%jec) :: tmp_for_sum
  integer :: i, j, k, nz

  nz = size(hthick,3)
  tmp_for_sum(:,:) = 0.0
  do k=1,nz ; do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
    tmp_for_sum(i,j) = tmp_for_sum(i,j) + hthick(i,j,k) * stuff(i,j,k) * areat(i,j)
  enddo ; enddo ; enddo
  totalstuff = reproducing_sum(tmp_for_sum)

end function totalstuff

!> This subroutine display the total thickness, temperature and salinity
!! as well as the change since the last call.
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 [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.
  ! Local variables
  logical, save :: firstcall = .true.
  real, dimension(HI%isc:HI%iec, HI%jsc:HI%jec) :: tmp_for_sum
  real :: thish, thist, thiss, delh, delt, dels
  integer :: i, j, k, nz

  nz = size(hthick,3)
  tmp_for_sum(:,:) = 0.0
  do k=1,nz ; do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
    tmp_for_sum(i,j) = tmp_for_sum(i,j) + hthick(i,j,k) * areat(i,j)
  enddo ; enddo ; enddo
  thish = reproducing_sum(tmp_for_sum)
  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