diff --git a/.github/actions/macos-setup/action.yml b/.github/actions/macos-setup/action.yml new file mode 100644 index 0000000000..7f38a81cc6 --- /dev/null +++ b/.github/actions/macos-setup/action.yml @@ -0,0 +1,23 @@ +# This file is part of MOM6, the Modular Ocean Model version 6. +# See the LICENSE file for licensing information. +# SPDX-License-Identifier: Apache-2. + +name: 'install-macos-prerequisites' + +description: 'Install prerequisites for Mac OS compilation' + +runs: + using: 'composite' + + steps: + - name: Install macOS packages + shell: bash + run: | + echo "::group::Install packages" + brew reinstall gcc + brew install automake + brew install netcdf + brew install netcdf-fortran + brew install openmpi + brew install libyaml + echo "::endgroup::" diff --git a/.github/workflows/github_cmake_macos.yml b/.github/workflows/github_cmake_macos.yml new file mode 100644 index 0000000000..721457ba94 --- /dev/null +++ b/.github/workflows/github_cmake_macos.yml @@ -0,0 +1,33 @@ +name: CMake build and unit tests on MacOS + +on: + push: + branches: + - main + +jobs: + build: + runs-on: macOS-latest + strategy: + matrix: + libyaml-flag: [ "", -DWITH_YAML=on ] + build-type: [ "-DCMAKE_BUILD_TYPE=Release", "-DCMAKE_BUILD_TYPE=Debug" ] + env: + CMAKE_FLAGS: "${{ matrix.build-type }} ${{ matrix.libyaml-flag }}" + EXCLUDE_TESTS: "test_mpp_nesting|test_mpp_clock_begin_end_id|test_time_*" + CC: gcc-15 + FC: gfortran + steps: + - uses: actions/checkout@v4 + - uses: ./.github/actions/macos-setup/ + - name: Generate makefiles with CMake + run: | + mkdir build + cd build + brew shellenv + cmake $CMAKE_FLAGS .. + - name: Build the library + run: | + brew shellenv + make -j -C build + diff --git a/CMakeLists.txt b/CMakeLists.txt index 05e941b181..7f9705443b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -314,6 +314,10 @@ foreach(kind ${kinds}) target_link_libraries(${libTgt}_c PRIVATE OpenMP::OpenMP_C) endif() + if(YAML_FOUND) + target_link_libraries(${libTgt}_c PRIVATE PkgConfig::YAML) + endif() + # Fortran add_library(${libTgt}_f OBJECT ${fms_fortran_src_files}) @@ -459,6 +463,10 @@ if(NOT kinds) target_link_libraries(${libTgt}_c PRIVATE OpenMP::OpenMP_C) endif() + if(YAML_FOUND) + target_link_libraries(${libTgt}_c PRIVATE PkgConfig::YAML) + endif() + # Fortran add_library(${libTgt}_f OBJECT ${fms_fortran_src_files}) @@ -717,6 +725,7 @@ if(UNIT_TESTS) test_fms/fms/test_fms.F90 test_fms/interpolator/test_interpolator.F90 test_fms/mpp/test_clock_init.F90 + test_fms/mpp/test_corner_mosaic.F90 test_fms/mpp/test_domains_simple.F90 test_fms/mpp/test_domains_utility_mod.F90 test_fms/mpp/test_global_arrays.F90 diff --git a/diag_manager/diag_output.F90 b/diag_manager/diag_output.F90 index 0336f86048..cf8cd4f85a 100644 --- a/diag_manager/diag_output.F90 +++ b/diag_manager/diag_output.F90 @@ -48,7 +48,6 @@ MODULE diag_output_mod use mpp_domains_mod, only: mpp_get_UG_io_domain use mpp_domains_mod, only: mpp_get_UG_domain_npes use mpp_domains_mod, only: mpp_get_UG_domain_pelist - use mpp_mod, only: mpp_gather use mpp_mod, only: uppercase,lowercase use fms2_io_mod diff --git a/diag_manager/fms_diag_axis_object.F90 b/diag_manager/fms_diag_axis_object.F90 index e13ae7551a..73904d0d40 100644 --- a/diag_manager/fms_diag_axis_object.F90 +++ b/diag_manager/fms_diag_axis_object.F90 @@ -458,7 +458,13 @@ subroutine write_axis_data(this, fms2io_fileobj, parent_axis) if (present(parent_axis)) then select type(parent_axis) type is (fmsDiagFullAxis_type) - call write_data(fms2io_fileobj, this%subaxis_name, parent_axis%axis_data(i:j)) + ! Added this select type so the the data is indexed correctly + select type (vardata => parent_axis%axis_data) + type is (real(kind=r8_kind)) + call write_data(fms2io_fileobj, this%subaxis_name, vardata(i:j)) + type is (real(kind=r4_kind)) + call write_data(fms2io_fileobj, this%subaxis_name, vardata(i:j)) + end select end select endif type is (fmsDiagDiurnalAxis_type) @@ -881,8 +887,10 @@ subroutine fill_subaxis(this, starting_index, ending_index, axis_id, parent_id, if (present(nz_subaxis)) nsubaxis = nz_subaxis this%axis_id = axis_id - this%starting_index = starting_index - this%ending_index = ending_index + + ! The min and max were added here to support axis that are both increasing and decreasing + this%starting_index = min(starting_index, ending_index) + this%ending_index = max(starting_index, ending_index) this%parent_axis_id = parent_id write(nsubaxis_char, '(i2.2)') nsubaxis this%subaxis_name = trim(parent_axis_name)//"_sub"//nsubaxis_char diff --git a/fms2_io/include/compressed_write.inc b/fms2_io/include/compressed_write.inc index 4f1ab51a7d..01a6db8ab1 100644 --- a/fms2_io/include/compressed_write.inc +++ b/fms2_io/include/compressed_write.inc @@ -137,7 +137,7 @@ subroutine compressed_write_1d(fileobj, variable_name, cdata, unlim_dim_level, & call error("unsupported variable type: "//trim(append_error_msg)) end select else - select type(cdata) + select type(cdata) type is (integer(kind=i4_kind)) allocate(buf_i4_kind(1)) type is (integer(kind=i8_kind)) @@ -148,7 +148,7 @@ subroutine compressed_write_1d(fileobj, variable_name, cdata, unlim_dim_level, & allocate(buf_r8_kind(1)) class default call error("unsupported variable type: "//trim(append_error_msg)) - end select + end select endif !Gather the data onto the I/O root and write it out. @@ -253,6 +253,19 @@ subroutine compressed_write_2d(fileobj, variable_name, cdata, unlim_dim_level, & class default call error("unsupported variable type: "//trim(append_error_msg)) end select + else + select type(cdata) + type is (integer(kind=i4_kind)) + allocate(buf_i4_kind(1, 1)) + type is (integer(kind=i8_kind)) + allocate(buf_i8_kind(1, 1)) + type is (real(kind=r4_kind)) + allocate(buf_r4_kind(1, 1)) + type is (real(kind=r8_kind)) + allocate(buf_r8_kind(1, 1)) + class default + call error("unsupported variable type: "//trim(append_error_msg)) + end select endif c(:) = 1 @@ -376,7 +389,20 @@ subroutine compressed_write_3d(fileobj, variable_name, cdata, unlim_dim_level, & call allocate_array(buf_r8_kind, e) class default call error("unsupported variable type: "//trim(append_error_msg)) - end select + end select + else + select type(cdata) + type is (integer(kind=i4_kind)) + allocate(buf_i4_kind(1, 1, 1)) + type is (integer(kind=i8_kind)) + allocate(buf_i8_kind(1, 1, 1)) + type is (real(kind=r4_kind)) + allocate(buf_r4_kind(1, 1, 1)) + type is (real(kind=r8_kind)) + allocate(buf_r8_kind(1, 1, 1)) + class default + call error("unsupported variable type: "//trim(append_error_msg)) + end select endif c(:) = 1 diff --git a/fms2_io/include/domain_read.inc b/fms2_io/include/domain_read.inc index fa404270f4..764d91a725 100644 --- a/fms2_io/include/domain_read.inc +++ b/fms2_io/include/domain_read.inc @@ -119,6 +119,7 @@ subroutine domain_read_2d(fileobj, variable_name, vdata, unlim_dim_level, & integer :: xgsize !< Size of global x io domain integer :: ygbegin !< Starting y index of global io domain integer :: ygsize !< Size of global y io domain + integer :: dim_order(2) !< Order of the dimensions type(domain2d), pointer :: io_domain !< pointer to the io_domain !< The global data is only allocated by the io root PEs @@ -164,14 +165,6 @@ subroutine domain_read_2d(fileobj, variable_name, vdata, unlim_dim_level, & return endif - if (xdim_index .ne. 1 .or. ydim_index .ne. 2) then - ! This is a KLUDGE - ! mpp_scatter assumes that the variable is (x,y), if that is not the case it remaps the data - ! to a 4D array and calls domain_read_4d which does not use mpp_scatter yet - vdata_dummy(1:size(vdata,1),1:size(vdata,2), 1:1, 1:1) => vdata(:,:) - call domain_read_4d(fileobj, variable_name, vdata_dummy, unlim_dim_level) - return - endif io_domain => mpp_get_io_domain(fileobj%domain) c(:) = 1 e(:) = shape(vdata) @@ -240,29 +233,32 @@ subroutine domain_read_2d(fileobj, variable_name, vdata, unlim_dim_level, & e(xdim_index) = xc_size e(ydim_index) = yc_size + dim_order(xdim_index) = 1 + dim_order(ydim_index) = 2 + select type(vdata) type is (integer(kind=i4_kind)) call allocate_array(buf_i4_kind_pe, e) call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, fileobj%pelist, & - buf_i4_kind_pe, buf_i4_kind, fileobj%is_root) + buf_i4_kind_pe, buf_i4_kind, dim_order, fileobj%is_root) call put_array_section(buf_i4_kind_pe, vdata, c, e) deallocate(buf_i4_kind_pe) type is (integer(kind=i8_kind)) call allocate_array(buf_i8_kind_pe, e) call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, fileobj%pelist, & - buf_i8_kind_pe, buf_i8_kind, fileobj%is_root) + buf_i8_kind_pe, buf_i8_kind, dim_order, fileobj%is_root) call put_array_section(buf_i8_kind_pe, vdata, c, e) deallocate(buf_i8_kind_pe) type is (real(kind=r4_kind)) call allocate_array(buf_r4_kind_pe, e) call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, fileobj%pelist, & - buf_r4_kind_pe, buf_r4_kind, fileobj%is_root) + buf_r4_kind_pe, buf_r4_kind, dim_order, fileobj%is_root) call put_array_section(buf_r4_kind_pe, vdata, c, e) deallocate(buf_r4_kind_pe) type is (real(kind=r8_kind)) call allocate_array(buf_r8_kind_pe, e) call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, fileobj%pelist, & - buf_r8_kind_pe, buf_r8_kind, fileobj%is_root) + buf_r8_kind_pe, buf_r8_kind, dim_order, fileobj%is_root) call put_array_section(buf_r8_kind_pe, vdata, c, e) deallocate(buf_r8_kind_pe) class default @@ -304,6 +300,7 @@ subroutine domain_read_3d(fileobj, variable_name, vdata, unlim_dim_level, & integer :: xdim_index !< The index of the variable that is the x dimension integer :: ydim_index !< The index of the variable that is the y dimension + integer :: zdim_index !< The index of the variable that is the z dimension integer :: xpos !< The position of the x axis integer :: ypos !< The position of the y axis integer :: i !< For do loops @@ -320,6 +317,7 @@ subroutine domain_read_3d(fileobj, variable_name, vdata, unlim_dim_level, & integer :: xgsize !< Size of global x io domain integer :: ygbegin !< Starting y index of global io domain integer :: ygsize !< Size of global y io domain + integer :: dim_order(3) !< Order of the dimensions type(domain2d), pointer :: io_domain !< pointer to the io_domain !< The global data is only allocated by the io root PEs @@ -365,14 +363,6 @@ subroutine domain_read_3d(fileobj, variable_name, vdata, unlim_dim_level, & return endif - if (xdim_index .ne. 1 .or. ydim_index .ne. 2) then - ! This is a KLUDGE - ! mpp_scatter assumes that the variable is (x,y), if that is not the case it remaps the data - ! to a 4D array and calls domain_read_4d which does not use mpp_scatter yet - vdata_dummy(1:size(vdata,1),1:size(vdata,2), 1:size(vdata,3), 1:1) => vdata(:,:,:) - call domain_read_4d(fileobj, variable_name, vdata_dummy, unlim_dim_level) - return - endif io_domain => mpp_get_io_domain(fileobj%domain) c(:) = 1 if (present(corner)) c = corner @@ -444,29 +434,36 @@ subroutine domain_read_3d(fileobj, variable_name, vdata, unlim_dim_level, & e(xdim_index) = xc_size e(ydim_index) = yc_size + ! Calculate the index of the z dimension + zdim_index = 6 - xdim_index - ydim_index + + dim_order(xdim_index) = 1 + dim_order(ydim_index) = 2 + dim_order(zdim_index) = 3 + select type(vdata) type is (integer(kind=i4_kind)) call allocate_array(buf_i4_kind_pe, e) - call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(3), fileobj%pelist, & - buf_i4_kind_pe, buf_i4_kind, fileobj%is_root) + call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(zdim_index), & + fileobj%pelist, buf_i4_kind_pe, buf_i4_kind, dim_order, fileobj%is_root) call put_array_section(buf_i4_kind_pe, vdata, c, e) deallocate(buf_i4_kind_pe) type is (integer(kind=i8_kind)) call allocate_array(buf_i8_kind_pe, e) - call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(3), fileobj%pelist, & - buf_i8_kind_pe, buf_i8_kind, fileobj%is_root) + call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(zdim_index), & + fileobj%pelist, buf_i8_kind_pe, buf_i8_kind, dim_order, fileobj%is_root) call put_array_section(buf_i8_kind_pe, vdata, c, e) deallocate(buf_i8_kind_pe) type is (real(kind=r4_kind)) call allocate_array(buf_r4_kind_pe, e) - call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(3), fileobj%pelist, & - buf_r4_kind_pe, buf_r4_kind, fileobj%is_root) + call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(zdim_index), & + fileobj%pelist, buf_r4_kind_pe, buf_r4_kind, dim_order, fileobj%is_root) call put_array_section(buf_r4_kind_pe, vdata, c, e) deallocate(buf_r4_kind_pe) type is (real(kind=r8_kind)) call allocate_array(buf_r8_kind_pe, e) - call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(3), fileobj%pelist, & - buf_r8_kind_pe, buf_r8_kind, fileobj%is_root) + call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(zdim_index), & + fileobj%pelist, buf_r8_kind_pe, buf_r8_kind, dim_order, fileobj%is_root) call put_array_section(buf_r8_kind_pe, vdata, c, e) deallocate(buf_r8_kind_pe) class default diff --git a/fms2_io/include/domain_write.inc b/fms2_io/include/domain_write.inc index 3d9f52c1bf..3e8e19570f 100644 --- a/fms2_io/include/domain_write.inc +++ b/fms2_io/include/domain_write.inc @@ -121,8 +121,9 @@ subroutine domain_write_2d(fileobj, variable_name, vdata, unlim_dim_level, & class(*), dimension(:,:,:,:), pointer :: vdata_dummy !< Vdata remapped as 4D integer :: xgmin !< Starting x index of the global io domain integer :: ygmin !< Ending y index of the global io domain - integer :: xgsize, ygsize - integer :: istart, jstart, iend, jend + integer :: gsize(2) !< Shape of global_buf + integer :: dim_order(2) !< Order of the dimensions + integer :: start(2), end(2) integer :: ioff, joff if (fileobj%use_netcdf_mpi) then @@ -157,8 +158,8 @@ subroutine domain_write_2d(fileobj, variable_name, vdata, unlim_dim_level, & msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name)) ! Get the global io domain: - call mpp_get_global_domain(io_domain, xbegin=xgmin, xsize=xgsize, position=xpos) - call mpp_get_global_domain(io_domain, ybegin=ygmin, ysize=ygsize, position=ypos) + call mpp_get_global_domain(io_domain, xbegin=xgmin, xsize=gsize(xdim_index), position=xpos) + call mpp_get_global_domain(io_domain, ybegin=ygmin, ysize=gsize(ydim_index), position=ypos) ! Root pe allocates room to gather the data and computes offsets for data placement if (fileobj%is_root) then @@ -169,25 +170,25 @@ subroutine domain_write_2d(fileobj, variable_name, vdata, unlim_dim_level, & ! Allocate recv buffer for gather select type(vdata) type is (integer(kind=i4_kind)) - allocate(global_buf_i4_kind(xgsize, ygsize)) + allocate(global_buf_i4_kind(gsize(1), gsize(2))) global_buf_i4_kind = 0 if (get_fill_value(fileobj, variable_name, fill_i4_kind, broadcast=.false.)) then global_buf_i4_kind = fill_i4_kind endif type is (integer(kind=i8_kind)) - allocate(global_buf_i8_kind(xgsize, ygsize)) + allocate(global_buf_i8_kind(gsize(1), gsize(2))) global_buf_i8_kind = 0 if (get_fill_value(fileobj, variable_name, fill_i8_kind, broadcast=.false.)) then global_buf_i8_kind = fill_i8_kind endif type is (real(kind=r4_kind)) - allocate(global_buf_r4_kind(xgsize, ygsize)) + allocate(global_buf_r4_kind(gsize(1), gsize(2))) global_buf_r4_kind = 0. if (get_fill_value(fileobj, variable_name, fill_r4_kind, broadcast=.false.)) then global_buf_r4_kind = fill_r4_kind endif type is (real(kind=r8_kind)) - allocate(global_buf_r8_kind(xgsize, ygsize)) + allocate(global_buf_r8_kind(gsize(1), gsize(2))) global_buf_r8_kind = 0. if (get_fill_value(fileobj, variable_name, fill_r8_kind, broadcast=.false.)) then global_buf_r8_kind = fill_r8_kind @@ -196,35 +197,55 @@ subroutine domain_write_2d(fileobj, variable_name, vdata, unlim_dim_level, & call error("unsupported variable type: domain_write_2d_mpp_gather: file: " & & //trim(fileobj%path)//" variable:"//trim(variable_name)) end select + else + select type(vdata) + type is (integer(kind=i4_kind)) + allocate(global_buf_i4_kind(1, 1)) + type is (integer(kind=i8_kind)) + allocate(global_buf_i8_kind(1, 1)) + type is (real(kind=r4_kind)) + allocate(global_buf_r4_kind(1, 1)) + type is (real(kind=r8_kind)) + allocate(global_buf_r8_kind(1, 1)) + class default + call error("unsupported variable type: domain_write_2d_mpp_gather: file: " & + & //trim(fileobj%path)//" variable:"//trim(variable_name)) + end select endif ! Get the starting and indices of the compute domain relative to vdata (note that vdata start indices at 1 #Fortran) - istart = 1 - jstart = 1 + start = 1 ! If the buffer contains halos, get the portion of vdata with only the compute domain if (buffer_includes_halos) then - istart = isc - isd + 1 - jstart = jsc - jsd + 1 + start(xdim_index) = isc - isd + 1 + start(ydim_index) = jsc - jsd + 1 endif - iend = istart + xc_size - 1 - jend = jstart + yc_size - 1 + end(xdim_index) = start(xdim_index) + xc_size - 1 + end(ydim_index) = start(ydim_index) + yc_size - 1 + + dim_order(xdim_index) = 1 + dim_order(ydim_index) = 2 ! Gather the data select type(vdata) type is (integer(kind=i4_kind)) - call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, fileobj%pelist, vdata(istart:iend, jstart:jend), & - & global_buf_i4_kind, fileobj%is_root, ishift=ioff, jshift=joff) + call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, fileobj%pelist, & + vdata(start(1):end(1), start(2):end(2)), global_buf_i4_kind, & + dim_order, fileobj%is_root, ishift=ioff, jshift=joff) type is (integer(kind=i8_kind)) - call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, fileobj%pelist, vdata(istart:iend, jstart:jend), & - & global_buf_i8_kind, fileobj%is_root, ishift=ioff, jshift=joff) + call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, fileobj%pelist, & + vdata(start(1):end(1), start(2):end(2)), global_buf_i8_kind, & + dim_order, fileobj%is_root, ishift=ioff, jshift=joff) type is (real(kind=r4_kind)) - call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, fileobj%pelist, vdata(istart:iend, jstart:jend), & - & global_buf_r4_kind, fileobj%is_root, ishift=ioff, jshift=joff) + call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, fileobj%pelist, & + vdata(start(1):end(1), start(2):end(2)), global_buf_r4_kind, & + dim_order, fileobj%is_root, ishift=ioff, jshift=joff) type is (real(kind=r8_kind)) - call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, fileobj%pelist, vdata(istart:iend, jstart:jend), & - & global_buf_r8_kind, fileobj%is_root, ishift=ioff, jshift=joff) + call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, fileobj%pelist, & + vdata(start(1):end(1), start(2):end(2)), global_buf_r8_kind, & + dim_order, fileobj%is_root, ishift=ioff, jshift=joff) class default call error("unsupported variable type: domain_write_2d_mpp_gather: file: " & & //trim(fileobj%path)//" variable:"// trim(variable_name)) @@ -293,6 +314,7 @@ subroutine domain_write_3d(fileobj, variable_name, vdata, unlim_dim_level, & integer :: ydim_index integer :: ypos integer :: yc_size + integer :: zdim_index integer(kind=i4_kind) :: fill_i4_kind !< Fill value of a i4_kind variable integer(kind=i8_kind) :: fill_i8_kind !< Fill value of a i8_kind variable real(kind=r4_kind) :: fill_r4_kind !< Fill value of a r4_kind variable @@ -300,8 +322,9 @@ subroutine domain_write_3d(fileobj, variable_name, vdata, unlim_dim_level, & class(*), dimension(:,:,:,:), pointer :: vdata_dummy !< Vdata remapped as 4D integer :: xgmin !< Starting x index of the global io domain integer :: ygmin !< Ending y index of the global io domain - integer :: xgsize, ygsize - integer :: istart, jstart, iend, jend + integer :: gsize(3) !< Shape of global_buf + integer :: dim_order(3) !< Order of the dimensions + integer :: start(3), end(3) integer :: ioff, joff if (fileobj%use_netcdf_mpi) then @@ -317,15 +340,6 @@ subroutine domain_write_3d(fileobj, variable_name, vdata, unlim_dim_level, & return endif - if (xdim_index .ne. 1 .or. ydim_index .ne. 2) then - ! This is a KLUDGE - ! mpp_scatter assumes that the variable is (x,y), if that is not the case it remaps the data - ! to a 4D array and calls domain_read_4d which does not use mpp_scatter yet - vdata_dummy(1:size(vdata,1),1:size(vdata,2), 1:size(vdata,3), 1:1) => vdata(:,:,:) - call domain_write_4d(fileobj, variable_name, vdata_dummy, unlim_dim_level) - return - endif - ! Get the io domain from the fileobj io_domain => mpp_get_io_domain(fileobj%domain) @@ -334,9 +348,13 @@ subroutine domain_write_3d(fileobj, variable_name, vdata, unlim_dim_level, & xpos, ypos, isd, isc, xc_size, jsd, jsc, yc_size, buffer_includes_halos, & msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name)) + ! Calculate the index of the z dimension + zdim_index = 6 - xdim_index - ydim_index + ! Get the global io domain: - call mpp_get_global_domain(io_domain, xbegin=xgmin, xsize=xgsize, position=xpos) - call mpp_get_global_domain(io_domain, ybegin=ygmin, ysize=ygsize, position=ypos) + call mpp_get_global_domain(io_domain, xbegin=xgmin, xsize=gsize(xdim_index), position=xpos) + call mpp_get_global_domain(io_domain, ybegin=ygmin, ysize=gsize(ydim_index), position=ypos) + gsize(zdim_index) = size(vdata, zdim_index) ! Root pe allocates room to gather the data and computes offsets for data placement if (fileobj%is_root) then @@ -347,25 +365,25 @@ subroutine domain_write_3d(fileobj, variable_name, vdata, unlim_dim_level, & ! Allocate recv buffer for gather select type(vdata) type is (integer(kind=i4_kind)) - allocate(global_buf_i4_kind(xgsize, ygsize, size(vdata,3))) + allocate(global_buf_i4_kind(gsize(1), gsize(2), gsize(3))) global_buf_i4_kind = 0 if (get_fill_value(fileobj, variable_name, fill_i4_kind, broadcast=.false.)) then global_buf_i4_kind = fill_i4_kind endif type is (integer(kind=i8_kind)) - allocate(global_buf_i8_kind(xgsize, ygsize, size(vdata,3))) + allocate(global_buf_i8_kind(gsize(1), gsize(2), gsize(3))) global_buf_i8_kind = 0 if (get_fill_value(fileobj, variable_name, fill_i8_kind, broadcast=.false.)) then global_buf_i8_kind = fill_i8_kind endif type is (real(kind=r4_kind)) - allocate(global_buf_r4_kind(xgsize, ygsize, size(vdata,3))) + allocate(global_buf_r4_kind(gsize(1), gsize(2), gsize(3))) global_buf_r4_kind = 0. if (get_fill_value(fileobj, variable_name, fill_r4_kind, broadcast=.false.)) then global_buf_r4_kind = fill_r4_kind endif type is (real(kind=r8_kind)) - allocate(global_buf_r8_kind(xgsize, ygsize, size(vdata,3))) + allocate(global_buf_r8_kind(gsize(1), gsize(2), gsize(3))) global_buf_r8_kind = 0. if (get_fill_value(fileobj, variable_name, fill_r8_kind, broadcast=.false.)) then global_buf_r8_kind = fill_r8_kind @@ -374,39 +392,61 @@ subroutine domain_write_3d(fileobj, variable_name, vdata, unlim_dim_level, & call error("unsupported variable type: domain_write_3d_mpp_gather: file: " & & //trim(fileobj%path)//" variable:"//trim(variable_name)) end select + else + select type(vdata) + type is (integer(kind=i4_kind)) + allocate(global_buf_i4_kind(1, 1, 1)) + type is (integer(kind=i8_kind)) + allocate(global_buf_i8_kind(1, 1, 1)) + type is (real(kind=r4_kind)) + allocate(global_buf_r4_kind(1, 1, 1)) + type is (real(kind=r8_kind)) + allocate(global_buf_r8_kind(1, 1, 1)) + class default + call error("unsupported variable type: domain_write_3d_mpp_gather: file: " & + & //trim(fileobj%path)//" variable:"//trim(variable_name)) + end select endif ! Get the starting and indices of the compute domain relative to vdata(note that vdata start indices at 1 #Fortran) - istart = 1 - jstart = 1 + start = 1 ! If the buffer contains halos, get the portion of vdata with only the compute domain if (buffer_includes_halos) then - istart = isc - isd + 1 - jstart = jsc - jsd + 1 + start(xdim_index) = isc - isd + 1 + start(ydim_index) = jsc - jsd + 1 endif - iend = istart + xc_size - 1 - jend = jstart + yc_size - 1 + end(xdim_index) = start(xdim_index) + xc_size - 1 + end(ydim_index) = start(ydim_index) + yc_size - 1 + end(zdim_index) = size(vdata, zdim_index) ! Get offsets for buffer ioff = 1-xgmin joff = 1-ygmin + dim_order(xdim_index) = 1 + dim_order(ydim_index) = 2 + dim_order(zdim_index) = 3 + ! Gather the data select type(vdata) type is (integer(kind=i4_kind)) - call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, size(vdata,3), fileobj%pelist, & - & vdata(istart:iend, jstart:jend,:), global_buf_i4_kind, fileobj%is_root, ishift=ioff, jshift=joff) + call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, size(vdata, zdim_index), fileobj%pelist, & + vdata(start(1):end(1), start(2):end(2), start(3):end(3)), global_buf_i4_kind, & + dim_order, fileobj%is_root, ishift=ioff, jshift=joff) type is (integer(kind=i8_kind)) - call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, size(vdata,3), fileobj%pelist, & - & vdata(istart:iend, jstart:jend,:), global_buf_i8_kind, fileobj%is_root, ishift=ioff, jshift=joff) + call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, size(vdata, zdim_index), fileobj%pelist, & + vdata(start(1):end(1), start(2):end(2), start(3):end(3)), global_buf_i8_kind, & + dim_order, fileobj%is_root, ishift=ioff, jshift=joff) type is (real(kind=r4_kind)) - call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, size(vdata,3), fileobj%pelist, & - & vdata(istart:iend, jstart:jend,:), global_buf_r4_kind, fileobj%is_root, ishift=ioff, jshift=joff) + call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, size(vdata, zdim_index), fileobj%pelist, & + vdata(start(1):end(1), start(2):end(2), start(3):end(3)), global_buf_r4_kind, & + dim_order, fileobj%is_root, ishift=ioff, jshift=joff) type is (real(kind=r8_kind)) - call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, size(vdata,3), fileobj%pelist, & - & vdata(istart:iend, jstart:jend,:), global_buf_r8_kind, fileobj%is_root, ishift=ioff, jshift=joff) + call mpp_gather(isc, isc+xc_size-1, jsc, jsc+yc_size-1, size(vdata, zdim_index), fileobj%pelist, & + vdata(start(1):end(1), start(2):end(2), start(3):end(3)), global_buf_r8_kind, & + dim_order, fileobj%is_root, ishift=ioff, jshift=joff) class default call error("unsupported variable type: domain_write_3d_mpp_gather: file: " & & //trim(fileobj%path)//" variable:"//trim(variable_name)) diff --git a/fms2_io/include/gather_data_bc.inc b/fms2_io/include/gather_data_bc.inc index 3734ca401c..c1deb03f0e 100644 --- a/fms2_io/include/gather_data_bc.inc +++ b/fms2_io/include/gather_data_bc.inc @@ -81,7 +81,11 @@ subroutine gather_data_bc_2d(fileobj, vdata, bc_info) i1=1; i2=1; j1=1; j2=1 else !! In this case there is data in fileobj's root, so there is no need for the dummy data - if(fileobj%is_root) allocate(global_buf_r4_kind(i_glob, j_glob)) + if(fileobj%is_root) then + allocate(global_buf_r4_kind(i_glob, j_glob)) + else + allocate(global_buf_r4_kind(1, 1)) + endif allocate(local_buf_r4_kind(size(vdata,1), size(vdata,2))) local_buf_r4_kind = vdata endif @@ -112,7 +116,11 @@ subroutine gather_data_bc_2d(fileobj, vdata, bc_info) i1=1; i2=1; j1=1; j2=1 else !! In this case there is data in fileobj's root, so there is no need for the dummy data - if(fileobj%is_root) allocate(global_buf_r8_kind(i_glob, j_glob)) + if(fileobj%is_root) then + allocate(global_buf_r8_kind(i_glob, j_glob)) + else + allocate(global_buf_r8_kind(1, 1)) + endif allocate(local_buf_r8_kind(size(vdata,1), size(vdata,2))) local_buf_r8_kind = vdata endif @@ -203,7 +211,11 @@ subroutine gather_data_bc_3d(fileobj, vdata, bc_info) i1=1; i2=1; j1=1; j2=1 else !! In this case there is data in fileobj's root, so there is no need for the dummy data - if(fileobj%is_root) allocate(global_buf_r4_kind(i_glob, j_glob, k_glob)) + if(fileobj%is_root) then + allocate(global_buf_r4_kind(i_glob, j_glob, k_glob)) + else + allocate(global_buf_r4_kind(1, 1, 1)) + endif allocate(local_buf_r4_kind(size(vdata,1), size(vdata,2), size(vdata,3))) local_buf_r4_kind = vdata endif @@ -233,7 +245,11 @@ subroutine gather_data_bc_3d(fileobj, vdata, bc_info) i1=1; i2=1; j1=1; j2=1 else !! In this case there is data in fileobj's root, so there is no need for the dummy data - if(fileobj%is_root) allocate(global_buf_r8_kind(i_glob, j_glob, k_glob)) + if(fileobj%is_root) then + allocate(global_buf_r8_kind(i_glob, j_glob, k_glob)) + else + allocate(global_buf_r8_kind(1, 1, 1)) + endif allocate(local_buf_r8_kind(size(vdata,1), size(vdata,2), size(vdata,3))) local_buf_r8_kind = vdata endif diff --git a/grid_utils/constant.h b/grid_utils/constant.h index dfadf9bd1c..6983b4784a 100644 --- a/grid_utils/constant.h +++ b/grid_utils/constant.h @@ -15,7 +15,7 @@ * PARTICULAR PURPOSE. See the License for the specific language * governing permissions and limitations under the License. ***********************************************************************/ -#define RADIUS (6371000.) + #define STRING 255 #define EPSLN8 (1.e-8) @@ -25,3 +25,11 @@ #define R2D (180/M_PI) #define TPI (2.0*M_PI) #define HPI (0.5*M_PI) + +/* This matches FMSconstants behavior, and sets Earth radius with gfs_constants.fh value, + or with gfdl_constants.fh == geos_constants.fh value */ +#ifdef GFS_CONSTANTS +#define RADIUS (6371200.) +#else +#define RADIUS (6371000.) +#endif \ No newline at end of file diff --git a/interpolator/include/interpolator.inc b/interpolator/include/interpolator.inc index fca7b66092..b2b31e44fe 100644 --- a/interpolator/include/interpolator.inc +++ b/interpolator/include/interpolator.inc @@ -1379,6 +1379,7 @@ integer, parameter :: lkind=FMS_INTP_KIND_ !! @param [in] OPTIONAL: Index for the physics window !! @param [out] The model fields with the interpolated climatology data !! @param [out] OPTIONAL: The units of field_name +!! @param [in] OPTIONAL: Order of the x (1), y (2), z (3), and field (4) dimensions !! !! @throw FATAL "interpolator_4D : You must call interpolator_init !! before calling interpolator" @@ -1397,7 +1398,7 @@ integer, parameter :: lkind=FMS_INTP_KIND_ !! @throw FATAL "Interpolator: the field name is not contained in this !! intepolate_type: " subroutine INTERPOLATOR_4D_(clim_type, Time, phalf, interp_data, & - field_name, is,js, clim_units) + field_name, is, js, clim_units, dim_order) ! ! Return 4-D field interpolated to model grid and time ! @@ -1412,6 +1413,7 @@ subroutine INTERPOLATOR_4D_(clim_type, Time, phalf, interp_data, & ! Time : The model time that you wish to interpolate to. ! phalf : The half level model pressure field. ! is, js : The indices of the physics window. +! dim_order : Order of the x (1), y (2), z (3), and field (4) dimensions ! ! INTENT OUT ! interp_data : The model fields with the interpolated climatology data. @@ -1421,15 +1423,13 @@ type(interpolate_type), intent(inout) :: clim_type character(len=*) , intent(in) :: field_name type(time_type) , intent(in) :: Time real(FMS_INTP_KIND_), dimension(:,:,:), intent(in) :: phalf -real(FMS_INTP_KIND_), dimension(:,:,:,:), intent(out) :: interp_data +real(FMS_INTP_KIND_), dimension(:,:,:,:), intent(out), target :: interp_data integer , intent(in) , optional :: is,js character(len=*) , intent(out), optional :: clim_units +integer , intent(in), optional :: dim_order(4) integer :: taum, taup, ilon -real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),size(interp_data,2),& - size(clim_type%FMS_INTP_TYPE_%levs(:)),size(clim_type%field_name(:))) -real(FMS_INTP_KIND_) :: p_fact(size(interp_data,1),size(interp_data,2)) -real(FMS_INTP_KIND_) :: col_data(size(interp_data,1),size(interp_data,2), & - size(clim_type%field_name(:))) +real(FMS_INTP_KIND_), allocatable :: hinterp_data(:,:,:,:), p_fact(:,:), col_data(:,:,:), phalf_diff(:,:,:) +real(FMS_INTP_KIND_), pointer :: interp_buf(:,:,:,:) real(FMS_INTP_KIND_) :: pclim(size(clim_type%FMS_INTP_TYPE_%halflevs(:))) integer :: istart,iend,jstart,jend logical :: result, found @@ -1442,14 +1442,45 @@ type(time_type) :: t_prev, t_next type(time_type), dimension(2) :: month integer :: indexm, indexp, yearm, yearp integer :: i, j, k, n +integer :: ni, nj, nk +integer :: nlevs, nfields +integer :: dim_map(4) !< x, y, z, field +logical :: use_buf character(len=256) :: err_msg +integer, parameter :: dim_map_default(4) = [1, 2, 3, 4] integer, parameter :: lkind=FMS_INTP_KIND_ if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) & call mpp_error(FATAL, "interpolator_4D : You must call interpolator_init before calling interpolator") - do n=2,size(clim_type%field_name(:)) +if (present(dim_order)) then + call inverse_permutation(dim_order, dim_map) +else + dim_map = dim_map_default +endif + +use_buf = any(dim_map.ne.dim_map_default) + +ni = size(interp_data, dim_map(1)) +nj = size(interp_data, dim_map(2)) +nk = size(interp_data, dim_map(3)) +nlevs = size(clim_type%FMS_INTP_TYPE_%levs) +nfields = size(clim_type%field_name) + +allocate(hinterp_data(ni,nj,nlevs,nfields)) +allocate(p_fact(ni,nj)) +allocate(col_data(ni,nj,nfields)) +allocate(phalf_diff(ni,nj,nk)) + +if (use_buf) then + allocate(interp_buf(ni,nj,nk,nfields)) + call mpp_error(NOTE, "interpolator_4D: Using temporary buffer") +else + interp_buf => interp_data +endif + + do n=2,nfields if (clim_type%vert_interp(n) /= clim_type%vert_interp(n-1) .or. & clim_type%out_of_bounds(n) /= clim_type%out_of_bounds(n-1)) then if (mpp_pe() == mpp_root_pe() ) then @@ -1460,18 +1491,15 @@ if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lo endif end do - - - istart = 1 if (present(is)) istart = is -iend = istart - 1 + size(interp_data,1) +iend = istart - 1 + ni jstart = 1 if (present(js)) jstart = js -jend = jstart - 1 + size(interp_data,2) +jend = jstart - 1 + nj - do i= 1,size(clim_type%field_name(:)) + do i= 1,nfields !!++lwh if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then !--lwh @@ -1610,7 +1638,7 @@ if ( .not. clim_type%separate_time_vary_calc) then clim_type%indexm(:) = indexm clim_type%indexp(:) = indexp clim_type%climatology(:) = climatology - do i=1, size(clim_type%field_name(:)) + do i=1, nfields call read_data(clim_type,clim_type%field_name(i), & clim_type%FMS_INTP_TYPE_%pmon_pyear(:,:,:,i), & clim_type%indexm(i)+(clim_type%climatology(i)-1)*12,i,Time) @@ -1631,7 +1659,7 @@ if ( .not. clim_type%separate_time_vary_calc) then else ! We are within a climatology data set - do i=1, size(clim_type%field_name(:)) + do i=1, nfields if (taum /= clim_type%time_init(i,1) .or. & taup /= clim_type%time_init(i,2) ) then @@ -1678,7 +1706,7 @@ if ( .not. clim_type%separate_time_vary_calc) then !Set up ! field(:,:,:,1) as the previous time slice. ! field(:,:,:,2) as the next time slice. - do i=1, size(clim_type%field_name(:)) + do i=1, nfields call read_data(clim_type,clim_type%field_name(i), clim_type%FMS_INTP_TYPE_%data5d(:,:,:,1,i), taum,i,Time) clim_type%time_init(i,1) = taum clim_type%itaum = 1 @@ -1696,7 +1724,7 @@ if ( .not. clim_type%separate_time_vary_calc) then !We have the previous time step but not the next time step data clim_type%itaup = 1 if (clim_type%itaum .eq. 1 ) clim_type%itaup = 2 - do i=1, size(clim_type%field_name(:)) + do i=1, nfields call read_data(clim_type,clim_type%field_name(i), & clim_type%FMS_INTP_TYPE_%data5d(:,:,:,clim_type%itaup,i), taup,i, Time) clim_type%time_init(i,clim_type%itaup)=taup @@ -1712,7 +1740,7 @@ endif ! (.not. separate_time_vary_calc) select case(clim_type%TIME_FLAG) case (LINEAR) - do n=1, size(clim_type%field_name(:)) + do n=1, nfields hinterp_data(:,:,:,n) = (1._lkind-clim_type%FMS_INTP_TYPE_%tweight)* & clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,clim_type%itaum,n) + & clim_type%FMS_INTP_TYPE_%tweight* & @@ -1721,7 +1749,7 @@ select case(clim_type%TIME_FLAG) ! case (SEASONAL) ! Do sine fit to data at this point case (BILINEAR) - do n=1, size(clim_type%field_name(:)) + do n=1, nfields hinterp_data(:,:,:,n) = (1._lkind-clim_type%FMS_INTP_TYPE_%tweight1)*& (1._lkind-clim_type%FMS_INTP_TYPE_%tweight3)* & clim_type%FMS_INTP_TYPE_%pmon_pyear(istart:iend,jstart:jend,:,n) + & @@ -1744,7 +1772,7 @@ select case(clim_type%level_type) end select col_data(:,:,:)=0.0_lkind - do i= 1, size(clim_type%field_name(:)) + do i= 1, nfields select case(clim_type%mr(i)) case(NO_CONV) @@ -1762,7 +1790,7 @@ select case(clim_type%mr(i)) end select enddo - do i= 1, size(clim_type%field_name(:)) + do i= 1, nfields found = .false. do j = 1,size(climo_diag_name(:)) if ( trim(adjustl(lowercase(climo_diag_name(j)))) .eq. trim(adjustl(lowercase(clim_type%field_name(i))))) then @@ -1813,28 +1841,31 @@ do j = 1, size(phalf,2) endif select case(clim_type%vert_interp(i)) case(INTERP_WEIGHTED_P) - call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,:),interp_data(ilon,j,:,:)) + call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,:),interp_buf(ilon,j,:,:)) case(INTERP_LINEAR_P) - do n=1, size(clim_type%field_name(:)) - call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,n),interp_data(ilon,j,:,n)) + do n=1, nfields + call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,n),interp_buf(ilon,j,:,n)) end do ! case(INTERP_LOG) end select enddo enddo -!--lwh - do i= 1, size(clim_type%field_name(:)) - -select case(clim_type%mr(i)) - case(KG_M2) - do k = 1,size(interp_data,3) - interp_data(:,:,k,i) = interp_data(:,:,k,i)*(phalf(:,:,k+1)-phalf(:,:,k)) - enddo -end select + phalf_diff = phalf(:,:,2:nk+1) - phalf(:,:,1:nk) +!--lwh + do i= 1, nfields + select case(clim_type%mr(i)) + case(KG_M2) + interp_buf(:,:,:,i) = interp_buf(:,:,:,i)*phalf_diff + end select end do + if (use_buf) then + interp_data = reshape(interp_buf, shape(interp_data), order=dim_map) + deallocate(interp_buf) + endif + if( .not. found_field) then !field name is not in interpolator file.ERROR. call mpp_error(FATAL,"Interpolator: the field name is not contained in this & &intepolate_type: "//trim(field_name)) @@ -1856,6 +1887,7 @@ end subroutine INTERPOLATOR_4D_ !! @param [in] OPTIONAL: Index for the physics window !! @param [out] The model fields with the interpolated climatology data !! @param [out] OPTIONAL: The units of field_name +!! @param [in] OPTIONAL: Order of the x (1), y (2), and z (3) dimensions !! !! @throw FATAL "interpolator_3D : You must call interpolator_init !! before calling interpolator" @@ -1872,7 +1904,7 @@ end subroutine INTERPOLATOR_4D_ !! climatology top pressure for " !! @throw FATAL "Interpolator: the field name is not contained in this !! intepolate_type: " -subroutine INTERPOLATOR_3D_(clim_type, Time, phalf, interp_data,field_name, is,js, clim_units) +subroutine INTERPOLATOR_3D_(clim_type, Time, phalf, interp_data, field_name, is, js, clim_units, dim_order) ! ! Return 3-D field interpolated to model grid and time ! @@ -1884,6 +1916,7 @@ subroutine INTERPOLATOR_3D_(clim_type, Time, phalf, interp_data,field_name, is,j ! Time : The model time that you wish to interpolate to. ! phalf : The half level model pressure field. ! is, js : The indices of the physics window. +! dim_order : Order of the x (1), y (2), and z (3) dimensions ! ! INTENT OUT ! interp_data : The model field with the interpolated climatology data. @@ -1893,14 +1926,13 @@ type(interpolate_type), intent(inout) :: clim_type character(len=*) , intent(in) :: field_name type(time_type) , intent(in) :: Time real(FMS_INTP_KIND_), dimension(:,:,:), intent(in) :: phalf -real(FMS_INTP_KIND_), dimension(:,:,:), intent(out) :: interp_data +real(FMS_INTP_KIND_), dimension(:,:,:), intent(out), target :: interp_data integer , intent(in) , optional :: is,js character(len=*) , intent(out), optional :: clim_units +integer , intent(in), optional :: dim_order(3) integer :: taum, taup, ilon -real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),& - size(interp_data,2),size(clim_type%FMS_INTP_TYPE_%levs(:))) -real(FMS_INTP_KIND_) :: p_fact(size(interp_data,1),size(interp_data,2)) -real(FMS_INTP_KIND_) :: col_data(size(interp_data,1),size(interp_data,2)) +real(FMS_INTP_KIND_), allocatable :: hinterp_data(:,:,:), p_fact(:,:), col_data(:,:), phalf_diff(:,:,:) +real(FMS_INTP_KIND_), pointer :: interp_buf(:,:,:) real(FMS_INTP_KIND_) :: pclim(size(clim_type%FMS_INTP_TYPE_%halflevs(:))) integer :: istart,iend,jstart,jend logical :: result, found @@ -1913,22 +1945,53 @@ type(time_type) :: t_prev, t_next type(time_type), dimension(2) :: month integer :: indexm, indexp, yearm, yearp integer :: i, j, k, n +integer :: ni, nj, nk +integer :: nlevs, nfields +integer :: dim_map(3) !< x, y, z +logical :: use_buf character(len=256) :: err_msg +integer, parameter :: dim_map_default(3) = [1, 2, 3] integer, parameter :: lkind=FMS_INTP_KIND_ if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) & call mpp_error(FATAL, "interpolator_3D : You must call interpolator_init before calling interpolator") +if (present(dim_order)) then + call inverse_permutation(dim_order, dim_map) +else + dim_map = dim_map_default +endif + +use_buf = any(dim_map.ne.dim_map_default) + +ni = size(interp_data, dim_map(1)) +nj = size(interp_data, dim_map(2)) +nk = size(interp_data, dim_map(3)) +nlevs = size(clim_type%FMS_INTP_TYPE_%levs) +nfields = size(clim_type%field_name) + +allocate(hinterp_data(ni,nj,nlevs)) +allocate(p_fact(ni,nj)) +allocate(col_data(ni,nj)) +allocate(phalf_diff(ni,nj,nk)) + +if (use_buf) then + allocate(interp_buf(ni,nj,nk)) + call mpp_error(NOTE, "interpolator_3D: Using temporary buffer") +else + interp_buf => interp_data +endif + istart = 1 if (present(is)) istart = is -iend = istart - 1 + size(interp_data,1) +iend = istart - 1 + ni jstart = 1 if (present(js)) jstart = js -jend = jstart - 1 + size(interp_data,2) +jend = jstart - 1 + nj -do i= 1,size(clim_type%field_name(:)) +do i= 1,nfields !++lwh if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then !--lwh @@ -2255,23 +2318,29 @@ do j = 1, size(phalf,2) endif select case(clim_type%vert_interp(i)) case(INTERP_WEIGHTED_P) - call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_data(ilon,j,:)) + call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_buf(ilon,j,:)) case(INTERP_LINEAR_P) - call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_data(ilon,j,:)) + call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_buf(ilon,j,:)) ! case(INTERP_LOG) end select enddo enddo !--lwh +phalf_diff = phalf(:,:,2:nk+1) - phalf(:,:,1:nk) + select case(clim_type%mr(i)) case(KG_M2) - do k = 1,size(interp_data,3) - interp_data(:,:,k) = interp_data(:,:,k)*(phalf(:,:,k+1)-phalf(:,:,k)) - enddo + interp_buf(:,:,:) = interp_buf(:,:,:)*phalf_diff end select endif !field_name enddo !End of i loop + +if (use_buf) then + interp_data = reshape(interp_buf, shape(interp_data), order=dim_map) + deallocate(interp_buf) +endif + if( .not. found_field) then !field name is not in interpolator file.ERROR. call mpp_error(FATAL,"Interpolator: the field name is not contained in this & &intepolate_type: "//trim(field_name)) @@ -2292,6 +2361,7 @@ end subroutine INTERPOLATOR_3D_ !! @param [in] OPTIONAL: Index for the physics window !! @param [out] The model fields with the interpolated climatology data !! @param [out] OPTIONAL: The units of field_name +!! @param [in] OPTIONAL: Order of the x (1) and y (2) dimensions !! !! @throw FATAL "interpolator_2D : You must call interpolator_init !! before calling interpolator" @@ -2304,7 +2374,7 @@ end subroutine INTERPOLATOR_3D_ !! time but we have the next time. How did this happen?" !! @throw FATAL "Interpolator: the field name is not contained in this !! intepolate_type: " -subroutine INTERPOLATOR_2D_(clim_type, Time, interp_data, field_name, is, js, clim_units) +subroutine INTERPOLATOR_2D_(clim_type, Time, interp_data, field_name, is, js, clim_units, dim_order) ! ! Return 2-D field interpolated to model grid and time ! @@ -2316,6 +2386,7 @@ subroutine INTERPOLATOR_2D_(clim_type, Time, interp_data, field_name, is, js, cl ! field_name : The name of the field that you wish to interpolate. ! Time : The model time that you wish to interpolate to. ! is, js : The indices of the physics window. +! dim_order : Order of the x (1) and y (2) dimensions ! ! INTENT OUT ! interp_data : The model field with the interpolated climatology data. @@ -2324,11 +2395,12 @@ subroutine INTERPOLATOR_2D_(clim_type, Time, interp_data, field_name, is, js, cl type(interpolate_type), intent(inout) :: clim_type character(len=*) , intent(in) :: field_name type(time_type) , intent(in) :: Time -real(FMS_INTP_KIND_), dimension(:,:), intent(out) :: interp_data +real(FMS_INTP_KIND_), dimension(:,:), intent(out), target :: interp_data integer , intent(in) , optional :: is,js character(len=*) , intent(out), optional :: clim_units +integer , intent(in), optional :: dim_order(2) integer :: taum, taup -real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),size(interp_data,2),size(clim_type%FMS_INTP_TYPE_%levs(:))) +real(FMS_INTP_KIND_), allocatable :: hinterp_data(:,:,:) integer :: istart,iend,jstart,jend logical :: result, found logical :: found_field=.false. @@ -2340,22 +2412,39 @@ type(time_type) :: t_prev, t_next type(time_type), dimension(2) :: month integer :: indexm, indexp, yearm, yearp integer :: j, i, n +integer :: ni, nj +integer :: nlevs, nfields +integer :: dim_map(2) !< x, y character(len=256) :: err_msg +integer, parameter :: dim_map_default(2) = [1, 2] integer, parameter :: lkind=FMS_INTP_KIND_ if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) & call mpp_error(FATAL, "interpolator_2D : You must call interpolator_init before calling interpolator") +if (present(dim_order)) then + call inverse_permutation(dim_order, dim_map) +else + dim_map = dim_map_default +endif + +ni = size(interp_data, dim_map(1)) +nj = size(interp_data, dim_map(2)) +nlevs = size(clim_type%FMS_INTP_TYPE_%levs) +nfields = size(clim_type%field_name) + +allocate(hinterp_data(ni,nj,nlevs)) + istart = 1 if (present(is)) istart = is -iend = istart - 1 + size(interp_data,1) +iend = istart - 1 + ni jstart = 1 if (present(js)) jstart = js -jend = jstart - 1 + size(interp_data,2) +jend = jstart - 1 + nj -do i= 1,size(clim_type%field_name(:)) +do i= 1,nfields !++lwh if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then !--lwh @@ -2539,9 +2628,6 @@ if ( .not. clim_type%separate_time_vary_calc) then clim_type%indexp(i)+clim_type%climatology(i)*12,i,Time) endif - - - else ! We are within a climatology data set if (taum /= clim_type%time_init(i,1) .or. & @@ -2614,8 +2700,6 @@ if ( .not. clim_type%separate_time_vary_calc) then endif ! (.not. separate_time_vary_calc) - - select case(clim_type%TIME_FLAG) case (LINEAR) hinterp_data = (1._lkind-clim_type%FMS_INTP_TYPE_%tweight)*& @@ -2652,7 +2736,7 @@ if (found) then endif endif - interp_data(:,:) = hinterp_data(:,:,1) + interp_data = reshape(hinterp_data(:,:,1), shape(interp_data), order=dim_map) endif !field_name enddo !End of i loop @@ -2677,6 +2761,7 @@ end subroutine INTERPOLATOR_2D_ !! @param [in] OPTIONAL: Index for the physics window !! @param [out] The model fields with the interpolated climatology data !! @param [out] OPTIONAL: The units of field_name +!! @param [in] OPTIONAL: Order of the x (1), y (2), z (3), and field (4) dimensions !! !! @throw FATAL "interpolator_4D_no_time_axis : You must call !! interpolator_init before calling interpolator" @@ -2688,7 +2773,7 @@ end subroutine INTERPOLATOR_2D_ !! pressure of input data for " !! @throw FATAL "Interpolator: the field name is not contained in this !! intepolate_type: " -subroutine INTERPOLATOR_4D_NO_TIME_AXIS_(clim_type, phalf, interp_data, field_name, is,js, clim_units) +subroutine INTERPOLATOR_4D_NO_TIME_AXIS_(clim_type, phalf, interp_data, field_name, is, js, clim_units, dim_order) ! Return 4-D field interpolated to model grid @@ -2702,6 +2787,7 @@ subroutine INTERPOLATOR_4D_NO_TIME_AXIS_(clim_type, phalf, interp_data, field_na ! be any one of the variables. ! phalf : The half level model pressure field. ! is, js : The indices of the physics window. +! dim_order : Order of the x (1), y (2), z (3), and field (4) dimensions ! INTENT OUT ! interp_data : The model fields @@ -2710,24 +2796,54 @@ subroutine INTERPOLATOR_4D_NO_TIME_AXIS_(clim_type, phalf, interp_data, field_na type(interpolate_type), intent(inout) :: clim_type character(len=*) , intent(in) :: field_name real(FMS_INTP_KIND_), dimension(:,:,:), intent(in) :: phalf -real(FMS_INTP_KIND_), dimension(:,:,:,:), intent(out) :: interp_data +real(FMS_INTP_KIND_), dimension(:,:,:,:), intent(out), target :: interp_data integer , intent(in) , optional :: is,js character(len=*) , intent(out), optional :: clim_units +integer , intent(in), optional :: dim_order(4) integer :: ilon -real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),& - size(interp_data,2),size(clim_type%FMS_INTP_TYPE_%levs(:)),size(clim_type%field_name(:))) -real(FMS_INTP_KIND_) :: p_fact(size(interp_data,1),size(interp_data,2)) +real(FMS_INTP_KIND_), allocatable :: hinterp_data(:,:,:,:), p_fact(:,:), phalf_diff(:,:,:) +real(FMS_INTP_KIND_), pointer :: interp_buf(:,:,:,:) real(FMS_INTP_KIND_) :: pclim(size(clim_type%FMS_INTP_TYPE_%halflevs(:))) integer :: istart,iend,jstart,jend logical :: found_field=.false. integer :: i, j, k, n +integer :: ni, nj, nk +integer :: nlevs, nfields +integer :: dim_map(4) !< x, y, z, field +logical :: use_buf +integer, parameter :: dim_map_default(4) = [1, 2, 3, 4] integer, parameter :: lkind=FMS_INTP_KIND_ if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) & call mpp_error(FATAL, "interpolator_4D_no_time_axis : You must call interpolator_init before calling interpolator") -do n=2,size(clim_type%field_name(:)) +if (present(dim_order)) then + call inverse_permutation(dim_order, dim_map) +else + dim_map = dim_map_default +endif + +use_buf = any(dim_map.ne.dim_map_default) + +ni = size(interp_data, dim_map(1)) +nj = size(interp_data, dim_map(2)) +nk = size(interp_data, dim_map(3)) +nlevs = size(clim_type%FMS_INTP_TYPE_%levs) +nfields = size(clim_type%field_name) + +allocate(hinterp_data(ni,nj,nlevs,nfields)) +allocate(p_fact(ni,nj)) +allocate(phalf_diff(ni,nj,nk)) + +if (use_buf) then + allocate(interp_buf(ni,nj,nk,nfields)) + call mpp_error(NOTE, "interpolator_4D_no_time_axis: Using temporary buffer") +else + interp_buf => interp_data +endif + +do n=2,nfields if (clim_type%vert_interp(n) /= clim_type%vert_interp(n-1) .or. & clim_type%out_of_bounds(n) /= clim_type%out_of_bounds(n-1)) then if (mpp_pe() == mpp_root_pe() ) then @@ -2740,13 +2856,13 @@ end do istart = 1 if (present(is)) istart = is -iend = istart - 1 + size(interp_data,1) +iend = istart - 1 + ni jstart = 1 if (present(js)) jstart = js -jend = jstart - 1 + size(interp_data,2) +jend = jstart - 1 + nj -do i= 1,size(clim_type%field_name(:)) +do i= 1,nfields if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then found_field=.true. exit @@ -2759,7 +2875,7 @@ if(present(clim_units)) then clim_units = chomp(clim_units) endif -do n=1, size(clim_type%field_name(:)) +do n=1, nfields hinterp_data(:,:,:,n) = clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,1,n) end do @@ -2770,7 +2886,7 @@ select case(clim_type%level_type) p_fact = maxval(phalf,3)! max pressure in the column !(:,:,size(phalf,3)) end select - do i= 1, size(clim_type%field_name(:)) + do i= 1, nfields select case(clim_type%mr(i)) case(KG_M2) do k = 1,size(hinterp_data,3) @@ -2809,26 +2925,29 @@ do j = 1, size(phalf,2) endif select case(clim_type%vert_interp(i)) case(INTERP_WEIGHTED_P) - call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,:),interp_data(ilon,j,:,:)) + call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,:),interp_buf(ilon,j,:,:)) case(INTERP_LINEAR_P) - do n=1, size(clim_type%field_name(:)) - call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,n),interp_data(ilon,j,:,n)) - end do + do n=1, nfields + call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,n),interp_buf(ilon,j,:,n)) + end do end select enddo enddo - do i= 1, size(clim_type%field_name(:)) - -select case(clim_type%mr(i)) - case(KG_M2) - do k = 1,size(interp_data,3) - interp_data(:,:,k,i) = interp_data(:,:,k,i)*(phalf(:,:,k+1)-phalf(:,:,k)) - enddo -end select + phalf_diff = phalf(:,:,2:nk+1) - phalf(:,:,1:nk) + do i= 1, nfields + select case(clim_type%mr(i)) + case(KG_M2) + interp_buf(:,:,:,i) = interp_buf(:,:,:,i)*phalf_diff + end select end do + if (use_buf) then + interp_data = reshape(interp_buf, shape(interp_data), order=dim_map) + deallocate(interp_buf) + endif + if( .not. found_field) then !field name is not in interpolator file.ERROR. call mpp_error(FATAL,"Interpolator: the field name is not contained in this & &intepolate_type: "//trim(field_name)) @@ -2848,6 +2967,7 @@ end subroutine INTERPOLATOR_4D_NO_TIME_AXIS_ !! @param [in] OPTIONAL: Index for the physics window !! @param [out] The model fields with the interpolated climatology data !! @param [out] OPTIONAL: The units of field_name +!! @param [in] OPTIONAL: Order of the x (1), y (2), and z (3) dimensions !! !! @throw FATAL "interpolator_3D_no_time_axis : You must call !! interpolator_init before calling interpolator" @@ -2859,7 +2979,7 @@ end subroutine INTERPOLATOR_4D_NO_TIME_AXIS_ !! climatology top pressure for " !! @throw FATAL "Interpolator: the field name is not contained in this !! intepolate_type: " -subroutine INTERPOLATOR_3D_NO_TIME_AXIS_(clim_type, phalf, interp_data, field_name, is,js, clim_units) +subroutine INTERPOLATOR_3D_NO_TIME_AXIS_(clim_type, phalf, interp_data, field_name, is, js, clim_units, dim_order) ! Return 3-D field interpolated to model grid @@ -2870,6 +2990,7 @@ subroutine INTERPOLATOR_3D_NO_TIME_AXIS_(clim_type, phalf, interp_data, field_na ! field_name : The name of the field that you wish to interpolate. ! phalf : The half level model pressure field. ! is, js : The indices of the physics window. +! dim_order : Order of the x (1), y (2), and z (3) dimensions ! INTENT OUT ! interp_data : The model field with the interpolated climatology data. @@ -2878,32 +2999,62 @@ subroutine INTERPOLATOR_3D_NO_TIME_AXIS_(clim_type, phalf, interp_data, field_na type(interpolate_type), intent(inout) :: clim_type character(len=*) , intent(in) :: field_name real(FMS_INTP_KIND_), dimension(:,:,:), intent(in) :: phalf -real(FMS_INTP_KIND_), dimension(:,:,:), intent(out) :: interp_data +real(FMS_INTP_KIND_), dimension(:,:,:), intent(out), target :: interp_data integer , intent(in) , optional :: is,js character(len=*) , intent(out), optional :: clim_units +integer , intent(in), optional :: dim_order(3) integer :: ilon !< No description -real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),& - size(interp_data,2),size(clim_type%FMS_INTP_TYPE_%levs(:))) !< No description -real(FMS_INTP_KIND_) :: p_fact(size(interp_data,1),size(interp_data,2)) !< No description +real(FMS_INTP_KIND_), allocatable :: hinterp_data(:,:,:), p_fact(:,:), phalf_diff(:,:,:) +real(FMS_INTP_KIND_), pointer :: interp_buf(:,:,:) real(FMS_INTP_KIND_) :: pclim(size(clim_type%FMS_INTP_TYPE_%halflevs(:))) !< No description integer :: istart,iend,jstart,jend !< No description logical :: found_field=.false. !< No description integer :: i, j, k !< No description +integer :: ni, nj, nk +integer :: nlevs, nfields +integer :: dim_map(3) !< x, y, z +logical :: use_buf +integer, parameter :: dim_map_default(3) = [1, 2, 3] integer, parameter :: lkind=FMS_INTP_KIND_ if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) & call mpp_error(FATAL, "interpolator_3D_no_time_axis : You must call interpolator_init before calling interpolator") +if (present(dim_order)) then + call inverse_permutation(dim_order, dim_map) +else + dim_map = dim_map_default +endif + +use_buf = any(dim_map.ne.dim_map_default) + +ni = size(interp_data, dim_map(1)) +nj = size(interp_data, dim_map(2)) +nk = size(interp_data, dim_map(3)) +nlevs = size(clim_type%FMS_INTP_TYPE_%levs) +nfields = size(clim_type%field_name) + +allocate(hinterp_data(ni,nj,nlevs)) +allocate(p_fact(ni,nj)) +allocate(phalf_diff(ni,nj,nk)) + +if (use_buf) then + allocate(interp_buf(ni,nj,nk)) + call mpp_error(NOTE, "interpolator_3D_no_time_axis: Using temporary buffer") +else + interp_buf => interp_data +endif + istart = 1 if (present(is)) istart = is -iend = istart - 1 + size(interp_data,1) +iend = istart - 1 + ni jstart = 1 if (present(js)) jstart = js -jend = jstart - 1 + size(interp_data,2) +jend = jstart - 1 + nj -do i= 1,size(clim_type%field_name(:)) +do i= 1,nfields if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then found_field=.true. if(present(clim_units)) then @@ -2955,20 +3106,25 @@ do j = 1, size(phalf,2) endif select case(clim_type%vert_interp(i)) case(INTERP_WEIGHTED_P) - call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_data(ilon,j,:)) + call interp_weighted_scalar(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_buf(ilon,j,:)) case(INTERP_LINEAR_P) - call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_data(ilon,j,:)) + call interp_linear(pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_buf(ilon,j,:)) end select enddo enddo +phalf_diff = phalf(:,:,2:nk+1) - phalf(:,:,1:nk) + select case(clim_type%mr(i)) case(KG_M2) - do k = 1,size(interp_data,3) - interp_data(:,:,k) = interp_data(:,:,k)*(phalf(:,:,k+1)-phalf(:,:,k)) - enddo + interp_buf(:,:,:) = interp_buf(:,:,:)*phalf_diff end select +if (use_buf) then + interp_data = reshape(interp_buf, shape(interp_data), order=dim_map) + deallocate(interp_buf) +endif + endif !field_name enddo !End of i loop if( .not. found_field) then !field name is not in interpolator file.ERROR. @@ -2989,12 +3145,13 @@ end subroutine INTERPOLATOR_3D_NO_TIME_AXIS_ !! @param [in] OPTIONAL: Index for the physics window !! @param [out] The model fields with the interpolated climatology data !! @param [out] OPTIONAL: The units of field_name +!! @param [in] OPTIONAL: Order of the x (1) and y (2) dimensions !! !! @throw FATAL "interpolator_2D_no_time_axis : You must call !! interpolator_init before calling interpolator" !! @throw FATAL "Interpolator: the field name is not contained in this !! intepolate_type: " -subroutine INTERPOLATOR_2D_NO_TIME_AXIS_(clim_type, interp_data, field_name, is, js, clim_units) +subroutine INTERPOLATOR_2D_NO_TIME_AXIS_(clim_type, interp_data, field_name, is, js, clim_units, dim_order) ! Return 2-D field interpolated to model grid @@ -3004,6 +3161,7 @@ subroutine INTERPOLATOR_2D_NO_TIME_AXIS_(clim_type, interp_data, field_name, is, ! INTENT IN ! field_name : The name of the field that you wish to interpolate. ! is, js : The indices of the physics window. +! dim_order : Order of the x (1) and y (2) dimensions ! INTENT OUT ! interp_data : The model field with the interpolated climatology data. @@ -3011,27 +3169,42 @@ subroutine INTERPOLATOR_2D_NO_TIME_AXIS_(clim_type, interp_data, field_name, is, type(interpolate_type), intent(inout) :: clim_type character(len=*) , intent(in) :: field_name -real(FMS_INTP_KIND_), dimension(:,:), intent(out) :: interp_data +real(FMS_INTP_KIND_), dimension(:,:), intent(out), target :: interp_data integer , intent(in) , optional :: is,js character(len=*) , intent(out), optional :: clim_units -real(FMS_INTP_KIND_) :: hinterp_data(size(interp_data,1),& - size(interp_data,2),size(clim_type%FMS_INTP_TYPE_%levs(:))) +integer , intent(in), optional :: dim_order(2) integer :: istart,iend,jstart,jend logical :: found_field=.false. integer :: i +integer :: ni, nj +integer :: nlevs, nfields +integer :: dim_map(2) !< x, y + +integer, parameter :: dim_map_default(2) = [1, 2] if (.not. module_is_initialized .or. .not. allocated(clim_type%FMS_INTP_TYPE_%lon)) & call mpp_error(FATAL, "interpolator_2D_no_time_axis : You must call interpolator_init before calling interpolator") +if (present(dim_order)) then + call inverse_permutation(dim_order, dim_map) +else + dim_map = dim_map_default +endif + +ni = size(interp_data, dim_map(1)) +nj = size(interp_data, dim_map(2)) +nlevs = size(clim_type%FMS_INTP_TYPE_%levs) +nfields = size(clim_type%field_name) + istart = 1 if (present(is)) istart = is -iend = istart - 1 + size(interp_data,1) +iend = istart - 1 + ni jstart = 1 if (present(js)) jstart = js -jend = jstart - 1 + size(interp_data,2) +jend = jstart - 1 + nj -do i= 1,size(clim_type%field_name(:)) +do i= 1,nfields if ( trim(adjustl(lowercase(field_name))) == trim(adjustl(lowercase(clim_type%field_name(i)))) ) then found_field=.true. @@ -3041,9 +3214,8 @@ do i= 1,size(clim_type%field_name(:)) clim_units = chomp(clim_units) endif - hinterp_data = clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,:,1,i) - - interp_data(:,:) = hinterp_data(:,:,1) + interp_data = reshape(clim_type%FMS_INTP_TYPE_%data5d(istart:iend,jstart:jend,1,1,i), & + shape(interp_data), order=dim_map) endif !field_name enddo !End of i loop diff --git a/interpolator/interpolator.F90 b/interpolator/interpolator.F90 index e37528e83d..817ac1b712 100644 --- a/interpolator/interpolator.F90 +++ b/interpolator/interpolator.F90 @@ -768,6 +768,33 @@ function chomp(string) chomp = string(:len) end function chomp + +!> @brief Produce an inverse permutation. For example, transform [2, 3, 1, 4] into [3, 1, 2, 4]. +!! +!! @param [in] The original permutation vector +!! @param [out] The inverted permutation vector +subroutine inverse_permutation(x, y) + use fms_string_utils_mod, only: string + integer, intent(in) :: x(:) + integer, intent(out) :: y(size(x)) + integer :: i, n + + y = 0 + + n = size(x) + do i=1,n + if (x(i).ge.1 .and. x(i).le.n) then + y(x(i)) = i + else + call mpp_error(FATAL, "interpolator_mod: Invalid dim_order. Values must be in the range from 1 through " & + // string(n) // ".") + endif + enddo + + if (any(y.eq.0)) then + call mpp_error(FATAL, "interpolator_mod: Invalid dim_order. Values must be non-repeating.") + endif +end subroutine inverse_permutation ! !######################################################################## diff --git a/mpp/include/mpp_gather.fh b/mpp/include/mpp_gather.fh index 54414da5fd..6c40dc4e8f 100644 --- a/mpp/include/mpp_gather.fh +++ b/mpp/include/mpp_gather.fh @@ -105,9 +105,9 @@ subroutine MPP_GATHER_PELIST_2D_(is, ie, js, je, pelist, array_seg, gather_data, logical, intent(in) :: is_root_pe integer, optional, intent(in) :: ishift, jshift - integer, dimension(3) :: dim_order + integer, dimension(2) :: dim_order - dim_order = (/1,2,3/) + dim_order = (/1,2/) call mpp_gather(is, ie, js, je, pelist, array_seg, gather_data, dim_order, is_root_pe, & ishift, jshift) @@ -121,7 +121,7 @@ subroutine MPP_GATHER_PELIST_GEN_2D_(is, ie, js, je, pelist, array_seg, gather_d integer, dimension(:), intent(in) :: pelist MPP_TYPE_, dimension(:,:), contiguous, target, intent(in) :: array_seg MPP_TYPE_, dimension(:,:), contiguous, target, intent(inout) :: gather_data - integer, dimension(3), intent(in) :: dim_order + integer, dimension(2), intent(in) :: dim_order logical, intent(in) :: is_root_pe integer, optional, intent(in) :: ishift, jshift @@ -135,7 +135,7 @@ subroutine MPP_GATHER_PELIST_GEN_2D_(is, ie, js, je, pelist, array_seg, gather_d data3D => null() endif - call mpp_gather(is, ie, js, je, 1, pelist, arr3D, data3D, dim_order, is_root_pe, & + call mpp_gather(is, ie, js, je, 1, pelist, arr3D, data3D, [dim_order, 3], is_root_pe, & ishift, jshift) return @@ -212,7 +212,11 @@ subroutine MPP_GATHER_PELIST_GEN_3D_(is, ie, js, je, nk, pelist, array_seg, gath if (present(jshift)) joff=jshift ! gather indices into global index on root_pe - if (is_root_pe) allocate(gind(4*size(pelist))) + if (is_root_pe) then + allocate(gind(4*size(pelist))) + else + allocate(gind(1)) + endif call mpp_gather((/is, ie, js, je/), gind, pelist) ! Compute recv counts and allocate 1d recv buffer (rbuf) diff --git a/mpp/include/mpp_scatter.fh b/mpp/include/mpp_scatter.fh index cca548ec16..b0e5c5d126 100644 --- a/mpp/include/mpp_scatter.fh +++ b/mpp/include/mpp_scatter.fh @@ -31,9 +31,9 @@ subroutine MPP_SCATTER_PELIST_2D_(is, ie, js, je, pelist, array_seg, input_data, MPP_TYPE_, dimension(:,:), contiguous, target, intent(in) :: input_data !< 2D array of input data logical, intent(in) :: is_root_pe !< operational root pe - integer, dimension(3) :: dim_order + integer, dimension(2) :: dim_order - dim_order = (/1,2,3/) + dim_order = (/1,2/) call mpp_scatter(is, ie, js, je, pelist, array_seg, input_data, dim_order, is_root_pe) @@ -47,7 +47,7 @@ subroutine MPP_SCATTER_PELIST_GEN_2D_(is, ie, js, je, pelist, array_seg, input_d !! must be in monotonic increasing order MPP_TYPE_, dimension(:,:), contiguous, target, intent(inout) :: array_seg !< 2D array of output data MPP_TYPE_, dimension(:,:), contiguous, target, intent(in) :: input_data !< 2D array of input data - integer, dimension(3), intent(in) :: dim_order + integer, dimension(2), intent(in) :: dim_order logical, intent(in) :: is_root_pe !< operational root pe MPP_TYPE_, pointer :: arr3D(:,:,:) @@ -60,7 +60,7 @@ subroutine MPP_SCATTER_PELIST_GEN_2D_(is, ie, js, je, pelist, array_seg, input_d data3D => null() endif - call mpp_scatter(is, ie, js, je, 1, pelist, arr3D, data3D, dim_order, is_root_pe) + call mpp_scatter(is, ie, js, je, 1, pelist, arr3D, data3D, [dim_order, 3], is_root_pe) return diff --git a/mpp/include/mpp_write_unlimited_axis.fh b/mpp/include/mpp_write_unlimited_axis.fh index 130c71d852..7626574298 100644 --- a/mpp/include/mpp_write_unlimited_axis.fh +++ b/mpp/include/mpp_write_unlimited_axis.fh @@ -43,7 +43,11 @@ nelems = sum(nelems_io(:)) - if(mpp_file(unit)%write_on_this_pe) allocate(rbuff(nelems)) + if(mpp_file(unit)%write_on_this_pe) then + allocate(rbuff(nelems)) + else + allocate(rbuff(1)) + endif ! Note that the gatherV implied here is asymmetric; only root needs to know the vector of recv size call mpp_gather(data,size(data),rbuff,nelems_io(:),pelist) diff --git a/test_fms/diag_manager/test_multiple_zbounds.F90 b/test_fms/diag_manager/test_multiple_zbounds.F90 index 1ccb6eaca8..03848d091a 100644 --- a/test_fms/diag_manager/test_multiple_zbounds.F90 +++ b/test_fms/diag_manager/test_multiple_zbounds.F90 @@ -29,9 +29,12 @@ program test_multiple_zbounds type(time_type) :: Time_step !< Time_step of the simulation integer :: id_z1 !< Axis id for the z dimension integer :: id_z2 !< Axis id for the z dimension + integer :: id_z_reverse !< Axis id for the z dimension (decreasing) integer :: id_var1 !< var id for the first variable integer :: id_var2 !< var_id for the second variable + integer :: id_var3 !< var_id for the third variable real, allocatable :: z(:) !< z axis data + real, allocatable :: z_reverse(:) !< z axis data (decreasing) integer :: nz !< Size of the z dimension integer :: i logical :: used !< Dummy argument to send_data @@ -42,8 +45,10 @@ program test_multiple_zbounds nz = 10 allocate(z(nz)) + allocate(z_reverse(nz)) do i=1, nz z(i) = i + z_reverse(i) = nz - i + 1 enddo Time = set_date(2,1,1,0,0,0) @@ -51,16 +56,17 @@ program test_multiple_zbounds id_z1 = diag_axis_init('zaxis1', z, 'z', 'z', long_name='Z1') id_z2 = diag_axis_init('zaxis2', z, 'z', 'z', long_name='Z2') + id_z_reverse = diag_axis_init('zaxis3', z_reverse, 'z_reverse', 'z', long_name='Z3') id_var1 = register_diag_field ('atmos', 'ua_1', (/id_z1/), Time) id_var2 = register_diag_field ('atmos', 'ua_2', (/id_z2/), Time) + id_var3 = register_diag_field ('atmos', 'ua_3', (/id_z_reverse/), Time) call diag_manager_set_time_end(set_date(2,1,2,0,0,0)) do i = 1, 24 Time = Time + Time_step - z = real(i) used = send_data(id_var1, z, Time) used = send_data(id_var2, z, Time) - + used = send_data(id_var3, z, Time) call diag_send_complete(Time_step) end do @@ -77,6 +83,9 @@ subroutine check_output() integer :: SUB1_SIZE = 3 integer :: SUB2_SIZE = 1 character(len=20) :: EXPECTED_DIM_NAMES(2) + real :: EXPECTED_ZSUBAXIS_1(3) + real :: EXPECTED_ZSUBAXIS_2(1) + real :: EXPECTED_ZSUBAXIS_3(3) EXPECTED_DIM_NAMES(2) = "time" if (.not. open_file(fileobj, "test_multiple_zbounds.nc", "read")) then @@ -84,8 +93,15 @@ subroutine check_output() endif call check_dimension(fileobj, "time", EXPECTED_NTIMES) - call check_dimension(fileobj, "zaxis1_sub01", SUB1_SIZE) - call check_dimension(fileobj, "zaxis2_sub02", SUB2_SIZE) + + EXPECTED_ZSUBAXIS_1 = (/3., 4., 5./) + call check_dimension(fileobj, "zaxis1_sub01", SUB1_SIZE, EXPECTED_ZSUBAXIS_1) + + EXPECTED_ZSUBAXIS_2 = (/1./) + call check_dimension(fileobj, "zaxis2_sub02", SUB2_SIZE, EXPECTED_ZSUBAXIS_2) + + EXPECTED_ZSUBAXIS_3 = (/5., 4., 3./) + call check_dimension(fileobj, "zaxis3_sub03", SUB1_SIZE, EXPECTED_ZSUBAXIS_3) EXPECTED_DIM_NAMES(1) = "zaxis1_sub01" call check_variable(fileobj, "ua_1", EXPECTED_DIM_NAMES) @@ -114,17 +130,40 @@ subroutine check_variable(fileobj, variable_name, expected_dimnames) end subroutine check_variable - subroutine check_dimension(fileobj, dimension_name, expected_size) + !> @brief Check dimension data + subroutine check_data(err_msg, actual_data, expected_data) + character(len=*), intent(in) :: err_msg !< Error message to append + real, intent(in) :: actual_data(:) !< Dimension data from file + real, intent(in) :: expected_data(:) !< Expected data + + integer :: i + + do i = 1, size(actual_data) + if (actual_data(i) .ne. expected_data(i)) & + call mpp_error(FATAL, "The data is not expected for "//trim(err_msg)) + enddo + end subroutine check_data + + subroutine check_dimension(fileobj, dimension_name, expected_size, expected_data) type(FmsNetcdfFile_t), intent(in) :: fileobj character(len=*), intent(in) :: dimension_name integer, intent(in) :: expected_size + real, optional, intent(in) :: expected_data(:) integer :: dim_size + real, allocatable :: z_data(:) call get_dimension_size(fileobj, dimension_name, dim_size) if (dim_size .ne. expected_size) then call mpp_error(FATAL, trim(dimension_name)//" is not the expected size!") endif + if (present(expected_data)) then + allocate(z_data(dim_size)) + call read_data(fileobj, dimension_name, z_data) + call check_data(dimension_name, z_data, expected_data) + endif + + end subroutine end program test_multiple_zbounds \ No newline at end of file diff --git a/test_fms/diag_manager/test_multiple_zbounds.sh b/test_fms/diag_manager/test_multiple_zbounds.sh index 44db78c036..29c79d06e8 100755 --- a/test_fms/diag_manager/test_multiple_zbounds.sh +++ b/test_fms/diag_manager/test_multiple_zbounds.sh @@ -47,6 +47,8 @@ diag_files: zbounds: 3 5 - var_name: ua_2 zbounds: 1 1 + - var_name: ua_3 + zbounds: 3 5 _EOF test_expect_success "Test with multiple zbounds limits (modern diag manager)" ' diff --git a/test_fms/interpolator/test_interpolator2.F90 b/test_fms/interpolator/test_interpolator2.F90 index 1cf450dd1b..354bf86136 100644 --- a/test_fms/interpolator/test_interpolator2.F90 +++ b/test_fms/interpolator/test_interpolator2.F90 @@ -40,7 +40,7 @@ program test_interpolator2 use constants_mod, only: PI use platform_mod, only: r4_kind, r8_kind use fms_test_mod, only: permutable_indices_2d, permutable_indices_3d, permutable_indices_4d, factorial, & - arr_compare_tol + permute_arr, arr_compare_tol use interpolator_mod implicit none @@ -156,26 +156,28 @@ subroutine test_interpolator_2d(clim_type, phalf_in, test_time, answer) real(TEST_FMS_KIND_), allocatable, dimension(:,:) :: interp_data type(permutable_indices_2d) :: dims - integer :: p + integer :: p, dim_order(2) do p=1,factorial(2) + dim_order = [1, 2] dims%lb = [1, 1] dims%ub = [nlonlat, nlonlat] + call permute_arr(dim_order, p) call dims%permute(p) allocate(interp_data(dims%ub(1), dims%ub(2))) interp_data = 0. !> test interpolator_2D_r4/8 - call interpolator(clim_type, test_time, interp_data, 'ozone') + call interpolator(clim_type, test_time, interp_data, 'ozone', dim_order=dim_order) call arr_compare_tol(interp_data, answer, tol, 'test interpolator_2D') interp_data = 0. !> Test obtain_interpolator_time_slices call obtain_interpolator_time_slices(clim_type,test_time) - call interpolator(clim_type, test_time, interp_data, 'ozone') + call interpolator(clim_type, test_time, interp_data, 'ozone', dim_order=dim_order) call unset_interpolator_time_flag(clim_type) call arr_compare_tol(interp_data, answer, tol, 'test interpolator_2D') @@ -191,19 +193,21 @@ subroutine test_interpolator_3d(clim_type, phalf_in, test_time, answer) real(TEST_FMS_KIND_), allocatable, dimension(:,:,:) :: interp_data type(permutable_indices_3d) :: dims - integer :: p + integer :: p, dim_order(3) do p=1,factorial(3) + dim_order = [1, 2, 3] dims%lb = [1, 1, 1] dims%ub = [nlonlat, nlonlat, npfull] + call permute_arr(dim_order, p) call dims%permute(p) allocate(interp_data(dims%ub(1), dims%ub(2), dims%ub(3))) interp_data = 0. !> test interpolator_3_r4/8 - call interpolator(clim_type, test_time, phalf_in, interp_data, 'ozone') + call interpolator(clim_type, test_time, phalf_in, interp_data, 'ozone', dim_order=dim_order) call arr_compare_tol(interp_data, answer, tol, 'test interpolator_3D') deallocate(interp_data) @@ -218,19 +222,21 @@ subroutine test_interpolator_4d(clim_type, phalf_in, test_time, answer) real(TEST_FMS_KIND_), allocatable, dimension(:,:,:,:) :: interp_data type(permutable_indices_4d) :: dims - integer :: p + integer :: p, dim_order(4) do p=1,factorial(4) + dim_order = [1, 2, 3, 4] dims%lb = [1, 1, 1, 1] dims%ub = [nlonlat, nlonlat, npfull, 1] + call permute_arr(dim_order, p) call dims%permute(p) allocate(interp_data(dims%ub(1), dims%ub(2), dims%ub(3), dims%ub(4))) interp_data = 0. !> test interpolator_4D_r4/8 - call interpolator(clim_type, test_time, phalf_in, interp_data, 'ozone') + call interpolator(clim_type, test_time, phalf_in, interp_data, 'ozone', dim_order=dim_order) call arr_compare_tol(interp_data, answer, tol, 'test interpolator_4D') deallocate(interp_data) diff --git a/test_fms/interpolator/test_interpolator2.sh b/test_fms/interpolator/test_interpolator2.sh index b43fd605c7..3178280190 100755 --- a/test_fms/interpolator/test_interpolator2.sh +++ b/test_fms/interpolator/test_interpolator2.sh @@ -26,16 +26,6 @@ # Set common test settings. . ../test-lib.sh -# TODO: Enable these tests once generalized indices work is complete -SKIP_TESTS="test_interpolator2.2 \ - test_interpolator2.3 \ - test_interpolator2.4 \ - test_interpolator2.5 \ - test_interpolator2.6 \ - test_interpolator2.7 \ - test_interpolator2.8 \ - test_interpolator2.9" - # Tests to skip if input files not present if test ! -z "$test_input_path" ; then rm -rf INPUT && mkdir INPUT diff --git a/test_fms/mpp/Makefile.am b/test_fms/mpp/Makefile.am index 4d12b04a18..c11829d4c0 100644 --- a/test_fms/mpp/Makefile.am +++ b/test_fms/mpp/Makefile.am @@ -71,7 +71,8 @@ check_PROGRAMS = test_mpp \ test_mpp_clock_begin_end_id \ test_mpp_nesting \ test_mpp_chksum \ - test_stdlog + test_stdlog \ + test_corner_mosaic # These are the sources for the tests. test_mpp_SOURCES = test_mpp.F90 @@ -140,6 +141,7 @@ test_mpp_clock_begin_end_id_SOURCES=test_mpp_clock_begin_end_id.F90 test_super_grid_SOURCES = test_super_grid.F90 test_mpp_chksum_SOURCES = test_mpp_chksum.F90 test_stdlog_SOURCES = test_stdlog.F90 +test_corner_mosaic_SOURCES = test_corner_mosaic.F90 test_mpp_global_field_r4_CPPFLAGS = $(AM_CPPFLAGS) -DFMS_TEST_TYPE_=real -DFMS_TEST_KIND_=r4_kind test_mpp_global_field_r8_CPPFLAGS = $(AM_CPPFLAGS) -DFMS_TEST_TYPE_=real -DFMS_TEST_KIND_=r8_kind diff --git a/test_fms/mpp/test_corner_mosaic.F90 b/test_fms/mpp/test_corner_mosaic.F90 new file mode 100644 index 0000000000..4355f1008a --- /dev/null +++ b/test_fms/mpp/test_corner_mosaic.F90 @@ -0,0 +1,446 @@ +!*********************************************************************** +!* 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. +!*********************************************************************** +! Corner mosaic test with resource reallocation +! This test constructs a single-tile domain, then a three-tile mosaic, then a single-tile domain once more. +! The number of PEs must be such that the single-tile domain uses all of the PEs (e.g., 32 PEs if the single-tile +! domain is 8x4). The layout of the corner mosaic is defined by ni(:) and nj(:). The corner mosaic does not need to +! use all of the available PEs. + !!!!!!!!! + ! --- ! + ! /|3| ! + ! ----- ! + ! |1|2| ! + ! ----- ! + !!!!!!!!! + +program test_corner_mosaic + use mpp_mod, only: mpp_error, FATAL + use mpp_mod, only: mpp_npes, mpp_pe, mpp_root_pe, mpp_declare_pelist, mpp_set_current_pelist + use fms_mod, only: fms_init, fms_end + use mpp_domains_mod, only: domain2d, mpp_define_mosaic, mpp_get_compute_domain, mpp_get_data_domain, & + mpp_get_tile_id, mpp_update_domains, mpp_domains_set_stack_size, CGRID_NE, & + DGRID_NE, mpp_define_domains, mpp_deallocate_domain + use platform_mod + use fms_string_utils_mod, only: string + use random_numbers_mod + implicit none + + integer, parameter :: ntiles = 3 + integer, parameter :: halo = 2 + + character(*), parameter :: name = "corner mosaic" + integer, dimension(ntiles) :: pe_start, pe_end + + integer, parameter :: n=16 !< Number of grid points along each axis on corner mosaic tiles + integer :: ni(3), nj(3) !< Layouts for the tiles comprising the corner mosaic + integer :: ni_st=8, nj_st=4 !< Layout for the single-tile domain + + real(r4_kind), parameter :: eps = 1e-3 + integer, parameter :: stackmax=10000000 + + integer :: isg, ieg, jsg, jeg + integer :: isd, ied, jsd, jed + integer :: isc, iec, jsc, jec + integer :: tid + + call fms_init + call mpp_domains_set_stack_size(stackmax) + + call define_single_tile(ni_st, nj_st) + + ! Corner mosaic layout: + ! - Tile 1: 4x4 + ! - Tile 2: 4x2 + ! - Tile 3: 4x2 + + ni(:) = [4, 4, 4] ! Corner mosaic layout (i axis) + nj(:) = [4, 2, 2] ! Corner mosaic layout (j axis) + call define_corner_mosaic(ni, nj) + call mpp_set_current_pelist ! Restore original pelist + + ni(:) = [2, 4, 2] ! Corner mosaic layout (i axis) + nj(:) = [2, 2, 2] ! Corner mosaic layout (j axis) + call define_corner_mosaic(ni, nj) + call mpp_set_current_pelist ! Restore original pelist + + ni(:) = [2, 6, 2] ! Corner mosaic layout (i axis) + nj(:) = [4, 2, 6] ! Corner mosaic layout (j axis) + call define_corner_mosaic(ni, nj) + call mpp_set_current_pelist ! Restore original pelist + + call define_single_tile(ni_st, nj_st) + + call fms_end + +contains + + subroutine define_single_tile(ni, nj) + integer, intent(in) :: ni, nj + type(domain2d) :: domain + integer, dimension(:), allocatable :: tile_id + integer :: npes + + npes = ni*nj + if (mpp_npes().ne.npes) then + call mpp_error(FATAL, "Cannot construct single-tile domain: Incorrect number of PEs") + endif + + isg=1 + ieg=2*n + jsg=1 + jeg=2*n + + call mpp_define_domains([isg,ieg,jsg,jeg], [ni, nj], domain, & + symmetry=.true., whalo=halo, ehalo=halo, shalo=halo, nhalo=halo) + + call mpp_get_compute_domain(domain, isc, iec, jsc, jec) + call mpp_get_data_domain(domain, isd, ied, jsd, jed) + + call single_tile_scalar_test(domain) + call mpp_deallocate_domain(domain) + end subroutine define_single_tile + + subroutine single_tile_scalar_test(domain) + type(domain2d), intent(inout) :: domain + real(r4_kind) :: global((isg-halo):(ieg+halo), (jsg-halo):(jeg+halo)) + real(r4_kind) :: x(isd:ied, jsd:jed) + integer :: i, j, h + + do concurrent (i=lbound(global,1):ubound(global,1), & + j=lbound(global,2):ubound(global,2)) + h = 0._r4_kind + if (i.lt.isg .or. i.gt.ieg) h = 999000._r4_kind + if (j.lt.jsg .or. j.gt.jeg) h = 999000._r4_kind + global(i,j) = i*1._r4_kind + j*.01_r4_kind + h + enddo + + x(isd:ied, jsd:jed) = global(isd:ied, jsd:jed) + call mpp_update_domains(x, domain) + + if (any(abs(x(isd:ied, jsd:jed) - global(isd:ied, jsd:jed)) .gt. eps)) then + call mpp_error(FATAL, "Error: A-grid scalar update (single tile) did not produce expected result") + else + print "(A)", "A-grid scalar update (single tile) succeeded on PE " // string(mpp_pe()) + endif + end subroutine single_tile_scalar_test + + subroutine define_corner_mosaic(ni, nj) + type contact + integer :: tile(2) + integer :: is(2) + integer :: ie(2) + integer :: js(2) + integer :: je(2) + end type contact + + integer, intent(in), dimension(3) :: ni, nj + type(domain2d) :: domain + integer, allocatable :: pelist(:) + integer, allocatable :: tile_id(:) + integer :: npes, i + + integer, parameter :: num_contacts = 3 + type(contact) :: contacts(num_contacts) + integer :: global_indices(4,ntiles), layout(2,ntiles) + + npes = dot_product(ni, nj) + + if (mpp_npes().lt.npes) then + call mpp_error(FATAL, "Cannot construct corner mosaic: Insufficient number of PEs") + endif + + if (mpp_pe().gt.npes - 1) then + ! This is an idle PE. + return + endif + + allocate (pelist(npes)) + do i=1,npes + pelist(i) = i-1 + enddo + call mpp_declare_pelist(pelist) + call mpp_set_current_pelist(pelist) + + isg=1 + ieg=n + jsg=1 + jeg=n + + contacts(:) = [ & + contact( & + tile = [1, 2], & + is = [ieg, isg], & + ie = [ieg, isg], & + js = [jsg, jsg], & + je = [jeg, jeg]), & + contact( & + tile = [2, 3], & + is = [isg, isg], & + ie = [ieg, ieg], & + js = [jeg, jsg], & + je = [jeg, jsg]), & + contact( & + tile = [1, 3], & + is = [isg, isg], & + ie = [ieg, isg], & + js = [jeg, jeg], & + je = [jeg, jsg]) & + ] + + do i=1,ntiles + global_indices(:, i) = [isg, ieg, jsg, jeg] + layout(:, i) = [ni(i), nj(i)] + + pe_start(i) = dot_product(ni(1:i-1), nj(1:i-1)) + pe_end(i) = dot_product(ni(1:i), nj(1:i)) - 1 + enddo + + call mpp_define_mosaic(global_indices, layout, domain, ntiles, num_contacts, & + contacts(:)%tile(1), contacts(:)%tile(2), & + contacts(:)%is(1), contacts(:)%ie(1), contacts(:)%js(1), contacts(:)%je(1), & + contacts(:)%is(2), contacts(:)%ie(2), contacts(:)%js(2), contacts(:)%je(2), & + pe_start, pe_end, symmetry=.true., whalo=halo, ehalo=halo, & + shalo=halo, nhalo=halo, name=name) + + call mpp_get_compute_domain(domain, isc, iec, jsc, jec) + call mpp_get_data_domain(domain, isd, ied, jsd, jed) + + tile_id = mpp_get_tile_id(domain) + if (size(tile_id).ne.1) then + call mpp_error(FATAL, "Too many tiles on current PE: Only one expected") + endif + tid = tile_id(1) + + call test_agrid_scalar_update(domain) + call test_cgrid_vector_update(domain) + call test_dgrid_vector_update(domain) + + call mpp_deallocate_domain(domain) + end subroutine define_corner_mosaic + + subroutine agrid_scalar_init(global, x) + real(r4_kind), intent(out) :: global((isg-halo):(ieg+halo), (jsg-halo):(jeg+halo), ntiles) + real(r4_kind), intent(out) :: x(isd:ied, jsd:jed) + integer :: i, j, t, h + + do concurrent (i=lbound(global,1):ubound(global,1), & + j=lbound(global,2):ubound(global,2), & + t=lbound(global,3):ubound(global,3)) + h = 0._r4_kind + if (i.lt.isg .or. i.gt.ieg) h = 999000._r4_kind + if (j.lt.jsg .or. j.gt.jeg) h = 999000._r4_kind + global(i,j,t) = t*100._r4_kind + i*1._r4_kind + j*.01_r4_kind + h + enddo + + x(isd:ied, jsd:jed) = global(isd:ied, jsd:jed, tid) + + ! 1-2 ! + global((n+1):(n+halo), 1:n, 1) = global(1:halo, 1:n, 2) + global((1-halo):0, 1:n, 2) = global((n+1-halo):n, 1:n, 1) + + ! 2-3 ! + global(1:n, (n+1):(n+halo), 2) = global(1:n, 1:halo, 3) + global(1:n, (1-halo):0, 3) = global(1:n, (n+1-halo):n, 2) + + ! 1-3 ! + do concurrent (i=1:n, j=1:halo) + global(i, n+j, 1) = global(j, n-i+1, 3) + global(j-halo, i, 3) = global(n-i+1, n-halo+j, 1) + enddo + end subroutine agrid_scalar_init + + subroutine cgrid_vector_init(gx, gy, x, y) + real(r4_kind) :: gx((isg-halo):(ieg+halo+1), (jsg-halo):(jeg+halo), ntiles), & + gy((isg-halo):(ieg+halo), (jsg-halo):(jeg+halo+1), ntiles) + real(r4_kind) :: x(isd:ied+1, jsd:jed), & + y(isd:ied, jsd:jed+1) + integer :: i, j, t, h, repeat + + do concurrent (i=lbound(gx,1):ubound(gx,1), & + j=lbound(gx,2):ubound(gx,2), & + t=lbound(gx,3):ubound(gx,3)) + h = 0._r4_kind + if (i.lt.isg .or. i.gt.ieg+1) h = 999000._r4_kind + if (j.lt.jsg .or. j.gt.jeg) h = 999000._r4_kind + gx(i,j,t) = t*100._r4_kind + i*1._r4_kind + j*.01_r4_kind + h + enddo + + do concurrent (i=lbound(gy,1):ubound(gy,1), & + j=lbound(gy,2):ubound(gy,2), & + t=lbound(gy,3):ubound(gy,3)) + h = 0._r4_kind + if (i.lt.isg .or. i.gt.ieg) h = 999000._r4_kind + if (j.lt.jsg .or. j.gt.jeg+1) h = 999000._r4_kind + gy(i,j,t) = -1._r4_kind * (t*100._r4_kind + i*1._r4_kind + j*.01_r4_kind + h) + enddo + + ! TODO: Copying the halos and the common rows/columns twice seems to make everything + ! consistent. A single round of copying would probably suffice if the copying + ! operations were more carefully ordered. + do repeat=1,2 + ! 1,2 (x) + gx((1-halo):1, 1:n, 2) = gx((n+1-halo):(n+1), 1:n, 1) + gx((n+1):(n+1+halo), 1:n, 1) = gx(1:(1+halo), 1:n, 2) + + ! 1,2 (y) + gy((1-halo):0, 1:(n+1), 2) = gy((n+1-halo):n, 1:(n+1), 1) + gy((n+1):(n+halo), 1:(n+1), 1) = gy(1:halo, 1:(n+1), 2) + + ! 2,3 (x) + gx(1:(n+1), (1-halo):0, 3) = gx(1:(n+1), (n+1-halo):n, 2) + gx(1:(n+1), (n+1):(n+halo), 2) = gx(1:(n+1), 1:halo, 3) + + ! 2,3 (y) + gy(1:n, (1-halo):1, 3) = gy(1:n, (n+1-halo):(n+1), 2) + gy(1:n, (n+1):(n+1+halo), 2) = gy(1:n, 1:(1+halo), 3) + + ! 3(x), 1(y) + gy(1:n, (n+1):(n+1+halo), 1) = transpose(gx(1:(1+halo), n:1:-1, 3)) + gx((1-halo):1, 1:n, 3) = transpose(gy(n:1:-1, (n+1-halo):(n+1), 1)) + + ! 3(y), 1(x) + gx(1:(n+1), (n+1):(n+halo), 1) = -transpose(gy(1:halo, (n+1):1:-1, 3)) + gy((1-halo):0, 1:(n+1), 3) = -transpose(gx((n+1):1:-1, (n+1-halo):n, 1)) + enddo + + x(isd:ied+1, jsd:jed) = gx(isd:ied+1, jsd:jed, tid) + y(isd:ied, jsd:jed+1) = gy(isd:ied, jsd:jed+1, tid) + end subroutine cgrid_vector_init + + subroutine dgrid_vector_init(gx, gy, x, y) + real(r4_kind) :: gx((isg-halo):(ieg+halo), (jsg-halo):(jeg+halo+1), ntiles), & + gy((isg-halo):(ieg+halo+1), (jsg-halo):(jeg+halo), ntiles) + real(r4_kind) :: x(isd:ied, jsd:jed+1), & + y(isd:ied+1, jsd:jed) + integer :: i, j, t, h, repeat + + do concurrent (i=lbound(gx,1):ubound(gx,1), & + j=lbound(gx,2):ubound(gx,2), & + t=lbound(gx,3):ubound(gx,3)) + h = 0._r4_kind + if (i.lt.isg .or. i.gt.ieg) h = 999000._r4_kind + if (j.lt.jsg .or. j.gt.jeg+1) h = 999000._r4_kind + gx(i,j,t) = t*100._r4_kind + i*1._r4_kind + j*.01_r4_kind + h + enddo + + do concurrent (i=lbound(gy,1):ubound(gy,1), & + j=lbound(gy,2):ubound(gy,2), & + t=lbound(gy,3):ubound(gy,3)) + h = 0._r4_kind + if (i.lt.isg .or. i.gt.ieg+1) h = 999000._r4_kind + if (j.lt.jsg .or. j.gt.jeg) h = 999000._r4_kind + gy(i,j,t) = -1._r4_kind * (t*100._r4_kind + i*1._r4_kind + j*.01_r4_kind + h) + enddo + + ! TODO: Copying the halos and the common rows/columns twice seems to make everything + ! consistent. A single round of copying would probably suffice if the copying + ! operations were more carefully ordered. + do repeat=1,2 + ! 1,2 (x) + gy((1-halo):1, 1:n, 2) = gy((n+1-halo):(n+1), 1:n, 1) + gy((n+1):(n+1+halo), 1:n, 1) = gy(1:(1+halo), 1:n, 2) + + ! 1,2 (y) + gx((1-halo):0, 1:(n+1), 2) = gx((n+1-halo):n, 1:(n+1), 1) + gx((n+1):(n+halo), 1:(n+1), 1) = gx(1:halo, 1:(n+1), 2) + + ! 2,3 (x) + gy(1:(n+1), (1-halo):0, 3) = gy(1:(n+1), (n+1-halo):n, 2) + gy(1:(n+1), (n+1):(n+halo), 2) = gy(1:(n+1), 1:halo, 3) + + ! 2,3 (y) + gx(1:n, (1-halo):1, 3) = gx(1:n, (n+1-halo):(n+1), 2) + gx(1:n, (n+1):(n+1+halo), 2) = gx(1:n, 1:(1+halo), 3) + + ! 3(x), 1(y) + gx(1:n, (n+1):(n+1+halo), 1) = -transpose(gy(1:(1+halo), n:1:-1, 3)) + gy((1-halo):1, 1:n, 3) = -transpose(gx(n:1:-1, (n+1-halo):(n+1), 1)) + + ! 3(y), 1(x) + gy(1:(n+1), (n+1):(n+halo), 1) = transpose(gx(1:halo, (n+1):1:-1, 3)) + gx((1-halo):0, 1:(n+1), 3) = transpose(gy((n+1):1:-1, (n+1-halo):n, 1)) + enddo + + x(isd:ied, jsd:jed+1) = gx(isd:ied, jsd:jed+1, tid) + y(isd:ied+1, jsd:jed) = gy(isd:ied+1, jsd:jed, tid) + end subroutine dgrid_vector_init + + subroutine test_agrid_scalar_update(domain) + type(domain2d), intent(inout) :: domain + real(r4_kind) :: global((isg-halo):(ieg+halo), (jsg-halo):(jeg+halo), ntiles) + real(r4_kind) :: x(isd:ied, jsd:jed) + + call agrid_scalar_init(global, x) + call mpp_update_domains(x, domain) + + if (any(abs(x(isd:ied, jsd:jed) - global(isd:ied, jsd:jed, tid)) .gt. eps)) then + call mpp_error(FATAL, "Error: A-grid scalar update did not produce expected result") + else + print "(A)", "A-grid scalar update succeeded on PE " // string(mpp_pe()) + endif + end subroutine test_agrid_scalar_update + + subroutine test_cgrid_vector_update(domain) + type(domain2d), intent(inout) :: domain + real(r4_kind) :: gx((isg-halo):(ieg+halo+1), (jsg-halo):(jeg+halo), ntiles), & + gy((isg-halo):(ieg+halo), (jsg-halo):(jeg+halo+1), ntiles) + real(r4_kind) :: x(isd:ied+1, jsd:jed), & + y(isd:ied, jsd:jed+1) + + call cgrid_vector_init(gx, gy, x, y) + call mpp_update_domains(x, y, domain, gridtype=CGRID_NE) + + if (any(abs(x(isd:ied+1, jsd:jed) - gx(isd:ied+1, jsd:jed, tid)) .gt. eps) .or. & + any(abs(y(isd:ied, jsd:jed+1) - gy(isd:ied, jsd:jed+1, tid)) .gt. eps)) then + call mpp_error(FATAL, "Error: C-grid vector update did not produce expected result") + else + print "(A)", "C-grid vector update succeeded on PE " // string(mpp_pe()) + endif + end subroutine test_cgrid_vector_update + + subroutine test_dgrid_vector_update(domain) + type(domain2d), intent(inout) :: domain + real(r4_kind) :: gx((isg-halo):(ieg+halo), (jsg-halo):(jeg+halo+1), ntiles), & + gy((isg-halo):(ieg+halo+1), (jsg-halo):(jeg+halo), ntiles) + real(r4_kind) :: x(isd:ied, jsd:jed+1), & + y(isd:ied+1, jsd:jed) + + call dgrid_vector_init(gx, gy, x, y) + call mpp_update_domains(x, y, domain, gridtype=DGRID_NE) + + if (any(abs(x(isd:ied, jsd:jed+1) - gx(isd:ied, jsd:jed+1, tid)) .gt. eps) .or. & + any(abs(y(isd:ied+1, jsd:jed) - gy(isd:ied+1, jsd:jed, tid)) .gt. eps)) then + call mpp_error(FATAL, "Error: D-grid vector update did not produce expected result") + else + print "(A)", "D-grid vector update succeeded on PE " // string(mpp_pe()) + endif + end subroutine test_dgrid_vector_update + + subroutine print_matrix(x) + real(r4_kind), dimension(:,:), intent(in) :: x + integer :: i,j + integer :: iostat + + do j=size(x,2),1,-1 + do i=1,size(x,1) + write(*,'(" ",F12.4)', advance="no", iostat=iostat) x(i,j) + enddo + write(*, "(A)") "" + enddo + end subroutine +end program test_corner_mosaic diff --git a/test_fms/mpp/test_mpp_update_domains.sh b/test_fms/mpp/test_mpp_update_domains.sh index 50debf96e2..68cb06c17e 100755 --- a/test_fms/mpp/test_mpp_update_domains.sh +++ b/test_fms/mpp/test_mpp_update_domains.sh @@ -34,4 +34,7 @@ test_expect_success "domain updates with single PE" ' test_expect_success "domain updates with multiple PEs" ' mpirun -n 2 ./test_mpp_update_domains ' +test_expect_success "corner mosaic with (32 PEs)" ' + mpirun -n 32 ./test_corner_mosaic +' test_done