MOM_checksums.F90

Go to the documentation of this file.

module mom_checksums

!***********************************************************************
!*                   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                      *
!***********************************************************************

use mom_coms, only : pe_here, root_pe, num_pes, sum_across_pes
use mom_coms, only : min_across_pes, max_across_pes
use mom_coms, only : reproducing_sum
use mom_error_handler, only : mom_error, fatal, is_root_pe
use mom_file_parser, only : log_version, param_file_type
use mom_hor_index, only : hor_index_type

implicit none ; private

public :: hchksum, bchksum, uchksum, vchksum, qchksum, is_nan, chksum
public :: hchksum_pair, uvchksum, bchksum_pair
public :: mom_checksums_init

interface hchksum_pair
  module procedure chksum_pair_h_2d, chksum_pair_h_3d
end interface

interface uvchksum
  module procedure chksum_uv_2d, chksum_uv_3d
end interface

interface uchksum
  module procedure chksum_u_2d, chksum_u_3d
end interface

interface vchksum
  module procedure chksum_v_2d, chksum_v_3d
end interface

interface bchksum_pair
  module procedure chksum_pair_b_2d, chksum_pair_b_3d
end interface

interface hchksum
  module procedure chksum_h_2d, chksum_h_3d
end interface

interface bchksum
  module procedure chksum_b_2d, chksum_b_3d
end interface

! This is an older interface that has been renamed Bchksum
interface qchksum
  module procedure chksum_b_2d, chksum_b_3d
end interface

interface chksum
  module procedure chksum1d, chksum2d, chksum3d
end interface

interface chk_sum_msg
  module procedure chk_sum_msg1, chk_sum_msg2, chk_sum_msg3, chk_sum_msg5
end interface

interface is_nan
  module procedure is_nan_0d, is_nan_1d, is_nan_2d, is_nan_3d
end interface

integer, parameter :: default_shift=0
logical :: calculatestatistics=.true. ! If true, report min, max and mean.
logical :: writechksums=.true. ! If true, report the bitcount checksum
logical :: checkfornans=.true. ! If true, checks array for NaNs and cause
                               ! FATAL error is any are found

contains

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

subroutine chksum_pair_h_2d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, scale)
  character(len=*),                 intent(in) :: mesg !< Identifying messages
  type(hor_index_type),             intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayA, arrayB !< The arrays to be checksummed
  integer,                optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,                optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                   optional, intent(in) :: scale     !< A scaling factor for this array.

  if (present(haloshift)) then
    call chksum_h_2d(arraya, 'x '//mesg, hi, haloshift, omit_corners, scale=scale)
    call chksum_h_2d(arrayb, 'y '//mesg, hi, haloshift, omit_corners, scale=scale)
  else
    call chksum_h_2d(arraya, 'x '//mesg, hi, scale=scale)
    call chksum_h_2d(arrayb, 'y '//mesg, hi, scale=scale)
  endif

end subroutine chksum_pair_h_2d

subroutine chksum_pair_h_3d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, scale)
  character(len=*),                    intent(in) :: mesg !< Identifying messages
  type(hor_index_type),                intent(in) :: HI   !< A horizontal index type
  real, dimension(HI%isd:,HI%jsd:, :), intent(in) :: arrayA, arrayB !< The arrays to be checksummed
  integer,                   optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,                   optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                      optional, intent(in) :: scale     !< A scaling factor for this array.

  if (present(haloshift)) then
    call chksum_h_3d(arraya, 'x '//mesg, hi, haloshift, omit_corners, scale=scale)
    call chksum_h_3d(arrayb, 'y '//mesg, hi, haloshift, omit_corners, scale=scale)
  else
    call chksum_h_3d(arraya, 'x '//mesg, hi, scale=scale)
    call chksum_h_3d(arrayb, 'y '//mesg, hi, scale=scale)
  endif

end subroutine chksum_pair_h_3d

!> chksum_h_2d performs checksums on a 2d array staggered at tracer points.
subroutine chksum_h_2d(array, mesg, HI, haloshift, omit_corners, scale)
  type(hor_index_type),            intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%isd:,HI%jsd:), intent(in) :: array !< The array to be checksummed
  character(len=*),                intent(in) :: mesg  !< An identifying message
  integer,               optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,               optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                  optional, intent(in) :: scale     !< A scaling factor for this array.

  real :: scaling
  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
  integer :: bcN, bcS, bcE, bcW
  logical :: do_corners

  if (checkfornans) then
    if (is_nan(array(hi%isc:hi%iec,hi%jsc:hi%jec))) &
      call chksum_error(fatal, 'NaN detected: '//trim(mesg))
!   if (is_NaN(array)) &
!     call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
  endif
  scaling = 1.0 ; if (present(scale)) scaling = scale

  if (calculatestatistics) call substats(hi, array, mesg, scaling)

  if (.not.writechksums) return

  hshift = default_shift
  if (present(haloshift)) hshift = haloshift
  if (hshift<0) hshift = hi%ied-hi%iec

  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
       hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
    write(0,*) 'chksum_h_2d: haloshift =',hshift
    write(0,*) 'chksum_h_2d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
    write(0,*) 'chksum_h_2d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
    call chksum_error(fatal,'Error in chksum_h_2d '//trim(mesg))
  endif

  bc0 = subchk(array, hi, 0, 0, scaling)

  if (hshift==0) then
    if (is_root_pe()) call chk_sum_msg("h-point:",bc0,mesg)
    return
  endif

  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners

  if (do_corners) then
    bcsw = subchk(array, hi, -hshift, -hshift, scaling)
    bcse = subchk(array, hi, hshift, -hshift, scaling)
    bcnw = subchk(array, hi, -hshift, hshift, scaling)
    bcne = subchk(array, hi, hshift, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg("h-point:",bc0,bcsw,bcse,bcnw,bcne,mesg)
  else
    bcs = subchk(array, hi, 0, -hshift, scaling)
    bce = subchk(array, hi, hshift, 0, scaling)
    bcw = subchk(array, hi, -hshift, 0, scaling)
    bcn = subchk(array, hi, 0, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg_nsew("h-point:",bc0,bcn,bcs,bce,bcw,mesg)
  endif

  contains

  integer function subchk(array, HI, di, dj, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%isd:,HI%jsd:), intent(in) :: array
    integer, intent(in) :: di, dj
    real, intent(in) :: scale
    integer :: bitcount, i, j, bc
    subchk = 0
    do j=hi%jsc+dj,hi%jec+dj; do i=hi%isc+di,hi%iec+di
      bc = bitcount(abs(scale*array(i,j)))
      subchk = subchk + bc
    enddo; enddo
    call sum_across_pes(subchk)
    subchk=mod(subchk,1000000000)
  end function subchk

  subroutine substats(HI, array, mesg, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%isd:,HI%jsd:), intent(in) :: array
    character(len=*), intent(in) :: mesg
    real, intent(in) :: scale
    integer :: i, j, n
    real :: aMean, aMin, aMax

    amin = array(hi%isc,hi%jsc)
    amax = array(hi%isc,hi%jsc)
    n = 0
    do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
      amin = min(amin, array(i,j))
      amax = max(amax, array(i,j))
      n = n + 1
    enddo ; enddo
    amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec))
    call sum_across_pes(n)
    call min_across_pes(amin)
    call max_across_pes(amax)
    amean = amean / real(n)
    if (is_root_pe()) call chk_sum_msg("h-point:",amean*scale,amin*scale,amax*scale,mesg)
  end subroutine substats

end subroutine chksum_h_2d

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

subroutine chksum_pair_b_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, omit_corners, scale)
  character(len=*),                 intent(in) :: mesg   !< Identifying messages
  type(hor_index_type),             intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayA, arrayB !< The arrays to be checksummed
  logical,                optional, intent(in) :: symmetric !< If true, do the checksums on the full symmetric computational domain.
  integer,                optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,                optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                   optional, intent(in) :: scale     !< A scaling factor for this array.

  logical :: sym

  sym = .false. ; if (present(symmetric)) sym = symmetric

  if (present(haloshift)) then
    call chksum_b_2d(arraya, 'x '//mesg, hi, haloshift, symmetric=sym, &
                     omit_corners=omit_corners, scale=scale)
    call chksum_b_2d(arrayb, 'y '//mesg, hi, haloshift, symmetric=sym, &
                     omit_corners=omit_corners, scale=scale)
  else
    call chksum_b_2d(arraya, 'x '//mesg, hi, symmetric=sym, scale=scale)
    call chksum_b_2d(arrayb, 'y '//mesg, hi, symmetric=sym, scale=scale)
  endif

end subroutine chksum_pair_b_2d

subroutine chksum_pair_b_3d(mesg, arrayA, arrayB, HI, haloshift, symmetric, omit_corners, scale)
  character(len=*),                    intent(in) :: mesg !< Identifying messages
  type(hor_index_type),                intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%IsdB:,HI%JsdB:, :), intent(in) :: arrayA, arrayB !< The arrays to be checksummed
  integer,                   optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,                   optional, intent(in) :: symmetric !< If true, do the checksums on the full symmetric computational domain.
  logical,                   optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                      optional, intent(in) :: scale     !< A scaling factor for this array.

  logical :: sym

  if (present(haloshift)) then
    call chksum_b_3d(arraya, 'x '//mesg, hi, haloshift, symmetric, &
                     omit_corners, scale=scale)
    call chksum_b_3d(arrayb, 'y '//mesg, hi, haloshift, symmetric, &
                     omit_corners, scale=scale)
  else
    call chksum_b_3d(arraya, 'x '//mesg, hi, symmetric=symmetric, scale=scale)
    call chksum_b_3d(arrayb, 'y '//mesg, hi, symmetric=symmetric, scale=scale)
  endif

end subroutine chksum_pair_b_3d

!> chksum_B_2d performs checksums on a 2d array staggered at corner points.
subroutine chksum_b_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scale)
  type(hor_index_type), intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%IsdB:,HI%JsdB:), &
                        intent(in) :: array !< The array to be checksummed
  character(len=*),     intent(in) :: mesg  !< An identifying message
  integer,    optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,    optional, intent(in) :: symmetric !< If true, do the checksums on the
                                                !! full symmetric computational domain.
  logical,    optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,       optional, intent(in) :: scale     !< A scaling factor for this array.

  real :: scaling
  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
  integer :: bcN, bcS, bcE, bcW
  logical :: do_corners, sym, sym_stats

  if (checkfornans) then
    if (is_nan(array(hi%IscB:hi%IecB,hi%JscB:hi%JecB))) &
      call chksum_error(fatal, 'NaN detected: '//trim(mesg))
!   if (is_NaN(array)) &
!     call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
  endif
  scaling = 1.0 ; if (present(scale)) scaling = scale

  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
  if (calculatestatistics) call substats(hi, array, mesg, sym_stats, scaling)

  if (.not.writechksums) return

  hshift = default_shift
  if (present(haloshift)) hshift = haloshift
  if (hshift<0) hshift = hi%ied-hi%iec

  if ( hi%iscB-hshift<hi%isdB .or. hi%iecB+hshift>hi%iedB .or. &
       hi%jscB-hshift<hi%jsdB .or. hi%jecB+hshift>hi%jedB ) then
    write(0,*) 'chksum_B_2d: haloshift =',hshift
    write(0,*) 'chksum_B_2d: isd,isc,iec,ied=',hi%isdB,hi%iscB,hi%iecB,hi%iedB
    write(0,*) 'chksum_B_2d: jsd,jsc,jec,jed=',hi%jsdB,hi%jscB,hi%jecB,hi%jedB
    call chksum_error(fatal,'Error in chksum_B_2d '//trim(mesg))
  endif

  bc0 = subchk(array, hi, 0, 0, scaling)

  sym = .false. ; if (present(symmetric)) sym = symmetric

  if ((hshift==0) .and. .not.sym) then
    if (is_root_pe()) call chk_sum_msg("B-point:",bc0,mesg)
    return
  endif

  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners

  if (do_corners) then
    if (sym) then
      bcsw = subchk(array, hi, -hshift-1, -hshift-1, scaling)
      bcse = subchk(array, hi, hshift, -hshift-1, scaling)
      bcnw = subchk(array, hi, -hshift-1, hshift, scaling)
    else
      bcsw = subchk(array, hi, -hshift, -hshift, scaling)
      bcse = subchk(array, hi, hshift, -hshift, scaling)
      bcnw = subchk(array, hi, -hshift, hshift, scaling)
    endif
    bcne = subchk(array, hi, hshift, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg("B-point:",bc0,bcsw,bcse,bcnw,bcne,mesg)
  else
    bcs = subchk(array, hi, 0, -hshift, scaling)
    bce = subchk(array, hi, hshift, 0, scaling)
    bcw = subchk(array, hi, -hshift, 0, scaling)
    bcn = subchk(array, hi, 0, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg_nsew("B-point:",bc0,bcn,bcs,bce,bcw,mesg)
  endif

  contains

  integer function subchk(array, HI, di, dj, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%IsdB:,HI%JsdB:), intent(in) :: array
    integer, intent(in) :: di, dj
    real, intent(in) :: scale
    integer :: bitcount, i, j, bc
    subchk = 0
    ! This line deliberately uses the h-point computational domain.
    do j=hi%jsc+dj,hi%jec+dj; do i=hi%isc+di,hi%iec+di
      bc = bitcount(abs(scale*array(i,j)))
      subchk = subchk + bc
    enddo; enddo
    call sum_across_pes(subchk)
    subchk=mod(subchk,1000000000)
  end function subchk

  subroutine substats(HI, array, mesg, sym_stats, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%IsdB:,HI%JsdB:), intent(in) :: array
    character(len=*), intent(in) :: mesg
    logical, intent(in) :: sym_stats
    real, intent(in) :: scale
    integer :: i, j, n, IsB, JsB
    real :: aMean, aMin, aMax

    isb = hi%isc ; if (sym_stats) isb = hi%isc-1
    jsb = hi%jsc ; if (sym_stats) jsb = hi%jsc-1

    amin = array(hi%isc,hi%jsc) ; amax = amin
    do j=jsb,hi%JecB ; do i=isb,hi%IecB
      amin = min(amin, array(i,j))
      amax = max(amax, array(i,j))
    enddo ; enddo
    ! This line deliberately uses the h-point computational domain.
    amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec))
    n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc)
    call sum_across_pes(n)
    call min_across_pes(amin)
    call max_across_pes(amax)
    amean = amean / real(n)
    if (is_root_pe()) call chk_sum_msg("B-point:",amean*scale,amin*scale,amax*scale,mesg)
  end subroutine substats

end subroutine chksum_b_2d

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

subroutine chksum_uv_2d(mesg, arrayU, arrayV, HI, haloshift, symmetric, omit_corners, scale)
  character(len=*),                  intent(in) :: mesg   !< Identifying messages
  type(hor_index_type),              intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: arrayU !< The u-component array to be checksummed
  real, dimension(HI%isd:,HI%JsdB:), intent(in) :: arrayV !< The v-component array to be checksummed
  integer,                 optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,                 optional, intent(in) :: symmetric !< If true, do the checksums on the full symmetric computational domain.
  logical,                 optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                    optional, intent(in) :: scale     !< A scaling factor for these arrays.

  if (present(haloshift)) then
    call chksum_u_2d(arrayu, 'u '//mesg, hi, haloshift, symmetric, omit_corners, scale)
    call chksum_v_2d(arrayv, 'v '//mesg, hi, haloshift, symmetric, omit_corners, scale)
  else
    call chksum_u_2d(arrayu, 'u '//mesg, hi, symmetric=symmetric)
    call chksum_v_2d(arrayv, 'v '//mesg, hi, symmetric=symmetric)
  endif

end subroutine chksum_uv_2d

subroutine chksum_uv_3d(mesg, arrayU, arrayV, HI, haloshift, symmetric, omit_corners, scale)
  character(len=*),                    intent(in) :: mesg   !< Identifying messages
  type(hor_index_type),                intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%IsdB:,HI%jsd:,:), intent(in) :: arrayU !< The u-component array to be checksummed
  real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: arrayV !< The v-component array to be checksummed
  integer,                   optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,                   optional, intent(in) :: symmetric !< If true, do the checksums on the full symmetric computational domain.
  logical,                   optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                      optional, intent(in) :: scale     !< A scaling factor for these arrays.

  if (present(haloshift)) then
    call chksum_u_3d(arrayu, 'u '//mesg, hi, haloshift, symmetric, omit_corners, scale)
    call chksum_v_3d(arrayv, 'v '//mesg, hi, haloshift, symmetric, omit_corners, scale)
  else
    call chksum_u_3d(arrayu, 'u '//mesg, hi, symmetric=symmetric)
    call chksum_v_3d(arrayv, 'v '//mesg, hi, symmetric=symmetric)
  endif

end subroutine chksum_uv_3d

!> chksum_u_2d performs checksums on a 2d array staggered at C-grid u points.
subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scale)
  type(hor_index_type),           intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: array !< The array to be checksummed
  character(len=*),                intent(in) :: mesg  !< An identifying message
  integer,               optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,               optional, intent(in) :: symmetric !< If true, do the checksums on the full symmetric computational domain.
  logical,               optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                    optional, intent(in) :: scale     !< A scaling factor for this array.

  real :: scaling
  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
  integer :: bcN, bcS, bcE, bcW
  logical :: do_corners, sym, sym_stats

  if (checkfornans) then
    if (is_nan(array(hi%IscB:hi%IecB,hi%jsc:hi%jec))) &
      call chksum_error(fatal, 'NaN detected: '//trim(mesg))
!   if (is_NaN(array)) &
!     call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
  endif
  scaling = 1.0 ; if (present(scale)) scaling = scale

  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
  if (calculatestatistics) call substats(hi, array, mesg, sym_stats, scaling)

  if (.not.writechksums) return

  hshift = default_shift
  if (present(haloshift)) hshift = haloshift
  if (hshift<0) hshift = hi%iedB-hi%iecB

  if ( hi%iscB-hshift<hi%isdB .or. hi%iecB+hshift>hi%iedB .or. &
       hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
    write(0,*) 'chksum_u_2d: haloshift =',hshift
    write(0,*) 'chksum_u_2d: isd,isc,iec,ied=',hi%isdB,hi%iscB,hi%iecB,hi%iedB
    write(0,*) 'chksum_u_2d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
    call chksum_error(fatal,'Error in chksum_u_2d '//trim(mesg))
  endif

  bc0 = subchk(array, hi, 0, 0, scaling)

  sym = .false. ; if (present(symmetric)) sym = symmetric

  if ((hshift==0) .and. .not.sym) then
    if (is_root_pe()) call chk_sum_msg("u-point:",bc0,mesg)
    return
  endif

  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners

  if (hshift==0) then
    bcw = subchk(array, hi, -hshift-1, 0, scaling)
    if (is_root_pe()) call chk_sum_msg_w("u-point:",bc0,bcw,mesg)
  elseif (do_corners) then
    if (sym) then
      bcsw = subchk(array, hi, -hshift-1, -hshift, scaling)
      bcnw = subchk(array, hi, -hshift-1, hshift, scaling)
    else
      bcsw = subchk(array, hi, -hshift, -hshift, scaling)
      bcnw = subchk(array, hi, -hshift, hshift, scaling)
    endif
    bcse = subchk(array, hi, hshift, -hshift, scaling)
    bcne = subchk(array, hi, hshift, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg("u-point:",bc0,bcsw,bcse,bcnw,bcne,mesg)
  else
    bcs = subchk(array, hi, 0, -hshift, scaling)
    bce = subchk(array, hi, hshift, 0, scaling)
    if (sym) then
      bcw = subchk(array, hi, -hshift-1, 0, scaling)
    else
      bcw = subchk(array, hi, -hshift, 0, scaling)
    endif
    bcn = subchk(array, hi, 0, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg_nsew("u-point:",bc0,bcn,bcs,bce,bcw,mesg)
  endif

  contains

  integer function subchk(array, HI, di, dj, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: array
    integer, intent(in) :: di, dj
    real, intent(in) :: scale
    integer :: bitcount, i, j, bc
    subchk = 0
    ! This line deliberately uses the h-point computational domain.
    do j=hi%jsc+dj,hi%jec+dj; do i=hi%isc+di,hi%iec+di
      bc = bitcount(abs(scale*array(i,j)))
      subchk = subchk + bc
    enddo; enddo
    call sum_across_pes(subchk)
    subchk=mod(subchk,1000000000)
  end function subchk

  subroutine substats(HI, array, mesg, sym_stats, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: array
    character(len=*), intent(in) :: mesg
    logical, intent(in) :: sym_stats
    real, intent(in) :: scale
    integer :: i, j, n, IsB
    real :: aMean, aMin, aMax

    isb = hi%isc ; if (sym_stats) isb = hi%isc-1

    amin = array(hi%isc,hi%jsc) ; amax = amin
    do j=hi%jsc,hi%jec ; do i=isb,hi%IecB
      amin = min(amin, array(i,j))
      amax = max(amax, array(i,j))
    enddo ; enddo
    ! This line deliberately uses the h-point computational domain.
    amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec))
    n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc)
    call sum_across_pes(n)
    call min_across_pes(amin)
    call max_across_pes(amax)
    amean = amean / real(n)
    if (is_root_pe()) call chk_sum_msg("u-point:",amean*scale,amin*scale,amax*scale,mesg)
  end subroutine substats

end subroutine chksum_u_2d

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

!> chksum_v_2d performs checksums on a 2d array staggered at C-grid v points.
subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scale)
  type(hor_index_type),           intent(in) :: HI     !< A horizontal index type
  real, dimension(HI%isd:,HI%JsdB:), intent(in) :: array !< The array to be checksummed
  character(len=*),                intent(in) :: mesg  !< An identifying message
  integer,               optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,               optional, intent(in) :: symmetric !< If true, do the checksums on the full symmetric computational domain.
  logical,               optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                  optional, intent(in) :: scale     !< A scaling factor for this array.

  real :: scaling
  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
  integer :: bcN, bcS, bcE, bcW
  logical :: do_corners, sym, sym_stats

  if (checkfornans) then
    if (is_nan(array(hi%isc:hi%iec,hi%JscB:hi%JecB))) &
      call chksum_error(fatal, 'NaN detected: '//trim(mesg))
!   if (is_NaN(array)) &
!     call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
  endif
  scaling = 1.0 ; if (present(scale)) scaling = scale

  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
  if (calculatestatistics) call substats(hi, array, mesg, sym_stats, scaling)

  if (.not.writechksums) return

  hshift = default_shift
  if (present(haloshift)) hshift = haloshift
  if (hshift<0) hshift = hi%ied-hi%iec

  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
       hi%jscB-hshift<hi%jsdB .or. hi%jecB+hshift>hi%jedB ) then
    write(0,*) 'chksum_v_2d: haloshift =',hshift
    write(0,*) 'chksum_v_2d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
    write(0,*) 'chksum_v_2d: jsd,jsc,jec,jed=',hi%jsdB,hi%jscB,hi%jecB,hi%jedB
    call chksum_error(fatal,'Error in chksum_v_2d '//trim(mesg))
  endif

  bc0 = subchk(array, hi, 0, 0, scaling)

  sym = .false. ; if (present(symmetric)) sym = symmetric

  if ((hshift==0) .and. .not.sym) then
    if (is_root_pe()) call chk_sum_msg("v-point:",bc0,mesg)
    return
  endif

  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners

  if (hshift==0) then
    bcs = subchk(array, hi, 0, -hshift-1, scaling)
    if (is_root_pe()) call chk_sum_msg_s("v-point:",bc0,bcs,mesg)
  elseif (do_corners) then
    if (sym) then
      bcsw = subchk(array, hi, -hshift, -hshift-1, scaling)
      bcse = subchk(array, hi, hshift, -hshift-1, scaling)
    else
      bcsw = subchk(array, hi, -hshift, -hshift, scaling)
      bcse = subchk(array, hi, hshift, -hshift, scaling)
    endif
    bcnw = subchk(array, hi, -hshift, hshift, scaling)
    bcne = subchk(array, hi, hshift, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg("v-point:",bc0,bcsw,bcse,bcnw,bcne,mesg)
  else
    if (sym) then
      bcs = subchk(array, hi, 0, -hshift-1, scaling)
    else
      bcs = subchk(array, hi, 0, -hshift, scaling)
    endif
    bce = subchk(array, hi, hshift, 0, scaling)
    bcw = subchk(array, hi, -hshift, 0, scaling)
    bcn = subchk(array, hi, 0, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg_nsew("v-point:",bc0,bcn,bcs,bce,bcw,mesg)
  endif

  contains

  integer function subchk(array, HI, di, dj, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%isd:,HI%JsdB:), intent(in) :: array
    integer, intent(in) :: di, dj
    real, intent(in) :: scale
    integer :: bitcount, i, j, bc
    subchk = 0
    ! This line deliberately uses the h-point computational domain.
    do j=hi%jsc+dj,hi%jec+dj; do i=hi%isc+di,hi%iec+di
      bc = bitcount(abs(scale*array(i,j)))
      subchk = subchk + bc
    enddo; enddo
    call sum_across_pes(subchk)
    subchk=mod(subchk,1000000000)
  end function subchk

  subroutine substats(HI, array, mesg, sym_stats, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%isd:,HI%JsdB:), intent(in) :: array
    character(len=*), intent(in) :: mesg
    logical, intent(in) :: sym_stats
    real, intent(in) :: scale
    integer :: i, j, n, JsB
    real :: aMean, aMin, aMax

    jsb = hi%jsc ; if (sym_stats) jsb = hi%jsc-1

    amin = array(hi%isc,hi%jsc) ; amax = amin
    do j=jsb,hi%JecB ; do i=hi%isc,hi%iec
      amin = min(amin, array(i,j))
      amax = max(amax, array(i,j))
    enddo ; enddo
    ! This line deliberately uses the h-computational domain.
    amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec))
    n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc)
    call sum_across_pes(n)
    call min_across_pes(amin)
    call max_across_pes(amax)
    amean = amean / real(n)
    if (is_root_pe()) call chk_sum_msg("v-point:",amean*scale,amin*scale,amax*scale,mesg)
  end subroutine substats

end subroutine chksum_v_2d

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

!> chksum_h_3d performs checksums on a 3d array staggered at tracer points.
subroutine chksum_h_3d(array, mesg, HI, haloshift, omit_corners, scale)
  type(hor_index_type),             intent(in) :: HI !< A horizontal index type
  real, dimension(HI%isd:,HI%jsd:,:),  intent(in) :: array !< The array to be checksummed
  character(len=*),                  intent(in) :: mesg  !< An identifying message
  integer,                 optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,                 optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                    optional, intent(in) :: scale     !< A scaling factor for this array.

  real :: scaling
  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
  integer :: bcN, bcS, bcE, bcW
  logical :: do_corners

  if (checkfornans) then
    if (is_nan(array(hi%isc:hi%iec,hi%jsc:hi%jec,:))) &
      call chksum_error(fatal, 'NaN detected: '//trim(mesg))
!   if (is_NaN(array)) &
!     call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
  endif
  scaling = 1.0 ; if (present(scale)) scaling = scale

  if (calculatestatistics) call substats(hi, array, mesg, scaling)

  if (.not.writechksums) return

  hshift = default_shift
  if (present(haloshift)) hshift = haloshift
  if (hshift<0) hshift = hi%ied-hi%iec

  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
       hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
    write(0,*) 'chksum_h_3d: haloshift =',hshift
    write(0,*) 'chksum_h_3d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
    write(0,*) 'chksum_h_3d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
    call chksum_error(fatal,'Error in chksum_h_3d '//trim(mesg))
  endif

  bc0 = subchk(array, hi, 0, 0, scaling)

  if (hshift==0) then
    if (is_root_pe()) call chk_sum_msg("h-point:",bc0,mesg)
    return
  endif

  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners

  if (do_corners) then
    bcsw = subchk(array, hi, -hshift, -hshift, scaling)
    bcse = subchk(array, hi, hshift, -hshift, scaling)
    bcnw = subchk(array, hi, -hshift, hshift, scaling)
    bcne = subchk(array, hi, hshift, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg("h-point:",bc0,bcsw,bcse,bcnw,bcne,mesg)
  else
    bcs = subchk(array, hi, 0, -hshift, scaling)
    bce = subchk(array, hi, hshift, 0, scaling)
    bcw = subchk(array, hi, -hshift, 0, scaling)
    bcn = subchk(array, hi, 0, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg_nsew("h-point:",bc0,bcn,bcs,bce,bcw,mesg)
  endif

  contains

  integer function subchk(array, HI, di, dj, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: array
    integer, intent(in) :: di, dj
    real, intent(in) :: scale
    integer :: bitcount, i, j, k, bc
    subchk = 0
    do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc+dj,hi%jec+dj ; do i=hi%isc+di,hi%iec+di
      bc = bitcount(abs(scale*array(i,j,k)))
      subchk = subchk + bc
    enddo ; enddo ; enddo
    call sum_across_pes(subchk)
    subchk=mod(subchk,1000000000)
  end function subchk

  subroutine substats(HI, array, mesg, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: array
    character(len=*), intent(in) :: mesg
    real, intent(in) :: scale
    integer :: i, j, k, n
    real :: aMean, aMin, aMax

    amin = array(hi%isc,hi%jsc,1)
    amax = array(hi%isc,hi%jsc,1)
    n = 0
    do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
      amin = min(amin, array(i,j,k))
      amax = max(amax, array(i,j,k))
      n = n + 1
    enddo ; enddo ; enddo
    amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec,:))
    call sum_across_pes(n)
    call min_across_pes(amin)
    call max_across_pes(amax)
    amean = amean / real(n)
    if (is_root_pe()) call chk_sum_msg("h-point:",amean*scale,amin*scale,amax*scale,mesg)
  end subroutine substats

end subroutine chksum_h_3d

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

!> chksum_B_3d performs checksums on a 3d array staggered at corner points.
subroutine chksum_b_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scale)
  type(hor_index_type),              intent(in) :: HI !< A horizontal index type
  real, dimension(HI%IsdB:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed
  character(len=*),                   intent(in) :: mesg  !< An identifying message
  integer,                  optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,                  optional, intent(in) :: symmetric !< If true, do the checksums on the full symmetric computational domain.
  logical,                  optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                     optional, intent(in) :: scale     !< A scaling factor for this array.

  real :: scaling
  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
  integer :: bcN, bcS, bcE, bcW
  logical :: do_corners, sym, sym_stats

  if (checkfornans) then
    if (is_nan(array(hi%IscB:hi%IecB,hi%JscB:hi%JecB,:))) &
      call chksum_error(fatal, 'NaN detected: '//trim(mesg))
!   if (is_NaN(array)) &
!     call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
  endif
  scaling = 1.0 ; if (present(scale)) scaling = scale

  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
  if (calculatestatistics) call substats(hi, array, mesg, sym_stats, scaling)

  if (.not.writechksums) return

  hshift = default_shift
  if (present(haloshift)) hshift = haloshift
  if (hshift<0) hshift = hi%ied-hi%iec

  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
       hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
    write(0,*) 'chksum_B_3d: haloshift =',hshift
    write(0,*) 'chksum_B_3d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
    write(0,*) 'chksum_B_3d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
    call chksum_error(fatal,'Error in chksum_B_3d '//trim(mesg))
  endif

  bc0 = subchk(array, hi, 0, 0, scaling)

  sym = .false. ; if (present(symmetric)) sym = symmetric

  if ((hshift==0) .and. .not.sym) then
    if (is_root_pe()) call chk_sum_msg("B-point:",bc0,mesg)
    return
  endif

  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners

  if (do_corners) then
    if (sym) then
      bcsw = subchk(array, hi, -hshift-1, -hshift-1, scaling)
      bcse = subchk(array, hi, hshift, -hshift-1, scaling)
      bcnw = subchk(array, hi, -hshift-1, hshift, scaling)
    else
      bcsw = subchk(array, hi, -hshift-1, -hshift-1, scaling)
      bcse = subchk(array, hi, hshift, -hshift-1, scaling)
      bcnw = subchk(array, hi, -hshift-1, hshift, scaling)
    endif
    bcne = subchk(array, hi, hshift, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg("B-point:",bc0,bcsw,bcse,bcnw,bcne,mesg)
  else
    if (sym) then
      bcs = subchk(array, hi, 0, -hshift-1, scaling)
      bcw = subchk(array, hi, -hshift-1, 0, scaling)
    else
      bcs = subchk(array, hi, 0, -hshift, scaling)
      bcw = subchk(array, hi, -hshift, 0, scaling)
    endif
    bce = subchk(array, hi, hshift, 0, scaling)
    bcn = subchk(array, hi, 0, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg_nsew("B-point:",bc0,bcn,bcs,bce,bcw,mesg)
  endif

  contains

  integer function subchk(array, HI, di, dj, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%IsdB:,HI%JsdB:,:), intent(in) :: array
    integer, intent(in) :: di, dj
    real, intent(in) :: scale
    integer :: bitcount, i, j, k, bc
    subchk = 0
    ! This line deliberately uses the h-point computational domain.
    do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc+dj,hi%jec+dj ; do i=hi%isc+di,hi%iec+di
      bc = bitcount(abs(scale*array(i,j,k)))
      subchk = subchk + bc
    enddo ; enddo ; enddo
    call sum_across_pes(subchk)
    subchk=mod(subchk,1000000000)
  end function subchk

  subroutine substats(HI, array, mesg, sym_stats, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%IsdB:,HI%JsdB:,:), intent(in) :: array
    character(len=*), intent(in) :: mesg
    logical, intent(in) :: sym_stats
    real, intent(in) :: scale
    integer :: i, j, k, n, IsB, JsB
    real :: aMean, aMin, aMax

    isb = hi%isc ; if (sym_stats) isb = hi%isc-1
    jsb = hi%jsc ; if (sym_stats) jsb = hi%jsc-1

    amin = array(hi%isc,hi%jsc,1) ; amax = amin
    do k=lbound(array,3),ubound(array,3) ; do j=jsb,hi%JecB ; do i=isb,hi%IecB
      amin = min(amin, array(i,j,k))
      amax = max(amax, array(i,j,k))
    enddo ; enddo ; enddo
    amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec,:))
    n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc) * size(array,3)
    call sum_across_pes(n)
    call min_across_pes(amin)
    call max_across_pes(amax)
    amean = amean / real(n)
    if (is_root_pe()) call chk_sum_msg("B-point:",amean*scale,amin*scale,amax*scale,mesg)
  end subroutine substats

end subroutine chksum_b_3d

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

!> chksum_u_3d performs checksums on a 3d array staggered at C-grid u points.
subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scale)
  type(hor_index_type),             intent(in) :: HI !< A horizontal index type
  real, dimension(HI%isdB:,HI%Jsd:,:), intent(in) :: array !< The array to be checksummed
  character(len=*),                  intent(in) :: mesg  !< An identifying message
  integer,                 optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,                 optional, intent(in) :: symmetric !< If true, do the checksums on the full symmetric computational domain.
  logical,                 optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                    optional, intent(in) :: scale     !< A scaling factor for this array.

  real :: scaling
  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
  integer :: bcN, bcS, bcE, bcW
  logical :: do_corners, sym, sym_stats

  if (checkfornans) then
    if (is_nan(array(hi%IscB:hi%IecB,hi%jsc:hi%jec,:))) &
      call chksum_error(fatal, 'NaN detected: '//trim(mesg))
!   if (is_NaN(array)) &
!     call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
  endif
  scaling = 1.0 ; if (present(scale)) scaling = scale

  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
  if (calculatestatistics) call substats(hi, array, mesg, sym_stats, scaling)

  if (.not.writechksums) return

  hshift = default_shift
  if (present(haloshift)) hshift = haloshift
  if (hshift<0) hshift = hi%ied-hi%iec

  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
       hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
    write(0,*) 'chksum_u_3d: haloshift =',hshift
    write(0,*) 'chksum_u_3d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
    write(0,*) 'chksum_u_3d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
    call chksum_error(fatal,'Error in chksum_u_3d '//trim(mesg))
  endif

  bc0 = subchk(array, hi, 0, 0, scaling)

  sym = .false. ; if (present(symmetric)) sym = symmetric

  if ((hshift==0) .and. .not.sym) then
    if (is_root_pe()) call chk_sum_msg("u-point:",bc0,mesg)
    return
  endif

  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners

  if (hshift==0) then
    bcw = subchk(array, hi, -hshift-1, 0, scaling)
    if (is_root_pe()) call chk_sum_msg_w("u-point:",bc0,bcw,mesg)
  elseif (do_corners) then
    if (sym) then
      bcsw = subchk(array, hi, -hshift-1, -hshift, scaling)
      bcnw = subchk(array, hi, -hshift-1, hshift, scaling)
    else
      bcsw = subchk(array, hi, -hshift, -hshift, scaling)
      bcnw = subchk(array, hi, -hshift, hshift, scaling)
    endif
    bcse = subchk(array, hi, hshift, -hshift, scaling)
    bcne = subchk(array, hi, hshift, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg("u-point:",bc0,bcsw,bcse,bcnw,bcne,mesg)
  else
    bcs = subchk(array, hi, 0, -hshift, scaling)
    bce = subchk(array, hi, hshift, 0, scaling)
    if (sym) then
      bcw = subchk(array, hi, -hshift-1, 0, scaling)
    else
      bcw = subchk(array, hi, -hshift, 0, scaling)
    endif
    bcn = subchk(array, hi, 0, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg_nsew("u-point:",bc0,bcn,bcs,bce,bcw,mesg)
  endif

  contains

  integer function subchk(array, HI, di, dj, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%IsdB:,HI%jsd:,:), intent(in) :: array
    integer, intent(in) :: di, dj
    real, intent(in) :: scale
    integer :: bitcount, i, j, k, bc
    subchk = 0
    ! This line deliberately uses the h-point computational domain.
    do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc+dj,hi%jec+dj ; do i=hi%isc+di,hi%iec+di
      bc = bitcount(abs(scale*array(i,j,k)))
      subchk = subchk + bc
    enddo ; enddo ; enddo
    call sum_across_pes(subchk)
    subchk=mod(subchk,1000000000)
  end function subchk

  subroutine substats(HI, array, mesg, sym_stats, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%IsdB:,HI%jsd:,:), intent(in) :: array
    character(len=*), intent(in) :: mesg
    logical, intent(in) :: sym_stats
    real, intent(in) :: scale
    integer :: i, j, k, n, IsB
    real :: aMean, aMin, aMax

    isb = hi%isc ; if (sym_stats) isb = hi%isc-1

    amin = array(hi%isc,hi%jsc,1) ; amax = amin
    do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc,hi%jec ; do i=isb,hi%IecB
      amin = min(amin, array(i,j,k))
      amax = max(amax, array(i,j,k))
    enddo ; enddo ; enddo
    ! This line deliberately uses the h-point computational domain.
    amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec,:))
    n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc) * size(array,3)
    call sum_across_pes(n)
    call min_across_pes(amin)
    call max_across_pes(amax)
    amean = amean / real(n)
    if (is_root_pe()) call chk_sum_msg("u-point:",amean*scale,amin*scale,amax*scale,mesg)
  end subroutine substats

end subroutine chksum_u_3d

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

!> chksum_v_3d performs checksums on a 3d array staggered at C-grid v points.
subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scale)
  type(hor_index_type),             intent(in) :: HI !< A horizontal index type
  real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed
  character(len=*),                  intent(in) :: mesg  !< An identifying message
  integer,                 optional, intent(in) :: haloshift !< The width of halos to check (default 0)
  logical,                 optional, intent(in) :: symmetric !< If true, do the checksums on the full symmetric computational domain.
  logical,                 optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
  real,                    optional, intent(in) :: scale     !< A scaling factor for this array.

  real :: scaling
  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
  integer :: bcN, bcS, bcE, bcW
  logical :: do_corners, sym, sym_stats

  if (checkfornans) then
    if (is_nan(array(hi%isc:hi%iec,hi%JscB:hi%JecB,:))) &
      call chksum_error(fatal, 'NaN detected: '//trim(mesg))
!   if (is_NaN(array)) &
!     call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
  endif
  scaling = 1.0 ; if (present(scale)) scaling = scale

  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
  if (calculatestatistics) call substats(hi, array, mesg, sym_stats, scaling)

  if (.not.writechksums) return

  hshift = default_shift
  if (present(haloshift)) hshift = haloshift
  if (hshift<0) hshift = hi%ied-hi%iec

  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
       hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
    write(0,*) 'chksum_v_3d: haloshift =',hshift
    write(0,*) 'chksum_v_3d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
    write(0,*) 'chksum_v_3d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
    call chksum_error(fatal,'Error in chksum_v_3d '//trim(mesg))
  endif

  bc0 = subchk(array, hi, 0, 0, scaling)

  sym = .false. ; if (present(symmetric)) sym = symmetric

  if ((hshift==0) .and. .not.sym) then
    if (is_root_pe()) call chk_sum_msg("v-point:",bc0,mesg)
    return
  endif

  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners

  if (hshift==0) then
    bcs = subchk(array, hi, 0, -hshift-1, scaling)
    if (is_root_pe()) call chk_sum_msg_s("v-point:",bc0,bcs,mesg)
  elseif (do_corners) then
    if (sym) then
      bcsw = subchk(array, hi, -hshift, -hshift-1, scaling)
      bcse = subchk(array, hi, hshift, -hshift-1, scaling)
    else
      bcsw = subchk(array, hi, -hshift, -hshift, scaling)
      bcse = subchk(array, hi, hshift, -hshift, scaling)
    endif
    bcnw = subchk(array, hi, -hshift, hshift, scaling)
    bcne = subchk(array, hi, hshift, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg("v-point:",bc0,bcsw,bcse,bcnw,bcne,mesg)
  else
    if (sym) then
      bcs = subchk(array, hi, 0, -hshift-1, scaling)
    else
      bcs = subchk(array, hi, 0, -hshift, scaling)
    endif
    bce = subchk(array, hi, hshift, 0, scaling)
    bcw = subchk(array, hi, -hshift, 0, scaling)
    bcn = subchk(array, hi, 0, hshift, scaling)

    if (is_root_pe()) call chk_sum_msg_nsew("v-point:",bc0,bcn,bcs,bce,bcw,mesg)
  endif

  contains

  integer function subchk(array, HI, di, dj, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: array
    integer, intent(in) :: di, dj
    real, intent(in) :: scale
    integer :: bitcount, i, j, k, bc
    subchk = 0
    ! This line deliberately uses the h-point computational domain.
    do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc+dj,hi%jec+dj ; do i=hi%isc+di,hi%iec+di
      bc = bitcount(abs(scale*array(i,j,k)))
      subchk = subchk + bc
    enddo ; enddo ; enddo
    call sum_across_pes(subchk)
    subchk=mod(subchk,1000000000)
  end function subchk

  subroutine substats(HI, array, mesg, sym_stats, scale)
    type(hor_index_type), intent(in) :: HI
    real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: array
    character(len=*), intent(in) :: mesg
    logical, intent(in) :: sym_stats
    real, intent(in) :: scale
    integer :: i, j, k, n, JsB
    real :: aMean, aMin, aMax

    jsb = hi%jsc ; if (sym_stats) jsb = hi%jsc-1

    amin = array(hi%isc,hi%jsc,1) ; amax = amin
    do k=lbound(array,3),ubound(array,3) ; do j=jsb,hi%JecB ; do i=hi%isc,hi%iec
      amin = min(amin, array(i,j,k))
      amax = max(amax, array(i,j,k))
    enddo ; enddo ; enddo
    ! This line deliberately uses the h-point computational domain.
    amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec,:))
    n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc) * size(array,3)
    call sum_across_pes(n)
    call min_across_pes(amin)
    call max_across_pes(amax)
    amean = amean / real(n)
    if (is_root_pe()) call chk_sum_msg("v-point:",amean*scale,amin*scale,amax*scale,mesg)
  end subroutine substats

end subroutine chksum_v_3d


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

!   These are the older version of chksum that do not take the grid staggering
! into account.

!> chksum1d does a checksum of a 1-dimensional array.
subroutine chksum1d(array, mesg, start_i, end_i, compare_PEs)
  real, dimension(:), intent(in) :: array   !< The array to be summed (index starts at 1).
  character(len=*),   intent(in) :: mesg    !< An identifying message.
  integer, optional,  intent(in) :: start_i !< The starting index for the sum (default 1)
  integer, optional,  intent(in) :: end_i   !< The ending index for the sum (default all)
  logical, optional,  intent(in) :: compare_PEs !< If true, compare across PEs instead of summing
                                                !! and list the root_PE value (default true)

  integer :: is, ie, i, bc, sum1, sum_bc
  integer :: bitcount
  real :: sum
  real, allocatable :: sum_here(:)
  logical :: compare
  integer :: pe_num   ! pe number of the data
  integer :: nPEs     ! Total number of processsors

  is = lbound(array,1) ; ie = ubound(array,1)
  if (present(start_i)) is = start_i
  if (present(end_i)) ie = end_i
  compare = .true. ; if (present(compare_pes)) compare = compare_pes

  sum = 0.0 ; sum_bc = 0
  do i=is,ie
    sum = sum + array(i)
    bc = bitcount(abs(array(i)))
    sum_bc = sum_bc + bc
  enddo

  pe_num = pe_here() + 1 - root_pe() ; npes = num_pes()
  allocate(sum_here(npes)) ; sum_here(:) = 0.0 ; sum_here(pe_num) = sum
  call sum_across_pes(sum_here,npes)

  sum1 = sum_bc
  call sum_across_pes(sum1)

  if (.not.compare) then
    sum = 0.0
    do i=1,npes ; sum = sum + sum_here(i) ; enddo
    sum_bc = sum1
  elseif (is_root_pe()) then
    if (sum1 /= npes*sum_bc) &
      write(0, '(A40," bitcounts do not match across PEs: ",I12,1X,I12)') &
            mesg, sum1, npes*sum_bc
    do i=1,npes ; if (sum /= sum_here(i)) then
      write(0, '(A40," PE ",i4," sum mismatches root_PE: ",3(ES22.13,1X))') &
            mesg, i, sum_here(i), sum, sum_here(i)-sum
    endif ; enddo
  endif
  deallocate(sum_here)

  if (is_root_pe()) &
    write(0,'(A50,1X,ES25.16,1X,I12)') mesg, sum, sum_bc

end subroutine chksum1d

! =====================================================================
!   These are the older version of chksum that do not take the grid staggering
! into account.

!> chksum2d does a checksum of all data in a 2-d array.
subroutine chksum2d(array, mesg)

  real, dimension(:,:) :: array
  character(len=*) :: mesg

  integer :: bitcount
  integer :: xs,xe,ys,ye,i,j,sum1,bc
  real :: sum

  xs = lbound(array,1) ; xe = ubound(array,1)
  ys = lbound(array,2) ; ye = ubound(array,2)

  sum = 0.0 ; sum1 = 0
  do i=xs,xe ; do j=ys,ye
    bc = bitcount(abs(array(i,j)))
    sum1 = sum1 + bc
  enddo ; enddo
  call sum_across_pes(sum1)

  sum = reproducing_sum(array(:,:))

  if (is_root_pe()) &
    write(0,'(A50,1X,ES25.16,1X,I12)') mesg, sum, sum1
!    write(0,'(A40,1X,Z16.16,1X,Z16.16,1X,ES25.16,1X,I12)') &
!      mesg, sum, sum1, sum, sum1

end subroutine chksum2d

!> chksum3d does a checksum of all data in a 2-d array.
subroutine chksum3d(array, mesg)

  real, dimension(:,:,:) :: array
  character(len=*) :: mesg

  integer :: bitcount
  integer :: xs,xe,ys,ye,zs,ze,i,j,k, bc,sum1
  real :: sum

  xs = lbound(array,1) ; xe = ubound(array,1)
  ys = lbound(array,2) ; ye = ubound(array,2)
  zs = lbound(array,3) ; ze = ubound(array,3)

  sum = 0.0 ; sum1 = 0
  do i=xs,xe ; do j=ys,ye ; do k=zs,ze
    bc = bitcount(abs(array(i,j,k)))
    sum1 = sum1 + bc
  enddo ; enddo ; enddo

  call sum_across_pes(sum1)
  sum = reproducing_sum(array(:,:,:))

  if (is_root_pe()) &
    write(0,'(A50,1X,ES25.16,1X,I12)') mesg, sum, sum1
!    write(0,'(A40,1X,Z16.16,1X,Z16.16,1X,ES25.16,1X,I12)') &
!      mesg, sum, sum1, sum, sum1

end subroutine chksum3d

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

!> This function returns .true. if x is a NaN, and .false. otherwise.
function is_nan_0d(x)
  real, intent(in) :: x !< The value to be checked for NaNs.
  logical :: is_NaN_0d

 !is_NaN_0d = (((x < 0.0) .and. (x >= 0.0)) .or. &
 !          (.not.(x < 0.0) .and. .not.(x >= 0.0)))
  if (((x < 0.0) .and. (x >= 0.0)) .or. &
            (.not.(x < 0.0) .and. .not.(x >= 0.0))) then
    is_nan_0d = .true.
  else
    is_nan_0d = .false.
  endif

end function is_nan_0d

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

!> This function returns .true. if any element of x is a NaN, and .false. otherwise.
function is_nan_1d(x, skip_mpp)
  real, dimension(:), intent(in) :: x !< The array to be checked for NaNs.
  logical :: is_NaN_1d
  logical, optional :: skip_mpp  !< If true, only check this array only on the local PE (default false).

  integer :: i, n
  logical :: call_mpp

  n = 0
  do i = lbound(x,1), ubound(x,1)
    if (is_nan_0d(x(i))) n = n + 1
  enddo
  call_mpp = .true.
  if (present(skip_mpp)) call_mpp = .not.skip_mpp

  if (call_mpp) call sum_across_pes(n)
  is_nan_1d = .false.
  if (n>0) is_nan_1d = .true.

end function is_nan_1d

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

!> This function returns .true. if any element of x is a NaN, and .false. otherwise.
function is_nan_2d(x)
  real, dimension(:,:), intent(in) :: x !< The array to be checked for NaNs.
  logical :: is_NaN_2d

  integer :: i, j, n

  n = 0
  do j = lbound(x,2), ubound(x,2) ; do i = lbound(x,1), ubound(x,1)
    if (is_nan_0d(x(i,j))) n = n + 1
  enddo ; enddo
  call sum_across_pes(n)
  is_nan_2d = .false.
  if (n>0) is_nan_2d = .true.

end function is_nan_2d

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

!> This function returns .true. if any element of x is a NaN, and .false. otherwise.
function is_nan_3d(x)
  real, dimension(:,:,:), intent(in) :: x !< The array to be checked for NaNs.
  logical :: is_NaN_3d

  integer :: i, j, k, n

  n = 0
  do k = lbound(x,3), ubound(x,3)
    do j = lbound(x,2), ubound(x,2) ; do i = lbound(x,1), ubound(x,1)
      if (is_nan_0d(x(i,j,k))) n = n + 1
    enddo ; enddo
  enddo
  call sum_across_pes(n)
  is_nan_3d = .false.
  if (n>0) is_nan_3d = .true.

end function is_nan_3d

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

subroutine chk_sum_msg1(fmsg,bc0,mesg)
  character(len=*), intent(in) :: fmsg, mesg
  integer,          intent(in) :: bc0
  if (is_root_pe()) write(0,'(A,1(A,I10,X),A)') fmsg," c=",bc0,trim(mesg)
end subroutine chk_sum_msg1

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

subroutine chk_sum_msg5(fmsg,bc0,bcSW,bcSE,bcNW,bcNE,mesg)
  character(len=*), intent(in) :: fmsg, mesg
  integer,          intent(in) :: bc0,bcSW,bcSE,bcNW,bcNE
  if (is_root_pe()) write(0,'(A,5(A,I10,1X),A)') &
     fmsg," c=",bc0,"sw=",bcsw,"se=",bcse,"nw=",bcnw,"ne=",bcne,trim(mesg)
end subroutine chk_sum_msg5

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

subroutine chk_sum_msg_nsew(fmsg,bc0,bcN,bcS,bcE,bcW,mesg)
  character(len=*), intent(in) :: fmsg, mesg
  integer,          intent(in) :: bc0, bcN, bcS, bcE, bcW
  if (is_root_pe()) write(0,'(A,5(A,I10,1X),A)') &
     fmsg," c=",bc0,"N=",bcn,"S=",bcs,"E=",bce,"W=",bcw,trim(mesg)
end subroutine chk_sum_msg_nsew

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

subroutine chk_sum_msg_s(fmsg,bc0,bcS,mesg)
  character(len=*), intent(in) :: fmsg, mesg
  integer,          intent(in) :: bc0, bcS
  if (is_root_pe()) write(0,'(A,2(A,I10,1X),A)') &
     fmsg," c=",bc0,"S=",bcs,trim(mesg)
end subroutine chk_sum_msg_s

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

subroutine chk_sum_msg_w(fmsg,bc0,bcW,mesg)
  character(len=*), intent(in) :: fmsg, mesg
  integer,          intent(in) :: bc0, bcW
  if (is_root_pe()) write(0,'(A,2(A,I10,1X),A)') &
     fmsg," c=",bc0,"W=",bcw,trim(mesg)
end subroutine chk_sum_msg_w

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

subroutine chk_sum_msg2(fmsg,bc0,bcSW,mesg)
  character(len=*), intent(in) :: fmsg, mesg
  integer,          intent(in) :: bc0,bcSW
  if (is_root_pe()) write(0,'(A,2(A,I9,1X),A)') &
     fmsg," c=",bc0,"s/w=",bcsw,trim(mesg)
end subroutine chk_sum_msg2

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

subroutine chk_sum_msg3(fmsg,aMean,aMin,aMax,mesg)
  character(len=*), intent(in) :: fmsg, mesg
  real,             intent(in) :: aMean,aMin,aMax
  if (is_root_pe()) write(0,'(A,3(A,ES25.16,1X),A)') &
     fmsg," mean=",amean,"min=",amin,"max=",amax,trim(mesg)
end subroutine chk_sum_msg3

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

!> MOM_checksums_init initializes the MOM_checksums module. As it happens, the
!! only thing that it does is to log the version of this module.
subroutine mom_checksums_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_checksums" ! This module's name.

  call log_version(param_file, mdl, version)

end subroutine mom_checksums_init

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

subroutine chksum_error(signal, message)
  ! Wrapper for MOM_error to help place specific break points in
  ! debuggers
  integer, intent(in) :: signal
  character(len=*), intent(in) :: message
  call mom_error(signal, message)
end subroutine chksum_error

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

end module mom_checksums