MOM_array_transform.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!> Module for supporting the rotation of a field's index map.
6!! The implementation of each angle is described below.
7!!
8!! +90deg: B(i,j) = A(n-j,i)
9!! = transpose, then row reverse
10!! 180deg: B(i,j) = A(m-i,n-j)
11!! = row reversal + column reversal
12!! -90deg: B(i,j) = A(j,m-i)
13!! = row reverse, then transpose
14!!
15!! 90 degree rotations change the shape of the field, and are handled
16!! separately from 180 degree rotations.
17!!
18!! It also provides the symmetric_sum functions to do a rotationally invariant
19!! sum of the contents of a 1d or 2d array.
20
22
23use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit
24
25implicit none ; private
26
27public rotate_array
28public rotate_array_pair
29public rotate_vector
30public allocate_rotated_array
31public symmetric_sum
33
34
35!> Rotate the elements of an array to the rotated set of indices.
36!! Rotation is applied across the first and second axes of the array.
37interface rotate_array
38 module procedure rotate_array_real_2d
39 module procedure rotate_array_real_3d
40 module procedure rotate_array_real_4d
41 module procedure rotate_array_integer
42 module procedure rotate_array_logical
43end interface rotate_array
44
45
46!> Rotate a pair of arrays which map to a rotated set of indices.
47!! Rotation is applied across the first and second axes of the array.
48!! This rotation should be applied when one field is mapped onto the other.
49!! For example, a tracer indexed along u or v face points will map from one
50!! to the other after a quarter turn, and back onto itself after a half turn.
51interface rotate_array_pair
52 module procedure rotate_array_pair_real_2d
53 module procedure rotate_array_pair_real_3d
54 module procedure rotate_array_pair_integer
55end interface rotate_array_pair
56
57
58!> Rotate an array pair representing the components of a vector.
59!! Rotation is applied across the first and second axes of the array.
60!! This rotation should be applied when the fields satisfy vector
61!! transformation rules. For example, the u and v components of a velocity
62!! will map from one to the other for quarter turns, with a sign change in one
63!! component. A half turn will map elements onto themselves with sign changes
64!! in both components.
65interface rotate_vector
66 module procedure rotate_vector_real_2d
67 module procedure rotate_vector_real_3d
68 module procedure rotate_vector_real_4d
69end interface rotate_vector
70
71
72!> Allocate an array based on the rotated index map of an unrotated reference array.
73interface allocate_rotated_array
74 module procedure allocate_rotated_array_real_2d
75 module procedure allocate_rotated_array_real_3d
76 module procedure allocate_rotated_array_real_4d
77 module procedure allocate_rotated_array_integer
78end interface allocate_rotated_array
79
80
81!> Return a rotationally symmetric sum of the elements of an array.
82interface symmetric_sum
83 module procedure symmetric_sum_1d, symmetric_sum_2d
84end interface symmetric_sum
85
86
87contains
88
89!> Rotate the elements of a 2d real array along first and second axes.
90subroutine rotate_array_real_2d(A_in, turns, A)
91 real, intent(in) :: A_in(:,:) !< Unrotated array [arbitrary]
92 integer, intent(in) :: turns !< Number of quarter turns
93 real, intent(out) :: A(:,:) !< Rotated array [arbitrary]
94
95 integer :: m, n
96
97 m = size(a_in, 1)
98 n = size(a_in, 2)
99
100 select case (modulo(turns, 4))
101 case(0)
102 a(:,:) = a_in(:,:)
103 case(1)
104 a(:,:) = transpose(a_in)
105 a(:,:) = a(n:1:-1, :)
106 case(2)
107 a(:,:) = a_in(m:1:-1, n:1:-1)
108 case(3)
109 a(:,:) = transpose(a_in(m:1:-1, :))
110 end select
111end subroutine rotate_array_real_2d
112
113
114!> Rotate the elements of a 3d real array along first and second axes.
115subroutine rotate_array_real_3d(A_in, turns, A)
116 real, intent(in) :: A_in(:,:,:) !< Unrotated array [arbitrary]
117 integer, intent(in) :: turns !< Number of quarter turns
118 real, intent(out) :: A(:,:,:) !< Rotated array [arbitrary]
119
120 integer :: k
121
122 do k = 1, size(a_in, 3)
123 call rotate_array(a_in(:,:,k), turns, a(:,:,k))
124 enddo
125end subroutine rotate_array_real_3d
126
127
128!> Rotate the elements of a 4d real array along first and second axes.
129subroutine rotate_array_real_4d(A_in, turns, A)
130 real, intent(in) :: A_in(:,:,:,:) !< Unrotated array [arbitrary]
131 integer, intent(in) :: turns !< Number of quarter turns
132 real, intent(out) :: A(:,:,:,:) !< Rotated array [arbitrary]
133
134 integer :: n
135
136 do n = 1, size(a_in, 4)
137 call rotate_array(a_in(:,:,:,n), turns, a(:,:,:,n))
138 enddo
139end subroutine rotate_array_real_4d
140
141
142!> Rotate the elements of a 2d integer array along first and second axes.
143subroutine rotate_array_integer(A_in, turns, A)
144 integer, intent(in) :: A_in(:,:) !< Unrotated array
145 integer, intent(in) :: turns !< Number of quarter turns
146 integer, intent(out) :: A(:,:) !< Rotated array
147
148 integer :: m, n
149
150 m = size(a_in, 1)
151 n = size(a_in, 2)
152
153 select case (modulo(turns, 4))
154 case(0)
155 a(:,:) = a_in(:,:)
156 case(1)
157 a(:,:) = transpose(a_in)
158 a(:,:) = a(n:1:-1, :)
159 case(2)
160 a(:,:) = a_in(m:1:-1, n:1:-1)
161 case(3)
162 a(:,:) = transpose(a_in(m:1:-1, :))
163 end select
164end subroutine rotate_array_integer
165
166
167!> Rotate the elements of a 2d logical array along first and second axes.
168subroutine rotate_array_logical(A_in, turns, A)
169 logical, intent(in) :: A_in(:,:) !< Unrotated array
170 integer, intent(in) :: turns !< Number of quarter turns
171 logical, intent(out) :: A(:,:) !< Rotated array
172
173 integer :: m, n
174
175 m = size(a_in, 1)
176 n = size(a_in, 2)
177
178 select case (modulo(turns, 4))
179 case(0)
180 a(:,:) = a_in(:,:)
181 case(1)
182 a(:,:) = transpose(a_in)
183 a(:,:) = a(n:1:-1, :)
184 case(2)
185 a(:,:) = a_in(m:1:-1, n:1:-1)
186 case(3)
187 a(:,:) = transpose(a_in(m:1:-1, :))
188 end select
189end subroutine rotate_array_logical
190
191
192!> Rotate the elements of a 2d real array pair along first and second axes.
193subroutine rotate_array_pair_real_2d(A_in, B_in, turns, A, B)
194 real, intent(in) :: A_in(:,:) !< Unrotated scalar array pair [arbitrary]
195 real, intent(in) :: B_in(:,:) !< Unrotated scalar array pair [arbitrary]
196 integer, intent(in) :: turns !< Number of quarter turns
197 real, intent(out) :: A(:,:) !< Rotated scalar array pair [arbitrary]
198 real, intent(out) :: B(:,:) !< Rotated scalar array pair [arbitrary]
199
200 if (modulo(turns, 2) /= 0) then
201 call rotate_array(b_in, turns, a)
202 call rotate_array(a_in, turns, b)
203 else
204 call rotate_array(a_in, turns, a)
205 call rotate_array(b_in, turns, b)
206 endif
207end subroutine rotate_array_pair_real_2d
208
209
210!> Rotate the elements of a 3d real array pair along first and second axes.
211subroutine rotate_array_pair_real_3d(A_in, B_in, turns, A, B)
212 real, intent(in) :: A_in(:,:,:) !< Unrotated scalar array pair [arbitrary]
213 real, intent(in) :: B_in(:,:,:) !< Unrotated scalar array pair [arbitrary]
214 integer, intent(in) :: turns !< Number of quarter turns
215 real, intent(out) :: A(:,:,:) !< Rotated scalar array pair [arbitrary]
216 real, intent(out) :: B(:,:,:) !< Rotated scalar array pair [arbitrary]
217
218 integer :: k
219
220 do k = 1, size(a_in, 3)
221 call rotate_array_pair(a_in(:,:,k), b_in(:,:,k), turns, &
222 a(:,:,k), b(:,:,k))
223 enddo
224end subroutine rotate_array_pair_real_3d
225
226
227!> Rotate the elements of a 4d real array pair along first and second axes.
228subroutine rotate_array_pair_integer(A_in, B_in, turns, A, B)
229 integer, intent(in) :: A_in(:,:) !< Unrotated scalar array pair
230 integer, intent(in) :: B_in(:,:) !< Unrotated scalar array pair
231 integer, intent(in) :: turns !< Number of quarter turns
232 integer, intent(out) :: A(:,:) !< Rotated scalar array pair
233 integer, intent(out) :: B(:,:) !< Rotated scalar array pair
234
235 if (modulo(turns, 2) /= 0) then
236 call rotate_array(b_in, turns, a)
237 call rotate_array(a_in, turns, b)
238 else
239 call rotate_array(a_in, turns, a)
240 call rotate_array(b_in, turns, b)
241 endif
242end subroutine rotate_array_pair_integer
243
244
245!> Rotate the elements of a 2d real vector along first and second axes.
246subroutine rotate_vector_real_2d(A_in, B_in, turns, A, B)
247 real, intent(in) :: A_in(:,:) !< First component of unrotated vector [arbitrary]
248 real, intent(in) :: B_in(:,:) !< Second component of unrotated vector [arbitrary]
249 integer, intent(in) :: turns !< Number of quarter turns
250 real, intent(out) :: A(:,:) !< First component of rotated vector [arbitrary]
251 real, intent(out) :: B(:,:) !< Second component of unrotated vector [arbitrary]
252
253 call rotate_array_pair(a_in, b_in, turns, a, b)
254
255 if (modulo(turns, 4) == 1 .or. modulo(turns, 4) == 2) &
256 a(:,:) = -a(:,:)
257
258 if (modulo(turns, 4) == 2 .or. modulo(turns, 4) == 3) &
259 b(:,:) = -b(:,:)
260end subroutine rotate_vector_real_2d
261
262
263!> Rotate the elements of a 3d real vector along first and second axes.
264subroutine rotate_vector_real_3d(A_in, B_in, turns, A, B)
265 real, intent(in) :: A_in(:,:,:) !< First component of unrotated vector [arbitrary]
266 real, intent(in) :: B_in(:,:,:) !< Second component of unrotated vector [arbitrary]
267 integer, intent(in) :: turns !< Number of quarter turns
268 real, intent(out) :: A(:,:,:) !< First component of rotated vector [arbitrary]
269 real, intent(out) :: B(:,:,:) !< Second component of unrotated vector [arbitrary]
270
271 integer :: k
272
273 do k = 1, size(a_in, 3)
274 call rotate_vector(a_in(:,:,k), b_in(:,:,k), turns, a(:,:,k), b(:,:,k))
275 enddo
276end subroutine rotate_vector_real_3d
277
278
279!> Rotate the elements of a 4d real vector along first and second axes.
280subroutine rotate_vector_real_4d(A_in, B_in, turns, A, B)
281 real, intent(in) :: A_in(:,:,:,:) !< First component of unrotated vector [arbitrary]
282 real, intent(in) :: B_in(:,:,:,:) !< Second component of unrotated vector [arbitrary]
283 integer, intent(in) :: turns !< Number of quarter turns
284 real, intent(out) :: A(:,:,:,:) !< First component of rotated vector [arbitrary]
285 real, intent(out) :: B(:,:,:,:) !< Second component of unrotated vector [arbitrary]
286
287 integer :: n
288
289 do n = 1, size(a_in, 4)
290 call rotate_vector(a_in(:,:,:,n), b_in(:,:,:,n), turns, &
291 a(:,:,:,n), b(:,:,:,n))
292 enddo
293end subroutine rotate_vector_real_4d
294
295
296!> Allocate a 2d real array on the rotated index map of a reference array.
297subroutine allocate_rotated_array_real_2d(A_in, lb, turns, A)
298 ! NOTE: lb must be declared before A_in
299 integer, intent(in) :: lb(2) !< Lower index bounds of A_in
300 real, intent(in) :: A_in(lb(1):, lb(2):) !< Reference array [arbitrary]
301 integer, intent(in) :: turns !< Number of quarter turns
302 real, allocatable, intent(inout) :: A(:,:) !< Array on rotated index [arbitrary]
303
304 integer :: ub(2)
305
306 ub(:) = ubound(a_in)
307
308 if (modulo(turns, 2) /= 0) then
309 allocate(a(lb(2):ub(2), lb(1):ub(1)))
310 else
311 allocate(a(lb(1):ub(1), lb(2):ub(2)))
312 endif
314
315
316!> Allocate a 3d real array on the rotated index map of a reference array.
317subroutine allocate_rotated_array_real_3d(A_in, lb, turns, A)
318 ! NOTE: lb must be declared before A_in
319 integer, intent(in) :: lb(3) !< Lower index bounds of A_in
320 real, intent(in) :: A_in(lb(1):, lb(2):, lb(3):) !< Reference array [arbitrary]
321 integer, intent(in) :: turns !< Number of quarter turns
322 real, allocatable, intent(inout) :: A(:,:,:) !< Array on rotated index [arbitrary]
323
324 integer :: ub(3)
325
326 ub(:) = ubound(a_in)
327
328 if (modulo(turns, 2) /= 0) then
329 allocate(a(lb(2):ub(2), lb(1):ub(1), lb(3):ub(3)))
330 else
331 allocate(a(lb(1):ub(1), lb(2):ub(2), lb(3):ub(3)))
332 endif
334
335
336!> Allocate a 4d real array on the rotated index map of a reference array.
337subroutine allocate_rotated_array_real_4d(A_in, lb, turns, A)
338 ! NOTE: lb must be declared before A_in
339 integer, intent(in) :: lb(4) !< Lower index bounds of A_in
340 real, intent(in) :: A_in(lb(1):,lb(2):,lb(3):,lb(4):) !< Reference array [arbitrary]
341 integer, intent(in) :: turns !< Number of quarter turns
342 real, allocatable, intent(inout) :: A(:,:,:,:) !< Array on rotated index [arbitrary]
343
344 integer:: ub(4)
345
346 ub(:) = ubound(a_in)
347
348 if (modulo(turns, 2) /= 0) then
349 allocate(a(lb(2):ub(2), lb(1):ub(1), lb(3):ub(3), lb(4):ub(4)))
350 else
351 allocate(a(lb(1):ub(1), lb(2):ub(2), lb(3):ub(3), lb(4):ub(4)))
352 endif
354
355
356!> Allocate a 2d integer array on the rotated index map of a reference array.
357subroutine allocate_rotated_array_integer(A_in, lb, turns, A)
358 integer, intent(in) :: lb(2) !< Lower index bounds of A_in
359 integer, intent(in) :: A_in(lb(1):,lb(2):) !< Reference array
360 integer, intent(in) :: turns !< Number of quarter turns
361 integer, allocatable, intent(inout) :: A(:,:) !< Array on rotated index
362
363 integer :: ub(2)
364
365 ub(:) = ubound(a_in)
366
367 if (modulo(turns, 2) /= 0) then
368 allocate(a(lb(2):ub(2), lb(1):ub(1)))
369 else
370 allocate(a(lb(1):ub(1), lb(2):ub(2)))
371 endif
373
374
375!> Do a rotationally symmetric sum of a 1-d array
376function symmetric_sum_1d(field) result(sum)
377 real, dimension(1:), intent(in) :: field !< The field to sum in arbitrary units [A ~> a]
378 real :: sum !< The rotationally symmetric sum of the entries in field [A ~> a]
379
380 ! Local variables
381 integer :: i, szi, szi_2
382
383 szi = size(field, 1)
384 szi_2 = szi / 2 ! Note that for an odd number szi_2 is rounded down.
385 sum = 0.0
386 if (2*szi_2 < szi) sum = field(szi_2+1)
387 ! Add pairs of values, working from the inside out.
388 do i=szi_2,1,-1
389 sum = sum + (field(i) + field(szi+1-i))
390 enddo
391end function symmetric_sum_1d
392
393
394!> Do a rotationally symmetric sum of a 2-d array using a recursive "Union-Jack" pattern of addition.
395recursive function symmetric_sum_2d(field) result(sum)
396 real, dimension(1:,1:), intent(in) :: field !< The field to sum in arbitrary units [A ~> a]
397 real :: sum !< The rotationally symmetric sum of the entries in field [A ~> a]
398
399 ! Local variables
400 real :: quad_sum(2,2) ! The sums in each of the quadrants [A ~> a]
401 logical :: odd_i, odd_j
402 integer :: ij, szi, szj, szi_2, szj_2, ic, jc
403
404 szi = size(field, 1) ; szj = size(field, 2)
405 ! These 5 special cases are equivalent to the general case, but they reduce the use
406 ! of complicated logic for common simple cases.
407 if ((szi == 1) .and. (szj == 1)) then
408 sum = field(1,1)
409 elseif ((szi == 2) .and. (szj == 2)) then
410 sum = (field(1,1) + field(2,2)) + (field(2,1) + field(1,2))
411 elseif ((szi == 3) .and. (szj == 3)) then
412 sum = (field(2,2) + ((field(1,2) + field(3,2)) + (field(2,1) + field(2,3)))) + &
413 ((field(1,1) + field(3,3)) + (field(3,1) + field(1,3)))
414 elseif (szi == 1) then
415 sum = symmetric_sum_1d(field(1,:))
416 elseif (szj == 1) then
417 sum = symmetric_sum_1d(field(:,1))
418 else
419 ! This is the general case.
420 ! Note that for odd numbers szi_2 and szj_2 are rounded down.
421 szi_2 = szi / 2
422 szj_2 = szj / 2
423
424 odd_i = (2*szi_2 < szi) ! This could be (modulo(szi,2) == 1)
425 odd_j = (2*szj_2 < szj)
426 ! Start by finding the sums along the central axes if there are an odd number of points.
427 if (odd_i .and. odd_j) then
428 ic = szi_2+1 ; jc = szj_2+1 ! The index of the central point
429 sum = field(ic,jc)
430 ! Add pairs of pairs of values, working from the inside out.
431 do ij=1,min(szi_2,szj_2)
432 sum = sum + ((field(ic-ij,jc) + field(ic+ij,jc)) + (field(ic,jc-ij) + field(ic,jc+ij)))
433 enddo
434 ! Add extra pairs of values, working from the inside out.
435 if (szi_2 > szj_2) then
436 do ij=szj_2+1,szi_2
437 sum = sum + (field(ic-ij,jc) + field(ic+ij,jc))
438 enddo
439 elseif (szj_2 > szi_2) then
440 do ij=szi_2+1,szj_2
441 sum = sum + (field(ic,jc-ij) + field(ic,jc+ij))
442 enddo
443 endif
444 elseif (odd_i) then
445 sum = symmetric_sum_1d(field(szi_2+1,1:szj))
446 elseif (odd_j) then
447 sum = symmetric_sum_1d(field(1:szi,szj_2+1))
448 else
449 sum = 0.0
450 endif
451
452 ! Find the sums in the four quadrants of the array.
453 if ((szi_2 > 1) .and. (szj_2 > 1)) then
454 ! Use a recursive call to symmetric_sum_2d to determine the sums in the corner quadrants.
455 quad_sum(1,1) = symmetric_sum_2d(field(1:szi_2,1:szj_2))
456 quad_sum(2,1) = symmetric_sum_2d(field(szi+1-szi_2:szi,1:szj_2))
457 quad_sum(1,2) = symmetric_sum_2d(field(1:szi_2,szj+1-szj_2:szj))
458 quad_sum(2,2) = symmetric_sum_2d(field(szi+1-szi_2:szi,szj+1-szj_2:szj))
459 elseif (szi_2 > 1) then
460 quad_sum(1,1) = symmetric_sum_1d(field(1:szi_2,1))
461 quad_sum(2,1) = symmetric_sum_1d(field(szi+1-szi_2:szi,1))
462 quad_sum(1,2) = symmetric_sum_1d(field(1:szi_2,szj))
463 quad_sum(2,2) = symmetric_sum_1d(field(szi+1-szi_2:szi,szj))
464 elseif (szj_2 > 1) then
465 quad_sum(1,1) = symmetric_sum_1d(field(1,1:szj_2))
466 quad_sum(2,1) = symmetric_sum_1d(field(szi,1:szj_2))
467 quad_sum(1,2) = symmetric_sum_1d(field(1,szj+1-szj_2:szj))
468 quad_sum(2,2) = symmetric_sum_1d(field(szi,szj+1-szj_2:szj))
469 else
470 quad_sum(1,1) = field(1,1)
471 quad_sum(2,1) = field(szi,1)
472 quad_sum(1,2) = field(1,szj)
473 quad_sum(2,2) = field(szi,szj)
474 endif
475
476 sum = sum + ((quad_sum(1,1) + quad_sum(2,2)) + (quad_sum(2,1) + quad_sum(1,2)))
477 endif
478end function symmetric_sum_2d
479
480
481!> Do a naive non-rotationally symmetric sum of a 2-d array. This function is only here for testing.
482function naive_sum_2d(field, abs_val) result(sum)
483 real, dimension(1:,1:), intent(in) :: field !< The field to sum in arbitrary units [A ~> a]
484 logical, optional, intent(in) :: abs_val !< If present and true, sum the absolute values
485 real :: sum !< The rotation dependent sum of the entries in field [A ~> a]
486
487 ! Local variables
488 logical :: sum_abs_val
489 integer :: i, j, szi, szj
490
491 szi = size(field, 1) ; szj = size(field, 2)
492 sum_abs_val = .false. ; if (present(abs_val)) sum_abs_val = abs_val
493 sum = 0.0
494 if (sum_abs_val) then
495 do j=1,szj ; do i=1,szi
496 sum = sum + abs(field(i,j))
497 enddo ; enddo
498 else
499 do j=1,szj ; do i=1,szi
500 sum = sum + field(i,j)
501 enddo ; enddo
502 endif
503end function naive_sum_2d
504
505
506!> Returns true if a unit test of the symmetric sums fails.
507logical function symmetric_sum_unit_tests(verbose)
508 ! Arguments
509 logical, intent(in) :: verbose !< If true, write results to stdout
510 ! Local variables
511 character(len=120) :: fail_message !< Blank or a description of the first failed test.
512 integer, parameter :: sz=13 ! The maximum size of the test arrays
513 real :: array(sz,sz) ! An array of inexact real values for testing in arbitrary units [A]
514 real :: ar_90(sz,sz) ! Array rotated by 90 degrees in arbitrary units [A]
515 real :: ar_180(sz,sz) ! Array rotated by 180 degrees in arbitrary units [A]
516 real :: ar_270(sz,sz) ! Array rotated by 270 degrees in arbitrary units [A]
517 real :: sum(5) ! Different versions of sums over a sub-array [A]
518 real :: abs_sum ! The sum of the absolute values of the array [A]
519 real :: tol ! The tolerance for an inexact test [A]
520
521 character(len=120) :: mesg
522 integer :: i, j, n, m, r
523 logical :: fail
524
525 fail = .false.
526 fail_message = ""
527
528 if (verbose) write(stdout,*) '==== MOM_array_transform: symmetric_sum_unit_tests ===='
529
530 ! Fill the array with real numbers that can not be represented exactly.
531 do j=1,sz ; do i=1,sz
532 array(i,j) = 1.0 / (2.0*(j*sz + i) + 1.0)
533 ! Combining positive and negative numbers amplifies differences from the order of arithmetic.
534 if (modulo(i+j, 2) == 0) array(i,j) = -array(i,j)
535 enddo ; enddo
536 call rotate_array_real_2d(array, 1, ar_90)
537 call rotate_array_real_2d(array, 2, ar_180)
538 call rotate_array_real_2d(array, 3, ar_270)
539
540 do n = 1, sz ; do m = 1, sz
541 sum(1) = symmetric_sum(array(1:n,1:m))
542 sum(2) = symmetric_sum(ar_90(sz+1-m:sz,1:n))
543 sum(3) = symmetric_sum(ar_180(sz+1-n:sz,sz+1-m:sz))
544 sum(4) = symmetric_sum(ar_270(1:m,sz+1-n:sz))
545 sum(5) = naive_sum_2d(array(1:n,1:m))
546 abs_sum = naive_sum_2d(array(1:n,1:m), abs_val=.true.)
547 tol = 2.0 * abs_sum * epsilon(abs_sum)
548 if (abs(sum(1) - sum(5)) > tol) then
549 write(mesg,'(i0," x ",i0," symmetric vs naive sum, sum=",ES13.5," diff=",ES13.5)') &
550 n, m, sum(1), sum(5) - sum(1)
551 write(stdout,*) "Symmetric_sum_failure: "//trim(mesg)
552 write(stderr,*) "Symmetric_sum_failure: "//trim(mesg)
553 if (.not.fail) fail_message = mesg ! This is the first failed test.
554 fail = .true.
555 endif
556 do r = 2, 4 ; if (abs(sum(1) - sum(r)) > 0.0) then
557 write(mesg,'(i0," x ",i0," with ",i0," degree rotation, sum=",ES13.5," diff=",ES13.5)') &
558 n, m, 90*(r-1), sum(1), sum(r) - sum(1)
559 write(stdout,*) "Symmetric_sum_failure: "//trim(mesg)
560 write(stderr,*) "Symmetric_sum_failure: "//trim(mesg)
561 if (.not.fail) fail_message = mesg ! This is the first failed test.
562 fail = .true.
563 endif ; enddo
564 enddo ; enddo
565
566 if (fail) then
567 write(stdout,*) "MOM_array_transform: One or more symmetric sum tests has failed."
568 write(stderr,*) "MOM_array_transform: One or more symmetric sum tests has failed."
569 else
570 if (verbose) write(stdout,*) ("MOM_array_transform: All symmetric sum tests have passed.")
571 endif
573
574end function symmetric_sum_unit_tests
575
576end module mom_array_transform