MOM6
MOM_checksums.F90
1 !> Routines to calculate checksums of various array and vector types
2 module mom_checksums
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_coms, only : pe_here, root_pe, num_pes, sum_across_pes
7 use mom_coms, only : min_across_pes, max_across_pes
8 use mom_coms, only : reproducing_sum
9 use mom_error_handler, only : mom_error, fatal, is_root_pe
10 use mom_file_parser, only : log_version, param_file_type
11 use mom_hor_index, only : hor_index_type
12 
13 use iso_fortran_env, only: error_unit
14 
15 implicit none ; private
16 
17 public :: chksum0, zchksum
18 public :: hchksum, bchksum, uchksum, vchksum, qchksum, is_nan, chksum
19 public :: hchksum_pair, uvchksum, bchksum_pair
20 public :: mom_checksums_init
21 
22 !> Checksums a pair of arrays (2d or 3d) staggered at tracer points
23 interface hchksum_pair
24  module procedure chksum_pair_h_2d, chksum_pair_h_3d
25 end interface
26 
27 !> Checksums a pair velocity arrays (2d or 3d) staggered at C-grid locations
28 interface uvchksum
29  module procedure chksum_uv_2d, chksum_uv_3d
30 end interface
31 
32 !> Checksums an array (2d or 3d) staggered at C-grid u points.
33 interface uchksum
34  module procedure chksum_u_2d, chksum_u_3d
35 end interface
36 
37 !> Checksums an array (2d or 3d) staggered at C-grid v points.
38 interface vchksum
39  module procedure chksum_v_2d, chksum_v_3d
40 end interface
41 
42 !> Checksums a pair of arrays (2d or 3d) staggered at corner points
43 interface bchksum_pair
44  module procedure chksum_pair_b_2d, chksum_pair_b_3d
45 end interface
46 
47 !> Checksums an array (2d or 3d) staggered at tracer points.
48 interface hchksum
49  module procedure chksum_h_2d, chksum_h_3d
50 end interface
51 
52 !> Checksums an array (2d or 3d) staggered at corner points.
53 interface bchksum
54  module procedure chksum_b_2d, chksum_b_3d
55 end interface
56 
57 !> This is an older interface that has been renamed Bchksum
58 interface qchksum
59  module procedure chksum_b_2d, chksum_b_3d
60 end interface
61 
62 !> This is an older interface for 1-, 2-, or 3-D checksums
63 interface chksum
64  module procedure chksum1d, chksum2d, chksum3d
65 end interface
66 
67 !> Write a message with either checksums or numerical statistics of arrays
68 interface chk_sum_msg
69  module procedure chk_sum_msg1, chk_sum_msg2, chk_sum_msg3, chk_sum_msg5
70 end interface
71 
72 !> Returns .true. if any element of x is a NaN, and .false. otherwise.
73 interface is_nan
74  module procedure is_nan_0d, is_nan_1d, is_nan_2d, is_nan_3d
75 end interface
76 
77 integer, parameter :: bc_modulus = 1000000000 !< Modulus of checksum bitcount
78 integer, parameter :: default_shift=0 !< The default array shift
79 logical :: calculatestatistics=.true. !< If true, report min, max and mean.
80 logical :: writechksums=.true. !< If true, report the bitcount checksum
81 logical :: checkfornans=.true. !< If true, checks array for NaNs and cause
82  !! FATAL error is any are found
83 
84 contains
85 
86 !> Checksum a scalar field (consistent with array checksums)
87 subroutine chksum0(scalar, mesg, scale, logunit)
88  real, intent(in) :: scalar !< The array to be checksummed
89  character(len=*), intent(in) :: mesg !< An identifying message
90  real, optional, intent(in) :: scale !< A scaling factor for this array.
91  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
92 
93  real :: scaling !< Explicit rescaling factor
94  integer :: iounit !< Log IO unit
95  real :: rs !< Rescaled scalar
96  integer :: bc !< Scalar bitcount
97 
98  if (checkfornans .and. is_nan(scalar)) &
99  call chksum_error(fatal, 'NaN detected: '//trim(mesg))
100 
101  scaling = 1.0 ; if (present(scale)) scaling = scale
102  iounit = error_unit; if(present(logunit)) iounit = logunit
103 
104  if (calculatestatistics) then
105  rs = scaling * scalar
106  if (is_root_pe()) &
107  call chk_sum_msg(" scalar:", rs, rs, rs, mesg, iounit)
108  endif
109 
110  if (.not. writechksums) return
111 
112  bc = mod(bitcount(abs(scaling * scalar)), bc_modulus)
113  if (is_root_pe()) &
114  call chk_sum_msg(" scalar:", bc, mesg, iounit)
115 
116 end subroutine chksum0
117 
118 
119 !> Checksum a 1d array (typically a column).
120 subroutine zchksum(array, mesg, scale, logunit)
121  real, dimension(:), intent(in) :: array !< The array to be checksummed
122  character(len=*), intent(in) :: mesg !< An identifying message
123  real, optional, intent(in) :: scale !< A scaling factor for this array.
124  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
125 
126  real, allocatable, dimension(:) :: rescaled_array
127  real :: scaling
128  integer :: iounit !< Log IO unit
129  integer :: k
130  real :: amean, amin, amax
131  integer :: bc0
132 
133  if (checkfornans) then
134  if (is_nan(array(:))) &
135  call chksum_error(fatal, 'NaN detected: '//trim(mesg))
136  endif
137 
138  scaling = 1.0 ; if (present(scale)) scaling = scale
139  iounit = error_unit; if(present(logunit)) iounit = logunit
140 
141  if (calculatestatistics) then
142  if (present(scale)) then
143  allocate(rescaled_array(lbound(array,1):ubound(array,1)))
144  rescaled_array(:) = 0.0
145  do k=1, size(array, 1)
146  rescaled_array(k) = scale * array(k)
147  enddo
148 
149  call substats(rescaled_array, amean, amin, amax)
150  deallocate(rescaled_array)
151  else
152  call substats(array, amean, amin, amax)
153  endif
154 
155  if (is_root_pe()) &
156  call chk_sum_msg(" column:", amean, amin, amax, mesg, iounit)
157  endif
158 
159  if (.not. writechksums) return
160 
161  bc0 = subchk(array, scaling)
162  if (is_root_pe()) call chk_sum_msg(" column:", bc0, mesg, iounit)
163 
164  contains
165 
166  integer function subchk(array, scale)
167  real, dimension(:), intent(in) :: array !< The array to be checksummed
168  real, intent(in) :: scale !< A scaling factor for this array.
169  integer :: k, bc
170  subchk = 0
171  do k=lbound(array, 1), ubound(array, 1)
172  bc = bitcount(abs(scale * array(k)))
173  subchk = subchk + bc
174  enddo
175  subchk=mod(subchk, bc_modulus)
176  end function subchk
177 
178  subroutine substats(array, aMean, aMin, aMax)
179  real, dimension(:), intent(in) :: array !< The array to be checksummed
180  real, intent(out) :: amean, amin, amax
181 
182  integer :: k, n
183 
184  amin = array(1)
185  amax = array(1)
186  n = 0
187  do k=lbound(array,1), ubound(array,1)
188  amin = min(amin, array(k))
189  amax = max(amax, array(k))
190  n = n + 1
191  enddo
192  amean = sum(array(:)) / real(n)
193  end subroutine substats
194 
195 end subroutine zchksum
196 
197 !> Checksums on a pair of 2d arrays staggered at tracer points.
198 subroutine chksum_pair_h_2d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, &
199  scale, logunit)
200  character(len=*), intent(in) :: mesg !< Identifying messages
201  type(hor_index_type), intent(in) :: HI !< A horizontal index type
202  real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayA !< The first array to be checksummed
203  real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayB !< The second array to be checksummed
204  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
205  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
206  real, optional, intent(in) :: scale !< A scaling factor for this array.
207  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
208 
209  if (present(haloshift)) then
210  call chksum_h_2d(arraya, 'x '//mesg, hi, haloshift, omit_corners, &
211  scale=scale, logunit=logunit)
212  call chksum_h_2d(arrayb, 'y '//mesg, hi, haloshift, omit_corners, &
213  scale=scale, logunit=logunit)
214  else
215  call chksum_h_2d(arraya, 'x '//mesg, hi, scale=scale, logunit=logunit)
216  call chksum_h_2d(arrayb, 'y '//mesg, hi, scale=scale, logunit=logunit)
217  endif
218 
219 end subroutine chksum_pair_h_2d
220 
221 !> Checksums on a pair of 3d arrays staggered at tracer points.
222 subroutine chksum_pair_h_3d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, &
223  scale, logunit)
224  character(len=*), intent(in) :: mesg !< Identifying messages
225  type(hor_index_type), intent(in) :: HI !< A horizontal index type
226  real, dimension(HI%isd:,HI%jsd:, :), intent(in) :: arrayA !< The first array to be checksummed
227  real, dimension(HI%isd:,HI%jsd:, :), intent(in) :: arrayB !< The second array to be checksummed
228  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
229  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
230  real, optional, intent(in) :: scale !< A scaling factor for this array.
231  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
232 
233  if (present(haloshift)) then
234  call chksum_h_3d(arraya, 'x '//mesg, hi, haloshift, omit_corners, &
235  scale=scale, logunit=logunit)
236  call chksum_h_3d(arrayb, 'y '//mesg, hi, haloshift, omit_corners, &
237  scale=scale, logunit=logunit)
238  else
239  call chksum_h_3d(arraya, 'x '//mesg, hi, scale=scale, logunit=logunit)
240  call chksum_h_3d(arrayb, 'y '//mesg, hi, scale=scale, logunit=logunit)
241  endif
242 
243 end subroutine chksum_pair_h_3d
244 
245 !> Checksums a 2d array staggered at tracer points.
246 subroutine chksum_h_2d(array, mesg, HI, haloshift, omit_corners, scale, logunit)
247  type(hor_index_type), intent(in) :: HI !< A horizontal index type
248  real, dimension(HI%isd:,HI%jsd:), intent(in) :: array !< The array to be checksummed
249  character(len=*), intent(in) :: mesg !< An identifying message
250  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
251  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
252  real, optional, intent(in) :: scale !< A scaling factor for this array.
253  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
254 
255  real, allocatable, dimension(:,:) :: rescaled_array
256  real :: scaling
257  integer :: iounit !< Log IO unit
258  integer :: i, j
259  real :: aMean, aMin, aMax
260  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
261  integer :: bcN, bcS, bcE, bcW
262  logical :: do_corners
263 
264  if (checkfornans) then
265  if (is_nan(array(hi%isc:hi%iec,hi%jsc:hi%jec))) &
266  call chksum_error(fatal, 'NaN detected: '//trim(mesg))
267 ! if (is_NaN(array)) &
268 ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
269  endif
270 
271  scaling = 1.0 ; if (present(scale)) scaling = scale
272  iounit = error_unit; if(present(logunit)) iounit = logunit
273 
274  if (calculatestatistics) then
275  if (present(scale)) then
276  allocate( rescaled_array(lbound(array,1):ubound(array,1), &
277  lbound(array,2):ubound(array,2)) )
278  rescaled_array(:,:) = 0.0
279  do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
280  rescaled_array(i,j) = scale*array(i,j)
281  enddo ; enddo
282  call substats(hi, rescaled_array, amean, amin, amax)
283  deallocate(rescaled_array)
284  else
285  call substats(hi, array, amean, amin, amax)
286  endif
287 
288  if (is_root_pe()) &
289  call chk_sum_msg("h-point:", amean, amin, amax, mesg, iounit)
290  endif
291 
292  if (.not.writechksums) return
293 
294  hshift = default_shift
295  if (present(haloshift)) hshift = haloshift
296  if (hshift<0) hshift = hi%ied-hi%iec
297 
298  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
299  hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
300  write(0,*) 'chksum_h_2d: haloshift =',hshift
301  write(0,*) 'chksum_h_2d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
302  write(0,*) 'chksum_h_2d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
303  call chksum_error(fatal,'Error in chksum_h_2d '//trim(mesg))
304  endif
305 
306  bc0 = subchk(array, hi, 0, 0, scaling)
307 
308  if (hshift==0) then
309  if (is_root_pe()) call chk_sum_msg("h-point:", bc0, mesg, iounit)
310  return
311  endif
312 
313  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners
314 
315  if (do_corners) then
316  bcsw = subchk(array, hi, -hshift, -hshift, scaling)
317  bcse = subchk(array, hi, hshift, -hshift, scaling)
318  bcnw = subchk(array, hi, -hshift, hshift, scaling)
319  bcne = subchk(array, hi, hshift, hshift, scaling)
320 
321  if (is_root_pe()) &
322  call chk_sum_msg("h-point:", bc0, bcsw, bcse, bcnw, bcne, mesg, iounit)
323  else
324  bcs = subchk(array, hi, 0, -hshift, scaling)
325  bce = subchk(array, hi, hshift, 0, scaling)
326  bcw = subchk(array, hi, -hshift, 0, scaling)
327  bcn = subchk(array, hi, 0, hshift, scaling)
328 
329  if (is_root_pe()) &
330  call chk_sum_msg_nsew("h-point:", bc0, bcn, bcs, bce, bcw, mesg, iounit)
331  endif
332 
333  contains
334  integer function subchk(array, HI, di, dj, scale)
335  type(hor_index_type), intent(in) :: HI !< A horizontal index type
336  real, dimension(HI%isd:,HI%jsd:), intent(in) :: array !< The array to be checksummed
337  integer, intent(in) :: di !< i- direction array shift for this checksum
338  integer, intent(in) :: dj !< j- direction array shift for this checksum
339  real, intent(in) :: scale !< A scaling factor for this array.
340  integer :: i, j, bc
341  subchk = 0
342  do j=hi%jsc+dj,hi%jec+dj; do i=hi%isc+di,hi%iec+di
343  bc = bitcount(abs(scale*array(i,j)))
344  subchk = subchk + bc
345  enddo ; enddo
346  call sum_across_pes(subchk)
347  subchk=mod(subchk, bc_modulus)
348  end function subchk
349 
350  subroutine substats(HI, array, aMean, aMin, aMax)
351  type(hor_index_type), intent(in) :: HI !< A horizontal index type
352  real, dimension(HI%isd:,HI%jsd:), intent(in) :: array !< The array to be checksummed
353  real, intent(out) :: aMean, aMin, aMax
354 
355  integer :: i, j, n
356 
357  amin = array(hi%isc,hi%jsc)
358  amax = array(hi%isc,hi%jsc)
359  n = 0
360  do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
361  amin = min(amin, array(i,j))
362  amax = max(amax, array(i,j))
363  n = n + 1
364  enddo ; enddo
365  amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec))
366  call sum_across_pes(n)
367  call min_across_pes(amin)
368  call max_across_pes(amax)
369  amean = amean / real(n)
370  end subroutine substats
371 
372 end subroutine chksum_h_2d
373 
374 !> Checksums on a pair of 2d arrays staggered at q-points.
375 subroutine chksum_pair_b_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, &
376  omit_corners, scale, logunit)
377  character(len=*), intent(in) :: mesg !< Identifying messages
378  type(hor_index_type), intent(in) :: HI !< A horizontal index type
379  real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayA !< The first array to be checksummed
380  real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayB !< The second array to be checksummed
381  logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full
382  !! symmetric computational domain.
383  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
384  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
385  real, optional, intent(in) :: scale !< A scaling factor for this array.
386  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
387 
388  logical :: sym
389 
390  sym = .false. ; if (present(symmetric)) sym = symmetric
391 
392  if (present(haloshift)) then
393  call chksum_b_2d(arraya, 'x '//mesg, hi, haloshift, symmetric=sym, &
394  omit_corners=omit_corners, scale=scale, logunit=logunit)
395  call chksum_b_2d(arrayb, 'y '//mesg, hi, haloshift, symmetric=sym, &
396  omit_corners=omit_corners, scale=scale, logunit=logunit)
397  else
398  call chksum_b_2d(arraya, 'x '//mesg, hi, symmetric=sym, scale=scale, &
399  logunit=logunit)
400  call chksum_b_2d(arrayb, 'y '//mesg, hi, symmetric=sym, scale=scale, &
401  logunit=logunit)
402  endif
403 
404 end subroutine chksum_pair_b_2d
405 
406 !> Checksums on a pair of 3d arrays staggered at q-points.
407 subroutine chksum_pair_b_3d(mesg, arrayA, arrayB, HI, haloshift, symmetric, &
408  omit_corners, scale, logunit)
409  character(len=*), intent(in) :: mesg !< Identifying messages
410  type(hor_index_type), intent(in) :: HI !< A horizontal index type
411  real, dimension(HI%IsdB:,HI%JsdB:, :), intent(in) :: arrayA !< The first array to be checksummed
412  real, dimension(HI%IsdB:,HI%JsdB:, :), intent(in) :: arrayB !< The second array to be checksummed
413  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
414  logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full
415  !! symmetric computational domain.
416  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
417  real, optional, intent(in) :: scale !< A scaling factor for this array.
418  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
419 
420  logical :: sym
421 
422  if (present(haloshift)) then
423  call chksum_b_3d(arraya, 'x '//mesg, hi, haloshift, symmetric, &
424  omit_corners, scale=scale, logunit=logunit)
425  call chksum_b_3d(arrayb, 'y '//mesg, hi, haloshift, symmetric, &
426  omit_corners, scale=scale, logunit=logunit)
427  else
428  call chksum_b_3d(arraya, 'x '//mesg, hi, symmetric=symmetric, scale=scale, &
429  logunit=logunit)
430  call chksum_b_3d(arrayb, 'y '//mesg, hi, symmetric=symmetric, scale=scale, &
431  logunit=logunit)
432  endif
433 
434 end subroutine chksum_pair_b_3d
435 
436 !> Checksums a 2d array staggered at corner points.
437 subroutine chksum_b_2d(array, mesg, HI, haloshift, symmetric, omit_corners, &
438  scale, logunit)
439  type(hor_index_type), intent(in) :: HI !< A horizontal index type
440  real, dimension(HI%IsdB:,HI%JsdB:), &
441  intent(in) :: array !< The array to be checksummed
442  character(len=*), intent(in) :: mesg !< An identifying message
443  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
444  logical, optional, intent(in) :: symmetric !< If true, do the checksums on the
445  !! full symmetric computational domain.
446  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
447  real, optional, intent(in) :: scale !< A scaling factor for this array.
448  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
449 
450  real, allocatable, dimension(:,:) :: rescaled_array
451  real :: scaling
452  integer :: iounit !< Log IO unit
453  integer :: i, j, Is, Js
454  real :: aMean, aMin, aMax
455  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
456  integer :: bcN, bcS, bcE, bcW
457  logical :: do_corners, sym, sym_stats
458 
459  if (checkfornans) then
460  if (is_nan(array(hi%IscB:hi%IecB,hi%JscB:hi%JecB))) &
461  call chksum_error(fatal, 'NaN detected: '//trim(mesg))
462 ! if (is_NaN(array)) &
463 ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
464  endif
465 
466  scaling = 1.0 ; if (present(scale)) scaling = scale
467  iounit = error_unit; if(present(logunit)) iounit = logunit
468  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
469  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
470 
471  if (calculatestatistics) then
472  if (present(scale)) then
473  allocate( rescaled_array(lbound(array,1):ubound(array,1), &
474  lbound(array,2):ubound(array,2)) )
475  rescaled_array(:,:) = 0.0
476  is = hi%isc ; if (sym_stats) is = hi%isc-1
477  js = hi%jsc ; if (sym_stats) js = hi%jsc-1
478  do j=js,hi%JecB ; do i=is,hi%IecB
479  rescaled_array(i,j) = scale*array(i,j)
480  enddo ; enddo
481  call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
482  deallocate(rescaled_array)
483  else
484  call substats(hi, array, sym_stats, amean, amin, amax)
485  endif
486  if (is_root_pe()) &
487  call chk_sum_msg("B-point:", amean, amin, amax, mesg, iounit)
488  endif
489 
490  if (.not.writechksums) return
491 
492  hshift = default_shift
493  if (present(haloshift)) hshift = haloshift
494  if (hshift<0) hshift = hi%ied-hi%iec
495 
496  if ( hi%iscB-hshift<hi%isdB .or. hi%iecB+hshift>hi%iedB .or. &
497  hi%jscB-hshift<hi%jsdB .or. hi%jecB+hshift>hi%jedB ) then
498  write(0,*) 'chksum_B_2d: haloshift =',hshift
499  write(0,*) 'chksum_B_2d: isd,isc,iec,ied=',hi%isdB,hi%iscB,hi%iecB,hi%iedB
500  write(0,*) 'chksum_B_2d: jsd,jsc,jec,jed=',hi%jsdB,hi%jscB,hi%jecB,hi%jedB
501  call chksum_error(fatal,'Error in chksum_B_2d '//trim(mesg))
502  endif
503 
504  bc0 = subchk(array, hi, 0, 0, scaling)
505 
506  sym = .false. ; if (present(symmetric)) sym = symmetric
507 
508  if ((hshift==0) .and. .not.sym) then
509  if (is_root_pe()) call chk_sum_msg("B-point:", bc0, mesg, iounit)
510  return
511  endif
512 
513  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners
514 
515  if (do_corners) then
516  if (sym) then
517  bcsw = subchk(array, hi, -hshift-1, -hshift-1, scaling)
518  bcse = subchk(array, hi, hshift, -hshift-1, scaling)
519  bcnw = subchk(array, hi, -hshift-1, hshift, scaling)
520  else
521  bcsw = subchk(array, hi, -hshift, -hshift, scaling)
522  bcse = subchk(array, hi, hshift, -hshift, scaling)
523  bcnw = subchk(array, hi, -hshift, hshift, scaling)
524  endif
525  bcne = subchk(array, hi, hshift, hshift, scaling)
526 
527  if (is_root_pe()) &
528  call chk_sum_msg("B-point:", bc0, bcsw, bcse, bcnw, bcne, mesg, iounit)
529  else
530  bcs = subchk(array, hi, 0, -hshift, scaling)
531  bce = subchk(array, hi, hshift, 0, scaling)
532  bcw = subchk(array, hi, -hshift, 0, scaling)
533  bcn = subchk(array, hi, 0, hshift, scaling)
534 
535  if (is_root_pe()) &
536  call chk_sum_msg_nsew("B-point:", bc0, bcn, bcs, bce, bcw, mesg, iounit)
537  endif
538 
539  contains
540 
541  integer function subchk(array, HI, di, dj, scale)
542  type(hor_index_type), intent(in) :: HI !< A horizontal index type
543  real, dimension(HI%IsdB:,HI%JsdB:), intent(in) :: array !< The array to be checksummed
544  integer, intent(in) :: di !< i- direction array shift for this checksum
545  integer, intent(in) :: dj !< j- direction array shift for this checksum
546  real, intent(in) :: scale !< A scaling factor for this array.
547  integer :: i, j, bc
548  subchk = 0
549  ! This line deliberately uses the h-point computational domain.
550  do j=hi%jsc+dj,hi%jec+dj; do i=hi%isc+di,hi%iec+di
551  bc = bitcount(abs(scale*array(i,j)))
552  subchk = subchk + bc
553  enddo ; enddo
554  call sum_across_pes(subchk)
555  subchk=mod(subchk, bc_modulus)
556  end function subchk
557 
558  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
559  type(hor_index_type), intent(in) :: HI !< A horizontal index type
560  real, dimension(HI%IsdB:,HI%JsdB:), intent(in) :: array !< The array to be checksummed
561  logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the
562  !! full symmetric computational domain.
563  real, intent(out) :: aMean, aMin, aMax
564 
565  integer :: i, j, n, IsB, JsB
566 
567  isb = hi%isc ; if (sym_stats) isb = hi%isc-1
568  jsb = hi%jsc ; if (sym_stats) jsb = hi%jsc-1
569 
570  amin = array(hi%isc,hi%jsc) ; amax = amin
571  do j=jsb,hi%JecB ; do i=isb,hi%IecB
572  amin = min(amin, array(i,j))
573  amax = max(amax, array(i,j))
574  enddo ; enddo
575  ! This line deliberately uses the h-point computational domain.
576  amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec))
577  n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc)
578  call sum_across_pes(n)
579  call min_across_pes(amin)
580  call max_across_pes(amax)
581  amean = amean / real(n)
582  end subroutine substats
583 
584 end subroutine chksum_b_2d
585 
586 !> Checksums a pair of 2d velocity arrays staggered at C-grid locations
587 subroutine chksum_uv_2d(mesg, arrayU, arrayV, HI, haloshift, symmetric, &
588  omit_corners, scale, logunit)
589  character(len=*), intent(in) :: mesg !< Identifying messages
590  type(hor_index_type), intent(in) :: HI !< A horizontal index type
591  real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: arrayU !< The u-component array to be checksummed
592  real, dimension(HI%isd:,HI%JsdB:), intent(in) :: arrayV !< The v-component array to be checksummed
593  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
594  logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full
595  !! symmetric computational domain.
596  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
597  real, optional, intent(in) :: scale !< A scaling factor for these arrays.
598  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
599 
600  if (present(haloshift)) then
601  call chksum_u_2d(arrayu, 'u '//mesg, hi, haloshift, symmetric, &
602  omit_corners, scale, logunit=logunit)
603  call chksum_v_2d(arrayv, 'v '//mesg, hi, haloshift, symmetric, &
604  omit_corners, scale, logunit=logunit)
605  else
606  call chksum_u_2d(arrayu, 'u '//mesg, hi, symmetric=symmetric, &
607  logunit=logunit)
608  call chksum_v_2d(arrayv, 'v '//mesg, hi, symmetric=symmetric, &
609  logunit=logunit)
610  endif
611 
612 end subroutine chksum_uv_2d
613 
614 !> Checksums a pair of 3d velocity arrays staggered at C-grid locations
615 subroutine chksum_uv_3d(mesg, arrayU, arrayV, HI, haloshift, symmetric, &
616  omit_corners, scale, logunit)
617  character(len=*), intent(in) :: mesg !< Identifying messages
618  type(hor_index_type), intent(in) :: HI !< A horizontal index type
619  real, dimension(HI%IsdB:,HI%jsd:,:), intent(in) :: arrayU !< The u-component array to be checksummed
620  real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: arrayV !< The v-component array to be checksummed
621  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
622  logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full
623  !! symmetric computational domain.
624  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
625  real, optional, intent(in) :: scale !< A scaling factor for these arrays.
626  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
627 
628  if (present(haloshift)) then
629  call chksum_u_3d(arrayu, 'u '//mesg, hi, haloshift, symmetric, &
630  omit_corners, scale, logunit=logunit)
631  call chksum_v_3d(arrayv, 'v '//mesg, hi, haloshift, symmetric, &
632  omit_corners, scale, logunit=logunit)
633  else
634  call chksum_u_3d(arrayu, 'u '//mesg, hi, symmetric=symmetric, &
635  logunit=logunit)
636  call chksum_v_3d(arrayv, 'v '//mesg, hi, symmetric=symmetric, &
637  logunit=logunit)
638  endif
639 
640 end subroutine chksum_uv_3d
641 
642 !> Checksums a 2d array staggered at C-grid u points.
643 subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, &
644  scale, logunit)
645  type(hor_index_type), intent(in) :: HI !< A horizontal index type
646  real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: array !< The array to be checksummed
647  character(len=*), intent(in) :: mesg !< An identifying message
648  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
649  logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full
650  !! symmetric computational domain.
651  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
652  real, optional, intent(in) :: scale !< A scaling factor for this array.
653  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
654 
655  real, allocatable, dimension(:,:) :: rescaled_array
656  real :: scaling
657  integer :: iounit !< Log IO unit
658  integer :: i, j, Is
659  real :: aMean, aMin, aMax
660  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
661  integer :: bcN, bcS, bcE, bcW
662  logical :: do_corners, sym, sym_stats
663 
664  if (checkfornans) then
665  if (is_nan(array(hi%IscB:hi%IecB,hi%jsc:hi%jec))) &
666  call chksum_error(fatal, 'NaN detected: '//trim(mesg))
667 ! if (is_NaN(array)) &
668 ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
669  endif
670 
671  scaling = 1.0 ; if (present(scale)) scaling = scale
672  iounit = error_unit; if(present(logunit)) iounit = logunit
673  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
674  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
675 
676  if (calculatestatistics) then
677  if (present(scale)) then
678  allocate( rescaled_array(lbound(array,1):ubound(array,1), &
679  lbound(array,2):ubound(array,2)) )
680  rescaled_array(:,:) = 0.0
681  is = hi%isc ; if (sym_stats) is = hi%isc-1
682  do j=hi%jsc,hi%jec ; do i=is,hi%IecB
683  rescaled_array(i,j) = scale*array(i,j)
684  enddo ; enddo
685  call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
686  deallocate(rescaled_array)
687  else
688  call substats(hi, array, sym_stats, amean, amin, amax)
689  endif
690 
691  if (is_root_pe()) &
692  call chk_sum_msg("u-point:", amean, amin, amax, mesg, iounit)
693  endif
694 
695  if (.not.writechksums) return
696 
697  hshift = default_shift
698  if (present(haloshift)) hshift = haloshift
699  if (hshift<0) hshift = hi%iedB-hi%iecB
700 
701  if ( hi%iscB-hshift<hi%isdB .or. hi%iecB+hshift>hi%iedB .or. &
702  hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
703  write(0,*) 'chksum_u_2d: haloshift =',hshift
704  write(0,*) 'chksum_u_2d: isd,isc,iec,ied=',hi%isdB,hi%iscB,hi%iecB,hi%iedB
705  write(0,*) 'chksum_u_2d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
706  call chksum_error(fatal,'Error in chksum_u_2d '//trim(mesg))
707  endif
708 
709  bc0 = subchk(array, hi, 0, 0, scaling)
710 
711  sym = .false. ; if (present(symmetric)) sym = symmetric
712 
713  if ((hshift==0) .and. .not.sym) then
714  if (is_root_pe()) call chk_sum_msg("u-point:", bc0, mesg, iounit)
715  return
716  endif
717 
718  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners
719 
720  if (hshift==0) then
721  bcw = subchk(array, hi, -hshift-1, 0, scaling)
722  if (is_root_pe()) call chk_sum_msg_w("u-point:", bc0, bcw, mesg, iounit)
723  elseif (do_corners) then
724  if (sym) then
725  bcsw = subchk(array, hi, -hshift-1, -hshift, scaling)
726  bcnw = subchk(array, hi, -hshift-1, hshift, scaling)
727  else
728  bcsw = subchk(array, hi, -hshift, -hshift, scaling)
729  bcnw = subchk(array, hi, -hshift, hshift, scaling)
730  endif
731  bcse = subchk(array, hi, hshift, -hshift, scaling)
732  bcne = subchk(array, hi, hshift, hshift, scaling)
733 
734  if (is_root_pe()) &
735  call chk_sum_msg("u-point:", bc0, bcsw, bcse, bcnw, bcne, mesg, iounit)
736  else
737  bcs = subchk(array, hi, 0, -hshift, scaling)
738  bce = subchk(array, hi, hshift, 0, scaling)
739  if (sym) then
740  bcw = subchk(array, hi, -hshift-1, 0, scaling)
741  else
742  bcw = subchk(array, hi, -hshift, 0, scaling)
743  endif
744  bcn = subchk(array, hi, 0, hshift, scaling)
745 
746  if (is_root_pe()) &
747  call chk_sum_msg_nsew("u-point:", bc0, bcn, bcs, bce, bcw, mesg, iounit)
748  endif
749 
750  contains
751 
752  integer function subchk(array, HI, di, dj, scale)
753  type(hor_index_type), intent(in) :: HI !< A horizontal index type
754  real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: array !< The array to be checksummed
755  integer, intent(in) :: di !< i- direction array shift for this checksum
756  integer, intent(in) :: dj !< j- direction array shift for this checksum
757  real, intent(in) :: scale !< A scaling factor for this array.
758  integer :: i, j, bc
759  subchk = 0
760  ! This line deliberately uses the h-point computational domain.
761  do j=hi%jsc+dj,hi%jec+dj; do i=hi%isc+di,hi%iec+di
762  bc = bitcount(abs(scale*array(i,j)))
763  subchk = subchk + bc
764  enddo ; enddo
765  call sum_across_pes(subchk)
766  subchk=mod(subchk, bc_modulus)
767  end function subchk
768 
769  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
770  type(hor_index_type), intent(in) :: HI !< A horizontal index type
771  real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: array !< The array to be checksummed
772  logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the
773  !! full symmetric computational domain.
774  real, intent(out) :: aMean, aMin, aMax
775 
776  integer :: i, j, n, IsB
777 
778  isb = hi%isc ; if (sym_stats) isb = hi%isc-1
779 
780  amin = array(hi%isc,hi%jsc) ; amax = amin
781  do j=hi%jsc,hi%jec ; do i=isb,hi%IecB
782  amin = min(amin, array(i,j))
783  amax = max(amax, array(i,j))
784  enddo ; enddo
785  ! This line deliberately uses the h-point computational domain.
786  amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec))
787  n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc)
788  call sum_across_pes(n)
789  call min_across_pes(amin)
790  call max_across_pes(amax)
791  amean = amean / real(n)
792  end subroutine substats
793 
794 end subroutine chksum_u_2d
795 
796 !> Checksums a 2d array staggered at C-grid v points.
797 subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, &
798  scale, logunit)
799  type(hor_index_type), intent(in) :: HI !< A horizontal index type
800  real, dimension(HI%isd:,HI%JsdB:), intent(in) :: array !< The array to be checksummed
801  character(len=*), intent(in) :: mesg !< An identifying message
802  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
803  logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full
804  !! symmetric computational domain.
805  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
806  real, optional, intent(in) :: scale !< A scaling factor for this array.
807  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
808 
809  real, allocatable, dimension(:,:) :: rescaled_array
810  real :: scaling
811  integer :: iounit !< Log IO unit
812  integer :: i, j, Js
813  real :: aMean, aMin, aMax
814  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
815  integer :: bcN, bcS, bcE, bcW
816  logical :: do_corners, sym, sym_stats
817 
818  if (checkfornans) then
819  if (is_nan(array(hi%isc:hi%iec,hi%JscB:hi%JecB))) &
820  call chksum_error(fatal, 'NaN detected: '//trim(mesg))
821 ! if (is_NaN(array)) &
822 ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
823  endif
824 
825  scaling = 1.0 ; if (present(scale)) scaling = scale
826  iounit = error_unit; if(present(logunit)) iounit = logunit
827  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
828  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
829 
830  if (calculatestatistics) then
831  if (present(scale)) then
832  allocate( rescaled_array(lbound(array,1):ubound(array,1), &
833  lbound(array,2):ubound(array,2)) )
834  rescaled_array(:,:) = 0.0
835  js = hi%jsc ; if (sym_stats) js = hi%jsc-1
836  do j=js,hi%JecB ; do i=hi%isc,hi%iec
837  rescaled_array(i,j) = scale*array(i,j)
838  enddo ; enddo
839  call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
840  deallocate(rescaled_array)
841  else
842  call substats(hi, array, sym_stats, amean, amin, amax)
843  endif
844 
845  if (is_root_pe()) &
846  call chk_sum_msg("v-point:", amean, amin, amax, mesg, iounit)
847  endif
848 
849  if (.not.writechksums) return
850 
851  hshift = default_shift
852  if (present(haloshift)) hshift = haloshift
853  if (hshift<0) hshift = hi%ied-hi%iec
854 
855  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
856  hi%jscB-hshift<hi%jsdB .or. hi%jecB+hshift>hi%jedB ) then
857  write(0,*) 'chksum_v_2d: haloshift =',hshift
858  write(0,*) 'chksum_v_2d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
859  write(0,*) 'chksum_v_2d: jsd,jsc,jec,jed=',hi%jsdB,hi%jscB,hi%jecB,hi%jedB
860  call chksum_error(fatal,'Error in chksum_v_2d '//trim(mesg))
861  endif
862 
863  bc0 = subchk(array, hi, 0, 0, scaling)
864 
865  sym = .false. ; if (present(symmetric)) sym = symmetric
866 
867  if ((hshift==0) .and. .not.sym) then
868  if (is_root_pe()) call chk_sum_msg("v-point:", bc0, mesg, iounit)
869  return
870  endif
871 
872  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners
873 
874  if (hshift==0) then
875  bcs = subchk(array, hi, 0, -hshift-1, scaling)
876  if (is_root_pe()) call chk_sum_msg_s("v-point:", bc0, bcs, mesg, iounit)
877  elseif (do_corners) then
878  if (sym) then
879  bcsw = subchk(array, hi, -hshift, -hshift-1, scaling)
880  bcse = subchk(array, hi, hshift, -hshift-1, scaling)
881  else
882  bcsw = subchk(array, hi, -hshift, -hshift, scaling)
883  bcse = subchk(array, hi, hshift, -hshift, scaling)
884  endif
885  bcnw = subchk(array, hi, -hshift, hshift, scaling)
886  bcne = subchk(array, hi, hshift, hshift, scaling)
887 
888  if (is_root_pe()) &
889  call chk_sum_msg("v-point:", bc0, bcsw, bcse, bcnw, bcne, mesg, iounit)
890  else
891  if (sym) then
892  bcs = subchk(array, hi, 0, -hshift-1, scaling)
893  else
894  bcs = subchk(array, hi, 0, -hshift, scaling)
895  endif
896  bce = subchk(array, hi, hshift, 0, scaling)
897  bcw = subchk(array, hi, -hshift, 0, scaling)
898  bcn = subchk(array, hi, 0, hshift, scaling)
899 
900  if (is_root_pe()) &
901  call chk_sum_msg_nsew("v-point:", bc0, bcn, bcs, bce, bcw, mesg, iounit)
902  endif
903 
904  contains
905 
906  integer function subchk(array, HI, di, dj, scale)
907  type(hor_index_type), intent(in) :: HI !< A horizontal index type
908  real, dimension(HI%isd:,HI%JsdB:), intent(in) :: array !< The array to be checksummed
909  integer, intent(in) :: di !< i- direction array shift for this checksum
910  integer, intent(in) :: dj !< j- direction array shift for this checksum
911  real, intent(in) :: scale !< A scaling factor for this array.
912  integer :: i, j, bc
913  subchk = 0
914  ! This line deliberately uses the h-point computational domain.
915  do j=hi%jsc+dj,hi%jec+dj; do i=hi%isc+di,hi%iec+di
916  bc = bitcount(abs(scale*array(i,j)))
917  subchk = subchk + bc
918  enddo ; enddo
919  call sum_across_pes(subchk)
920  subchk=mod(subchk, bc_modulus)
921  end function subchk
922 
923  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
924  type(hor_index_type), intent(in) :: HI !< A horizontal index type
925  real, dimension(HI%isd:,HI%JsdB:), intent(in) :: array !< The array to be checksummed
926  logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the
927  !! full symmetric computational domain.
928  real, intent(out) :: aMean, aMin, aMax
929 
930  integer :: i, j, n, JsB
931 
932  jsb = hi%jsc ; if (sym_stats) jsb = hi%jsc-1
933 
934  amin = array(hi%isc,hi%jsc) ; amax = amin
935  do j=jsb,hi%JecB ; do i=hi%isc,hi%iec
936  amin = min(amin, array(i,j))
937  amax = max(amax, array(i,j))
938  enddo ; enddo
939  ! This line deliberately uses the h-computational domain.
940  amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec))
941  n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc)
942  call sum_across_pes(n)
943  call min_across_pes(amin)
944  call max_across_pes(amax)
945  amean = amean / real(n)
946  end subroutine substats
947 
948 end subroutine chksum_v_2d
949 
950 !> Checksums a 3d array staggered at tracer points.
951 subroutine chksum_h_3d(array, mesg, HI, haloshift, omit_corners, scale, logunit)
952  type(hor_index_type), intent(in) :: HI !< A horizontal index type
953  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed
954  character(len=*), intent(in) :: mesg !< An identifying message
955  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
956  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
957  real, optional, intent(in) :: scale !< A scaling factor for this array.
958  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
959 
960  real, allocatable, dimension(:,:,:) :: rescaled_array
961  real :: scaling
962  integer :: iounit !< Log IO unit
963  integer :: i, j, k
964  real :: aMean, aMin, aMax
965  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
966  integer :: bcN, bcS, bcE, bcW
967  logical :: do_corners
968 
969  if (checkfornans) then
970  if (is_nan(array(hi%isc:hi%iec,hi%jsc:hi%jec,:))) &
971  call chksum_error(fatal, 'NaN detected: '//trim(mesg))
972 ! if (is_NaN(array)) &
973 ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
974  endif
975 
976  scaling = 1.0 ; if (present(scale)) scaling = scale
977  iounit = error_unit; if(present(logunit)) iounit = logunit
978 
979  if (calculatestatistics) then
980  if (present(scale)) then
981  allocate( rescaled_array(lbound(array,1):ubound(array,1), &
982  lbound(array,2):ubound(array,2), &
983  lbound(array,3):ubound(array,3)) )
984  rescaled_array(:,:,:) = 0.0
985  do k=1,size(array,3) ; do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
986  rescaled_array(i,j,k) = scale*array(i,j,k)
987  enddo ; enddo ; enddo
988 
989  call substats(hi, rescaled_array, amean, amin, amax)
990  deallocate(rescaled_array)
991  else
992  call substats(hi, array, amean, amin, amax)
993  endif
994 
995  if (is_root_pe()) &
996  call chk_sum_msg("h-point:", amean, amin, amax, mesg, iounit)
997  endif
998 
999  if (.not.writechksums) return
1000 
1001  hshift = default_shift
1002  if (present(haloshift)) hshift = haloshift
1003  if (hshift<0) hshift = hi%ied-hi%iec
1004 
1005  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
1006  hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
1007  write(0,*) 'chksum_h_3d: haloshift =',hshift
1008  write(0,*) 'chksum_h_3d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
1009  write(0,*) 'chksum_h_3d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
1010  call chksum_error(fatal,'Error in chksum_h_3d '//trim(mesg))
1011  endif
1012 
1013  bc0 = subchk(array, hi, 0, 0, scaling)
1014 
1015  if (hshift==0) then
1016  if (is_root_pe()) call chk_sum_msg("h-point:", bc0, mesg, iounit)
1017  return
1018  endif
1019 
1020  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners
1021 
1022  if (do_corners) then
1023  bcsw = subchk(array, hi, -hshift, -hshift, scaling)
1024  bcse = subchk(array, hi, hshift, -hshift, scaling)
1025  bcnw = subchk(array, hi, -hshift, hshift, scaling)
1026  bcne = subchk(array, hi, hshift, hshift, scaling)
1027 
1028  if (is_root_pe()) &
1029  call chk_sum_msg("h-point:", bc0, bcsw, bcse, bcnw, bcne, mesg, iounit)
1030  else
1031  bcs = subchk(array, hi, 0, -hshift, scaling)
1032  bce = subchk(array, hi, hshift, 0, scaling)
1033  bcw = subchk(array, hi, -hshift, 0, scaling)
1034  bcn = subchk(array, hi, 0, hshift, scaling)
1035 
1036  if (is_root_pe()) &
1037  call chk_sum_msg_nsew("h-point:", bc0, bcn, bcs, bce, bcw, mesg, iounit)
1038  endif
1039 
1040  contains
1041 
1042  integer function subchk(array, HI, di, dj, scale)
1043  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1044  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed
1045  integer, intent(in) :: di !< i- direction array shift for this checksum
1046  integer, intent(in) :: dj !< j- direction array shift for this checksum
1047  real, intent(in) :: scale !< A scaling factor for this array.
1048  integer :: i, j, k, bc
1049  subchk = 0
1050  do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc+dj,hi%jec+dj ; do i=hi%isc+di,hi%iec+di
1051  bc = bitcount(abs(scale*array(i,j,k)))
1052  subchk = subchk + bc
1053  enddo ; enddo ; enddo
1054  call sum_across_pes(subchk)
1055  subchk=mod(subchk, bc_modulus)
1056  end function subchk
1057 
1058  subroutine substats(HI, array, aMean, aMin, aMax)
1059  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1060  real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed
1061  real, intent(out) :: aMean, aMin, aMax
1062 
1063  integer :: i, j, k, n
1064 
1065  amin = array(hi%isc,hi%jsc,1)
1066  amax = array(hi%isc,hi%jsc,1)
1067  n = 0
1068  do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
1069  amin = min(amin, array(i,j,k))
1070  amax = max(amax, array(i,j,k))
1071  n = n + 1
1072  enddo ; enddo ; enddo
1073  amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec,:))
1074  call sum_across_pes(n)
1075  call min_across_pes(amin)
1076  call max_across_pes(amax)
1077  amean = amean / real(n)
1078  end subroutine substats
1079 
1080 end subroutine chksum_h_3d
1081 
1082 !> Checksums a 3d array staggered at corner points.
1083 subroutine chksum_b_3d(array, mesg, HI, haloshift, symmetric, omit_corners, &
1084  scale, logunit)
1085  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1086  real, dimension(HI%IsdB:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed
1087  character(len=*), intent(in) :: mesg !< An identifying message
1088  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
1089  logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full
1090  !! symmetric computational domain.
1091  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
1092  real, optional, intent(in) :: scale !< A scaling factor for this array.
1093  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
1094 
1095  real, allocatable, dimension(:,:,:) :: rescaled_array
1096  real :: scaling
1097  integer :: iounit !< Log IO unit
1098  integer :: i, j, k, Is, Js
1099  real :: aMean, aMin, aMax
1100  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
1101  integer :: bcN, bcS, bcE, bcW
1102  logical :: do_corners, sym, sym_stats
1103 
1104  if (checkfornans) then
1105  if (is_nan(array(hi%IscB:hi%IecB,hi%JscB:hi%JecB,:))) &
1106  call chksum_error(fatal, 'NaN detected: '//trim(mesg))
1107 ! if (is_NaN(array)) &
1108 ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
1109  endif
1110 
1111  scaling = 1.0 ; if (present(scale)) scaling = scale
1112  iounit = error_unit; if(present(logunit)) iounit = logunit
1113  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
1114  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
1115 
1116  if (calculatestatistics) then
1117  if (present(scale)) then
1118  allocate( rescaled_array(lbound(array,1):ubound(array,1), &
1119  lbound(array,2):ubound(array,2), &
1120  lbound(array,3):ubound(array,3)) )
1121  rescaled_array(:,:,:) = 0.0
1122  is = hi%isc ; if (sym_stats) is = hi%isc-1
1123  js = hi%jsc ; if (sym_stats) js = hi%jsc-1
1124  do k=1,size(array,3) ; do j=js,hi%JecB ; do i=is,hi%IecB
1125  rescaled_array(i,j,k) = scale*array(i,j,k)
1126  enddo ; enddo ; enddo
1127  call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
1128  deallocate(rescaled_array)
1129  else
1130  call substats(hi, array, sym_stats, amean, amin, amax)
1131  endif
1132 
1133  if (is_root_pe()) &
1134  call chk_sum_msg("B-point:", amean, amin, amax, mesg, iounit)
1135  endif
1136 
1137  if (.not.writechksums) return
1138 
1139  hshift = default_shift
1140  if (present(haloshift)) hshift = haloshift
1141  if (hshift<0) hshift = hi%ied-hi%iec
1142 
1143  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
1144  hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
1145  write(0,*) 'chksum_B_3d: haloshift =',hshift
1146  write(0,*) 'chksum_B_3d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
1147  write(0,*) 'chksum_B_3d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
1148  call chksum_error(fatal,'Error in chksum_B_3d '//trim(mesg))
1149  endif
1150 
1151  bc0 = subchk(array, hi, 0, 0, scaling)
1152 
1153  sym = .false. ; if (present(symmetric)) sym = symmetric
1154 
1155  if ((hshift==0) .and. .not.sym) then
1156  if (is_root_pe()) call chk_sum_msg("B-point:", bc0, mesg, iounit)
1157  return
1158  endif
1159 
1160  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners
1161 
1162  if (do_corners) then
1163  if (sym) then
1164  bcsw = subchk(array, hi, -hshift-1, -hshift-1, scaling)
1165  bcse = subchk(array, hi, hshift, -hshift-1, scaling)
1166  bcnw = subchk(array, hi, -hshift-1, hshift, scaling)
1167  else
1168  bcsw = subchk(array, hi, -hshift-1, -hshift-1, scaling)
1169  bcse = subchk(array, hi, hshift, -hshift-1, scaling)
1170  bcnw = subchk(array, hi, -hshift-1, hshift, scaling)
1171  endif
1172  bcne = subchk(array, hi, hshift, hshift, scaling)
1173 
1174  if (is_root_pe()) &
1175  call chk_sum_msg("B-point:", bc0, bcsw, bcse, bcnw, bcne, mesg, iounit)
1176  else
1177  if (sym) then
1178  bcs = subchk(array, hi, 0, -hshift-1, scaling)
1179  bcw = subchk(array, hi, -hshift-1, 0, scaling)
1180  else
1181  bcs = subchk(array, hi, 0, -hshift, scaling)
1182  bcw = subchk(array, hi, -hshift, 0, scaling)
1183  endif
1184  bce = subchk(array, hi, hshift, 0, scaling)
1185  bcn = subchk(array, hi, 0, hshift, scaling)
1186 
1187  if (is_root_pe()) &
1188  call chk_sum_msg_nsew("B-point:", bc0, bcn, bcs, bce, bcw, mesg, iounit)
1189  endif
1190 
1191  contains
1192 
1193  integer function subchk(array, HI, di, dj, scale)
1194  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1195  real, dimension(HI%IsdB:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed
1196  integer, intent(in) :: di !< i- direction array shift for this checksum
1197  integer, intent(in) :: dj !< j- direction array shift for this checksum
1198  real, intent(in) :: scale !< A scaling factor for this array.
1199  integer :: i, j, k, bc
1200  subchk = 0
1201  ! This line deliberately uses the h-point computational domain.
1202  do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc+dj,hi%jec+dj ; do i=hi%isc+di,hi%iec+di
1203  bc = bitcount(abs(scale*array(i,j,k)))
1204  subchk = subchk + bc
1205  enddo ; enddo ; enddo
1206  call sum_across_pes(subchk)
1207  subchk=mod(subchk, bc_modulus)
1208  end function subchk
1209 
1210  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
1211  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1212  real, dimension(HI%IsdB:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed
1213  logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the
1214  !! full symmetric computational domain.
1215  real, intent(out) :: aMean, aMin, aMax
1216 
1217  integer :: i, j, k, n, IsB, JsB
1218 
1219  isb = hi%isc ; if (sym_stats) isb = hi%isc-1
1220  jsb = hi%jsc ; if (sym_stats) jsb = hi%jsc-1
1221 
1222  amin = array(hi%isc,hi%jsc,1) ; amax = amin
1223  do k=lbound(array,3),ubound(array,3) ; do j=jsb,hi%JecB ; do i=isb,hi%IecB
1224  amin = min(amin, array(i,j,k))
1225  amax = max(amax, array(i,j,k))
1226  enddo ; enddo ; enddo
1227  amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec,:))
1228  n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc) * size(array,3)
1229  call sum_across_pes(n)
1230  call min_across_pes(amin)
1231  call max_across_pes(amax)
1232  amean = amean / real(n)
1233  end subroutine substats
1234 
1235 end subroutine chksum_b_3d
1236 
1237 !> Checksums a 3d array staggered at C-grid u points.
1238 subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, &
1239  scale, logunit)
1240  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1241  real, dimension(HI%isdB:,HI%Jsd:,:), intent(in) :: array !< The array to be checksummed
1242  character(len=*), intent(in) :: mesg !< An identifying message
1243  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
1244  logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full
1245  !! symmetric computational domain.
1246  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
1247  real, optional, intent(in) :: scale !< A scaling factor for this array.
1248  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
1249 
1250  real, allocatable, dimension(:,:,:) :: rescaled_array
1251  real :: scaling
1252  integer :: iounit !< Log IO unit
1253  integer :: i, j, k, Is
1254  real :: aMean, aMin, aMax
1255  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
1256  integer :: bcN, bcS, bcE, bcW
1257  logical :: do_corners, sym, sym_stats
1258 
1259  if (checkfornans) then
1260  if (is_nan(array(hi%IscB:hi%IecB,hi%jsc:hi%jec,:))) &
1261  call chksum_error(fatal, 'NaN detected: '//trim(mesg))
1262 ! if (is_NaN(array)) &
1263 ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
1264  endif
1265 
1266  scaling = 1.0 ; if (present(scale)) scaling = scale
1267  iounit = error_unit; if(present(logunit)) iounit = logunit
1268  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
1269  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
1270 
1271  if (calculatestatistics) then
1272  if (present(scale)) then
1273  allocate( rescaled_array(lbound(array,1):ubound(array,1), &
1274  lbound(array,2):ubound(array,2), &
1275  lbound(array,3):ubound(array,3)) )
1276  rescaled_array(:,:,:) = 0.0
1277  is = hi%isc ; if (sym_stats) is = hi%isc-1
1278  do k=1,size(array,3) ; do j=hi%jsc,hi%jec ; do i=is,hi%IecB
1279  rescaled_array(i,j,k) = scale*array(i,j,k)
1280  enddo ; enddo ; enddo
1281  call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
1282  deallocate(rescaled_array)
1283  else
1284  call substats(hi, array, sym_stats, amean, amin, amax)
1285  endif
1286  if (is_root_pe()) &
1287  call chk_sum_msg("u-point:", amean, amin, amax, mesg, iounit)
1288  endif
1289 
1290  if (.not.writechksums) return
1291 
1292  hshift = default_shift
1293  if (present(haloshift)) hshift = haloshift
1294  if (hshift<0) hshift = hi%ied-hi%iec
1295 
1296  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
1297  hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
1298  write(0,*) 'chksum_u_3d: haloshift =',hshift
1299  write(0,*) 'chksum_u_3d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
1300  write(0,*) 'chksum_u_3d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
1301  call chksum_error(fatal,'Error in chksum_u_3d '//trim(mesg))
1302  endif
1303 
1304  bc0 = subchk(array, hi, 0, 0, scaling)
1305 
1306  sym = .false. ; if (present(symmetric)) sym = symmetric
1307 
1308  if ((hshift==0) .and. .not.sym) then
1309  if (is_root_pe()) call chk_sum_msg("u-point:", bc0, mesg, iounit)
1310  return
1311  endif
1312 
1313  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners
1314 
1315  if (hshift==0) then
1316  bcw = subchk(array, hi, -hshift-1, 0, scaling)
1317  if (is_root_pe()) call chk_sum_msg_w("u-point:", bc0, bcw, mesg, iounit)
1318  elseif (do_corners) then
1319  if (sym) then
1320  bcsw = subchk(array, hi, -hshift-1, -hshift, scaling)
1321  bcnw = subchk(array, hi, -hshift-1, hshift, scaling)
1322  else
1323  bcsw = subchk(array, hi, -hshift, -hshift, scaling)
1324  bcnw = subchk(array, hi, -hshift, hshift, scaling)
1325  endif
1326  bcse = subchk(array, hi, hshift, -hshift, scaling)
1327  bcne = subchk(array, hi, hshift, hshift, scaling)
1328 
1329  if (is_root_pe()) &
1330  call chk_sum_msg("u-point:", bc0, bcsw, bcse, bcnw, bcne, mesg, iounit)
1331  else
1332  bcs = subchk(array, hi, 0, -hshift, scaling)
1333  bce = subchk(array, hi, hshift, 0, scaling)
1334  if (sym) then
1335  bcw = subchk(array, hi, -hshift-1, 0, scaling)
1336  else
1337  bcw = subchk(array, hi, -hshift, 0, scaling)
1338  endif
1339  bcn = subchk(array, hi, 0, hshift, scaling)
1340 
1341  if (is_root_pe()) &
1342  call chk_sum_msg_nsew("u-point:", bc0, bcn, bcs, bce, bcw, mesg, iounit)
1343  endif
1344 
1345  contains
1346 
1347  integer function subchk(array, HI, di, dj, scale)
1348  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1349  real, dimension(HI%IsdB:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed
1350  integer, intent(in) :: di !< i- direction array shift for this checksum
1351  integer, intent(in) :: dj !< j- direction array shift for this checksum
1352  real, intent(in) :: scale !< A scaling factor for this array.
1353  integer :: i, j, k, bc
1354  subchk = 0
1355  ! This line deliberately uses the h-point computational domain.
1356  do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc+dj,hi%jec+dj ; do i=hi%isc+di,hi%iec+di
1357  bc = bitcount(abs(scale*array(i,j,k)))
1358  subchk = subchk + bc
1359  enddo ; enddo ; enddo
1360  call sum_across_pes(subchk)
1361  subchk=mod(subchk, bc_modulus)
1362  end function subchk
1363 
1364  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
1365  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1366  real, dimension(HI%IsdB:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed
1367  logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the
1368  !! full symmetric computational domain.
1369  real, intent(out) :: aMean, aMin, aMax
1370 
1371  integer :: i, j, k, n, IsB
1372 
1373  isb = hi%isc ; if (sym_stats) isb = hi%isc-1
1374 
1375  amin = array(hi%isc,hi%jsc,1) ; amax = amin
1376  do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc,hi%jec ; do i=isb,hi%IecB
1377  amin = min(amin, array(i,j,k))
1378  amax = max(amax, array(i,j,k))
1379  enddo ; enddo ; enddo
1380  ! This line deliberately uses the h-point computational domain.
1381  amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec,:))
1382  n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc) * size(array,3)
1383  call sum_across_pes(n)
1384  call min_across_pes(amin)
1385  call max_across_pes(amax)
1386  amean = amean / real(n)
1387  end subroutine substats
1388 
1389 end subroutine chksum_u_3d
1390 
1391 !> Checksums a 3d array staggered at C-grid v points.
1392 subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, &
1393  scale, logunit)
1394  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1395  real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed
1396  character(len=*), intent(in) :: mesg !< An identifying message
1397  integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0)
1398  logical, optional, intent(in) :: symmetric !< If true, do the checksums on the full
1399  !! symmetric computational domain.
1400  logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts
1401  real, optional, intent(in) :: scale !< A scaling factor for this array.
1402  integer, optional, intent(in) :: logunit !< IO unit for checksum logging
1403 
1404  real, allocatable, dimension(:,:,:) :: rescaled_array
1405  real :: scaling
1406  integer :: iounit !< Log IO unit
1407  integer :: i, j, k, Js
1408  integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift
1409  integer :: bcN, bcS, bcE, bcW
1410  real :: aMean, aMin, aMax
1411  logical :: do_corners, sym, sym_stats
1412 
1413  if (checkfornans) then
1414  if (is_nan(array(hi%isc:hi%iec,hi%JscB:hi%JecB,:))) &
1415  call chksum_error(fatal, 'NaN detected: '//trim(mesg))
1416 ! if (is_NaN(array)) &
1417 ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg))
1418  endif
1419 
1420  scaling = 1.0 ; if (present(scale)) scaling = scale
1421  iounit = error_unit; if(present(logunit)) iounit = logunit
1422  sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric
1423  if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif
1424 
1425  if (calculatestatistics) then
1426  if (present(scale)) then
1427  allocate( rescaled_array(lbound(array,1):ubound(array,1), &
1428  lbound(array,2):ubound(array,2), &
1429  lbound(array,3):ubound(array,3)) )
1430  rescaled_array(:,:,:) = 0.0
1431  js = hi%jsc ; if (sym_stats) js = hi%jsc-1
1432  do k=1,size(array,3) ; do j=js,hi%JecB ; do i=hi%isc,hi%iec
1433  rescaled_array(i,j,k) = scale*array(i,j,k)
1434  enddo ; enddo ; enddo
1435  call substats(hi, rescaled_array, sym_stats, amean, amin, amax)
1436  deallocate(rescaled_array)
1437  else
1438  call substats(hi, array, sym_stats, amean, amin, amax)
1439  endif
1440  if (is_root_pe()) &
1441  call chk_sum_msg("v-point:", amean, amin, amax, mesg, iounit)
1442  endif
1443 
1444  if (.not.writechksums) return
1445 
1446  hshift = default_shift
1447  if (present(haloshift)) hshift = haloshift
1448  if (hshift<0) hshift = hi%ied-hi%iec
1449 
1450  if ( hi%isc-hshift<hi%isd .or. hi%iec+hshift>hi%ied .or. &
1451  hi%jsc-hshift<hi%jsd .or. hi%jec+hshift>hi%jed ) then
1452  write(0,*) 'chksum_v_3d: haloshift =',hshift
1453  write(0,*) 'chksum_v_3d: isd,isc,iec,ied=',hi%isd,hi%isc,hi%iec,hi%ied
1454  write(0,*) 'chksum_v_3d: jsd,jsc,jec,jed=',hi%jsd,hi%jsc,hi%jec,hi%jed
1455  call chksum_error(fatal,'Error in chksum_v_3d '//trim(mesg))
1456  endif
1457 
1458  bc0 = subchk(array, hi, 0, 0, scaling)
1459 
1460  sym = .false. ; if (present(symmetric)) sym = symmetric
1461 
1462  if ((hshift==0) .and. .not.sym) then
1463  if (is_root_pe()) call chk_sum_msg("v-point:", bc0, mesg, iounit)
1464  return
1465  endif
1466 
1467  do_corners = .true. ; if (present(omit_corners)) do_corners = .not.omit_corners
1468 
1469  if (hshift==0) then
1470  bcs = subchk(array, hi, 0, -hshift-1, scaling)
1471  if (is_root_pe()) call chk_sum_msg_s("v-point:", bc0, bcs, mesg, iounit)
1472  elseif (do_corners) then
1473  if (sym) then
1474  bcsw = subchk(array, hi, -hshift, -hshift-1, scaling)
1475  bcse = subchk(array, hi, hshift, -hshift-1, scaling)
1476  else
1477  bcsw = subchk(array, hi, -hshift, -hshift, scaling)
1478  bcse = subchk(array, hi, hshift, -hshift, scaling)
1479  endif
1480  bcnw = subchk(array, hi, -hshift, hshift, scaling)
1481  bcne = subchk(array, hi, hshift, hshift, scaling)
1482 
1483  if (is_root_pe()) &
1484  call chk_sum_msg("v-point:", bc0, bcsw, bcse, bcnw, bcne, mesg, iounit)
1485  else
1486  if (sym) then
1487  bcs = subchk(array, hi, 0, -hshift-1, scaling)
1488  else
1489  bcs = subchk(array, hi, 0, -hshift, scaling)
1490  endif
1491  bce = subchk(array, hi, hshift, 0, scaling)
1492  bcw = subchk(array, hi, -hshift, 0, scaling)
1493  bcn = subchk(array, hi, 0, hshift, scaling)
1494 
1495  if (is_root_pe()) &
1496  call chk_sum_msg_nsew("v-point:", bc0, bcn, bcs, bce, bcw, mesg, iounit)
1497  endif
1498 
1499  contains
1500 
1501  integer function subchk(array, HI, di, dj, scale)
1502  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1503  real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed
1504  integer, intent(in) :: di !< i- direction array shift for this checksum
1505  integer, intent(in) :: dj !< j- direction array shift for this checksum
1506  real, intent(in) :: scale !< A scaling factor for this array.
1507  integer :: i, j, k, bc
1508  subchk = 0
1509  ! This line deliberately uses the h-point computational domain.
1510  do k=lbound(array,3),ubound(array,3) ; do j=hi%jsc+dj,hi%jec+dj ; do i=hi%isc+di,hi%iec+di
1511  bc = bitcount(abs(scale*array(i,j,k)))
1512  subchk = subchk + bc
1513  enddo ; enddo ; enddo
1514  call sum_across_pes(subchk)
1515  subchk=mod(subchk, bc_modulus)
1516  end function subchk
1517 
1518  !subroutine subStats(HI, array, mesg, sym_stats)
1519  subroutine substats(HI, array, sym_stats, aMean, aMin, aMax)
1520  type(hor_index_type), intent(in) :: HI !< A horizontal index type
1521  real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed
1522  logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the
1523  !! full symmetric computational domain.
1524  real, intent(out) :: aMean, aMin, aMax !< Mean/min/max of array over domain
1525 
1526  integer :: i, j, k, n, JsB
1527 
1528  jsb = hi%jsc ; if (sym_stats) jsb = hi%jsc-1
1529 
1530  amin = array(hi%isc,hi%jsc,1) ; amax = amin
1531  do k=lbound(array,3),ubound(array,3) ; do j=jsb,hi%JecB ; do i=hi%isc,hi%iec
1532  amin = min(amin, array(i,j,k))
1533  amax = max(amax, array(i,j,k))
1534  enddo ; enddo ; enddo
1535  ! This line deliberately uses the h-point computational domain.
1536  amean = reproducing_sum(array(hi%isc:hi%iec,hi%jsc:hi%jec,:))
1537  n = (1 + hi%jec - hi%jsc) * (1 + hi%iec - hi%isc) * size(array,3)
1538  call sum_across_pes(n)
1539  call min_across_pes(amin)
1540  call max_across_pes(amax)
1541  amean = amean / real(n)
1542  end subroutine substats
1543 
1544 end subroutine chksum_v_3d
1545 
1546 ! These are the older version of chksum that do not take the grid staggering
1547 ! into account.
1548 
1549 !> chksum1d does a checksum of a 1-dimensional array.
1550 subroutine chksum1d(array, mesg, start_i, end_i, compare_PEs)
1551  real, dimension(:), intent(in) :: array !< The array to be summed (index starts at 1).
1552  character(len=*), intent(in) :: mesg !< An identifying message.
1553  integer, optional, intent(in) :: start_i !< The starting index for the sum (default 1)
1554  integer, optional, intent(in) :: end_i !< The ending index for the sum (default all)
1555  logical, optional, intent(in) :: compare_PEs !< If true, compare across PEs instead of summing
1556  !! and list the root_PE value (default true)
1557 
1558  integer :: is, ie, i, bc, sum1, sum_bc
1559  real :: sum
1560  real, allocatable :: sum_here(:)
1561  logical :: compare
1562  integer :: pe_num ! pe number of the data
1563  integer :: nPEs ! Total number of processsors
1564 
1565  is = lbound(array,1) ; ie = ubound(array,1)
1566  if (present(start_i)) is = start_i
1567  if (present(end_i)) ie = end_i
1568  compare = .true. ; if (present(compare_pes)) compare = compare_pes
1569 
1570  sum = 0.0 ; sum_bc = 0
1571  do i=is,ie
1572  sum = sum + array(i)
1573  bc = bitcount(abs(array(i)))
1574  sum_bc = sum_bc + bc
1575  enddo
1576 
1577  pe_num = pe_here() + 1 - root_pe() ; npes = num_pes()
1578  allocate(sum_here(npes)) ; sum_here(:) = 0.0 ; sum_here(pe_num) = sum
1579  call sum_across_pes(sum_here,npes)
1580 
1581  sum1 = sum_bc
1582  call sum_across_pes(sum1)
1583 
1584  if (.not.compare) then
1585  sum = 0.0
1586  do i=1,npes ; sum = sum + sum_here(i) ; enddo
1587  sum_bc = sum1
1588  elseif (is_root_pe()) then
1589  if (sum1 /= npes*sum_bc) &
1590  write(0, '(A40," bitcounts do not match across PEs: ",I12,1X,I12)') &
1591  mesg, sum1, npes*sum_bc
1592  do i=1,npes ; if (sum /= sum_here(i)) then
1593  write(0, '(A40," PE ",i4," sum mismatches root_PE: ",3(ES22.13,1X))') &
1594  mesg, i, sum_here(i), sum, sum_here(i)-sum
1595  endif ; enddo
1596  endif
1597  deallocate(sum_here)
1598 
1599  if (is_root_pe()) &
1600  write(0,'(A50,1X,ES25.16,1X,I12)') mesg, sum, sum_bc
1601 
1602 end subroutine chksum1d
1603 
1604 ! These are the older version of chksum that do not take the grid staggering
1605 ! into account.
1606 
1607 !> chksum2d does a checksum of all data in a 2-d array.
1608 subroutine chksum2d(array, mesg)
1609 
1610  real, dimension(:,:) :: array !< The array to be checksummed
1611  character(len=*) :: mesg !< An identifying message
1612 
1613  integer :: xs,xe,ys,ye,i,j,sum1,bc
1614  real :: sum
1615 
1616  xs = lbound(array,1) ; xe = ubound(array,1)
1617  ys = lbound(array,2) ; ye = ubound(array,2)
1618 
1619  sum = 0.0 ; sum1 = 0
1620  do i=xs,xe ; do j=ys,ye
1621  bc = bitcount(abs(array(i,j)))
1622  sum1 = sum1 + bc
1623  enddo ; enddo
1624  call sum_across_pes(sum1)
1625 
1626  sum = reproducing_sum(array(:,:))
1627 
1628  if (is_root_pe()) &
1629  write(0,'(A50,1X,ES25.16,1X,I12)') mesg, sum, sum1
1630 ! write(0,'(A40,1X,Z16.16,1X,Z16.16,1X,ES25.16,1X,I12)') &
1631 ! mesg, sum, sum1, sum, sum1
1632 
1633 end subroutine chksum2d
1634 
1635 !> chksum3d does a checksum of all data in a 2-d array.
1636 subroutine chksum3d(array, mesg)
1637 
1638  real, dimension(:,:,:) :: array !< The array to be checksummed
1639  character(len=*) :: mesg !< An identifying message
1640 
1641  integer :: xs,xe,ys,ye,zs,ze,i,j,k, bc,sum1
1642  real :: sum
1643 
1644  xs = lbound(array,1) ; xe = ubound(array,1)
1645  ys = lbound(array,2) ; ye = ubound(array,2)
1646  zs = lbound(array,3) ; ze = ubound(array,3)
1647 
1648  sum = 0.0 ; sum1 = 0
1649  do i=xs,xe ; do j=ys,ye ; do k=zs,ze
1650  bc = bitcount(abs(array(i,j,k)))
1651  sum1 = sum1 + bc
1652  enddo ; enddo ; enddo
1653 
1654  call sum_across_pes(sum1)
1655  sum = reproducing_sum(array(:,:,:))
1656 
1657  if (is_root_pe()) &
1658  write(0,'(A50,1X,ES25.16,1X,I12)') mesg, sum, sum1
1659 ! write(0,'(A40,1X,Z16.16,1X,Z16.16,1X,ES25.16,1X,I12)') &
1660 ! mesg, sum, sum1, sum, sum1
1661 
1662 end subroutine chksum3d
1663 
1664 !> This function returns .true. if x is a NaN, and .false. otherwise.
1665 function is_nan_0d(x)
1666  real, intent(in) :: x !< The value to be checked for NaNs.
1667  logical :: is_nan_0d
1668 
1669  !is_NaN_0d = (((x < 0.0) .and. (x >= 0.0)) .or. &
1670  ! (.not.(x < 0.0) .and. .not.(x >= 0.0)))
1671  if (((x < 0.0) .and. (x >= 0.0)) .or. &
1672  (.not.(x < 0.0) .and. .not.(x >= 0.0))) then
1673  is_nan_0d = .true.
1674  else
1675  is_nan_0d = .false.
1676  endif
1677 
1678 end function is_nan_0d
1679 
1680 !> Returns .true. if any element of x is a NaN, and .false. otherwise.
1681 function is_nan_1d(x, skip_mpp)
1682  real, dimension(:), intent(in) :: x !< The array to be checked for NaNs.
1683  logical, optional, intent(in) :: skip_mpp !< If true, only check this array only
1684  !! on the local PE (default false).
1685  logical :: is_nan_1d
1686 
1687  integer :: i, n
1688  logical :: call_mpp
1689 
1690  n = 0
1691  do i = lbound(x,1), ubound(x,1)
1692  if (is_nan_0d(x(i))) n = n + 1
1693  enddo
1694  call_mpp = .true.
1695  if (present(skip_mpp)) call_mpp = .not.skip_mpp
1696 
1697  if (call_mpp) call sum_across_pes(n)
1698  is_nan_1d = .false.
1699  if (n>0) is_nan_1d = .true.
1700 
1701 end function is_nan_1d
1702 
1703 !> Returns .true. if any element of x is a NaN, and .false. otherwise.
1704 function is_nan_2d(x)
1705  real, dimension(:,:), intent(in) :: x !< The array to be checked for NaNs.
1706  logical :: is_nan_2d
1707 
1708  integer :: i, j, n
1709 
1710  n = 0
1711  do j = lbound(x,2), ubound(x,2) ; do i = lbound(x,1), ubound(x,1)
1712  if (is_nan_0d(x(i,j))) n = n + 1
1713  enddo ; enddo
1714  call sum_across_pes(n)
1715  is_nan_2d = .false.
1716  if (n>0) is_nan_2d = .true.
1717 
1718 end function is_nan_2d
1719 
1720 !> Returns .true. if any element of x is a NaN, and .false. otherwise.
1721 function is_nan_3d(x)
1722  real, dimension(:,:,:), intent(in) :: x !< The array to be checked for NaNs.
1723  logical :: is_nan_3d
1724 
1725  integer :: i, j, k, n
1726 
1727  n = 0
1728  do k = lbound(x,3), ubound(x,3)
1729  do j = lbound(x,2), ubound(x,2) ; do i = lbound(x,1), ubound(x,1)
1730  if (is_nan_0d(x(i,j,k))) n = n + 1
1731  enddo ; enddo
1732  enddo
1733  call sum_across_pes(n)
1734  is_nan_3d = .false.
1735  if (n>0) is_nan_3d = .true.
1736 
1737 end function is_nan_3d
1738 
1739 !> Write a message including the checksum of the non-shifted array
1740 subroutine chk_sum_msg1(fmsg, bc0, mesg, iounit)
1741  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
1742  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
1743  integer, intent(in) :: bc0 !< The bitcount of the non-shifted array
1744  integer, intent(in) :: iounit !< Checksum logger IO unit
1745 
1746  if (is_root_pe()) &
1747  write(iounit, '(A,1(A,I10,X),A)') fmsg, " c=", bc0, trim(mesg)
1748 end subroutine chk_sum_msg1
1749 
1750 !> Write a message including checksums of non-shifted and diagonally shifted arrays
1751 subroutine chk_sum_msg5(fmsg, bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit)
1752  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
1753  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
1754  integer, intent(in) :: bc0 !< The bitcount of the non-shifted array
1755  integer, intent(in) :: bcSW !< The bitcount for SW shifted array
1756  integer, intent(in) :: bcSE !< The bitcount for SE shifted array
1757  integer, intent(in) :: bcNW !< The bitcount for NW shifted array
1758  integer, intent(in) :: bcNE !< The bitcount for NE shifted array
1759  integer, intent(in) :: iounit !< Checksum logger IO unit
1760 
1761  if (is_root_pe()) write(iounit, '(A,5(A,I10,1X),A)') &
1762  fmsg, " c=", bc0, "sw=", bcsw, "se=", bcse, "nw=", bcnw, "ne=", bcne, trim(mesg)
1763 end subroutine chk_sum_msg5
1764 
1765 !> Write a message including checksums of non-shifted and laterally shifted arrays
1766 subroutine chk_sum_msg_nsew(fmsg, bc0, bcN, bcS, bcE, bcW, mesg, iounit)
1767  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
1768  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
1769  integer, intent(in) :: bc0 !< The bitcount of the non-shifted array
1770  integer, intent(in) :: bcN !< The bitcount for N shifted array
1771  integer, intent(in) :: bcS !< The bitcount for S shifted array
1772  integer, intent(in) :: bcE !< The bitcount for E shifted array
1773  integer, intent(in) :: bcW !< The bitcount for W shifted array
1774  integer, intent(in) :: iounit !< Checksum logger IO unit
1775 
1776  if (is_root_pe()) write(iounit, '(A,5(A,I10,1X),A)') &
1777  fmsg, " c=", bc0, "N=", bcn, "S=", bcs, "E=", bce, "W=", bcw, trim(mesg)
1778 end subroutine chk_sum_msg_nsew
1779 
1780 !> Write a message including checksums of non-shifted and southward shifted arrays
1781 subroutine chk_sum_msg_s(fmsg, bc0, bcS, mesg, iounit)
1782  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
1783  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
1784  integer, intent(in) :: bc0 !< The bitcount of the non-shifted array
1785  integer, intent(in) :: bcS !< The bitcount of the south-shifted array
1786  integer, intent(in) :: iounit !< Checksum logger IO unit
1787 
1788  if (is_root_pe()) write(iounit, '(A,2(A,I10,1X),A)') &
1789  fmsg, " c=", bc0, "S=", bcs, trim(mesg)
1790 end subroutine chk_sum_msg_s
1791 
1792 !> Write a message including checksums of non-shifted and westward shifted arrays
1793 subroutine chk_sum_msg_w(fmsg, bc0, bcW, mesg, iounit)
1794  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
1795  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
1796  integer, intent(in) :: bc0 !< The bitcount of the non-shifted array
1797  integer, intent(in) :: bcW !< The bitcount of the west-shifted array
1798  integer, intent(in) :: iounit !< Checksum logger IO unit
1799 
1800  if (is_root_pe()) write(iounit, '(A,2(A,I10,1X),A)') &
1801  fmsg, " c=", bc0, "W=", bcw, trim(mesg)
1802 end subroutine chk_sum_msg_w
1803 
1804 !> Write a message including checksums of non-shifted and southwestward shifted arrays
1805 subroutine chk_sum_msg2(fmsg, bc0, bcSW, mesg, iounit)
1806  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
1807  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
1808  integer, intent(in) :: bc0 !< The bitcount of the non-shifted array
1809  integer, intent(in) :: bcSW !< The bitcount of the southwest-shifted array
1810  integer, intent(in) :: iounit !< Checksum logger IO unit
1811 
1812  if (is_root_pe()) write(iounit, '(A,2(A,I9,1X),A)') &
1813  fmsg, " c=", bc0, "s/w=", bcsw, trim(mesg)
1814 end subroutine chk_sum_msg2
1815 
1816 !> Write a message including the global mean, maximum and minimum of an array
1817 subroutine chk_sum_msg3(fmsg, aMean, aMin, aMax, mesg, iounit)
1818  character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble
1819  character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller
1820  real, intent(in) :: aMean !< The mean value of the array
1821  real, intent(in) :: aMin !< The minimum value of the array
1822  real, intent(in) :: aMax !< The maximum value of the array
1823  integer, intent(in) :: iounit !< Checksum logger IO unit
1824 
1825  if (is_root_pe()) write(iounit, '(A,3(A,ES25.16,1X),A)') &
1826  fmsg, " mean=", amean, "min=", amin, "max=", amax, trim(mesg)
1827 end subroutine chk_sum_msg3
1828 
1829 !> MOM_checksums_init initializes the MOM_checksums module. As it happens, the
1830 !! only thing that it does is to log the version of this module.
1831 subroutine mom_checksums_init(param_file)
1832  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1833 ! This include declares and sets the variable "version".
1834 #include "version_variable.h"
1835  character(len=40) :: mdl = "MOM_checksums" ! This module's name.
1836 
1837  call log_version(param_file, mdl, version)
1838 
1839 end subroutine mom_checksums_init
1840 
1841 !> A wrapper for MOM_error used in the checksum code
1842 subroutine chksum_error(signal, message)
1843  ! Wrapper for MOM_error to help place specific break points in debuggers
1844  integer, intent(in) :: signal !< An error severity level, such as FATAL or WARNING
1845  character(len=*), intent(in) :: message !< An error message
1846  call mom_error(signal, message)
1847 end subroutine chksum_error
1848 
1849 !> Does a bitcount of a number by first casting to an integer and then using BTEST
1850 !! to check bit by bit
1851 integer function bitcount(x)
1852  real :: x !< Number to be bitcount
1853 
1854  ! Local variables
1855  integer(kind(x)) :: y !< Store the integer representation of the memory used by x
1856  integer :: bit
1857 
1858  bitcount = 0
1859  y = transfer(x,y)
1860 
1861  ! Fortran standard says that bit indexing start at 0
1862  do bit = 0, bit_size(y)-1
1863  if (btest(y,bit)) bitcount = bitcount+1
1864  enddo
1865 
1866 end function bitcount
1867 
1868 end module mom_checksums
mom_checksums::bchksum
Checksums an array (2d or 3d) staggered at corner points.
Definition: MOM_checksums.F90:53
mom_checksums::vchksum
Checksums an array (2d or 3d) staggered at C-grid v points.
Definition: MOM_checksums.F90:38
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_checksums::is_nan
Returns .true. if any element of x is a NaN, and .false. otherwise.
Definition: MOM_checksums.F90:73
mom_checksums::qchksum
This is an older interface that has been renamed Bchksum.
Definition: MOM_checksums.F90:58
mom_checksums::uvchksum
Checksums a pair velocity arrays (2d or 3d) staggered at C-grid locations.
Definition: MOM_checksums.F90:28
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_checksums::chksum
This is an older interface for 1-, 2-, or 3-D checksums.
Definition: MOM_checksums.F90:63
mom_hor_index
Defines the horizontal index type (hor_index_type) used for providing index ranges.
Definition: MOM_hor_index.F90:2
mom_coms
Interfaces to non-domain-oriented communication subroutines, including the MOM6 reproducing sums faci...
Definition: MOM_coms.F90:3
mom_checksums
Routines to calculate checksums of various array and vector types.
Definition: MOM_checksums.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_hor_index::hor_index_type
Container for horizontal index ranges for data, computational and global domains.
Definition: MOM_hor_index.F90:15
mom_checksums::bchksum_pair
Checksums a pair of arrays (2d or 3d) staggered at corner points.
Definition: MOM_checksums.F90:43
mom_checksums::hchksum_pair
Checksums a pair of arrays (2d or 3d) staggered at tracer points.
Definition: MOM_checksums.F90:23
mom_checksums::uchksum
Checksums an array (2d or 3d) staggered at C-grid u points.
Definition: MOM_checksums.F90:33
mom_checksums::chk_sum_msg
Write a message with either checksums or numerical statistics of arrays.
Definition: MOM_checksums.F90:68
mom_coms::reproducing_sum
Find an accurate and order-invariant sum of distributed 2d or 3d fields.
Definition: MOM_coms.F90:54
mom_checksums::hchksum
Checksums an array (2d or 3d) staggered at tracer points.
Definition: MOM_checksums.F90:48
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2