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