portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 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 | ||
| 21 | module MOM_array_transform | |
| 22 | ||
| 23 | use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit | |
| 24 | ||
| 25 | implicit none ; private | |
| 26 | ||
| 27 | public rotate_array | |
| 28 | public rotate_array_pair | |
| 29 | public rotate_vector | |
| 30 | public allocate_rotated_array | |
| 31 | public symmetric_sum | |
| 32 | public symmetric_sum_unit_tests | |
| 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. | |
| 37 | interface 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 | |
| 43 | end 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. | |
| 51 | interface 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 | |
| 55 | end 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. | |
| 65 | interface rotate_vector | |
| 66 | module procedure rotate_vector_real_2d | |
| 67 | module procedure rotate_vector_real_3d | |
| 68 | module procedure rotate_vector_real_4d | |
| 69 | end interface rotate_vector | |
| 70 | ||
| 71 | ||
| 72 | !> Allocate an array based on the rotated index map of an unrotated reference array. | |
| 73 | interface 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 | |
| 78 | end interface allocate_rotated_array | |
| 79 | ||
| 80 | ||
| 81 | !> Return a rotationally symmetric sum of the elements of an array. | |
| 82 | interface symmetric_sum | |
| 83 | module procedure symmetric_sum_1d, symmetric_sum_2d | |
| 84 | end interface symmetric_sum | |
| 85 | ||
| 86 | ||
| 87 | contains | |
| 88 | ||
| 89 | !> Rotate the elements of a 2d real array along first and second axes. | |
| 90 | 2360 | subroutine 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 | 2360 | m = size(A_in, 1) |
| 98 | 2360 | n = size(A_in, 2) |
| 99 | ||
| 100 | 4720 | select case (modulo(turns, 4)) |
| 101 | case(0) | |
| 102 | 20852970 | A(:,:) = A_in(:,:) |
| 103 | case(1) | |
| 104 | 0 | A(:,:) = transpose(A_in) |
| 105 | 0 | A(:,:) = A(n:1:-1, :) |
| 106 | case(2) | |
| 107 | 0 | A(:,:) = A_in(m:1:-1, n:1:-1) |
| 108 | case(3) | |
| 109 | 2360 | A(:,:) = transpose(A_in(m:1:-1, :)) |
| 110 | end select | |
| 111 | 2360 | end subroutine rotate_array_real_2d |
| 112 | ||
| 113 | ||
| 114 | !> Rotate the elements of a 3d real array along first and second axes. | |
| 115 | 31 | subroutine 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 | 2363 | do k = 1, size(A_in, 3) |
| 123 | 2363 | call rotate_array(A_in(:,:,k), turns, A(:,:,k)) |
| 124 | enddo | |
| 125 | 31 | end subroutine rotate_array_real_3d |
| 126 | ||
| 127 | ||
| 128 | !> Rotate the elements of a 4d real array along first and second axes. | |
| 129 | 0 | subroutine 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 | 0 | do n = 1, size(A_in, 4) |
| 137 | 0 | call rotate_array(A_in(:,:,:,n), turns, A(:,:,:,n)) |
| 138 | enddo | |
| 139 | 0 | end subroutine rotate_array_real_4d |
| 140 | ||
| 141 | ||
| 142 | !> Rotate the elements of a 2d integer array along first and second axes. | |
| 143 | 0 | subroutine 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 | 0 | m = size(A_in, 1) |
| 151 | 0 | n = size(A_in, 2) |
| 152 | ||
| 153 | 0 | select case (modulo(turns, 4)) |
| 154 | case(0) | |
| 155 | 0 | A(:,:) = A_in(:,:) |
| 156 | case(1) | |
| 157 | 0 | A(:,:) = transpose(A_in) |
| 158 | 0 | A(:,:) = A(n:1:-1, :) |
| 159 | case(2) | |
| 160 | 0 | A(:,:) = A_in(m:1:-1, n:1:-1) |
| 161 | case(3) | |
| 162 | 0 | A(:,:) = transpose(A_in(m:1:-1, :)) |
| 163 | end select | |
| 164 | 0 | end subroutine rotate_array_integer |
| 165 | ||
| 166 | ||
| 167 | !> Rotate the elements of a 2d logical array along first and second axes. | |
| 168 | 0 | subroutine 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 | 0 | m = size(A_in, 1) |
| 176 | 0 | n = size(A_in, 2) |
| 177 | ||
| 178 | 0 | select case (modulo(turns, 4)) |
| 179 | case(0) | |
| 180 | 0 | A(:,:) = A_in(:,:) |
| 181 | case(1) | |
| 182 | 0 | A(:,:) = transpose(A_in) |
| 183 | 0 | A(:,:) = A(n:1:-1, :) |
| 184 | case(2) | |
| 185 | 0 | A(:,:) = A_in(m:1:-1, n:1:-1) |
| 186 | case(3) | |
| 187 | 0 | A(:,:) = transpose(A_in(m:1:-1, :)) |
| 188 | end select | |
| 189 | 0 | end subroutine rotate_array_logical |
| 190 | ||
| 191 | ||
| 192 | !> Rotate the elements of a 2d real array pair along first and second axes. | |
| 193 | 0 | subroutine 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 | 0 | if (modulo(turns, 2) /= 0) then |
| 201 | 0 | call rotate_array(B_in, turns, A) |
| 202 | 0 | call rotate_array(A_in, turns, B) |
| 203 | else | |
| 204 | 0 | call rotate_array(A_in, turns, A) |
| 205 | 0 | call rotate_array(B_in, turns, B) |
| 206 | endif | |
| 207 | 0 | end subroutine rotate_array_pair_real_2d |
| 208 | ||
| 209 | ||
| 210 | !> Rotate the elements of a 3d real array pair along first and second axes. | |
| 211 | 0 | subroutine 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 | 0 | do k = 1, size(A_in, 3) |
| 221 | call rotate_array_pair(A_in(:,:,k), B_in(:,:,k), turns, & | |
| 222 | 0 | A(:,:,k), B(:,:,k)) |
| 223 | enddo | |
| 224 | 0 | end subroutine rotate_array_pair_real_3d |
| 225 | ||
| 226 | ||
| 227 | !> Rotate the elements of a 4d real array pair along first and second axes. | |
| 228 | 0 | subroutine 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 | 0 | if (modulo(turns, 2) /= 0) then |
| 236 | 0 | call rotate_array(B_in, turns, A) |
| 237 | 0 | call rotate_array(A_in, turns, B) |
| 238 | else | |
| 239 | 0 | call rotate_array(A_in, turns, A) |
| 240 | 0 | call rotate_array(B_in, turns, B) |
| 241 | endif | |
| 242 | 0 | end subroutine rotate_array_pair_integer |
| 243 | ||
| 244 | ||
| 245 | !> Rotate the elements of a 2d real vector along first and second axes. | |
| 246 | 0 | subroutine 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 | 0 | call rotate_array_pair(A_in, B_in, turns, A, B) |
| 254 | ||
| 255 | 0 | if (modulo(turns, 4) == 1 .or. modulo(turns, 4) == 2) & |
| 256 | 0 | A(:,:) = -A(:,:) |
| 257 | ||
| 258 | 0 | if (modulo(turns, 4) == 2 .or. modulo(turns, 4) == 3) & |
| 259 | 0 | B(:,:) = -B(:,:) |
| 260 | 0 | end subroutine rotate_vector_real_2d |
| 261 | ||
| 262 | ||
| 263 | !> Rotate the elements of a 3d real vector along first and second axes. | |
| 264 | 0 | subroutine 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 | 0 | do k = 1, size(A_in, 3) |
| 274 | 0 | call rotate_vector(A_in(:,:,k), B_in(:,:,k), turns, A(:,:,k), B(:,:,k)) |
| 275 | enddo | |
| 276 | 0 | end subroutine rotate_vector_real_3d |
| 277 | ||
| 278 | ||
| 279 | !> Rotate the elements of a 4d real vector along first and second axes. | |
| 280 | 0 | subroutine 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 | 0 | do n = 1, size(A_in, 4) |
| 290 | call rotate_vector(A_in(:,:,:,n), B_in(:,:,:,n), turns, & | |
| 291 | 0 | A(:,:,:,n), B(:,:,:,n)) |
| 292 | enddo | |
| 293 | 0 | end subroutine rotate_vector_real_4d |
| 294 | ||
| 295 | ||
| 296 | !> Allocate a 2d real array on the rotated index map of a reference array. | |
| 297 | 28 | subroutine 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 | 84 | ub(:) = ubound(A_in) |
| 307 | ||
| 308 | 28 | if (modulo(turns, 2) /= 0) then |
| 309 | 0 | allocate(A(lb(2):ub(2), lb(1):ub(1))) |
| 310 | else | |
| 311 | 28 | allocate(A(lb(1):ub(1), lb(2):ub(2))) |
| 312 | endif | |
| 313 | 28 | end subroutine allocate_rotated_array_real_2d |
| 314 | ||
| 315 | ||
| 316 | !> Allocate a 3d real array on the rotated index map of a reference array. | |
| 317 | 31 | subroutine 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 | 124 | ub(:) = ubound(A_in) |
| 327 | ||
| 328 | 31 | if (modulo(turns, 2) /= 0) then |
| 329 | 0 | allocate(A(lb(2):ub(2), lb(1):ub(1), lb(3):ub(3))) |
| 330 | else | |
| 331 | 31 | allocate(A(lb(1):ub(1), lb(2):ub(2), lb(3):ub(3))) |
| 332 | endif | |
| 333 | 31 | end subroutine allocate_rotated_array_real_3d |
| 334 | ||
| 335 | ||
| 336 | !> Allocate a 4d real array on the rotated index map of a reference array. | |
| 337 | 0 | subroutine 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 | 0 | ub(:) = ubound(A_in) |
| 347 | ||
| 348 | 0 | if (modulo(turns, 2) /= 0) then |
| 349 | 0 | allocate(A(lb(2):ub(2), lb(1):ub(1), lb(3):ub(3), lb(4):ub(4))) |
| 350 | else | |
| 351 | 0 | allocate(A(lb(1):ub(1), lb(2):ub(2), lb(3):ub(3), lb(4):ub(4))) |
| 352 | endif | |
| 353 | 0 | end subroutine allocate_rotated_array_real_4d |
| 354 | ||
| 355 | ||
| 356 | !> Allocate a 2d integer array on the rotated index map of a reference array. | |
| 357 | 0 | subroutine 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 | 0 | ub(:) = ubound(A_in) |
| 366 | ||
| 367 | 0 | if (modulo(turns, 2) /= 0) then |
| 368 | 0 | allocate(A(lb(2):ub(2), lb(1):ub(1))) |
| 369 | else | |
| 370 | 0 | allocate(A(lb(1):ub(1), lb(2):ub(2))) |
| 371 | endif | |
| 372 | 0 | end subroutine allocate_rotated_array_integer |
| 373 | ||
| 374 | ||
| 375 | !> Do a rotationally symmetric sum of a 1-d array | |
| 376 | 0 | function 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 | 0 | szi = size(field, 1) |
| 384 | 0 | szi_2 = szi / 2 ! Note that for an odd number szi_2 is rounded down. |
| 385 | 0 | sum = 0.0 |
| 386 | 0 | if (2*szi_2 < szi) sum = field(szi_2+1) |
| 387 | ! Add pairs of values, working from the inside out. | |
| 388 | 0 | do i=szi_2,1,-1 |
| 389 | 0 | sum = sum + (field(i) + field(szi+1-i)) |
| 390 | enddo | |
| 391 | 0 | end 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. | |
| 395 | 0 | recursive 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 | 0 | 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 | 0 | if ((szi == 1) .and. (szj == 1)) then |
| 408 | 0 | sum = field(1,1) |
| 409 | 0 | elseif ((szi == 2) .and. (szj == 2)) then |
| 410 | 0 | sum = (field(1,1) + field(2,2)) + (field(2,1) + field(1,2)) |
| 411 | 0 | 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 | 0 | ((field(1,1) + field(3,3)) + (field(3,1) + field(1,3))) |
| 414 | 0 | elseif (szi == 1) then |
| 415 | 0 | sum = symmetric_sum_1d(field(1,:)) |
| 416 | 0 | elseif (szj == 1) then |
| 417 | 0 | 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 | 0 | szi_2 = szi / 2 |
| 422 | 0 | szj_2 = szj / 2 |
| 423 | ||
| 424 | 0 | odd_i = (2*szi_2 < szi) ! This could be (modulo(szi,2) == 1) |
| 425 | 0 | 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 | 0 | if (odd_i .and. odd_j) then |
| 428 | 0 | ic = szi_2+1 ; jc = szj_2+1 ! The index of the central point |
| 429 | 0 | sum = field(ic,jc) |
| 430 | ! Add pairs of pairs of values, working from the inside out. | |
| 431 | 0 | do ij=1,min(szi_2,szj_2) |
| 432 | 0 | 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 | 0 | if (szi_2 > szj_2) then |
| 436 | 0 | do ij=szj_2+1,szi_2 |
| 437 | 0 | sum = sum + (field(ic-ij,jc) + field(ic+ij,jc)) |
| 438 | enddo | |
| 439 | 0 | elseif (szj_2 > szi_2) then |
| 440 | 0 | do ij=szi_2+1,szj_2 |
| 441 | 0 | sum = sum + (field(ic,jc-ij) + field(ic,jc+ij)) |
| 442 | enddo | |
| 443 | endif | |
| 444 | 0 | elseif (odd_i) then |
| 445 | 0 | sum = symmetric_sum_1d(field(szi_2+1,1:szj)) |
| 446 | 0 | elseif (odd_j) then |
| 447 | 0 | sum = symmetric_sum_1d(field(1:szi,szj_2+1)) |
| 448 | else | |
| 449 | 0 | sum = 0.0 |
| 450 | endif | |
| 451 | ||
| 452 | ! Find the sums in the four quadrants of the array. | |
| 453 | 0 | 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 | 0 | quad_sum(1,1) = symmetric_sum_2d(field(1:szi_2,1:szj_2)) |
| 456 | 0 | quad_sum(2,1) = symmetric_sum_2d(field(szi+1-szi_2:szi,1:szj_2)) |
| 457 | 0 | quad_sum(1,2) = symmetric_sum_2d(field(1:szi_2,szj+1-szj_2:szj)) |
| 458 | 0 | quad_sum(2,2) = symmetric_sum_2d(field(szi+1-szi_2:szi,szj+1-szj_2:szj)) |
| 459 | 0 | elseif (szi_2 > 1) then |
| 460 | 0 | quad_sum(1,1) = symmetric_sum_1d(field(1:szi_2,1)) |
| 461 | 0 | quad_sum(2,1) = symmetric_sum_1d(field(szi+1-szi_2:szi,1)) |
| 462 | 0 | quad_sum(1,2) = symmetric_sum_1d(field(1:szi_2,szj)) |
| 463 | 0 | quad_sum(2,2) = symmetric_sum_1d(field(szi+1-szi_2:szi,szj)) |
| 464 | 0 | elseif (szj_2 > 1) then |
| 465 | 0 | quad_sum(1,1) = symmetric_sum_1d(field(1,1:szj_2)) |
| 466 | 0 | quad_sum(2,1) = symmetric_sum_1d(field(szi,1:szj_2)) |
| 467 | 0 | quad_sum(1,2) = symmetric_sum_1d(field(1,szj+1-szj_2:szj)) |
| 468 | 0 | quad_sum(2,2) = symmetric_sum_1d(field(szi,szj+1-szj_2:szj)) |
| 469 | else | |
| 470 | 0 | quad_sum(1,1) = field(1,1) |
| 471 | 0 | quad_sum(2,1) = field(szi,1) |
| 472 | 0 | quad_sum(1,2) = field(1,szj) |
| 473 | 0 | quad_sum(2,2) = field(szi,szj) |
| 474 | endif | |
| 475 | ||
| 476 | 0 | sum = sum + ((quad_sum(1,1) + quad_sum(2,2)) + (quad_sum(2,1) + quad_sum(1,2))) |
| 477 | endif | |
| 478 | 0 | end 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. | |
| 482 | 0 | function 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 | 0 | szi = size(field, 1) ; szj = size(field, 2) |
| 492 | 0 | sum_abs_val = .false. ; if (present(abs_val)) sum_abs_val = abs_val |
| 493 | 0 | sum = 0.0 |
| 494 | 0 | if (sum_abs_val) then |
| 495 | 0 | do j=1,szj ; do i=1,szi |
| 496 | 0 | sum = sum + abs(field(i,j)) |
| 497 | enddo ; enddo | |
| 498 | else | |
| 499 | 0 | do j=1,szj ; do i=1,szi |
| 500 | 0 | sum = sum + field(i,j) |
| 501 | enddo ; enddo | |
| 502 | endif | |
| 503 | 0 | end function naive_sum_2d |
| 504 | ||
| 505 | ||
| 506 | !> Returns true if a unit test of the symmetric sums fails. | |
| 507 | 0 | logical 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 | 0 | fail = .false. |
| 526 | 0 | fail_message = "" |
| 527 | ||
| 528 | 0 | 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 | 0 | do j=1,sz ; do i=1,sz |
| 532 | 0 | 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 | 0 | if (modulo(i+j, 2) == 0) array(i,j) = -array(i,j) |
| 535 | enddo ; enddo | |
| 536 | 0 | call rotate_array_real_2d(array, 1, ar_90) |
| 537 | 0 | call rotate_array_real_2d(array, 2, ar_180) |
| 538 | 0 | call rotate_array_real_2d(array, 3, ar_270) |
| 539 | ||
| 540 | 0 | do n = 1, sz ; do m = 1, sz |
| 541 | 0 | sum(1) = symmetric_sum(array(1:n,1:m)) |
| 542 | 0 | sum(2) = symmetric_sum(ar_90(sz+1-m:sz,1:n)) |
| 543 | 0 | sum(3) = symmetric_sum(ar_180(sz+1-n:sz,sz+1-m:sz)) |
| 544 | 0 | sum(4) = symmetric_sum(ar_270(1:m,sz+1-n:sz)) |
| 545 | 0 | sum(5) = naive_sum_2d(array(1:n,1:m)) |
| 546 | 0 | abs_sum = naive_sum_2d(array(1:n,1:m), abs_val=.true.) |
| 547 | 0 | tol = 2.0 * abs_sum * epsilon(abs_sum) |
| 548 | 0 | 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 | 0 | n, m, sum(1), sum(5) - sum(1) |
| 551 | 0 | write(stdout,*) "Symmetric_sum_failure: "//trim(mesg) |
| 552 | 0 | write(stderr,*) "Symmetric_sum_failure: "//trim(mesg) |
| 553 | 0 | if (.not.fail) fail_message = mesg ! This is the first failed test. |
| 554 | 0 | fail = .true. |
| 555 | endif | |
| 556 | 0 | 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 | 0 | n, m, 90*(r-1), sum(1), sum(r) - sum(1) |
| 559 | 0 | write(stdout,*) "Symmetric_sum_failure: "//trim(mesg) |
| 560 | 0 | write(stderr,*) "Symmetric_sum_failure: "//trim(mesg) |
| 561 | 0 | if (.not.fail) fail_message = mesg ! This is the first failed test. |
| 562 | 0 | fail = .true. |
| 563 | endif ; enddo | |
| 564 | enddo ; enddo | |
| 565 | ||
| 566 | 0 | if (fail) then |
| 567 | 0 | write(stdout,*) "MOM_array_transform: One or more symmetric sum tests has failed." |
| 568 | 0 | write(stderr,*) "MOM_array_transform: One or more symmetric sum tests has failed." |
| 569 | else | |
| 570 | 0 | if (verbose) write(stdout,*) ("MOM_array_transform: All symmetric sum tests have passed.") |
| 571 | endif | |
| 572 | 0 | symmetric_sum_unit_tests = fail |
| 573 | ||
| 574 | 0 | end function symmetric_sum_unit_tests |
| 575 | ||
| 576 | end module MOM_array_transform |