MOM_checksums.F90

Go to the documentation of this file.

!> Routines to calculate checksums of various array and vector types
module mom_checksums
 
! This file is part of MOM6. See LICENSE.md for the license.
 
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
 
use iso_fortran_env, only: error_unit
 
implicit none ; private
 
public :: chksum0, zchksum
public :: hchksum, bchksum, uchksum, vchksum, qchksum, is_nan, chksum
public :: hchksum_pair, uvchksum, bchksum_pair
public :: mom_checksums_init
 
!> Checksums a pair of arrays (2d or 3d) staggered at tracer points
interface hchksum_pair
  module procedure chksum_pair_h_2d, chksum_pair_h_3d
end interface
 
!> Checksums a pair velocity arrays (2d or 3d) staggered at C-grid locations
interface uvchksum
  module procedure chksum_uv_2d, chksum_uv_3d
end interface
 
!> Checksums an array (2d or 3d) staggered at C-grid u points.
interface uchksum
  module procedure chksum_u_2d, chksum_u_3d
end interface
 
!> Checksums an array (2d or 3d) staggered at C-grid v points.
interface vchksum
  module procedure chksum_v_2d, chksum_v_3d
end interface
 
!> Checksums a pair of arrays (2d or 3d) staggered at corner points
interface bchksum_pair
  module procedure chksum_pair_b_2d, chksum_pair_b_3d
end interface
 
!> Checksums an array (2d or 3d) staggered at tracer points.
interface hchksum
  module procedure chksum_h_2d, chksum_h_3d
end interface
 
!> Checksums an array (2d or 3d) staggered at corner points.
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
 
!> This is an older interface for 1-, 2-, or 3-D checksums
interface chksum
  module procedure chksum1d, chksum2d, chksum3d
end interface
 
!> Write a message with either checksums or numerical statistics of arrays
interface chk_sum_msg
  module procedure chk_sum_msg1, chk_sum_msg2, chk_sum_msg3, chk_sum_msg5
end interface
 
!> Returns .true. if any element of x is a NaN, and .false. otherwise.
interface is_nan
  module procedure is_nan_0d, is_nan_1d, is_nan_2d, is_nan_3d
end interface
 
integer, parameter :: bc_modulus = 1000000000 !< Modulus of checksum bitcount
integer, parameter :: default_shift=0 !< The default array shift
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
 
!> Checksum a scalar field (consistent with array checksums)
subroutine chksum0(scalar, mesg, scale, logunit)
  real, intent(in) :: scalar                !< The array to be checksummed
  character(len=*), intent(in) :: mesg     !< An identifying message
  real, optional, intent(in) :: scale      !< A scaling factor for this array.
  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
 
  real :: scaling   !< Explicit rescaling factor
  integer :: iounit !< Log IO unit
  real :: rs        !< Rescaled scalar
  integer :: bc     !< Scalar bitcount
 
  if (checkfornans .and. is_nan(scalar)) &
    call chksum_error(fatal, 'NaN detected: '//trim(mesg))
 
  scaling = 1.0 ; if (present(scale)) scaling = scale
  iounit = error_unit; if(present(logunit)) iounit = logunit
 
  if (calculatestatistics) then
    rs = scaling * scalar
    if (is_root_pe()) &
      call chk_sum_msg(" scalar:", rs, rs, rs, mesg, iounit)
  endif
 
  if (.not. writechksums) return
 
  bc = mod(bitcount(abs(scaling * scalar)), bc_modulus)
  if (is_root_pe()) &
    call chk_sum_msg(" scalar:", bc, mesg, iounit)
 
end subroutine chksum0
 
 
!> Checksum a 1d array (typically a column).
subroutine zchksum(array, mesg, scale, logunit)
  real, dimension(:), intent(in) :: array  !< The array to be checksummed
  character(len=*), intent(in) :: mesg     !< An identifying message
  real, optional, intent(in) :: scale      !< A scaling factor for this array.
  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
 
  real, allocatable, dimension(:) :: rescaled_array
  real :: scaling
  integer :: iounit !< Log IO unit
  integer :: k
  real :: amean, amin, amax
  integer :: bc0
 
  if (checkfornans) then
    if (is_nan(array(:))) &
      call chksum_error(fatal, 'NaN detected: '//trim(mesg))
  endif
 
  scaling = 1.0 ; if (present(scale)) scaling = scale
  iounit = error_unit; if(present(logunit)) iounit = logunit
 
  if (calculatestatistics) then
    if (present(scale)) then
      allocate(rescaled_array(lbound(array,1):ubound(array,1)))
      rescaled_array(:) = 0.0
      do k=1, size(array, 1)
        rescaled_array(k) = scale * array(k)
      enddo
 
      call substats(rescaled_array, amean, amin, amax)
      deallocate(rescaled_array)
    else
      call substats(array, amean, amin, amax)
    endif
 
    if (is_root_pe()) &
      call chk_sum_msg(" column:", amean, amin, amax, mesg, iounit)
  endif
 
  if (.not. writechksums) return
 
  bc0 = subchk(array, scaling)
  if (is_root_pe()) call chk_sum_msg(" column:", bc0, mesg, iounit)
 
  contains
 
  integer function subchk(array, scale)
    real, dimension(:), intent(in) :: array !< The array to be checksummed
    real, intent(in) :: scale !< A scaling factor for this array.
    integer :: k, bc
    subchk = 0
    do k=lbound(array, 1), ubound(array, 1)
      bc = bitcount(abs(scale * array(k)))
      subchk = subchk + bc
    enddo
    subchk=mod(subchk, bc_modulus)
  end function subchk
 
  subroutine substats(array, aMean, aMin, aMax)
    real, dimension(:), intent(in) :: array !< The array to be checksummed
    real, intent(out) :: aMean, aMin, aMax
 
    integer :: k, n
 
    amin = array(1)
    amax = array(1)
    n = 0
    do k=lbound(array,1), ubound(array,1)
      amin = min(amin, array(k))
      amax = max(amax, array(k))
      n = n + 1
    enddo
    amean = sum(array(:)) / real(n)
  end subroutine substats
 
end subroutine zchksum
 
!> Checksums on a pair of 2d arrays staggered at tracer points.
subroutine chksum_pair_h_2d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, &
                            scale, logunit)
  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 !< The first array to be checksummed
  real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayB !< The second array 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.
  integer,                optional, intent(in) :: logunit !< IO unit for checksum logging
 
  if (present(haloshift)) then
    call chksum_h_2d(arraya, 'x '//mesg, hi, haloshift, omit_corners, &
                     scale=scale, logunit=logunit)
    call chksum_h_2d(arrayb, 'y '//mesg, hi, haloshift, omit_corners, &
                     scale=scale, logunit=logunit)
  else
    call chksum_h_2d(arraya, 'x '//mesg, hi, scale=scale, logunit=logunit)
    call chksum_h_2d(arrayb, 'y '//mesg, hi, scale=scale, logunit=logunit)
  endif
 
end subroutine chksum_pair_h_2d
 
!> Checksums on a pair of 3d arrays staggered at tracer points.
subroutine chksum_pair_h_3d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, &
                            scale, logunit)
  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 !< The first array to be checksummed
  real, dimension(HI%isd:,HI%jsd:, :), intent(in) :: arrayB !< The second array 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.
  integer,                   optional, intent(in) :: logunit !< IO unit for checksum logging
 
  if (present(haloshift)) then
    call chksum_h_3d(arraya, 'x '//mesg, hi, haloshift, omit_corners, &
                     scale=scale, logunit=logunit)
    call chksum_h_3d(arrayb, 'y '//mesg, hi, haloshift, omit_corners, &
                     scale=scale, logunit=logunit)
  else
    call chksum_h_3d(arraya, 'x '//mesg, hi, scale=scale, logunit=logunit)
    call chksum_h_3d(arrayb, 'y '//mesg, hi, scale=scale, logunit=logunit)
  endif
 
end subroutine chksum_pair_h_3d
 
!> Checksums a 2d array staggered at tracer points.
subroutine chksum_h_2d(array, mesg, HI, haloshift, omit_corners, scale, logunit)
  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.
  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
 
  real, allocatable, dimension(:,:) :: rescaled_array
  real :: scaling
  integer :: iounit !< Log IO unit
  integer :: i, j
  real :: aMean, aMin, aMax
  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
  iounit = error_unit; if(present(logunit)) iounit = logunit
 
  if (calculatestatistics) then
    if (present(scale)) then
      allocate( rescaled_array(lbound(array,1):ubound(array,1), &
                               lbound(array,2):ubound(array,2)) )
      rescaled_array(:,:) = 0.0
      do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
        rescaled_array(i,j) = scale*array(i,j)
      enddo ; enddo
      call substats(hi, rescaled_array, amean, amin, amax)
      deallocate(rescaled_array)
    else
      call substats(hi, array, amean, amin, amax)
    endif
 
    if (is_root_pe()) &
      call chk_sum_msg("h-point:", amean, amin, amax, mesg, iounit)
  endif
 
  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, iounit)
    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, iounit)
  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, iounit)
  endif
 
  contains
  integer function subchk(array, HI, di, dj, 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
    integer, intent(in) :: di    !< i- direction array shift for this checksum
    integer, intent(in) :: dj    !< j- direction array shift for this checksum
    real, intent(in)    :: scale !< A scaling factor for this array.
    integer :: 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, bc_modulus)
  end function subchk
 
  subroutine substats(HI, array, aMean, aMin, aMax)
    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
    real, intent(out) :: aMean, aMin, aMax
 
    integer :: i, j, n
 
    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)
  end subroutine substats
 
end subroutine chksum_h_2d
 
!> Checksums on a pair of 2d arrays staggered at q-points.
subroutine chksum_pair_b_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, &
                            omit_corners, scale, logunit)
  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 !< The first array to be checksummed
  real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayB !< The second array 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.
  integer,                optional, intent(in) :: logunit !< IO unit for checksum logging
 
  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, logunit=logunit)
    call chksum_b_2d(arrayb, 'y '//mesg, hi, haloshift, symmetric=sym, &
                     omit_corners=omit_corners, scale=scale, logunit=logunit)
  else
    call chksum_b_2d(arraya, 'x '//mesg, hi, symmetric=sym, scale=scale, &
                     logunit=logunit)
    call chksum_b_2d(arrayb, 'y '//mesg, hi, symmetric=sym, scale=scale, &
                     logunit=logunit)
  endif
 
end subroutine chksum_pair_b_2d
 
!> Checksums on a pair of 3d arrays staggered at q-points.
subroutine chksum_pair_b_3d(mesg, arrayA, arrayB, HI, haloshift, symmetric, &
                            omit_corners, scale, logunit)
  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 !< The first array to be checksummed
  real, dimension(HI%IsdB:,HI%JsdB:, :), intent(in) :: arrayB !< The second 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 this array.
  integer,                   optional, intent(in) :: logunit !< IO unit for checksum logging
 
  logical :: sym
 
  if (present(haloshift)) then
    call chksum_b_3d(arraya, 'x '//mesg, hi, haloshift, symmetric, &
                     omit_corners, scale=scale, logunit=logunit)
    call chksum_b_3d(arrayb, 'y '//mesg, hi, haloshift, symmetric, &
                     omit_corners, scale=scale, logunit=logunit)
  else
    call chksum_b_3d(arraya, 'x '//mesg, hi, symmetric=symmetric, scale=scale, &
                     logunit=logunit)
    call chksum_b_3d(arrayb, 'y '//mesg, hi, symmetric=symmetric, scale=scale, &
                     logunit=logunit)
  endif
 
end subroutine chksum_pair_b_3d
 
!> Checksums a 2d array staggered at corner points.
subroutine chksum_b_2d(array, mesg, HI, haloshift, symmetric, omit_corners, &
                       scale, logunit)
  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.
  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
 
  real, allocatable, dimension(:,:) :: rescaled_array
  real :: scaling
  integer :: iounit !< Log IO unit
  integer :: i, j, Is, Js
  real :: aMean, aMin, aMax
  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
  iounit = error_unit; if(present(logunit)) iounit = logunit
  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
  if (calculatestatistics) then
    if (present(scale)) then
      allocate( rescaled_array(lbound(array,1):ubound(array,1), &
                               lbound(array,2):ubound(array,2)) )
      rescaled_array(:,:) = 0.0
      is = hi%isc ; if (sym_stats) is = hi%isc-1
      js = hi%jsc ; if (sym_stats) js = hi%jsc-1
      do j=js,hi%JecB ; do i=is,hi%IecB
        rescaled_array(i,j) = scale*array(i,j)
      enddo ; enddo
      call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
      deallocate(rescaled_array)
    else
      call substats(hi, array, sym_stats, amean, amin, amax)
    endif
    if (is_root_pe()) &
      call chk_sum_msg("B-point:", amean, amin, amax, mesg, iounit)
  endif
 
  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, iounit)
    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, iounit)
  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, iounit)
  endif
 
  contains
 
  integer function subchk(array, HI, di, dj, 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
    integer, intent(in) :: di    !< i- direction array shift for this checksum
    integer, intent(in) :: dj    !< j- direction array shift for this checksum
    real, intent(in)    :: scale !< A scaling factor for this array.
    integer :: 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, bc_modulus)
  end function subchk
 
  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
    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
    logical,          intent(in) :: sym_stats !< If true, evaluate the statistics on the
                                              !! full symmetric computational domain.
    real, intent(out) :: aMean, aMin, aMax
 
    integer :: i, j, n, IsB, JsB
 
    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)
  end subroutine substats
 
end subroutine chksum_b_2d
 
!> Checksums a pair of 2d velocity arrays staggered at C-grid locations
subroutine chksum_uv_2d(mesg, arrayU, arrayV, HI, haloshift, symmetric, &
                        omit_corners, scale, logunit)
  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.
  integer,                 optional, intent(in) :: logunit !< IO unit for checksum logging
 
  if (present(haloshift)) then
    call chksum_u_2d(arrayu, 'u '//mesg, hi, haloshift, symmetric, &
                     omit_corners, scale, logunit=logunit)
    call chksum_v_2d(arrayv, 'v '//mesg, hi, haloshift, symmetric, &
                     omit_corners, scale, logunit=logunit)
  else
    call chksum_u_2d(arrayu, 'u '//mesg, hi, symmetric=symmetric, &
                     logunit=logunit)
    call chksum_v_2d(arrayv, 'v '//mesg, hi, symmetric=symmetric, &
                     logunit=logunit)
  endif
 
end subroutine chksum_uv_2d
 
!> Checksums a pair of 3d velocity arrays staggered at C-grid locations
subroutine chksum_uv_3d(mesg, arrayU, arrayV, HI, haloshift, symmetric, &
                        omit_corners, scale, logunit)
  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.
  integer,                   optional, intent(in) :: logunit !< IO unit for checksum logging
 
  if (present(haloshift)) then
    call chksum_u_3d(arrayu, 'u '//mesg, hi, haloshift, symmetric, &
                     omit_corners, scale, logunit=logunit)
    call chksum_v_3d(arrayv, 'v '//mesg, hi, haloshift, symmetric, &
                     omit_corners, scale, logunit=logunit)
  else
    call chksum_u_3d(arrayu, 'u '//mesg, hi, symmetric=symmetric, &
                     logunit=logunit)
    call chksum_v_3d(arrayv, 'v '//mesg, hi, symmetric=symmetric, &
                     logunit=logunit)
  endif
 
end subroutine chksum_uv_3d
 
!> Checksums a 2d array staggered at C-grid u points.
subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, &
                       scale, logunit)
  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.
  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
 
  real, allocatable, dimension(:,:) :: rescaled_array
  real :: scaling
  integer :: iounit !< Log IO unit
  integer :: i, j, Is
  real :: aMean, aMin, aMax
  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
  iounit = error_unit; if(present(logunit)) iounit = logunit
  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
  if (calculatestatistics) then
    if (present(scale)) then
      allocate( rescaled_array(lbound(array,1):ubound(array,1), &
                               lbound(array,2):ubound(array,2)) )
      rescaled_array(:,:) = 0.0
      is = hi%isc ; if (sym_stats) is = hi%isc-1
      do j=hi%jsc,hi%jec ; do i=is,hi%IecB
        rescaled_array(i,j) = scale*array(i,j)
      enddo ; enddo
      call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
      deallocate(rescaled_array)
    else
      call substats(hi, array, sym_stats, amean, amin, amax)
    endif
 
    if (is_root_pe()) &
      call chk_sum_msg("u-point:", amean, amin, amax, mesg, iounit)
  endif
 
  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, iounit)
    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, iounit)
  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, iounit)
  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, iounit)
  endif
 
  contains
 
  integer function subchk(array, HI, di, dj, 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
    integer, intent(in) :: di    !< i- direction array shift for this checksum
    integer, intent(in) :: dj    !< j- direction array shift for this checksum
    real, intent(in)    :: scale !< A scaling factor for this array.
    integer :: 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, bc_modulus)
  end function subchk
 
  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
    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
    logical,          intent(in) :: sym_stats !< If true, evaluate the statistics on the
                                              !! full symmetric computational domain.
    real, intent(out) :: aMean, aMin, aMax
 
    integer :: i, j, n, IsB
 
    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)
  end subroutine substats
 
end subroutine chksum_u_2d
 
!> Checksums a 2d array staggered at C-grid v points.
subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, &
                       scale, logunit)
  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.
  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
 
  real, allocatable, dimension(:,:) :: rescaled_array
  real :: scaling
  integer :: iounit !< Log IO unit
  integer :: i, j, Js
  real :: aMean, aMin, aMax
  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
  iounit = error_unit; if(present(logunit)) iounit = logunit
  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
  if (calculatestatistics) then
    if (present(scale)) then
      allocate( rescaled_array(lbound(array,1):ubound(array,1), &
                               lbound(array,2):ubound(array,2)) )
      rescaled_array(:,:) = 0.0
      js = hi%jsc ; if (sym_stats) js = hi%jsc-1
      do j=js,hi%JecB ; do i=hi%isc,hi%iec
        rescaled_array(i,j) = scale*array(i,j)
      enddo ; enddo
      call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
      deallocate(rescaled_array)
    else
      call substats(hi, array, sym_stats, amean, amin, amax)
    endif
 
    if (is_root_pe()) &
      call chk_sum_msg("v-point:", amean, amin, amax, mesg, iounit)
  endif
 
  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, iounit)
    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, iounit)
  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, iounit)
  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, iounit)
  endif
 
  contains
 
  integer function subchk(array, HI, di, dj, 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
    integer, intent(in) :: di    !< i- direction array shift for this checksum
    integer, intent(in) :: dj    !< j- direction array shift for this checksum
    real, intent(in)    :: scale !< A scaling factor for this array.
    integer :: 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, bc_modulus)
  end function subchk
 
  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
    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
    logical,          intent(in) :: sym_stats !< If true, evaluate the statistics on the
                                              !! full symmetric computational domain.
    real, intent(out) :: aMean, aMin, aMax
 
    integer :: i, j, n, JsB
 
    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)
  end subroutine substats
 
end subroutine chksum_v_2d
 
!> Checksums a 3d array staggered at tracer points.
subroutine chksum_h_3d(array, mesg, HI, haloshift, omit_corners, scale, logunit)
  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.
  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
 
  real, allocatable, dimension(:,:,:) :: rescaled_array
  real :: scaling
  integer :: iounit !< Log IO unit
  integer :: i, j, k
  real :: aMean, aMin, aMax
  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
  iounit = error_unit; if(present(logunit)) iounit = logunit
 
  if (calculatestatistics) then
    if (present(scale)) then
      allocate( rescaled_array(lbound(array,1):ubound(array,1), &
                               lbound(array,2):ubound(array,2), &
                               lbound(array,3):ubound(array,3)) )
      rescaled_array(:,:,:) = 0.0
      do k=1,size(array,3) ; do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
        rescaled_array(i,j,k) = scale*array(i,j,k)
      enddo ; enddo ; enddo
 
      call substats(hi, rescaled_array, amean, amin, amax)
      deallocate(rescaled_array)
    else
      call substats(hi, array, amean, amin, amax)
    endif
 
    if (is_root_pe()) &
      call chk_sum_msg("h-point:", amean, amin, amax, mesg, iounit)
  endif
 
  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, iounit)
    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, iounit)
  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, iounit)
  endif
 
  contains
 
  integer function subchk(array, HI, di, dj, 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
    integer, intent(in) :: di    !< i- direction array shift for this checksum
    integer, intent(in) :: dj    !< j- direction array shift for this checksum
    real, intent(in)    :: scale !< A scaling factor for this array.
    integer :: 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, bc_modulus)
  end function subchk
 
  subroutine substats(HI, array, aMean, aMin, aMax)
    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
    real, intent(out) :: aMean, aMin, aMax
 
    integer :: i, j, k, n
 
    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)
  end subroutine substats
 
end subroutine chksum_h_3d
 
!> Checksums a 3d array staggered at corner points.
subroutine chksum_b_3d(array, mesg, HI, haloshift, symmetric, omit_corners, &
                       scale, logunit)
  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.
  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
 
  real, allocatable, dimension(:,:,:) :: rescaled_array
  real :: scaling
  integer :: iounit !< Log IO unit
  integer :: i, j, k, Is, Js
  real :: aMean, aMin, aMax
  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
  iounit = error_unit; if(present(logunit)) iounit = logunit
  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
  if (calculatestatistics) then
    if (present(scale)) then
      allocate( rescaled_array(lbound(array,1):ubound(array,1), &
                               lbound(array,2):ubound(array,2), &
                               lbound(array,3):ubound(array,3)) )
      rescaled_array(:,:,:) = 0.0
      is = hi%isc ; if (sym_stats) is = hi%isc-1
      js = hi%jsc ; if (sym_stats) js = hi%jsc-1
      do k=1,size(array,3) ; do j=js,hi%JecB ; do i=is,hi%IecB
        rescaled_array(i,j,k) = scale*array(i,j,k)
      enddo ; enddo ; enddo
      call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
      deallocate(rescaled_array)
    else
      call substats(hi, array, sym_stats, amean, amin, amax)
    endif
 
    if (is_root_pe()) &
      call chk_sum_msg("B-point:", amean, amin, amax, mesg, iounit)
  endif
 
  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, iounit)
    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, iounit)
  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, iounit)
  endif
 
  contains
 
  integer function subchk(array, HI, di, dj, 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
    integer, intent(in) :: di    !< i- direction array shift for this checksum
    integer, intent(in) :: dj    !< j- direction array shift for this checksum
    real, intent(in)    :: scale !< A scaling factor for this array.
    integer :: 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, bc_modulus)
  end function subchk
 
  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
    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
    logical,          intent(in) :: sym_stats !< If true, evaluate the statistics on the
                                              !! full symmetric computational domain.
    real, intent(out) :: aMean, aMin, aMax
 
    integer :: i, j, k, n, IsB, JsB
 
    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)
  end subroutine substats
 
end subroutine chksum_b_3d
 
!> Checksums a 3d array staggered at C-grid u points.
subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, &
                       scale, logunit)
  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.
  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
 
  real, allocatable, dimension(:,:,:) :: rescaled_array
  real :: scaling
  integer :: iounit !< Log IO unit
  integer :: i, j, k, Is
  real :: aMean, aMin, aMax
  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
  iounit = error_unit; if(present(logunit)) iounit = logunit
  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
  if (calculatestatistics) then
    if (present(scale)) then
      allocate( rescaled_array(lbound(array,1):ubound(array,1), &
                               lbound(array,2):ubound(array,2), &
                               lbound(array,3):ubound(array,3)) )
      rescaled_array(:,:,:) = 0.0
      is = hi%isc ; if (sym_stats) is = hi%isc-1
      do k=1,size(array,3) ; do j=hi%jsc,hi%jec ; do i=is,hi%IecB
        rescaled_array(i,j,k) = scale*array(i,j,k)
      enddo ; enddo ; enddo
      call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
      deallocate(rescaled_array)
    else
      call substats(hi, array, sym_stats, amean, amin, amax)
    endif
    if (is_root_pe()) &
      call chk_sum_msg("u-point:", amean, amin, amax, mesg, iounit)
  endif
 
  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, iounit)
    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, iounit)
  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, iounit)
  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, iounit)
  endif
 
  contains
 
  integer function subchk(array, HI, di, dj, 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
    integer, intent(in) :: di    !< i- direction array shift for this checksum
    integer, intent(in) :: dj    !< j- direction array shift for this checksum
    real, intent(in)    :: scale !< A scaling factor for this array.
    integer :: 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, bc_modulus)
  end function subchk
 
  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
    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
    logical,          intent(in) :: sym_stats !< If true, evaluate the statistics on the
                                              !! full symmetric computational domain.
    real, intent(out) :: aMean, aMin, aMax
 
    integer :: i, j, k, n, IsB
 
    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)
  end subroutine substats
 
end subroutine chksum_u_3d
 
!> Checksums a 3d array staggered at C-grid v points.
subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, &
                       scale, logunit)
  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.
  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
 
  real, allocatable, dimension(:,:,:) :: rescaled_array
  real :: scaling
  integer :: iounit !< Log IO unit
  integer :: i, j, k, Js
  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
  integer :: bcN, bcS, bcE, bcW
  real :: aMean, aMin, aMax
  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
  iounit = error_unit; if(present(logunit)) iounit = logunit
  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
 
  if (calculatestatistics) then
    if (present(scale)) then
      allocate( rescaled_array(lbound(array,1):ubound(array,1), &
                               lbound(array,2):ubound(array,2), &
                               lbound(array,3):ubound(array,3)) )
      rescaled_array(:,:,:) = 0.0
      js = hi%jsc ; if (sym_stats) js = hi%jsc-1
      do k=1,size(array,3) ; do j=js,hi%JecB ; do i=hi%isc,hi%iec
        rescaled_array(i,j,k) = scale*array(i,j,k)
      enddo ; enddo ; enddo
      call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
      deallocate(rescaled_array)
    else
      call substats(hi, array, sym_stats, amean, amin, amax)
    endif
    if (is_root_pe()) &
      call chk_sum_msg("v-point:", amean, amin, amax, mesg, iounit)
  endif
 
  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, iounit)
    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, iounit)
  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, iounit)
  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, iounit)
  endif
 
  contains
 
  integer function subchk(array, HI, di, dj, 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
    integer, intent(in) :: di    !< i- direction array shift for this checksum
    integer, intent(in) :: dj    !< j- direction array shift for this checksum
    real, intent(in)    :: scale !< A scaling factor for this array.
    integer :: 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, bc_modulus)
  end function subchk
 
  !subroutine subStats(HI, array, mesg, sym_stats)
  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
    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
    logical,          intent(in) :: sym_stats !< If true, evaluate the statistics on the
                                              !! full symmetric computational domain.
    real, intent(out) :: aMean, aMin, aMax    !< Mean/min/max of array over domain
 
    integer :: i, j, k, n, JsB
 
    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)
  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
  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 !< The array to be checksummed
  character(len=*) :: mesg  !< An identifying message
 
  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 !< The array to be checksummed
  character(len=*) :: mesg  !< An identifying message
 
  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
 
!> 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,  optional, intent(in) :: skip_mpp  !< If true, only check this array only
                                              !! on the local PE (default false).
  logical :: is_nan_1d
 
  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
 
!> 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
 
!> 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
 
!> Write a message including the checksum of the non-shifted array
subroutine chk_sum_msg1(fmsg, bc0, mesg, iounit)
  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
  integer,          intent(in) :: bc0  !< The bitcount of the non-shifted array
  integer,          intent(in) :: iounit !< Checksum logger IO unit
 
  if (is_root_pe()) &
    write(iounit, '(A,1(A,I10,X),A)') fmsg, " c=", bc0, trim(mesg)
end subroutine chk_sum_msg1
 
!> Write a message including checksums of non-shifted and diagonally shifted arrays
subroutine chk_sum_msg5(fmsg, bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit)
  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
  integer,          intent(in) :: bc0  !< The bitcount of the non-shifted array
  integer,          intent(in) :: bcSW !< The bitcount for SW shifted array
  integer,          intent(in) :: bcSE !< The bitcount for SE shifted array
  integer,          intent(in) :: bcNW !< The bitcount for NW shifted array
  integer,          intent(in) :: bcNE !< The bitcount for NE shifted array
  integer,          intent(in) :: iounit !< Checksum logger IO unit
 
  if (is_root_pe()) write(iounit, '(A,5(A,I10,1X),A)') &
    fmsg, " c=", bc0, "sw=", bcsw, "se=", bcse, "nw=", bcnw, "ne=", bcne, trim(mesg)
end subroutine chk_sum_msg5
 
!> Write a message including checksums of non-shifted and laterally shifted arrays
subroutine chk_sum_msg_nsew(fmsg, bc0, bcN, bcS, bcE, bcW, mesg, iounit)
  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
  integer,          intent(in) :: bc0  !< The bitcount of the non-shifted array
  integer,          intent(in) :: bcN !< The bitcount for N shifted array
  integer,          intent(in) :: bcS !< The bitcount for S shifted array
  integer,          intent(in) :: bcE !< The bitcount for E shifted array
  integer,          intent(in) :: bcW !< The bitcount for W shifted array
  integer,          intent(in) :: iounit !< Checksum logger IO unit
 
  if (is_root_pe()) write(iounit, '(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
 
!> Write a message including checksums of non-shifted and southward shifted arrays
subroutine chk_sum_msg_s(fmsg, bc0, bcS, mesg, iounit)
  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
  integer,          intent(in) :: bc0  !< The bitcount of the non-shifted array
  integer,          intent(in) :: bcS  !< The bitcount of the south-shifted array
  integer,          intent(in) :: iounit !< Checksum logger IO unit
 
  if (is_root_pe()) write(iounit, '(A,2(A,I10,1X),A)') &
    fmsg, " c=", bc0, "S=", bcs, trim(mesg)
end subroutine chk_sum_msg_s
 
!> Write a message including checksums of non-shifted and westward shifted arrays
subroutine chk_sum_msg_w(fmsg, bc0, bcW, mesg, iounit)
  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
  integer,          intent(in) :: bc0  !< The bitcount of the non-shifted array
  integer,          intent(in) :: bcW  !< The bitcount of the west-shifted array
  integer,          intent(in) :: iounit !< Checksum logger IO unit
 
  if (is_root_pe()) write(iounit, '(A,2(A,I10,1X),A)') &
    fmsg, " c=", bc0, "W=", bcw, trim(mesg)
end subroutine chk_sum_msg_w
 
!> Write a message including checksums of non-shifted and southwestward shifted arrays
subroutine chk_sum_msg2(fmsg, bc0, bcSW, mesg, iounit)
  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
  integer,          intent(in) :: bc0  !< The bitcount of the non-shifted array
  integer,          intent(in) :: bcSW !< The bitcount of the southwest-shifted array
  integer,          intent(in) :: iounit !< Checksum logger IO unit
 
  if (is_root_pe()) write(iounit, '(A,2(A,I9,1X),A)') &
    fmsg, " c=", bc0, "s/w=", bcsw, trim(mesg)
end subroutine chk_sum_msg2
 
!> Write a message including the global mean, maximum and minimum of an array
subroutine chk_sum_msg3(fmsg, aMean, aMin, aMax, mesg, iounit)
  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
  real,             intent(in) :: aMean !< The mean value of the array
  real,             intent(in) :: aMin !< The minimum value of the array
  real,             intent(in) :: aMax !< The maximum value of the array
  integer,          intent(in) :: iounit !< Checksum logger IO unit
 
  ! NOTE: We add zero to aMin and aMax to remove any negative zeros.
  ! This is due to inconsistencies of signed zero in local vs MPI calculations.
 
  if (is_root_pe()) write(iounit, '(A,3(A,ES25.16,1X),A)') &
    fmsg, " mean=", amean, "min=", (0. + amin), "max=", (0. + 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
 
!> A wrapper for MOM_error used in the checksum code
subroutine chksum_error(signal, message)
  ! Wrapper for MOM_error to help place specific break points in debuggers
  integer, intent(in) :: signal !< An error severity level, such as FATAL or WARNING
  character(len=*), intent(in) :: message !< An error message
  call mom_error(signal, message)
end subroutine chksum_error
 
!> Does a bitcount of a number by first casting to an integer and then using BTEST
!! to check bit by bit
integer function bitcount(x)
  real, intent(in) :: x !< Number to be bitcount
 
  integer, parameter :: xk = kind(x)  !< Kind type of x
 
  ! NOTE: Assumes that reals and integers of kind=xk are the same size
  bitcount = popcnt(transfer(x, 1_xk))
end function bitcount
 
end module mom_checksums