diff --git a/CMakeLists.txt b/CMakeLists.txt index ad7564cc7b..4a880284dd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -663,8 +663,10 @@ if(UNIT_TESTS) test_fms/diag_manager/test_flush_nc_file.F90 test_fms/diag_manager/test_diag_out_yaml.F90 test_fms/diag_manager/test_reduction_methods.F90 + test_fms/diag_manager/test_diag_generalized_indices.F90 test_fms/diag_manager/testing_utils.F90 test_fms/diag_manager/test_diag_diurnal.F90 + test_fms/diag_manager/check_diag_generalized_indices.F90 test_fms/diag_manager/check_time_none.F90 test_fms/diag_manager/check_time_min.F90 test_fms/diag_manager/check_time_max.F90 diff --git a/diag_manager/diag_manager.F90 b/diag_manager/diag_manager.F90 index 40617ec19c..1594cc481e 100644 --- a/diag_manager/diag_manager.F90 +++ b/diag_manager/diag_manager.F90 @@ -1683,12 +1683,12 @@ END FUNCTION send_data_3d LOGICAL FUNCTION diag_send_data(diag_field_id, field, time, is_in, js_in, ks_in, & & mask, rmask, ie_in, je_in, ke_in, weight, err_msg) INTEGER, INTENT(in) :: diag_field_id - CLASS(*), DIMENSION(:,:,:), INTENT(in),TARGET,CONTIGUOUS :: field + CLASS(*), DIMENSION(:,:,:), INTENT(in),TARGET :: field CLASS(*), INTENT(in), OPTIONAL :: weight TYPE (time_type), INTENT(in), OPTIONAL :: time INTEGER, INTENT(in), OPTIONAL :: is_in, js_in, ks_in,ie_in,je_in, ke_in - LOGICAL, DIMENSION(:,:,:), INTENT(in), OPTIONAL, contiguous, target :: mask - CLASS(*), DIMENSION(:,:,:), INTENT(in), OPTIONAL, contiguous, target :: rmask + LOGICAL, DIMENSION(:,:,:), INTENT(in), OPTIONAL, target :: mask + CLASS(*), DIMENSION(:,:,:), INTENT(in), OPTIONAL, target :: rmask CHARACTER(len=*), INTENT(out), OPTIONAL :: err_msg REAL :: weight1 diff --git a/test_fms/diag_manager/Makefile.am b/test_fms/diag_manager/Makefile.am index c3b8d7d944..a0ca6427bc 100644 --- a/test_fms/diag_manager/Makefile.am +++ b/test_fms/diag_manager/Makefile.am @@ -34,7 +34,7 @@ check_PROGRAMS = test_diag_manager test_diag_manager_time \ check_time_pow check_time_rms check_subregional test_cell_measures test_var_masks \ check_var_masks test_multiple_send_data test_diag_out_yaml test_output_every_freq \ test_dm_weights test_prepend_date test_ens_runs test_diag_multi_file test_diag_attribute_add \ - check_new_file_freq test_zbounds_limits test_multiple_zbounds + check_new_file_freq test_zbounds_limits test_multiple_zbounds test_diag_generalized_indices check_diag_generalized_indices # This is the source code for the test. test_output_every_freq_SOURCES = test_output_every_freq.F90 @@ -71,6 +71,8 @@ test_diag_attribute_add_SOURCES = test_diag_attribute_add.F90 check_new_file_freq_SOURCES = check_new_file_freq.F90 test_zbounds_limits_SOURCES = test_zbounds_limits.F90 test_multiple_zbounds_SOURCES = test_multiple_zbounds.F90 +test_diag_generalized_indices_SOURCES = testing_utils.F90 test_diag_generalized_indices.F90 +check_diag_generalized_indices_SOURCES = testing_utils.F90 check_diag_generalized_indices.F90 TEST_EXTENSIONS = .sh SH_LOG_DRIVER = env AM_TAP_AWK='$(AWK)' $(SHELL) \ diff --git a/test_fms/diag_manager/check_diag_generalized_indices.F90 b/test_fms/diag_manager/check_diag_generalized_indices.F90 new file mode 100644 index 0000000000..8fc947ac49 --- /dev/null +++ b/test_fms/diag_manager/check_diag_generalized_indices.F90 @@ -0,0 +1,99 @@ +!*********************************************************************** +!* Apache License 2.0 +!* +!* This file is part of the GFDL Flexible Modeling System (FMS). +!* +!* Licensed under the Apache License, Version 2.0 (the "License"); +!* you may not use this file except in compliance with the License. +!* You may obtain a copy of the License at +!* +!* http://www.apache.org/licenses/LICENSE-2.0 +!* +!* FMS is distributed in the hope that it will be useful, but WITHOUT +!* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied; +!* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A +!* PARTICULAR PURPOSE. See the License for the specific language +!* governing permissions and limitations under the License. +!*********************************************************************** + +!> @brief Checker for test_generalized_indices output. +!! Verifies permuted-axis variables match identity variables under axis permutations +program check_diag_generalized_indices + use fms_mod, only: fms_init, fms_end, string + use testing_utils, only: check_perm + use fms2_io_mod, only: FmsNetcdfFile_t, read_data, open_file, close_file, get_global_attribute + use mpp_mod, only: mpp_error, FATAL, mpp_pe + use platform_mod, only: r8_kind + + implicit none + + type(FmsNetcdfFile_t) :: fileobj + integer :: nx, ny, nz + integer :: i, ierr + + real(kind=r8_kind), allocatable :: var2_id(:,:) ! (x,y) + real(kind=r8_kind), allocatable :: var2_yx(:,:) ! (y,x) + real(kind=r8_kind), allocatable :: var3_id(:,:,:) ! (x,y,z) + real(kind=r8_kind), allocatable :: var3_zx(:,:,:) ! (z,y,x) + real(kind=r8_kind), allocatable :: var3_yzx(:,:,:) ! (y,z,x) + real(kind=r8_kind), allocatable :: var3_zxy(:,:,:) ! (z,x,y) + + call fms_init() + + nx = 96 + ny = 96 + nz = 5 + + if (.not. open_file(fileobj, "test_gen.nc", "read")) & + call mpp_error(FATAL, "unable to open test_gen.nc") + + call check_global_attribute(fileobj, "test_generalized_indices") + + allocate(var2_id(nx,ny), var2_yx(ny,nx)) + allocate(var3_id(nx,ny,nz), var3_zx(nz,ny,nx), var3_yzx(ny,nz,nx), var3_zxy(nz,nx,ny)) + + ! Output every 6 hours over 48 hours => 8 records + do i = 1, 8 + var2_id = -999._r8_kind + var2_yx = -999._r8_kind + var3_id = -999._r8_kind + var3_zx = -999._r8_kind + var3_yzx = -999._r8_kind + var3_zxy = -999._r8_kind + + print *, "Checking var2_yx vs var2_id - time_level:", i + call read_data(fileobj, "var2_id", var2_id, unlim_dim_level=i) + call read_data(fileobj, "var2_yx", var2_yx, unlim_dim_level=i) + call check_perm(var2_id, var2_yx, [2,1], ierr) + + print *, "Checking var3_zx vs var3_id - time_level:", i + call read_data(fileobj, "var3_id", var3_id, unlim_dim_level=i) + call read_data(fileobj, "var3_zx", var3_zx, unlim_dim_level=i) + call check_perm(var3_id, var3_zx, [3,2,1], ierr) + + print *, "Checking var3_yzx vs var3_id - time_level:", i + call read_data(fileobj, "var3_yzx", var3_yzx, unlim_dim_level=i) + call check_perm(var3_id, var3_yzx, [2,3,1], ierr) + + print *, "Checking var3_zxy vs var3_id - time_level:", i + call read_data(fileobj, "var3_zxy", var3_zxy, unlim_dim_level=i) + call check_perm(var3_id, var3_zxy, [3,1,2], ierr) + enddo + + call close_file(fileobj) + call fms_end() + +contains + + subroutine check_global_attribute(fileobj, expected_title) + type(FmsNetcdfFile_t), intent(in) :: fileobj + character(len=*), intent(in) :: expected_title + + character(len=100) :: attribute_value + + call get_global_attribute(fileobj, "title", attribute_value) + if (trim(attribute_value) .ne. trim(expected_title)) then + call mpp_error(FATAL, "Global attribute 'title' not expected value.") + endif + end subroutine check_global_attribute +end program check_diag_generalized_indices diff --git a/test_fms/diag_manager/test_diag_generalized_indices.F90 b/test_fms/diag_manager/test_diag_generalized_indices.F90 new file mode 100644 index 0000000000..91b692ce8c --- /dev/null +++ b/test_fms/diag_manager/test_diag_generalized_indices.F90 @@ -0,0 +1,255 @@ +!*********************************************************************** +!* Apache License 2.0 +!* +!* This file is part of the GFDL Flexible Modeling System (FMS). +!* +!* Licensed under the Apache License, Version 2.0 (the "License"); +!* you may not use this file except in compliance with the License. +!* You may obtain a copy of the License at +!* +!* http://www.apache.org/licenses/LICENSE-2.0 +!* +!* FMS is distributed in the hope that it will be useful, but WITHOUT +!* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied; +!* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A +!* PARTICULAR PURPOSE. See the License for the specific language +!* governing permissions and limitations under the License. +!*********************************************************************** + +!> @brief Test generalized axis permutations for send_data. +!! Applies predefined permutations to canonical (x,y,z,w) storage +!! and verifies consistency between data layout and axis metadata. +!! Assumes default configuration parameters: test_normal + no_mask. +program test_diag_generalized_indices + use fms_mod, only: fms_init, fms_end + use testing_utils, only: allocate_buffer, permute + use platform_mod, only: r8_kind + use mpp_mod, only: FATAL, mpp_error, mpp_npes, mpp_pe, mpp_root_pe + use time_manager_mod, only: time_type, set_calendar_type, set_date, JULIAN, set_time, OPERATOR(+) + use diag_manager_mod, only: diag_manager_init, diag_manager_end, diag_axis_init, register_diag_field, & + diag_send_complete, diag_manager_set_time_end, send_data + use mpp_domains_mod, only: domain2d, mpp_define_domains, mpp_define_io_domain, mpp_get_compute_domain + + implicit none + + integer :: nx, ny, nz, nw + integer :: layout(2), io_layout(2) + type(domain2d) :: Domain + integer :: isc, iec, jsc, jec + integer :: nhalox, nhaloy + integer :: ntimes, i + type(time_type) :: Time, Time_step + real(r8_kind) :: missing_value + integer :: ierr + + ! Axes + integer :: id_x, id_y, id_z, id_w + integer :: axis(4) + + ! Data + real(r8_kind), allocatable :: cdata(:,:,:,:) ! canonical storage: (x,y,z,w) + + ! Permutation test + integer, parameter :: LAYOUT_XY = 1 + integer, parameter :: LAYOUT_YX = 2 + integer, parameter :: LAYOUT_ZX = 3 + integer, parameter :: LAYOUT_YZX = 4 + integer, parameter :: LAYOUT_ZXY = 5 + integer, parameter :: PERM_TABLE(3,5) = reshape([ & + 1,2,3, & ! XY (identity) + 2,1,3, & ! YX + 3,2,1, & ! ZX + 2,3,1, & ! YZX + 3,1,2 & ! ZXY + ], [3,5]) + + integer :: id_var2_id, id_var2_yx + integer :: id_var3_id, id_var3_zx, id_var3_yzx, id_var3_zxy + + call fms_init + call set_calendar_type(JULIAN) + call diag_manager_init + + Time = set_date(2,1,1,0,0,0) + Time_step = set_time(3600,0) ! 1 hour + + nx = 96 + ny = 96 + nz = 5 + nw = 2 + ntimes = 48 + + nhalox = 2 + nhaloy = 2 + layout = (/1, mpp_npes()/) + io_layout = (/1, 1/) + + ! Domain + call mpp_define_domains((/1,nx,1,ny/), layout, Domain, name='2D domain', symmetry=.true., & + xhalo=nhalox, yhalo=nhaloy) + call mpp_define_io_domain(Domain, io_layout) + call mpp_get_compute_domain(Domain, isc, iec, jsc, jec) + + ! Data allocation + init (canonical x,y storage) + cdata = allocate_buffer(isc, iec, jsc, jec, nz, nw) + call init_buffer(cdata, isc, iec, jsc, jec, 0) + + ! Axes + id_x = diag_axis_init('x', real((/(i, i=1,nx)/), kind=r8_kind), 'point_E', 'x', long_name='point_E', Domain2=Domain) + id_y = diag_axis_init('y', real((/(i, i=1,ny)/), kind=r8_kind), 'point_N', 'y', long_name='point_N', Domain2=Domain) + id_z = diag_axis_init('z', real((/(i, i=1,nz)/), kind=r8_kind), 'point_Z', 'z', long_name='point_Z') + id_w = diag_axis_init('w', real((/(i, i=1,nw)/), kind=r8_kind), 'point_W', 'n', long_name='point_W') + axis = [id_x, id_y, id_z, id_w] + + missing_value = -666._r8_kind + + ! Register permuted diagnostic fields ONCE + id_var2_id = register_perm_diag_field('var2_id', 'Var2d id', axis(1:2), LAYOUT_XY) + id_var3_id = register_perm_diag_field('var3_id', 'Var3d id', axis(1:3), LAYOUT_XY) + + id_var2_yx = register_perm_diag_field('var2_yx', 'Var2d yx', axis(1:2), LAYOUT_YX) + id_var3_zx = register_perm_diag_field('var3_zx', 'Var3d zx', axis(1:3), LAYOUT_ZX) + + id_var3_yzx = register_perm_diag_field('var3_yzx', 'Var3d yzx', axis(1:3), LAYOUT_YZX) + id_var3_zxy = register_perm_diag_field('var3_zxy', 'Var3d zxy', axis(1:3), LAYOUT_ZXY) + + if (mpp_pe() == mpp_root_pe()) then + print *, "Testing generalized indices in default mode (test_normal + no_mask)" + print *, " canonical storage is (x,y,z,w)" + print *, " sending:" + print *, " var2_id with axes (x,y)" + print *, " var2_yx with axes (y,x)" + print *, " var3_id with axes (x,y,z)" + print *, " var3_zx with axes (z,y,x)" + print *, " var3_yzx with axes (y,z,x)" + print *, " var3_zxy with axes (z,x,y)" + end if + + call diag_manager_set_time_end(set_date(2,1,3,0,0,0)) + + do i = 1, ntimes + Time = Time + Time_step + call set_buffer(cdata, i) + + ! Identity: axes (x,y) / (x,y,z) with canonical storage + call send_perm_data(id_var2_id, cdata, PERM_TABLE(1:2, LAYOUT_XY), Time) + call send_perm_data(id_var3_id, cdata, PERM_TABLE(1:3, LAYOUT_XY), Time) + + ! Swap: axes (y,x) / (z,y,x) while canonical storage remains (x,y,...) -> pack to temp and send + call send_perm_data(id_var2_yx, cdata, PERM_TABLE(1:2, LAYOUT_YX), Time) + call send_perm_data(id_var3_zx, cdata, PERM_TABLE(1:3, LAYOUT_ZX), Time) + + ! Cyclic: axes (y,z,x) / (z,x,y) while canonical storage remains (x, y, ...) -> pack to temp and send + call send_perm_data(id_var3_yzx, cdata, PERM_TABLE(1:3, LAYOUT_YZX), Time) + call send_perm_data(id_var3_zxy, cdata, PERM_TABLE(1:3, LAYOUT_ZXY), Time) + + call diag_send_complete(Time_step) + call diag_send_complete(Time_step) + end do + + call diag_manager_end(Time) + call fms_end + +contains + + !> @brief Apply a predefined permutation to an axis array. + !> Maps canonical axis ordering to a permuted layout using PERM_TABLE. + !> Supports rank-2 and rank-3 axis subsets. + subroutine permute_axis(axis_in, perm_id, axis_out) + integer, intent(in) :: axis_in(:) + integer, intent(in) :: perm_id + integer, intent(out) :: axis_out(:) + + integer :: order(3) + + order = PERM_TABLE(:, perm_id) + + if (any(order(1:size(axis_out)) > size(axis_in))) then + call mpp_error(FATAL, "permute_axis: invalid permutation for given rank") + endif + + axis_out = axis_in(order(1:size(axis_out))) + end subroutine permute_axis + + !> @brief Register a diagnostic field with permuted axes. + !> Applies axis permutation before calling register_diag_field. + function register_perm_diag_field(var_name, long_name, axis, perm_id) result(id_var) + character(len=*), intent(in) :: var_name, long_name + integer, intent(in) :: axis(:) + integer, intent(in) :: perm_id + + integer :: id_var + integer, allocatable :: axis_perm(:) + + allocate(axis_perm(size(axis))) + call permute_axis(axis, perm_id, axis_perm) + id_var = register_diag_field('ocn_mod', var_name, axis_perm, Time, long_name, & + 'mullions', missing_value=missing_value) + deallocate(axis_perm) + end function register_perm_diag_field + + !> @brief Send data with optional axis permutation. + !> Applies 2D or 3D permutation to canonical (x,y,z,w) buffers before send_data. + !> Skips permutation for identity mappings. + subroutine send_perm_data(id_field, buf, order, Time_in) + integer, intent(in) :: id_field + real(r8_kind), intent(in) :: buf(:,:,:,:) ! canonical (x,y,z,w) + integer, intent(in) :: order(:) + type(time_type), intent(in) :: Time_in + + logical :: used_local + real(r8_kind), allocatable :: tmp2(:,:), tmp3(:,:,:) + + if (size(order) == 2) then + if (all(order == [1,2])) then + used_local = send_data(id_field, buf(:,:,1,1), Time_in) + else + tmp2 = permute(buf(:,:,1,1), order) + used_local = send_data(id_field, tmp2, Time_in) + endif + + else if (size(order) == 3) then + if (all(order == [1,2,3])) then + used_local = send_data(id_field, buf(:,:,:,1), Time_in) + else + tmp3 = permute(buf(:,:,:,1), order) + used_local = send_data(id_field, tmp3, Time_in) + endif + + else + call mpp_error(FATAL, "send_var_perm: unsupported permutation rank") + + endif + end subroutine send_perm_data + + !> @brief initialized the buffer based on the starting/ending indices + subroutine init_buffer(buffer, is, ie, js, je, nhalo) + real(r8_kind), intent(inout) :: buffer(:,:,:,:) + integer, intent(in) :: is, ie, js, je + integer, intent(in) :: nhalo + + integer :: ii, j, k, l + + do ii = is, ie + do j = js, je + do k = 1, size(buffer, 3) + do l = 1, size(buffer, 4) + buffer(ii-is+1+nhalo, j-js+1+nhalo, k, l) = real(ii, kind=r8_kind)*1000._r8_kind + & + real(j, kind=r8_kind)* 10._r8_kind + & + real(k, kind=r8_kind) + end do + end do + end do + end do + end subroutine init_buffer + + !> @brief Set the buffer based on the time_index + subroutine set_buffer(buffer, time_index) + real(r8_kind), intent(inout) :: buffer(:,:,:,:) + integer, intent(in) :: time_index + + buffer = nint(buffer) + real(time_index, kind=r8_kind)/100._r8_kind + end subroutine set_buffer + +end program test_diag_generalized_indices + diff --git a/test_fms/diag_manager/test_time_none.sh b/test_fms/diag_manager/test_time_none.sh index baa1c9954c..7b75b386da 100755 --- a/test_fms/diag_manager/test_time_none.sh +++ b/test_fms/diag_manager/test_time_none.sh @@ -27,6 +27,56 @@ if [ -z "${parser_skip}" ]; then # create and enter directory for in/output files output_dir +cat <<_EOF > diag_table.yaml +title: test_generalized_indices +base_date: 2 1 1 0 0 0 +diag_files: +- file_name: test_gen + freq: 6 hours + time_units: hours + unlimdim: time + varlist: + - module: ocn_mod + var_name: var2_id + output_name: var2_id + reduction: none + kind: r8 + - module: ocn_mod + var_name: var2_yx + output_name: var2_yx + reduction: none + kind: r8 + - module: ocn_mod + var_name: var3_id + output_name: var3_id + reduction: none + kind: r8 + - module: ocn_mod + var_name: var3_zx + output_name: var3_zx + reduction: none + kind: r8 + - module: ocn_mod + var_name: var3_yzx + output_name: var3_yzx + reduction: none + kind: r8 + - module: ocn_mod + var_name: var3_zxy + output_name: var3_zxy + reduction: none + kind: r8 +_EOF + +touch input.nml +printf "&diag_manager_nml \n use_modern_diag=.true. \n / \n&test_reduction_methods_nml \n test_case = 0 \n/" | cat > input.nml +test_expect_success "Write and read domain-decomposed data with generalized indices" ' + mpirun -n 6 ../test_diag_generalized_indices +' +test_expect_success "Checking answers for the generalized indices" ' + mpirun -n 1 ../check_diag_generalized_indices +' + cat <<_EOF > diag_table.yaml title: test_none base_date: 2 1 1 0 0 0 @@ -175,6 +225,8 @@ test_expect_success "Checking answers for the "none" reduction method with halo mpirun -n 1 ../check_time_none ' + + cat <<_EOF > diag_table.yaml title: test_none base_date: 2 1 1 0 0 0 diff --git a/test_fms/diag_manager/testing_utils.F90 b/test_fms/diag_manager/testing_utils.F90 index ba4616fa1f..579c2f4596 100644 --- a/test_fms/diag_manager/testing_utils.F90 +++ b/test_fms/diag_manager/testing_utils.F90 @@ -19,9 +19,10 @@ !> @brief Utilities used in multiple test module testing_utils use platform_mod, only: r8_kind + private - public :: allocate_buffer + public :: allocate_buffer, permute, check_perm public :: test_normal, test_openmp, test_halos public :: no_mask, logical_mask, real_mask @@ -32,6 +33,16 @@ module testing_utils integer, parameter :: logical_mask = 1 !< Using a logical mask integer, parameter :: real_mask = 2 !< Using a real mask + interface permute + module procedure permute_2d + module procedure permute_3d + end interface permute + + interface check_perm + module procedure check_perm_2d + module procedure check_perm_3d + end interface check_perm + contains !> @brief Allocate the output buffer based on the starting/ending indices @@ -49,4 +60,119 @@ function allocate_buffer(is, ie, js, je, k, l) & allocate(buffer(is:ie, js:je, 1:k, 1:l)) buffer = -999_r8_kind end function allocate_buffer + + !> @brief Apply a 2D axis permutation to an array. + !! @return permuted array with axes reordered according to order + function permute_2d(src, order) result(dst) + real(r8_kind), intent(in) :: src(:,:) + integer, intent(in) :: order(2) + real(r8_kind), allocatable :: dst(:,:) + + integer :: i, j + integer :: idx(2) + + allocate(dst(size(src,order(1)), size(src,order(2)))) + + do j = 1, size(src,2) + do i = 1, size(src,1) + idx = [i,j] + dst(idx(order(1)), idx(order(2))) = src(i,j) + end do + end do + end function permute_2d + + !> @brief Apply a 3D axis permutation to an array. + !! @return permuted array with axes reordered according to order + function permute_3d(src, order) result(dst) + real(r8_kind), intent(in) :: src(:,:,:) + integer, intent(in) :: order(3) + real(r8_kind), allocatable :: dst(:,:,:) + + integer :: i, j, k + integer :: idx(3) + + allocate(dst( size(src,order(1)), size(src,order(2)), size(src,order(3)) )) + + do k = 1, size(src,3) + do j = 1, size(src,2) + do i = 1, size(src,1) + idx = [i,j,k] + dst(idx(order(1)), idx(order(2)), idx(order(3))) = src(i,j,k) + enddo + enddo + enddo + end function permute_3d + + !> @brief Verify correctness of a 2D axis permutation. + !! Aborts if shape or values do not match expected permutation. + subroutine check_perm_2d(var, var_perm, order, ierr) + real(kind=r8_kind), intent(in) :: var(:,:) ! canonical (x,y) + real(kind=r8_kind), intent(in) :: var_perm(:,:) ! permuted + integer, intent(in) :: order(2) + integer, intent(out) :: ierr + + integer :: i, j + integer :: idx(2) + + ierr = 0 + + ! Check shape consistency + if ( size(var,order(1)) /= size(var_perm,1) .or. & + size(var,order(2)) /= size(var_perm,2) ) then + print *, "check_perm_2d: dimension mismatch" + ierr = 1 + endif + + do j = 1, size(var,2) + do i = 1, size(var,1) + idx = [i,j] + + if (abs(var(i,j) - var_perm(idx(order(1)), idx(order(2)))) > 0) then + print *, "perm mismatch at (x,y)=", i, j, "order=", order, & + " var =", var(i,j), & + " perm =", var_perm(idx(order(1)), idx(order(2))) + ierr = 1 + endif + + end do + end do + end subroutine check_perm_2d + + !> @brief Verify correctness of a 3D axis permutation. + !! Aborts if shape or values do not match expected permutation. + subroutine check_perm_3d(var, var_perm, order, ierr) + real(kind=r8_kind), intent(in) :: var(:,:,:) ! canonical (x,y,z) + real(kind=r8_kind), intent(in) :: var_perm(:,:,:) ! permuted + integer, intent(in) :: order(3) + integer, intent(out) :: ierr + + integer :: i, j, k + integer :: idx(3) + + ierr = 0 + + ! Check shape consistency + if ( size(var,order(1)) /= size(var_perm,1) .or. & + size(var,order(2)) /= size(var_perm,2) .or. & + size(var,order(3)) /= size(var_perm,3) ) then + print *, "check_perm_3d: dimension mismatch" + ierr = 1 + endif + + do k = 1, size(var,3) + do j = 1, size(var,2) + do i = 1, size(var,1) + idx = [i,j,k] + + if (abs(var(i,j,k) - var_perm( idx(order(1)), idx(order(2)), idx(order(3)) )) > 0) then + + print *, "perm mismatch at (x,y,z)=", i, j, k, "order=", order, & + " var =", var(i,j,k), & + " perm =", var_perm(idx(order(1)), idx(order(2)), idx(order(3))) + ierr = 1 + endif + enddo + enddo + enddo + end subroutine check_perm_3d end module