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