diff --git a/axis_utils/include/axis_utils2.inc b/axis_utils/include/axis_utils2.inc index a1055dfd41..bac9251e07 100644 --- a/axis_utils/include/axis_utils2.inc +++ b/axis_utils/include/axis_utils2.inc @@ -46,6 +46,9 @@ !! indicating to reproduce the mpp_io bug where !! the null characters were not removed after reading a string attribute integer, parameter :: lkind = FMS_AU_KIND_ + integer :: edge_index(2) !< Index to use when reading the edges from the file + !! (/1, 2/) if the axis data is monotonically increasing + !! (/2, 1/) if the axis data is monotonically decreasing ndims = get_variable_num_dimensions(fileobj, name) allocate(dim_sizes(ndims)) @@ -108,8 +111,10 @@ allocate(r2d(dim_sizes(1), dim_sizes(2))) call read_data(fileobj, buffer, r2d) - edge_data(1:dim_sizes(2)) = r2d(1,:) - edge_data(dim_sizes(2)+1) = r2d(2,dim_sizes(2)) + edge_index = (/1, 2/) + if (r2d(1,1) .gt. r2d(1,2)) edge_index = (/2, 1 /) + edge_data(1:dim_sizes(2)) = r2d(edge_index(1),:) + edge_data(dim_sizes(2)+1) = r2d(edge_index(2),dim_sizes(2)) deallocate(r2d) endif deallocate(dim_sizes) @@ -265,7 +270,7 @@ !! inputs: !! !! rval = arbitrary data...same units as elements in "array" - !! array = array of data points (must be monotonically increasing) + !! array = array of data points (must be monotonically) !! ia = dimension of "array" !! !! output: @@ -293,48 +298,71 @@ !! z(k1) would be the nearest data point to 12.5 and z(k2) would !! be the nearest data point to 0.0 !! @return integer nearest_index - - - function NEAREST_INDEX_(rval, array) + real(kind=FMS_AU_KIND_), intent(in) :: rval !< arbitrary data...same units as elements in "array" + real(kind=FMS_AU_KIND_), intent(in), dimension(:) :: array !< array of data points (must be monotonic) integer :: NEAREST_INDEX_ integer :: ia !< dimension of "array" - integer :: i, ii, iunit - real(kind=FMS_AU_KIND_) :: rval !< arbitrary data...same units as elements in "array" - real(kind=FMS_AU_KIND_), dimension(:) :: array !< array of data points (must be monotonically increasing) - logical :: keep_going - - ia = size(array(:)) - - do i = 2, ia - if (array(i) < array(i-1)) then - iunit = stdout() - write (iunit,*) '=> Error: "nearest_index" array must be monotonically increasing' & - & // 'when searching for nearest value to ', rval - write (iunit,*) ' array(i) < array(i-1) for i=',i - write (iunit,*) ' array(i) for i=1..ia follows:' - do ii = 1, ia - write (iunit,*) 'i=',ii, ' array(i)=', array(ii) - enddo - call mpp_error(FATAL,' "nearest_index" array must be monotonically increasing.') - endif - enddo + integer :: i !< For looping through "array" + + ia = SIZE(array(:)) + + ! check if array is monotonous + DO i = 2, ia-1 + IF( (array(i-1)array(i+1)) .OR. (array(i-1)>array(i).AND.array(i)array NOT monotonously ordered + CALL mpp_error(FATAL, 'axis_utils2::nearest_index array is NOT monotonously ordered') + END IF + END DO + + if (array(1) < array(ia)) then + !< increasing array + + !< Check if the rval is outside the range of the array + if (rval < array(1)) then + NEAREST_INDEX_ = 1 + return + elseif (rval > array(ia)) then + NEAREST_INDEX_ = ia + return + endif - if (rval < array(1) .or. rval > array(ia)) then - if (rval < array(1)) NEAREST_INDEX_ = 1 - if (rval > array(ia)) NEAREST_INDEX_ = ia + DO i = 1, ia-1 + IF ( (array(i)<=rval).AND.(array(i+1)>= rval) ) THEN + IF( rval - array(i) <= array(i+1) - rval ) THEN + NEAREST_INDEX_ = i + return + ELSE + NEAREST_INDEX_ = i+1 + return + ENDIF + EXIT + END IF + END DO else - i = 1 - keep_going = .true. - do while (i <= ia .and. keep_going) - i = i+1 - if (rval <= array(i)) then - NEAREST_INDEX_ = i - if (array(i)-rval > rval-array(i-1)) NEAREST_INDEX_ = i-1 - keep_going = .false. - endif - enddo + !< Decreasing Array + + !< Check if the rval is outside the range of the array + if (rval < array(ia)) then + NEAREST_INDEX_ = ia + return + elseif (rval > array(1)) then + NEAREST_INDEX_ = 1 + return + endif + + DO i = 1, ia-1 + IF ( (array(i)>=rval).AND.(array(i+1)<= rval) ) THEN + IF ( array(i)-rval <= rval-array(i+1) ) THEN + NEAREST_INDEX_ = i + return + ELSE + NEAREST_INDEX_ = i+1 + return + END IF + END IF + END DO endif end function NEAREST_INDEX_ diff --git a/data_override/include/data_override.inc b/data_override/include/data_override.inc index aa2df0d0c2..06d1d9379b 100644 --- a/data_override/include/data_override.inc +++ b/data_override/include/data_override.inc @@ -1031,14 +1031,14 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d ie_src = axis_sizes(1) js_src = 1 je_src = axis_sizes(2) - do j = 1, axis_sizes(2)+1 - if( override_array(curr_position)%lat_in(j) > lat_min ) exit - js_src = j - enddo - do j = 1, axis_sizes(2)+1 - je_src = j - if( override_array(curr_position)%lat_in(j) > lat_max ) exit - enddo + ! do j = 1, axis_sizes(2)+1 + ! if( override_array(curr_position)%lat_in(j) > lat_min ) exit + ! js_src = j + ! enddo + ! do j = 1, axis_sizes(2)+1 + ! je_src = j + ! if( override_array(curr_position)%lat_in(j) > lat_max ) exit + ! enddo !--- bicubic interpolation need one extra point in each direction. Also add !--- one more point for because lat_in is in the corner but the interpolation diff --git a/horiz_interp/horiz_interp_bilinear.F90 b/horiz_interp/horiz_interp_bilinear.F90 index 64abf15263..318d2c039b 100644 --- a/horiz_interp/horiz_interp_bilinear.F90 +++ b/horiz_interp/horiz_interp_bilinear.F90 @@ -34,6 +34,7 @@ module horiz_interp_bilinear_mod use constants_mod, only: PI use horiz_interp_type_mod, only: horiz_interp_type, stats use platform_mod, only: r4_kind, r8_kind + use axis_utils2_mod, only: nearest_index implicit none private @@ -63,12 +64,6 @@ module horiz_interp_bilinear_mod integer, parameter :: DUMMY = -999 !! Private helper routines, interfaces for mixed real precision support - - interface indp - module procedure indp_r4 - module procedure indp_r8 - end interface - interface intersect module procedure intersect_r4 module procedure intersect_r8 diff --git a/horiz_interp/include/horiz_interp_bilinear.inc b/horiz_interp/include/horiz_interp_bilinear.inc index 8f81e86f79..9e352d9c31 100644 --- a/horiz_interp/include/horiz_interp_bilinear.inc +++ b/horiz_interp/include/horiz_interp_bilinear.inc @@ -34,12 +34,20 @@ real(FMS_HI_KIND_) :: wtw, wte, wts, wtn, lon, lat, tpi, hpi real(FMS_HI_KIND_) :: glt_min, glt_max, gln_min, gln_max, min_lon, max_lon integer,parameter :: kindl = FMS_HI_KIND_ + logical :: decreasing_lat !< .True. if latitude is monotically decreasing + logical :: decreasing_lon !< .True. if longitude is monotically decreasing warns = 0 if(present(verbose)) warns = verbose src_is_modulo = .true. if (present(src_modulo)) src_is_modulo = src_modulo + decreasing_lat = .false. + if (lat_in(1) > lat_in(2)) decreasing_lat = .true. + + decreasing_lon = .false. + if (lon_in(1) > lon_in(2)) decreasing_lon = .true. + hpi = 0.5_kindl * real(pi, FMS_HI_KIND_) tpi = 4.0_kindl * hpi glt_min = hpi @@ -95,11 +103,23 @@ glt_min = min(lat,glt_min); glt_max = max(lat,glt_max) gln_min = min(lon,gln_min); gln_max = max(lon,gln_max) - is = indp(lon, lon_in ) - if( lon_in(is) .gt. lon ) is = max(is-1,1) + is = nearest_index(lon, lon_in ) + if (decreasing_lon) then + ! Lon_in is increasing + ! This is so that is is the lower bound. + ! For example, if the array is [50 40 30 20 10] and lon is 11, `is` is going to be 5 from `nearest_index` + ! but it needs to be 4 so that it can use the data at lon_in(4) and lon_in(5) + if( lon_in(is) .lt. lon ) is = max(is-1,1) + else + ! Lon_in is increasing + ! This is so that is is the lower bound. + ! For example, if the array is [10 20 30 40 50] and lon is 49, `is` is going to be 5 from `nearest_index` + ! but it needs to be 4 so that it can use the data at lon_in(4) and lon_in(5) + if( lon_in(is) .gt. lon ) is = max(is-1,1) + endif if( lon_in(is) .eq. lon .and. is .eq. nlon_in) is = max(is - 1,1) ie = min(is+1,nlon_in) - if(lon_in(is) .ne. lon_in(ie) .and. lon_in(is) .le. lon) then + if(lon_in(is) .ne. lon_in(ie) .and. (decreasing_lon .or. lon_in(is) .le. lon)) then wtw = ( lon_in(ie) - lon) / (lon_in(ie) - lon_in(is) ) else ! east or west of the last data value. this could be because a @@ -115,13 +135,18 @@ endif wte = 1.0_kindl - wtw - js = indp(lat, lat_in ) - - if( lat_in(js) .gt. lat ) js = max(js - 1, 1) + js = nearest_index(lat, lat_in ) + if (decreasing_lat) then + ! Lat_in is decreasing + if( lat_in(js) .lt. lat ) js = max(js - 1, 1) + else + ! Lat_in is increasing + if( lat_in(js) .gt. lat ) js = max(js - 1, 1) + endif if( lat_in(js) .eq. lat .and. js .eq. nlat_in) js = max(js - 1, 1) je = min(js + 1, nlat_in) - if ( lat_in(js) .ne. lat_in(je) .and. lat_in(js) .le. lat) then + if ( lat_in(js) .ne. lat_in(je) .and. (decreasing_lat .or. lat_in(js) .le. lat)) then wts = ( lat_in(je) - lat )/(lat_in(je)-lat_in(js)) else ! north or south of the last data value. this could be because a @@ -1180,51 +1205,4 @@ return end subroutine - - - !####################################################################### - !> @returns index of nearest data point to "value" - !! if "value" is outside the domain of "array" then INDP_ = 1 - !! or "ia" depending on whether array(1) or array(ia) is - !! closest to "value" - function INDP_ (rval, array) - integer :: INDP_ !< index of nearest data point within "array" - !! corresponding to "value". - real(FMS_HI_KIND_), dimension(:), intent(in) :: array !< array of data points (must be monotonically increasing) - real(FMS_HI_KIND_), intent(in) :: rval !< arbitrary data, same units as elements in 'array' - - !======================================================================= - - integer i, ia, iunit - logical keep_going - ! - ia = size(array(:)) - do i=2,ia - if (array(i) .lt. array(i-1)) then - iunit = stdout() - write (iunit,*) & - ' => Error: array must be monotonically increasing in "INDP_"' , & - ' when searching for nearest element to value=',rval - write (iunit,*) ' array(i) < array(i-1) for i=',i - write (iunit,*) ' array(i) for i=1..ia follows:' - call mpp_error() - endif - enddo - if (rval .lt. array(1) .or. rval .gt. array(ia)) then - if (rval .lt. array(1)) INDP_ = 1 - if (rval .gt. array(ia)) INDP_ = ia - else - i=1 - keep_going = .true. - do while (i .le. ia .and. keep_going) - i = i+1 - if (rval .le. array(i)) then - INDP_ = i - if (array(i)-rval .gt. rval-array(i-1)) INDP_ = i-1 - keep_going = .false. - endif - enddo - endif - return - end function INDP_ !> @} diff --git a/horiz_interp/include/horiz_interp_bilinear_r4.fh b/horiz_interp/include/horiz_interp_bilinear_r4.fh index 07cf041ccc..8880914e43 100644 --- a/horiz_interp/include/horiz_interp_bilinear_r4.fh +++ b/horiz_interp/include/horiz_interp_bilinear_r4.fh @@ -45,8 +45,5 @@ #undef INTERSECT_ #define INTERSECT_ intersect_r4 -#undef INDP_ -#define INDP_ indp_r4 - #include "horiz_interp_bilinear.inc" !> @} diff --git a/horiz_interp/include/horiz_interp_bilinear_r8.fh b/horiz_interp/include/horiz_interp_bilinear_r8.fh index 526fecfb69..37a2e6920b 100644 --- a/horiz_interp/include/horiz_interp_bilinear_r8.fh +++ b/horiz_interp/include/horiz_interp_bilinear_r8.fh @@ -45,8 +45,5 @@ #undef INTERSECT_ #define INTERSECT_ intersect_r8 -#undef INDP_ -#define INDP_ indp_r8 - #include "horiz_interp_bilinear.inc" !> @} diff --git a/test_fms/axis_utils/test_axis_utils.F90 b/test_fms/axis_utils/test_axis_utils.F90 index aac74de010..c86d2ee4c4 100644 --- a/test_fms/axis_utils/test_axis_utils.F90 +++ b/test_fms/axis_utils/test_axis_utils.F90 @@ -76,17 +76,25 @@ program test_axis_utils print "(A)", "Testing frac_index (FAILURE)" call test_frac_index_fail - case ('--nearest-index') - print "(A)", "Testing nearest_index" - call test_nearest_index + case ('--nearest-index-increasing') + print "(A)", "Testing nearest_index with a monotonically increasing array" + call test_nearest_index(.true.) + + case ('--nearest-index-decreasing') + print "(A)", "Testing nearest_index with a monotonically decreasing array" + call test_nearest_index(.false.) case ('--nearest-index-fail') print "(A)", "Testing nearest_index (FAILURE)" call test_nearest_index_fail - case ('--axis-edges') - print "(A)", "Testing axis_edges" - call test_axis_edges + case ('--axis-edges-increasing') + print "(A)", "Testing axis_edges-increasing" + call test_axis_edges(.true.) + + case ('--axis-edges-decreasing') + print "(A)", "Testing axis_edges-decreasing" + call test_axis_edges(.false.) case ('--tranlon') print "(A)", "Testing tranlon" @@ -385,28 +393,36 @@ subroutine test_frac_index_fail ret_test = frac_index(1.5_k, values) end subroutine -subroutine test_nearest_index +subroutine test_nearest_index(increasing_array) + logical, intent(in) :: increasing_array !< .True. if test using an increasing array real(k) :: arr(5) + integer :: ans(12) - arr = [5._k, 12._k, 20._k, 40._k, 100._k] + if (increasing_array) then + arr = [5._k, 12._k, 20._k, 40._k, 100._k] + ans=(/1, 5, 1, 2, 3, 4, 5, 1, 2, 2, 3, 3/) + else + arr = [100._k, 40._k, 20._k, 12._k, 5._k] + ans=(/5, 1, 5, 4, 3, 2, 1, 5, 4, 4, 3, 3/) + endif ! Test values beyond array boundaries - call nearest_index_assert(4._k, arr, 1) - call nearest_index_assert(1000._k, arr, size(arr)) + call nearest_index_assert(4._k, arr, ans(1)) + call nearest_index_assert(1000._k, arr, ans(2)) ! Test values actually in the array - call nearest_index_assert(5._k, arr, 1) - call nearest_index_assert(12._k, arr, 2) - call nearest_index_assert(20._k, arr, 3) - call nearest_index_assert(40._k, arr, 4) - call nearest_index_assert(100._k, arr, 5) + call nearest_index_assert(5._k, arr, ans(3)) + call nearest_index_assert(12._k, arr, ans(4)) + call nearest_index_assert(20._k, arr, ans(5)) + call nearest_index_assert(40._k, arr, ans(6)) + call nearest_index_assert(100._k, arr, ans(7)) ! Test the intervals between array values - call nearest_index_assert(6._k, arr, 1) - call nearest_index_assert(11._k, arr, 2) - call nearest_index_assert(15._k, arr, 2) - call nearest_index_assert(18._k, arr, 3) - call nearest_index_assert(29._k, arr, 3) + call nearest_index_assert(6._k, arr, ans(8)) + call nearest_index_assert(11._k, arr, ans(9)) + call nearest_index_assert(15._k, arr, ans(10)) + call nearest_index_assert(18._k, arr, ans(11)) + call nearest_index_assert(29._k, arr, ans(12)) end subroutine subroutine nearest_index_assert(val, arr, ret_expected) @@ -433,24 +449,43 @@ subroutine test_nearest_index_fail ret_test = nearest_index(5._k, arr) end subroutine -subroutine test_axis_edges +subroutine test_axis_edges(increasing_array) + logical, intent(in) :: increasing_array !< .True. if test using an increasing array real(k) :: data_in_var(10) real(k) :: data_in_var_edges(2,10) real(k) :: data_in_answers(11) type(FmsNetcdfFile_t) :: fileobj real(k) :: answers(11) - integer :: i + integer :: count + integer :: count_factor + integer :: factor + integer :: index !< For looping through the data + + if (increasing_array) then + count = 0 + factor = 1 + count_factor = -1 + else + count = 11 + factor = -1 + count_factor = 0 + endif - do i=1,10 - data_in_var(i) = real(i, k) - 0.5_k + do index=1,10 + count = count + factor + data_in_var(index) = real(count, k) - 0.5_k - data_in_var_edges(1,i) = real(i-1, k) - data_in_var_edges(2,i) = real(i, k) + data_in_var_edges(1,index) = real(count-1, k) + data_in_var_edges(2,index) = real(count, k) - data_in_answers(i) = real(i-1, k) + data_in_answers(index) = real(count + count_factor, k) enddo - data_in_answers(11) = 10._k + if (increasing_array) then + data_in_answers(11) = real(count, k) + else + data_in_answers(11) = real(count + factor, k) + endif call open_netcdf_w(fileobj) diff --git a/test_fms/axis_utils/test_axis_utils2.sh b/test_fms/axis_utils/test_axis_utils2.sh index 1288822481..8b49c6bdd5 100755 --- a/test_fms/axis_utils/test_axis_utils2.sh +++ b/test_fms/axis_utils/test_axis_utils2.sh @@ -27,11 +27,11 @@ # Prepare the directory to run the tests. touch input.nml -TESTS_SUCCESS='--get-axis-modulo --get-axis-modulo-times --get-axis-cart --lon-in-range --frac-index --nearest-index --axis-edges --tranlon --interp-1d-1d --interp-1d-2d --interp-1d-3d' +TESTS_SUCCESS='--get-axis-modulo --get-axis-modulo-times --get-axis-cart --lon-in-range --frac-index --nearest-index-increasing --nearest-index-decreasing --axis-edges-increasing --axis_edges-decreasing --tranlon --interp-1d-1d --interp-1d-2d --interp-1d-3d' TESTS_FAIL='--frac-index-fail --nearest-index-fail' # TODO: Enable these tests after tranlon's memory corruption bug is fixed. -SKIP_TESTS="test_axis_utils2.15 test_axis_utils2.16" +SKIP_TESTS="test_axis_utils2.19 test_axis_utils2.20" # Run the tests diff --git a/test_fms/data_override/Makefile.am b/test_fms/data_override/Makefile.am index 567f787434..d71512a364 100644 --- a/test_fms/data_override/Makefile.am +++ b/test_fms/data_override/Makefile.am @@ -71,10 +71,11 @@ TESTS_ENVIRONMENT= test_input_path="@TEST_INPUT_PATH@" \ parser_skip=${skipflag} # Run the test program. -TESTS = test_data_override2.sh test_data_override_init.sh + +TESTS = test_data_override2.sh test_data_override_init.sh test_data_override2_mono.sh test_data_override2_ongrid.sh # Include these files with the distribution. -EXTRA_DIST = test_data_override2.sh test_data_override_init.sh +EXTRA_DIST = test_data_override2.sh test_data_override_init.sh test_data_override2_mono.sh test_data_override2_ongrid.sh # Clean up CLEANFILES = input.nml *.nc* *.out diag_table data_table data_table.yaml INPUT/* *.dpi *.spi *.dyn *.spl *-files/* diff --git a/test_fms/data_override/test_data_override2.sh b/test_fms/data_override/test_data_override2.sh index ffa42c05c6..48940debf7 100755 --- a/test_fms/data_override/test_data_override2.sh +++ b/test_fms/data_override/test_data_override2.sh @@ -27,62 +27,14 @@ output_dir rm -rf data_table data_table.yaml input.nml input_base.nml -if [ ! -z $parser_skip ]; then - cat <<_EOF > input_base.nml -&data_override_nml -use_data_table_yaml=.False. -/ - -&test_data_override_ongrid_nml - nhalox=halo_size - nhaloy=halo_size -/ -_EOF - printf '"OCN", "runoff", "runoff", "./INPUT/runoff.daitren.clim.1440x1080.v20180328.nc", "none" , 1.0' | cat > data_table -else -cat <<_EOF > input_base.nml -&data_override_nml -use_data_table_yaml=.True. -/ - -&test_data_override_ongrid_nml - nhalox=halo_size - nhaloy=halo_size -/ -_EOF -cat <<_EOF > data_table.yaml -data_table: - - gridname : OCN - fieldname_code : runoff - fieldname_file : runoff - file_name : INPUT/runoff.daitren.clim.1440x1080.v20180328.nc - interpol_method : none - factor : 1.0 -_EOF -fi - -[ ! -d "INPUT" ] && mkdir -p "INPUT" -for KIND in r4 r8 -do -sed 's/halo_size/2/g' input_base.nml > input.nml -test_expect_success "data_override on grid with 2 halos in x and y (${KIND})" ' - mpirun -n 6 ../test_data_override_ongrid_${KIND} -' -sed 's/halo_size/0/g' input_base.nml > input.nml -test_expect_success "data_override on grid with 2 halos in x and y (${KIND})" ' - mpirun -n 6 ../test_data_override_ongrid_${KIND} -' - -test_expect_success "data_override get_grid_v1 (${KIND})" ' - mpirun -n 1 ../test_get_grid_v1_${KIND} -' -done - for KIND in r4 r8 do # Run tests with input if enabled # skips if built with yaml parser(tests older behavior) if test ! -z "$test_input_path" && test ! -z "$parser_skip" ; then + cat <<_EOF > input.nml +_EOF + cp -r $test_input_path/data_override/INPUT . cat <<_EOF > diag_table test_data_override diff --git a/test_fms/data_override/test_data_override2_mono.sh b/test_fms/data_override/test_data_override2_mono.sh new file mode 100755 index 0000000000..cf47a152f9 --- /dev/null +++ b/test_fms/data_override/test_data_override2_mono.sh @@ -0,0 +1,88 @@ +#!/bin/sh + +#*********************************************************************** +#* GNU Lesser General Public License +#* +#* This file is part of the GFDL Flexible Modeling System (FMS). +#* +#* FMS is free software: you can redistribute it and/or modify it under +#* the terms of the GNU Lesser General Public License as published by +#* the Free Software Foundation, either version 3 of the License, or (at +#* your option) any later version. +#* +#* FMS is distributed in the hope that it will be useful, but WITHOUT +#* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +#* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +#* for more details. +#* +#* You should have received a copy of the GNU Lesser General Public +#* License along with FMS. If not, see . +#*********************************************************************** +# +# Copyright (c) 2019-2021 Ed Hartnett, Uriel Ramirez, Seth Underwood + +# Set common test settings. +. ../test-lib.sh + +output_dir +[ ! -d "INPUT" ] && mkdir -p "INPUT" + +cat <<_EOF > input.nml +&test_data_override_ongrid_nml + test_case = 2 +/ +_EOF + +cat <<_EOF > data_table +"OCN", "runoff_increasing", "runoff", "./INPUT/bilinear_increasing.nc", "bilinear" , 1.0 +"OCN", "runoff_decreasing", "runoff", "./INPUT/bilinear_decreasing.nc", "bilinear" , 1.0 +_EOF + +for KIND in r4 r8 +do + rm -rf INPUT/* + test_expect_success "test_data_override with monotonically increasing and decreasing data sets (${KIND})" ' + mpirun -n 6 ../test_data_override_ongrid_${KIND} + ' +done + +rm -rf data_table + +cat <<_EOF > input.nml +&test_data_override_ongrid_nml + test_case = 2 +/ +&data_override_nml + use_data_table_yaml = .True. +/ +_EOF + +cat <<_EOF > data_table.yaml +data_table: +- gridname: OCN + fieldname_code: runoff_increasing + fieldname_file: runoff + file_name: ./INPUT/bilinear_increasing.nc + interpol_method: bilinear + factor: 1.0 +- gridname: OCN + fieldname_code: runoff_decreasing + fieldname_file: runoff + file_name: ./INPUT/bilinear_decreasing.nc + interpol_method: bilinear + factor: 1.0 +_EOF + +#Repeat the test with yaml if needed +if [ -z $parser_skip ]; then + for KIND in r4 r8 + do + rm -rf INPUT/* + test_expect_success "test_data_override with monotonically increasing and decreasing data sets -yaml (${KIND})" ' + mpirun -n 6 ../test_data_override_ongrid_${KIND} + ' + done +fi + +rm -rf INPUT +test_done diff --git a/test_fms/data_override/test_data_override2_ongrid.sh b/test_fms/data_override/test_data_override2_ongrid.sh new file mode 100755 index 0000000000..2e1d7a1b03 --- /dev/null +++ b/test_fms/data_override/test_data_override2_ongrid.sh @@ -0,0 +1,86 @@ +#!/bin/sh + +#*********************************************************************** +#* GNU Lesser General Public License +#* +#* This file is part of the GFDL Flexible Modeling System (FMS). +#* +#* FMS is free software: you can redistribute it and/or modify it under +#* the terms of the GNU Lesser General Public License as published by +#* the Free Software Foundation, either version 3 of the License, or (at +#* your option) any later version. +#* +#* FMS is distributed in the hope that it will be useful, but WITHOUT +#* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +#* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +#* for more details. +#* +#* You should have received a copy of the GNU Lesser General Public +#* License along with FMS. If not, see . +#*********************************************************************** +# +# Copyright (c) 2019-2021 Ed Hartnett, Uriel Ramirez, Seth Underwood + +# Set common test settings. +. ../test-lib.sh + +output_dir +rm -rf data_table data_table.yaml input.nml input_base.nml + +if [ ! -z $parser_skip ]; then + cat <<_EOF > input_base.nml +&data_override_nml +use_data_table_yaml=.False. +/ + +&test_data_override_ongrid_nml + nhalox=halo_size + nhaloy=halo_size +/ +_EOF + printf '"OCN", "runoff", "runoff", "./INPUT/runoff.daitren.clim.1440x1080.v20180328.nc", "none" , 1.0' | cat > data_table +else +cat <<_EOF > input_base.nml +&data_override_nml +use_data_table_yaml=.True. +/ + +&test_data_override_ongrid_nml + nhalox=halo_size + nhaloy=halo_size +/ +_EOF +cat <<_EOF > data_table.yaml +data_table: + - gridname : OCN + fieldname_code : runoff + fieldname_file : runoff + file_name : INPUT/runoff.daitren.clim.1440x1080.v20180328.nc + interpol_method : none + factor : 1.0 +_EOF +fi + +[ ! -d "INPUT" ] && mkdir -p "INPUT" +for KIND in r4 r8 +do +rm -rf INPUT/* +sed 's/halo_size/2/g' input_base.nml > input.nml +test_expect_success "data_override on grid with 2 halos in x and y (${KIND})" ' + mpirun -n 6 ../test_data_override_ongrid_${KIND} +' + +rm -rf INPUT/* +sed 's/halo_size/0/g' input_base.nml > input.nml +test_expect_success "data_override on grid with 0 halos in x and y (${KIND})" ' + mpirun -n 6 ../test_data_override_ongrid_${KIND} +' + +rm -rf INPUT/* +test_expect_success "data_override get_grid_v1 (${KIND})" ' + mpirun -n 1 ../test_get_grid_v1_${KIND} +' +done +rm -rf INPUT *.nc # remove any leftover files to reduce size + +test_done \ No newline at end of file diff --git a/test_fms/data_override/test_data_override_ongrid.F90 b/test_fms/data_override/test_data_override_ongrid.F90 index 3f031547fa..ac97f03a11 100644 --- a/test_fms/data_override/test_data_override_ongrid.F90 +++ b/test_fms/data_override/test_data_override_ongrid.F90 @@ -20,19 +20,20 @@ program test_data_override_ongrid !> @brief This programs tests data_override ability to override data for an -!! on grid case +!! on grid case and when using bilinear interpolation use platform_mod, only: r4_kind, r8_kind use mpp_domains_mod, only: mpp_define_domains, mpp_define_io_domain, mpp_get_data_domain, & mpp_domains_set_stack_size, mpp_get_compute_domain, domain2d use mpp_mod, only: mpp_init, mpp_exit, mpp_pe, mpp_root_pe, mpp_error, FATAL, & - input_nml_file, mpp_sync + input_nml_file, mpp_sync, NOTE use data_override_mod, only: data_override_init, data_override -use fms2_io_mod, only: fms2_io_init +use fms2_io_mod use time_manager_mod, only: set_calendar_type, time_type, set_date, NOLEAP use netcdf, only: nf90_create, nf90_def_dim, nf90_def_var, nf90_enddef, nf90_put_var, & nf90_close, nf90_put_att, nf90_clobber, nf90_64bit_offset, nf90_char, & nf90_double, nf90_unlimited +use fms_mod, only: string implicit none @@ -41,23 +42,17 @@ program test_data_override_ongrid integer :: nlon !< Number of points in x axis integer :: nlat !< Number of points in y axis type(domain2d) :: Domain !< Domain with mask table -real(lkind), allocatable, dimension(:,:) :: runoff !< Data to be written integer :: is !< Starting x index integer :: ie !< Ending x index integer :: js !< Starting y index integer :: je !< Ending y index -type(time_type) :: Time !< Time -integer :: i !< Helper indices -integer :: ncid !< Netcdf file id -integer :: err !< Error Code -integer :: dim1d, dim2d, dim3d, dim4d !< Dimension ids -integer :: varid, varid2, varid3, varid4 !< Variable ids -real(lkind), allocatable, dimension(:,:,:) :: runoff_in !< Data to be written to file -real(lkind) :: expected_result !< Expected result from data_override integer :: nhalox=2, nhaloy=2 integer :: io_status +integer, parameter :: ongrid = 1 +integer, parameter :: bilinear = 2 +integer :: test_case = ongrid -namelist / test_data_override_ongrid_nml / nhalox, nhaloy +namelist / test_data_override_ongrid_nml / nhalox, nhaloy, test_case call mpp_init call fms2_io_init @@ -65,64 +60,7 @@ program test_data_override_ongrid read (input_nml_file, test_data_override_ongrid_nml, iostat=io_status) if (io_status > 0) call mpp_error(FATAL,'=>test_data_override_ongrid: Error reading input.nml') -!< Create some files needed by data_override! -if (mpp_pe() .eq. mpp_root_pe()) then - allocate(runoff_in(1440, 1080, 10)) - do i = 1, 10 - runoff_in(:,:,i) = real(i, lkind) - enddo - err = nf90_create('INPUT/grid_spec.nc', ior(nf90_clobber, nf90_64bit_offset), ncid) - err = nf90_def_dim(ncid, 'str', 255, dim1d) - err = nf90_def_var(ncid, 'ocn_mosaic_file', nf90_char, (/dim1d/), varid) - err = nf90_enddef(ncid) - err = nf90_put_var(ncid, varid, "ocean_mosaic.nc") - err = nf90_close(ncid) - - err = nf90_create('INPUT/ocean_mosaic.nc', ior(nf90_clobber, nf90_64bit_offset), ncid) - err = nf90_def_dim(ncid, 'str', 255, dim1d) - err = nf90_def_dim(ncid, 'ntiles', 1, dim2d) - err = nf90_def_var(ncid, 'gridfiles', nf90_char, (/dim1d, dim2d/), varid) - err = nf90_enddef(ncid) - err = nf90_put_var(ncid, varid, "ocean_hgrid.nc") - err = nf90_close(ncid) - - err = nf90_create('INPUT/ocean_hgrid.nc', ior(nf90_clobber, nf90_64bit_offset), ncid) - err = nf90_def_dim(ncid, 'nx', 2880, dim1d) - err = nf90_def_dim(ncid, 'ny', 2160, dim2d) - err = nf90_def_dim(ncid, 'nxp', 2881, dim3d) - err = nf90_def_dim(ncid, 'nyp', 2161, dim4d) - err = nf90_def_var(ncid, 'x', nf90_double, (/dim3d, dim4d/), varid) - err = nf90_def_var(ncid, 'y', nf90_double, (/dim3d, dim4d/), varid2) - err = nf90_def_var(ncid, 'area', nf90_double, (/dim1d, dim2d/), varid3) - err = nf90_enddef(ncid) - err = nf90_close(ncid) - - err = nf90_create('INPUT/runoff.daitren.clim.1440x1080.v20180328.nc', ior(nf90_clobber, nf90_64bit_offset), ncid) - err = nf90_def_dim(ncid, 'i', 1440, dim1d) - err = nf90_def_dim(ncid, 'j', 1080, dim2d) - err = nf90_def_dim(ncid, 'time', nf90_unlimited, dim3d) - - err = nf90_def_var(ncid, 'i', nf90_double, (/dim1d/), varid3) - err = nf90_put_att(ncid, varid3, "cartesian_axis", "x") - - err = nf90_def_var(ncid, 'j', nf90_double, (/dim2d/), varid4) - err = nf90_put_att(ncid, varid4, "cartesian_axis", "y") - - err = nf90_def_var(ncid, 'time', nf90_double, (/dim3d/), varid2) - err = nf90_put_att(ncid, varid2, "cartesian_axis", "T") - err = nf90_put_att(ncid, varid2, "calendar", "noleap") - err = nf90_put_att(ncid, varid2, "units", "days since 0001-01-01 00:00:00") - err = nf90_def_var(ncid, 'runoff', nf90_double, (/dim1d, dim2d, dim3d/), varid) - - err = nf90_enddef(ncid) - err = nf90_put_var(ncid, varid, runoff_in) - err = nf90_put_var(ncid, varid2, (/1., 2., 3., 5., 6., 7., 8., 9., 10., 11./)) - err = nf90_close(ncid) - - deallocate(runoff_in) - -endif !< Wait for the root PE to catch up call mpp_sync @@ -131,8 +69,8 @@ program test_data_override_ongrid call set_calendar_type(NOLEAP) -nlon = 1440 -nlat = 1080 +nlon = 360 +nlat = 180 !< Create a domain nlonXnlat with mask call mpp_domains_set_stack_size(17280000) @@ -140,70 +78,61 @@ program test_data_override_ongrid call mpp_define_io_domain(Domain, (/1,1/)) call mpp_get_data_domain(Domain, is, ie, js, je) -print *, nhalox, nhaloy - -!< Set up the data -allocate(runoff(is:ie,js:je)) +select case (test_case) +case (ongrid) + call generate_ongrid_input_file () +case (bilinear) + call generate_bilinear_input_file () +end select -runoff = 999. +call mpp_sync() +call mpp_error(NOTE, "Finished creating INPUT Files") !< Initiliaze data_override call data_override_init(Ocean_domain_in=Domain, mode=lkind) -!< Run it when time=3 -Time = set_date(1,1,4,0,0,0) -call data_override('OCN','runoff',runoff, Time) -!< Because you are getting the data when time=3, and this is an "ongrid" case, the expected result is just -!! equal to the data at time=3, which is 3. -expected_result = 3._lkind -call compare_data(Domain, runoff, expected_result) - -!< Run it when time=4 -runoff = 999. -Time = set_date(1,1,5,0,0,0) -call data_override('OCN','runoff',runoff, Time) -!< You are getting the data when time=4, the data at time=3 is 3. and at time=5 is 4., so the expected result -!! is the average of the 2 (because this is is an "ongrid" case and there is no horizontal interpolation). -expected_result = (3._lkind + 4._lkind) / 2._lkind -call compare_data(Domain, runoff, expected_result) - -deallocate(runoff) +select case (test_case) +case (ongrid) + call ongrid_test() +case (bilinear) + call bilinear_test() +end select call mpp_exit contains -subroutine compare_data(Domain, actual_result, expected_result) -type(domain2d), intent(in) :: Domain !< Domain with mask table +subroutine compare_data(Domain_in, actual_result, expected_result) +type(domain2d), intent(in) :: Domain_in !< Domain with mask table real(lkind), intent(in) :: expected_result !< Expected result from data_override real(lkind), dimension(:,:), intent(in) :: actual_result !< Result from data_override integer :: xsizec, ysizec !< Size of the compute domain integer :: xsized, ysized !< Size of the data domain integer :: nx, ny !< Size of acual_result -integer :: nhalox, nhaloy !< Size of the halos +integer :: nhx, nhy !< Size of the halos integer :: i, j !< Helper indices !< Data is only expected to be overriden for the compute domain -not at the halos. -call mpp_get_compute_domain(Domain, xsize=xsizec, ysize=ysizec) -call mpp_get_data_domain(Domain, xsize=xsized, ysize=ysized) +call mpp_get_compute_domain(Domain_in, xsize=xsizec, ysize=ysizec) +call mpp_get_data_domain(Domain_in, xsize=xsized, ysize=ysized) !< Note that actual_result has indices at (1:nx,1:ny) not (is:ie,js:je) -nhalox= (xsized-xsizec)/2 -nhaloy = (ysized-ysizec)/2 +nhx= (xsized-xsizec)/2 +nhy = (ysized-ysizec)/2 nx = size(actual_result, 1) ny = size(actual_result, 2) do i = 1, nx do j = 1, ny - if (i <= nhalox .or. i > (nx-nhalox) .or. j <= nhaloy .or. j > (ny-nhaloy)) then + if (i <= nhx .or. i > (nx-nhx) .or. j <= nhy .or. j > (ny-nhy)) then !< This is the result at the halos it should 999. - if (actual_result(i,j) .ne. 999.) then + if (actual_result(i,j) .ne. 999._lkind) then print *, "for i=", i, " and j=", j, " result=", actual_result(i,j) call mpp_error(FATAL, "test_data_override_ongrid: Data was overriden in the halos!!") endif else if (actual_result(i,j) .ne. expected_result) then - print *, "for i=", i, " and j=", j, " result=", actual_result(i,j) + print *, "for i=", i, " and j=", j, " result=", actual_result(i,j), " expected=", expected_result call mpp_error(FATAL, "test_data_override_ongrid: Result is different from expected answer!") endif endif @@ -212,4 +141,300 @@ subroutine compare_data(Domain, actual_result, expected_result) end subroutine +subroutine create_grid_spec_file() + type(FmsNetcdfFile_t) :: fileobj + + if (open_file(fileobj, 'INPUT/grid_spec.nc', 'overwrite')) then + call register_axis(fileobj, 'str', 255) + call register_field(fileobj, 'ocn_mosaic_file', 'char', (/'str'/)) + call write_data(fileobj, 'ocn_mosaic_file', "ocean_mosaic.nc") + call close_file(fileobj) + else + call mpp_error(FATAL, "Error opening the file: 'INPUT/grid_spec.nc' to write") + endif +end subroutine create_grid_spec_file + +subroutine create_ocean_mosaic_file() + type(FmsNetcdfFile_t) :: fileobj + character(len=10) :: dimnames(2) + + dimnames(1) = 'str' + dimnames(2) = 'ntiles' + if (open_file(fileobj, 'INPUT/ocean_mosaic.nc', 'overwrite')) then + call register_axis(fileobj, dimnames(1) , 255) + call register_axis(fileobj, dimnames(2), 1) + call register_field(fileobj, 'gridfiles', 'char', dimnames) + call write_data(fileobj, 'gridfiles', (/"ocean_hgrid.nc"/)) + call close_file(fileobj) + else + call mpp_error(FATAL, "Error opening the file: 'INPUT/ocean_mosaic.nc' to write") + endif +end subroutine create_ocean_mosaic_file + +subroutine create_ocean_hgrid_file() + type(FmsNetcdfFile_t) :: fileobj + + real(lkind), allocatable, dimension(:,:) :: xdata, ydata + integer :: nx, nxp, ny, nyp, i, j + nx = nlon*2 + nxp = nx+1 + ny = nlat*2 + nyp = ny+1 + + allocate(xdata(nxp, nyp)) + xdata(1,:) = 0.0_lkind + do i = 2, nxp + xdata(i,:) = xdata(i-1,:) + 0.5 + enddo + + allocate(ydata(nxp, nyp)) + ydata(:,1) = -90.0_lkind + do i = 2, nyp + ydata(:,i) = ydata(:, i-1) + 0.5 + enddo + + if (open_file(fileobj, 'INPUT/ocean_hgrid.nc', 'overwrite')) then + call register_axis(fileobj, "nx", nx) + call register_axis(fileobj, "ny", ny) + call register_axis(fileobj, "nxp", nxp) + call register_axis(fileobj, "nyp", nyp) + call register_field(fileobj, 'x', 'float', (/'nxp', 'nyp'/)) + call register_field(fileobj, 'y', 'float', (/'nxp', 'nyp'/)) + call register_field(fileobj, 'area', 'float', (/'nx', 'ny'/)) + call write_data(fileobj, "x", xdata) + call write_data(fileobj, "y", ydata) + call close_file(fileobj) + else + call mpp_error(FATAL, "Error opening the file: 'INPUT/ocean_hgrid.nc' to write") + endif +end subroutine create_ocean_hgrid_file + +subroutine create_ongrid_data_file() + type(FmsNetcdfFile_t) :: fileobj + character(len=10) :: dimnames(3) + real(lkind), allocatable, dimension(:,:,:) :: runoff_in + real(lkind), allocatable, dimension(:) :: time_data + integer :: i + + allocate(runoff_in(nlon, nlat, 10)) + allocate(time_data(10)) + do i = 1, 10 + runoff_in(:,:,i) = real(i, lkind) + enddo + time_data = (/1., 2., 3., 5., 6., 7., 8., 9., 10., 11./) + + dimnames(1) = 'i' + dimnames(2) = 'j' + dimnames(3) = 'time' + + if (open_file(fileobj, 'INPUT/runoff.daitren.clim.1440x1080.v20180328.nc', 'overwrite')) then + call register_axis(fileobj, "i", nlon) + call register_axis(fileobj, "j", nlat) + call register_axis(fileobj, "time", unlimited) + + call register_field(fileobj, "i", "float", (/"i"/)) + call register_variable_attribute(fileobj, "i", "cartesian_axis", "x", str_len=1) + + call register_field(fileobj, "j", "float", (/"j"/)) + call register_variable_attribute(fileobj, "j", "cartesian_axis", "y", str_len=1) + + call register_field(fileobj, "time", "float", (/"time"/)) + call register_variable_attribute(fileobj, "time", "cartesian_axis", "T", str_len=1) + call register_variable_attribute(fileobj, "time", "calendar", "noleap", str_len=6) + call register_variable_attribute(fileobj, "time", "units", "days since 0001-01-01 00:00:00", str_len=30) + + call register_field(fileobj, "runoff", "float", dimnames) + call write_data(fileobj, "runoff", runoff_in) + call write_data(fileobj, "time", time_data) + call close_file(fileobj) + else + call mpp_error(FATAL, "Error opening the file: 'INPUT/runoff.daitren.clim.1440x1080.v20180328.nc' to write") + endif + deallocate(runoff_in) +end subroutine create_ongrid_data_file + +subroutine generate_ongrid_input_file() + !< Create some files needed by data_override! + if (mpp_pe() .eq. mpp_root_pe()) then + call create_grid_spec_file() + call create_ocean_mosaic_file() + call create_ocean_hgrid_file() + call create_ongrid_data_file() + endif + call mpp_sync() +end subroutine + +!> @brief Tests ongrid data overrides. +!! In the first case there is no time interpolation +!! In the second case there is time interpolation +subroutine ongrid_test() + real(lkind) :: expected_result !< Expected result from data_override + type(time_type) :: Time !< Time + real(lkind), allocatable, dimension(:,:) :: runoff !< Data to be written + + allocate(runoff(is:ie,js:je)) + + runoff = 999._lkind + !< Run it when time=3 + Time = set_date(1,1,4,0,0,0) + call data_override('OCN','runoff',runoff, Time) + !< Because you are getting the data when time=3, and this is an "ongrid" case, the expected result is just + !! equal to the data at time=3, which is 3. + expected_result = 3._lkind + call compare_data(Domain, runoff, expected_result) + + !< Run it when time=4 + runoff = 999._lkind + Time = set_date(1,1,5,0,0,0) + call data_override('OCN','runoff',runoff, Time) + !< You are getting the data when time=4, the data at time=3 is 3. and at time=5 is 4., so the expected result + !! is the average of the 2 (because this is is an "ongrid" case and there is no horizontal interpolation). + expected_result = (3._lkind + 4._lkind) / 2._lkind + call compare_data(Domain, runoff, expected_result) + + deallocate(runoff) +end subroutine ongrid_test + +!> @brief Creates an input netcdf data file to use for the ongrid data_override test case +!! with either an increasing or decreasing lat, lon grid +subroutine create_bilinear_data_file(increasing_grid) + logical, intent(in) :: increasing_grid !< .true. if increasing a file with an increasing lat/lon + + type(FmsNetcdfFile_t) :: fileobj !< Fms2_io fileobj + character(len=10) :: dimnames(3) !< dimension names for the variable + real(lkind), allocatable :: runoff_in(:,:,:) !< Data to write + real(lkind), allocatable :: time_data(:) !< Time dimension data + real(lkind), allocatable :: lat_data(:) !< Lat dimension data + real(lkind), allocatable :: lon_data(:) !< Lon dimension data + character(len=:), allocatable :: filename !< Name of the file + integer :: factor !< This is used when creating the grid data + !! -1 if the grid is decreasing + !! +1 if the grid is increasing + integer :: i, j, k !< For looping through variables + + integer :: nlon_data + integer :: nlat_data + + nlon_data = nlon + 1 + nlat_data = nlat - 1 + allocate(runoff_in(nlon_data, nlat_data, 10)) + allocate(time_data(10)) + allocate(lat_data(nlat_data)) + allocate(lon_data(nlon_data)) + + if (.not. increasing_grid) then + filename = 'INPUT/bilinear_decreasing.nc' + lon_data(1) = 360.0_lkind + lat_data(1) = 89.0_lkind + factor = -1 + do i = 1, nlon_data + do j = 1, nlat_data + do k = 1, 10 + runoff_in(i, j, k) = real(362-i, kind=lkind) * 1000._lkind + & + real(180-j, kind=lkind) + real(k, kind=lkind)/100._lkind + enddo + enddo + enddo + else + filename = 'INPUT/bilinear_increasing.nc' + lon_data(1) = 0.0_lkind + lat_data(1) = -89.0_lkind + factor = 1 + + do i = 1, nlon_data + do j = 1, nlat_data + do k = 1, 10 + runoff_in(i, j, k) = real(i, kind=lkind) * 1000._lkind + real(j, kind=lkind) + & + real(k, kind=lkind)/100._lkind + enddo + enddo + enddo + endif + + do i = 2, nlon_data + lon_data(i) = real(lon_data(i-1) + 1*factor, lkind) + enddo + + do i = 2, nlat_data + lat_data(i) =real(lat_data(i-1) + 1*factor, lkind) + enddo + + time_data = (/1., 2., 3., 5., 6., 7., 8., 9., 10., 11./) + + dimnames(1) = 'i' + dimnames(2) = 'j' + dimnames(3) = 'time' + + if (open_file(fileobj, filename, 'overwrite')) then + call register_axis(fileobj, "i", nlon_data) + call register_axis(fileobj, "j", nlat_data) + call register_axis(fileobj, "time", unlimited) + + call register_field(fileobj, "i", "float", (/"i"/)) + call register_variable_attribute(fileobj, "i", "cartesian_axis", "x", str_len=1) + + call register_field(fileobj, "j", "float", (/"j"/)) + call register_variable_attribute(fileobj, "j", "cartesian_axis", "y", str_len=1) + + call register_field(fileobj, "time", "float", (/"time"/)) + call register_variable_attribute(fileobj, "time", "cartesian_axis", "T", str_len=1) + call register_variable_attribute(fileobj, "time", "calendar", "noleap", str_len=6) + call register_variable_attribute(fileobj, "time", "units", "days since 0001-01-01 00:00:00", str_len=30) + + call register_field(fileobj, "runoff", "float", dimnames) + call write_data(fileobj, "runoff", runoff_in) + call write_data(fileobj, "i", lon_data) + call write_data(fileobj, "j", lat_data) + call write_data(fileobj, "time", time_data) + call close_file(fileobj) + else + call mpp_error(FATAL, "Error opening the file: 'INPUT/bilinear_increasing.nc' to write") + endif + deallocate(runoff_in) +end subroutine create_bilinear_data_file + +!> @brief Generates the input for the bilinear data_override test_case +subroutine generate_bilinear_input_file() + if (mpp_pe() .eq. mpp_root_pe()) then + call create_grid_spec_file () + call create_ocean_mosaic_file() + call create_ocean_hgrid_file() + call create_bilinear_data_file(.true.) + call create_bilinear_data_file(.false.) + endif + call mpp_sync() +end subroutine generate_bilinear_input_file + +!> @brief Tests bilinear data_override with and increasing and decreasing grid case +!! and comares the output betweeen the cases to ensure it is correct +subroutine bilinear_test() + type(time_type) :: Time !< Time + real(lkind), allocatable, dimension(:,:) :: runoff_decreasing !< Data to be written + real(lkind), allocatable, dimension(:,:) :: runoff_increasing !< Data to be written + + integer :: i, j, k + logical :: success + + allocate(runoff_decreasing(is:ie,js:je)) + allocate(runoff_increasing(is:ie,js:je)) + + runoff_decreasing = 999_lkind + runoff_increasing = 999_lkind + Time = set_date(1,1,4,0,0,0) + call data_override('OCN','runoff_increasing',runoff_increasing, Time, override=success) + if (.not. success) call mpp_error(FATAL, "Data override failed") + call data_override('OCN','runoff_decreasing',runoff_decreasing, Time, override=success) + if (.not. success) call mpp_error(FATAL, "Data override failed") + + do i = is, ie + do j = js, je + if (abs(runoff_decreasing(i,j) - runoff_increasing(i,j)) .gt. 1) then + call mpp_error(FATAL, "The data is not the same: "// & + string(i)//","//string(j)//":"// & + string(runoff_decreasing(i,j))//" vs "//string(runoff_increasing(i,j))) + endif + enddo + enddo + deallocate(runoff_decreasing, runoff_increasing) +end subroutine bilinear_test end program test_data_override_ongrid diff --git a/test_fms/horiz_interp/test_horiz_interp.F90 b/test_fms/horiz_interp/test_horiz_interp.F90 index c4e333b786..fd0d077a91 100644 --- a/test_fms/horiz_interp/test_horiz_interp.F90 +++ b/test_fms/horiz_interp/test_horiz_interp.F90 @@ -52,10 +52,11 @@ program horiz_interp_test integer :: ni_dst = 144, nj_dst = 72 integer, parameter :: max_neighbors = 400 !! took this from spherical mod !! max amount found neighbors to loop through in spherical search - + logical :: decreasing_lat = .False. !< latitude and longitude are decreasing instead of increasing + !! This is only for the bilinear test case namelist /test_horiz_interp_nml/ test_conserve, test_bicubic, test_spherical, test_bilinear, test_assign, test_solo,& - ni_src, nj_src, ni_dst,nj_dst + ni_src, nj_src, ni_dst,nj_dst, decreasing_lat type(domain2d) :: domain @@ -196,12 +197,20 @@ subroutine test_horiz_interp_bilinear real(HI_TEST_KIND_), allocatable, dimension(:) :: lon1D_src, lat1D_src, lon1D_dst, lat1D_dst real(HI_TEST_KIND_), allocatable, dimension(:,:) :: lon2D_src, lat2d_src, lon2D_dst, lat2D_dst real(HI_TEST_KIND_), allocatable, dimension(:,:) :: data_src, data_dst - real(HI_TEST_KIND_), parameter :: lon_src_beg = 0._lkind, lon_src_end = 360.0_lkind - real(HI_TEST_KIND_), parameter :: lat_src_beg = -90._lkind, lat_src_end = 90._lkind + real(HI_TEST_KIND_) :: lon_src_beg = 0._lkind, lon_src_end = 360.0_lkind + real(HI_TEST_KIND_) :: lat_src_beg = -90._lkind, lat_src_end = 90._lkind real(HI_TEST_KIND_), parameter :: D2R = real(PI,lkind)/180._lkind + real(HI_TEST_KIND_), parameter :: R2D = 180._lkind/real(PI,lkind) type(horiz_interp_type) :: interp + if (decreasing_lat) then + lon_src_beg = 360.0_lkind + lon_src_end = 0._lkind + lat_src_beg = 90._lkind + lat_src_end = -90._lkind + endif + allocate( lon1D_src(ni_src+1), lat1D_src(nj_src+1) ) allocate( lon1D_dst(ni_src+1), lat1D_dst(nj_src+1) ) allocate( lon2d_src(ni_src,nj_src), lat2d_src(ni_src,nj_src) ) @@ -377,6 +386,7 @@ subroutine test_horiz_interp_bilinear call check_dealloc(interp) endif + if (decreasing_lat) return ! --- 2dx1d version bilinear interpolation data_dst = 0.0_lkind lon1d_dst = lon1d_src diff --git a/test_fms/horiz_interp/test_horiz_interp2.sh b/test_fms/horiz_interp/test_horiz_interp2.sh index 485be882fd..cb8868d4e6 100755 --- a/test_fms/horiz_interp/test_horiz_interp2.sh +++ b/test_fms/horiz_interp/test_horiz_interp2.sh @@ -111,10 +111,45 @@ cat <<_EOF > input.nml / _EOF -test_expect_success "bilinear method with real kind=4" ' +test_expect_success "bilinear method with real kind=4 increasing latitude/longitude" ' mpirun -n 2 ./test_horiz_interp_r4 ' -test_expect_success "bilinear method with real kind=8" ' +test_expect_success "bilinear method with real kind=8 increasing latitude/longitude" ' + mpirun -n 2 ./test_horiz_interp_r8 +' + +cat <<_EOF > input.nml +&test_horiz_interp_nml + test_bilinear= .true. + ni_src = 360 + nj_src = 180 + ni_dst = 144 + nj_dst = 72 + decreasing_lat = .true. +/ +_EOF + +test_expect_success "bilinear method with real kind=4 decreasing latitude/longitude" ' + mpirun -n 2 ./test_horiz_interp_r4 +' +test_expect_success "bilinear method with real kind=8 decreasing latitude/longitude" ' + mpirun -n 2 ./test_horiz_interp_r8 +' +cat <<_EOF > input.nml +&test_horiz_interp_nml + test_bilinear= .true. + test_solo = .true. + ni_src = 360 + nj_src = 180 + ni_dst = 144 + nj_dst = 72 +/ +_EOF + +test_expect_success "bilinear method solo wrapper with real kind=4 increasing latitude/longitude" ' + mpirun -n 2 ./test_horiz_interp_r4 +' +test_expect_success "bilinear method solo wrapper with real kind=8 increasing latitude/longitude" ' mpirun -n 2 ./test_horiz_interp_r8 ' @@ -126,13 +161,14 @@ cat <<_EOF > input.nml nj_src = 180 ni_dst = 144 nj_dst = 72 + decreasing_lat = .true. / _EOF -test_expect_success "bilinear method solo wrapper with real kind=4" ' +test_expect_success "bilinear method solo wrapper with real kind=4 decreasing latitude/longitude" ' mpirun -n 2 ./test_horiz_interp_r4 ' -test_expect_success "bilinear method solo wrapper with real kind=8" ' +test_expect_success "bilinear method solo wrapper with real kind=8 decreasing latitude/longitude" ' mpirun -n 2 ./test_horiz_interp_r8 '