test_fms/axis_utils/test_axis_utils.F90

   1 !***********************************************************************
   2 !*                   GNU Lesser General Public License
   3 !*
   4 !* This file is part of the GFDL Flexible Modeling System (FMS).
   5 !*
   6 !* FMS is free software: you can redistribute it and/or modify it under
   7 !* the terms of the GNU Lesser General Public License as published by
   8 !* the Free Software Foundation, either version 3 of the License, or (at
   9 !* your option) any later version.
  10 !*
  11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
  12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 !* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 !* for more details.
  15 !*
  16 !* You should have received a copy of the GNU Lesser General Public
  17 !* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
  18 !***********************************************************************
  19
  20 program test_axis_utils
  21
  22 use fms_mod,         only : fms_init, fms_end, lowercase
  23 use fms2_io_mod, only: FmsNetcdfFile_t, open_file, close_file, register_axis, register_field, &
  24                      & register_variable_attribute, write_data
  25 use platform_mod, only: r4_kind, r8_kind
  26 use mpp_mod, only: mpp_error, fatal, stderr
  27 use fms_string_utils_mod, only: string, stringify
  28 use axis_utils2_mod
  29
  30 implicit none
  31
  32 type GetAxisCartTest_t
  33   type(FmsNetcdfFile_t) :: fileobj
  34   type(GetAxisCartTestCase_t), pointer :: test0, test1
  35 end type
  36
  37 type GetAxisCartTestCase_t
  38   character(:), allocatable :: var
  39   character(1) :: cart
  40   type(GetAxisCartTestCase_t), pointer :: next => NULL()
  41 end type
  42
  43 integer, parameter :: k = AU_TEST_KIND_
  44 real(k), parameter :: pi = 4._k * atan(1._k)
  45
  46 integer :: i
  47 character(100) :: arg
  48
  49 call fms_init
  50
  51 do i=1,command_argument_count()
  52   call get_command_argument(i, arg)
  53
  54   select case (arg)
  55     case ('--get-axis-modulo')
  56       print "(A)", "Testing get_axis_modulo"
  57       call test_get_axis_modulo
  58
  59     case ('--get-axis-modulo-times')
  60       print "(A)", "Testing get_axis_modulo_times"
  61       call test_get_axis_modulo_times
  62
  63     case ('--get-axis-cart')
  64       print "(A)", "Testing get_axis_cart"
  65       call test_get_axis_cart
  66
  67     case ('--lon-in-range')
  68       print "(A)", "Testing lon_in_range"
  69       call test_lon_in_range
  70
  71     case ('--frac-index')
  72       print "(A)", "Testing frac_index"
  73       call test_frac_index
  74
  75     case ('--frac-index-fail')
  76       print "(A)", "Testing frac_index (FAILURE)"
  77       call test_frac_index_fail
  78
  79     case ('--nearest-index-increasing')
  80       print "(A)", "Testing nearest_index with a monotonically increasing array"
  81       call test_nearest_index(.true.)
  82
  83     case ('--nearest-index-decreasing')
  84       print "(A)", "Testing nearest_index with a monotonically decreasing array"
  85       call test_nearest_index(.false.)
  86
  87     case ('--nearest-index-fail')
  88       print "(A)", "Testing nearest_index (FAILURE)"
  89       call test_nearest_index_fail
  90
  91     case ('--axis-edges-increasing')
  92       print "(A)", "Testing axis_edges-increasing"
  93       call test_axis_edges(.true.)
  94
  95     case ('--axis-edges-decreasing')
  96       print "(A)", "Testing axis_edges-decreasing"
  97       call test_axis_edges(.false.)
  98
  99     case ('--tranlon')
 100       print "(A)", "Testing tranlon"
 101       call test_tranlon
 102
 103     case ('--interp-1d-1d')
 104       print "(A)", "Testing interp_1d_1d"
 105       call test_interp_1d_1d
 106
 107     case ('--interp-1d-2d')
 108       print "(A)", "Testing interp_1d_2d"
 109       call test_interp_1d_2d
 110
 111     case ('--interp-1d-3d')
 112       print "(A)", "Testing interp_1d_3d"
 113       call test_interp_1d_3d
 114
 115     case default
 116       write(stderr(),"(A)") "Unrecognized command line option: " // trim(arg)
 117   end select
 118 enddo
 119
 120 call fms_end
 121
 122 contains
 123
 124 ! Status: TODO
 125 ! function get_axis_modulo(fileobj, axisname)
 126 subroutine test_get_axis_modulo
 127   type(FmsNetcdfFile_t) :: fileobj
 128
 129   write(stderr(), "(A)") "Warning: get_axis_modulo unit test not yet implemented"
 130 end subroutine
 131
 132 ! Status: TODO
 133 ! function get_axis_modulo_times(fileobj, axisname, tbeg, tend)
 134 subroutine test_get_axis_modulo_times
 135   type(FmsNetcdfFile_t) :: fileobj
 136
 137   write(stderr(), "(A)") "Warning: get_axis_modulo_times unit test not yet implemented"
 138 end subroutine
 139
 140 subroutine test_get_axis_cart
 141   type(GetAxisCartTest_t) :: test
 142   type(GetAxisCartTestCase_t), pointer :: test_nonexistent_var
 143   character(:), allocatable :: var_name, attr_name, attr_value
 144   integer :: i, j
 145
 146   character(*), parameter, dimension(*) :: &
 147     & special_axis_names_x = [character(12) :: "lon", "x", "degrees_e", "degrees_east", "degreese"], &
 148     & special_axis_names_y = [character(13) :: "lat", "y", "degrees_n", "degrees_north", "degreesn"], &
 149     & special_axis_names_z = [character(6) :: "depth", "height", "z", "cm", "m", "pa", "hpa"], &
 150     & special_axis_names_t = [character(4) :: "time", "t", "sec", "min", "hou", "day", "mon", "yea"], &
 151     & attr_names           = [character(14) :: "cartesian_axis", "axis"], &
 152     & xyzt_uc              = ["X", "Y", "Z", "T"]
 153
 154   call open_netcdf_w(test%fileobj)
 155   call register_axis(test%fileobj, "dim1", 1)
 156
 157   ! Check a variable which does not exist
 158
 159   allocate(test_nonexistent_var)
 160   test_nonexistent_var%var = "does_not_exist"
 161   test_nonexistent_var%cart = "N"
 162
 163   test%test0 => test_nonexistent_var
 164   test%test1 => test_nonexistent_var
 165
 166   ! Check a variable which exists, but which has neither a "cartesian_axis" nor an "axis" attribute.
 167   var_name = "exists_no_attributes"
 168   call get_axis_cart_test_add(test, var_name, "N")
 169
 170   do i=1,size(attr_names)
 171     attr_name = trim(attr_names(i))
 172
 173     ! Check an unknown value on a "cartesian_axis" or "axis" attribute.
 174     ! TODO: This test fails. It should be uncommented if/when get_axis_cart's behavior is fixed.
 175
 176     !attr_value = "unexpected"
 177     !var_name = attr_name // "_attr_value_" // attr_value
 178     !call get_axis_cart_test_add(test, var_name, "N")
 179     !call register_variable_attribute(test%fileobj, var_name, attr_name, attr_value, str_len=len(attr_value))
 180
 181     do j=1,size(xyzt_uc)
 182       ! Check upper-case "axis" attributes"
 183       attr_value = xyzt_uc(j)
 184       var_name = attr_name // "_attr_value_" // attr_value
 185       call get_axis_cart_test_add(test, var_name, xyzt_uc(j))
 186       call register_variable_attribute(test%fileobj, var_name, attr_name, attr_value, str_len=len(attr_value))
 187
 188       ! Check lower-case "axis" attributes"
 189       attr_value = lowercase(xyzt_uc(j))
 190       var_name = attr_name // "_attr_value_" // attr_value
 191       call get_axis_cart_test_add(test, var_name, xyzt_uc(j))
 192       call register_variable_attribute(test%fileobj, var_name, attr_name, attr_value, str_len=len(attr_value))
 193     enddo
 194   enddo
 195
 196   call test_special_axis_names(test, special_axis_names_x, "X")
 197   call test_special_axis_names(test, special_axis_names_y, "Y")
 198   call test_special_axis_names(test, special_axis_names_z, "Z")
 199   call test_special_axis_names(test, special_axis_names_t, "T")
 200
 201   call close_file(test%fileobj)
 202
 203   call get_axis_cart_tests_run(test)
 204 end subroutine
 205
 206 subroutine get_axis_cart_test_add(test, var_name, cart)
 207   type(GetAxisCartTest_t), intent(inout) :: test
 208   type(GetAxisCartTestCase_t), pointer :: test_case
 209   character(*), intent(in) :: var_name
 210   character(1), intent(in) :: cart
 211   character(:), allocatable :: kind_str
 212
 213   if (k .eq. r4_kind) then
 214     kind_str = "float"
 215   else
 216     kind_str = "double"
 217   endif
 218
 219   call register_field(test%fileobj, var_name, kind_str, dimensions=["dim1"])
 220
 221   allocate(test_case)
 222   test_case%var = var_name
 223   test_case%cart = cart
 224
 225   test%test1%next => test_case
 226   test%test1 => test_case
 227 end subroutine
 228
 229 subroutine get_axis_cart_tests_run(test)
 230   type(GetAxisCartTest_t), intent(inout) :: test
 231   type(GetAxisCartTestCase_t), pointer :: test_case, next
 232   character(1) :: cart_test
 233   integer :: i
 234
 235   call open_netcdf_r(test%fileobj)
 236
 237   test_case => test%test0
 238
 239   do while (associated(test_case))
 240     cart_test = " "
 241     call get_axis_cart(test%fileobj, test_case%var, cart_test)
 242
 243     if (cart_test .ne. test_case%cart) then
 244       write(stderr(), "(A)") "get_axis_cart result for variable '" // test_case%var // "': " // cart_test
 245       write(stderr(), "(A)") "Expected result: " // test_case%cart
 246       call mpp_error(FATAL, "get_axis_cart unit test failed")
 247     endif
 248
 249     next => test_case%next
 250     deallocate(test_case)
 251     test_case => next
 252   enddo
 253
 254   call close_file(test%fileobj)
 255 end subroutine
 256
 257 subroutine test_special_axis_names(test, special_axis_names, ret_expected)
 258   type(GetAxisCartTest_t), intent(inout) :: test
 259   character(*), intent(in) :: special_axis_names(:), ret_expected
 260   character(:), allocatable :: var_name
 261   integer :: i
 262
 263   do i=1,size(special_axis_names)
 264     var_name = trim(special_axis_names(i))
 265     call get_axis_cart_test_add(test, var_name, ret_expected)
 266   enddo
 267 end subroutine
 268
 269 subroutine test_lon_in_range
 270   real(k), parameter :: eps_big = 1e-3_k, eps_tiny = 1e-5_k
 271   real(k), parameter :: pi_plus_360 = 360._k + pi
 272
 273   ! Test some cases where no translation is needed
 274   call lon_in_range_assert(0._k,              0._k,  0._k)
 275   call lon_in_range_assert(1._k,              0._k,  1._k)
 276   call lon_in_range_assert(350._k,            0._k,  350._k)
 277   call lon_in_range_assert(1._k,              1._k,  1._k)
 278   call lon_in_range_assert(350._k,            1._k,  350._k)
 279   call lon_in_range_assert(359._k,            0._k,  359._k)
 280   call lon_in_range_assert(359._k,            1._k,  359._k)
 281   call lon_in_range_assert(pi,                0._k,  pi)
 282
 283   ! Test up-translation
 284   call lon_in_range_assert(-2._k,             -1._k, 358._k)
 285   call lon_in_range_assert(-2._k,             0._k,  358._k)
 286   call lon_in_range_assert(-2._k,             5._k,  358._k)
 287   call lon_in_range_assert(-1._k,             0._k,  359._k)
 288   call lon_in_range_assert(-1._k,             5._k,  359._k)
 289   call lon_in_range_assert(0._k,              5._k,  360._k)
 290   call lon_in_range_assert(1._k,              5._k,  361._k)
 291   call lon_in_range_assert(-pi,               0._k,  360._k - pi)
 292
 293   ! Test down-translation
 294   call lon_in_range_assert(359._k,            -1._k, -1._k)
 295   call lon_in_range_assert(360._k,            -1._k, 0._k)
 296   call lon_in_range_assert(360._k,            0._k,  0._k)
 297   call lon_in_range_assert(361._k,            -1._k, 1._k)
 298   call lon_in_range_assert(361._k,            0._k,  1._k)
 299   call lon_in_range_assert(362._k,            -1._k, 2._k)
 300   call lon_in_range_assert(362._k,            0._k,  2._k)
 301   call lon_in_range_assert(pi_plus_360,       0._k,  pi_plus_360 - 360._k)
 302
 303   ! Test rounding behavior
 304   call lon_in_range_assert(eps_tiny,          0._k,  0._k)
 305   call lon_in_range_assert(eps_big,           0._k,  eps_big)
 306   call lon_in_range_assert(360._k - eps_tiny, 0._k,  0._k)
 307   call lon_in_range_assert(360._k - eps_big,  0._k,  360._k - eps_big)
 308 end subroutine
 309
 310 subroutine lon_in_range_assert(lon, l_start, ret_expected)
 311   real(k), intent(in) :: lon, l_start, ret_expected
 312   real(k) :: ret_test
 313
 314   ret_test = lon_in_range(lon, l_start)
 315
 316   if (ret_test /= ret_expected) then
 317     write(stderr(), "(A)") "lon_in_range(" // string(lon) // ", " // string(l_start) // &
 318                          & ") returned erroneous value: " // string(ret_test)
 319     write(stderr(), "(A)") "Expected return value: " // string(ret_expected)
 320     call mpp_error(FATAL, "lon_in_range unit test failed")
 321   endif
 322 end subroutine
 323
 324 #define CALC_FRAC_INDEX_(i, v, values) real(i, k) + (v - values(i)) / (values(i + 1) - values(i))
 325
 326 subroutine test_frac_index
 327   real(k) :: values(6), v, fi
 328   integer :: i, n
 329   real(k), parameter :: f10=.1_k, f25=.25_k, f50=.5_k, f99=.99_k
 330
 331   values = [1._k, 2._k, 3._k, 5._k, 10._k, 11._k]
 332   n = size(values)
 333
 334   ! Test values outside of the input array
 335   call frac_index_assert(real(values(1), k) - f50, values, -1._k)
 336   call frac_index_assert(real(values(n), k) + f50, values, -1._k)
 337
 338   ! Test the actual indices
 339   do i=1,n
 340     v = values(i)
 341     call frac_index_assert(v, values, real(i, k))
 342   enddo
 343
 344   ! Test the 10% point
 345   do i=1,n-1
 346     v = values(i) + f10*(values(i+1) - values(i))
 347     fi = CALC_FRAC_INDEX_(i, v, values)
 348     call frac_index_assert(v, values, fi)
 349   enddo
 350
 351   ! Test the 25% point
 352   do i=1,n-1
 353     v = values(i) + f25*(values(i+1) - values(i))
 354     fi = CALC_FRAC_INDEX_(i, v, values)
 355     call frac_index_assert(v, values, fi)
 356   enddo
 357
 358   ! Test the mid-point
 359   do i=1,n-1
 360     v = values(i) + f50*(values(i+1) - values(i))
 361     fi = CALC_FRAC_INDEX_(i, v, values)
 362     call frac_index_assert(v, values, fi)
 363   enddo
 364
 365   ! Test the 99% point
 366   do i=1,n-1
 367     v = values(i) + f99*(values(i+1) - values(i))
 368     fi = CALC_FRAC_INDEX_(i, v, values)
 369     call frac_index_assert(v, values, fi)
 370   enddo
 371 end subroutine
 372
 373 subroutine frac_index_assert(fval, arr, ret_expected)
 374   real(k), intent(in) :: fval, arr(:), ret_expected
 375   real(k) :: ret_test
 376
 377   ret_test = frac_index(fval, arr)
 378
 379   if (ret_test /= ret_expected) then
 380     write(stderr(), "(A)") "frac_index(" // string(fval) // ", " // stringify(arr) // &
 381                          & ") returned erroneous value: " // string(ret_test)
 382     write(stderr(), "(A)") "Expected return value: " // string(ret_expected)
 383     call mpp_error(FATAL, "frac_index unit test failed")
 384   endif
 385 end subroutine
 386
 387 ! Test that frac_index fails with a non-monotonic array
 388 subroutine test_frac_index_fail
 389   real(k) :: values(5)
 390   real(k) :: ret_test
 391
 392   values = [1._k, 2._k, 4._k, 3._k, 5._k]
 393   ret_test = frac_index(1.5_k, values)
 394 end subroutine
 395
 396 subroutine test_nearest_index(increasing_array)
 397   logical, intent(in) :: increasing_array !< .True. if test using an increasing array
 398   real(k) :: arr(7)
 399   integer :: ans(16)
 400
 401   if (increasing_array) then
 402     arr = [-6._k, -3._k, 5._k, 12._k, 20._k, 40._k, 100._k]
 403     ans=(/1, 7, 3, 4, 5, 6, 7, 3, 4, 4, 5, 5, 1, 2, 1, 2/)
 404   else
 405     arr = [100._k, 40._k, 20._k, 12._k, 5._k, -3._k, -6._k]
 406     ans=(/7, 1, 5, 4, 3, 2, 1, 5, 4, 4, 3, 3, 7, 6, 7, 6/)
 407   endif
 408
 409   ! Test values beyond array boundaries
 410   call nearest_index_assert(-7._k,    arr, ans(1))
 411   call nearest_index_assert(1000._k, arr, ans(2))
 412
 413   ! Test values actually in the array
 414   call nearest_index_assert(5._k,    arr, ans(3))
 415   call nearest_index_assert(12._k,   arr, ans(4))
 416   call nearest_index_assert(20._k,   arr, ans(5))
 417   call nearest_index_assert(40._k,   arr, ans(6))
 418   call nearest_index_assert(100._k,  arr, ans(7))
 419
 420   ! Test the intervals between array values
 421   call nearest_index_assert(6._k,    arr, ans(8))
 422   call nearest_index_assert(11._k,   arr, ans(9))
 423   call nearest_index_assert(15._k,   arr, ans(10))
 424   call nearest_index_assert(18._k,   arr, ans(11))
 425   call nearest_index_assert(29._k,   arr, ans(12))
 426
 427   ! Test the negative numbers
 428   call nearest_index_assert(-6._k,    arr, ans(13))
 429   call nearest_index_assert(-3._k,    arr, ans(14))
 430   call nearest_index_assert(-5._k,    arr, ans(15))
 431   call nearest_index_assert(-1._k,    arr, ans(16))
 432
 433 end subroutine
 434
 435 subroutine nearest_index_assert(val, arr, ret_expected)
 436   real(k), intent(in) :: val, arr(:)
 437   integer, intent(in) :: ret_expected
 438   integer :: ret_test
 439
 440   ret_test = nearest_index(val, arr)
 441
 442   if (ret_test /= ret_expected) then
 443     write(stderr(), "(A)") "nearest_index(" // string(val) // ", " // stringify(arr) // &
 444                          & ") returned erroneous value: " // string(ret_test)
 445     write(stderr(), "(A)") "Expected return value: " // string(ret_expected)
 446     call mpp_error(FATAL, "nearest_index unit test failed")
 447   endif
 448 end subroutine
 449
 450 ! Test that nearest_index fails with a non-monotonic array
 451 subroutine test_nearest_index_fail
 452   real(k) :: arr(5)
 453   integer :: ret_test
 454
 455   arr=[5._k, 12._k, 40._k, 20._k, 100._k]
 456   ret_test = nearest_index(5._k, arr)
 457 end subroutine
 458
 459 subroutine test_axis_edges(increasing_array)
 460   logical, intent(in) :: increasing_array !< .True. if test using an increasing array
 461   real(k) :: data_in_var(10)
 462   real(k) :: data_in_var_edges(2,10)
 463   real(k) :: data_in_answers(11)
 464   type(FmsNetcdfFile_t) :: fileobj
 465   real(k)    :: answers(11)
 466   integer :: count
 467   integer :: count_factor
 468   integer :: factor
 469   integer :: index !< For looping through the data
 470
 471   if (increasing_array) then
 472     count = 0
 473     factor = 1
 474     count_factor = -1
 475   else
 476     count = 11
 477     factor = -1
 478     count_factor = 0
 479   endif
 480
 481   do index=1,10
 482      count = count + factor
 483      data_in_var(index) = real(count, k) - 0.5_k
 484
 485      data_in_var_edges(1,index) = real(count-1, k)
 486      data_in_var_edges(2,index) = real(count, k)
 487
 488      data_in_answers(index) = real(count + count_factor, k)
 489   enddo
 490
 491   if (increasing_array) then
 492     data_in_answers(11) = real(count, k)
 493   else
 494     data_in_answers(11) = real(count + factor, k)
 495   endif
 496
 497   call open_netcdf_w(fileobj)
 498
 499   call register_axis(fileobj, "dim1", 10)
 500   call register_axis(fileobj, "dim2", 2)
 501
 502   call register_field(fileobj, "axis", "double", dimensions=["dim1"])
 503
 504   call register_field(fileobj, "axis_with_bounds", "double", dimensions=["dim1"])
 505   call register_variable_attribute(fileobj, "axis_with_bounds", "bounds", "bounds", str_len=6)
 506   call register_field(fileobj, "bounds", "double", dimensions=["dim2", "dim1"])
 507
 508   call register_field(fileobj, "axis_with_edges", "double", dimensions=["dim1"])
 509   call register_variable_attribute(fileobj, "axis_with_edges", "edges", "edges"//char(0), str_len=6)
 510   call register_field(fileobj, "edges", "double", dimensions=["dim2", "dim1"])
 511
 512   call write_data(fileobj, "axis", data_in_var)
 513   call write_data(fileobj, "axis_with_bounds", data_in_var)
 514   call write_data(fileobj, "axis_with_edges", data_in_var)
 515   call write_data(fileobj, "bounds", data_in_var_edges)
 516   call write_data(fileobj, "edges", data_in_var_edges)
 517
 518   call close_file(fileobj)
 519
 520   call open_netcdf_r(fileobj)
 521
 522   !< Case 1: Here the variable "axis" in the file does not have the attribute "bounds" or "edges", so
 523   !! it calculates them from the data in "axis"
 524   answers = 0._k
 525   call axis_edges(fileobj, "axis", answers)
 526   call array_compare_1d(answers, data_in_answers, "axis_edges unit test failed (case 1)")
 527
 528   !< Case 2: Here the variable "axis_with_bounds" in the file has the attribute
 529   !! "bounds", so the data is read from the variable "bounds"
 530   answers = 0._k
 531   call axis_edges(fileobj, "axis_with_bounds", answers)
 532   call array_compare_1d(answers, data_in_answers, "axis_edges unit test failed (case 2)")
 533
 534   !< Case 3: Here the variable "axis_with_edges" in the file has the attribute
 535   !"edges", so the data is read from the variable "edges"
 536   answers = 0._k
 537   call axis_edges(fileobj, "axis_with_edges", answers)
 538   call array_compare_1d(answers, data_in_answers, "axis_edges unit test failed (case 3)")
 539
 540   !< Case 4: Here the flag "reproduce_null_char_bug_flag" is turned on, so the
 541   !! edges are calculated from the data in axis because edges has a null character
 542   !! in the end
 543   answers = 0._k
 544   call axis_edges(fileobj, "axis_with_edges", answers, reproduce_null_char_bug_flag=.true.)
 545   call array_compare_1d(answers, data_in_answers, "axis_edges unit test failed (case 4)")
 546
 547   call close_file(fileobj)
 548 end subroutine
 549
 550 subroutine test_tranlon
 551   real(k), dimension(5) :: lon1, lon2, lon3
 552
 553   lon1 = [1._k, 2._k, 3._k, 4._k,   5._k]
 554   lon2 = [2._k, 3._k, 4._k, 5._k,   361._k]
 555   lon3 = [3._k, 4._k, 5._k, 361._k, 362._k]
 556
 557   ! The first two cases fail due to tranlon's unexpected behavior when no elements are translated.
 558   ! TODO: Uncomment these tests if/when tranlon's behavior is fixed.
 559
 560   !call tranlon_assert(lon1, lon1, 0.0_k,    1)
 561   !call tranlon_assert(lon1, lon1, 1.0_k,    1)
 562
 563   call tranlon_assert(lon1, lon2, 1.5_k,    2)
 564   call tranlon_assert(lon1, lon2, 2.0_k,    2)
 565   call tranlon_assert(lon1, lon3, 2.001_k,  3)
 566 end subroutine
 567
 568 subroutine tranlon_assert(lon0, lon_expected, lon_start, istrt_expected)
 569   real(k), intent(in) :: lon0(:), lon_expected(:), lon_start
 570   integer, intent(in) :: istrt_expected
 571   integer :: istrt_test, i
 572   real(k) :: lon_test(size(lon0))
 573   character(:), allocatable :: test_name
 574
 575   test_name = "tranlon(" // stringify(lon0) // ", " // string(lon_start) // ", istrt)"
 576
 577   lon_test = lon0
 578   call tranlon(lon_test, lon_start, istrt_test)
 579   call array_compare_1d(lon_test, lon_expected, test_name // " unit test failed")
 580
 581   if (istrt_test.ne.istrt_expected) then
 582     write(stderr(), "(A)") test_name // " returned erroneous istrt value: " // string(istrt_test)
 583     write(stderr(), "(A)") "Expected istrt value: " // string(istrt_expected)
 584     call mpp_error(FATAL, "tranlon unit test failed")
 585   endif
 586 end subroutine
 587
 588 ! Status: SKELETAL
 589 ! TODO: More comprehensive interp_1d_1d test
 590 subroutine test_interp_1d_1d
 591   real(k) :: grid1(8), grid2(5), data1(8), data2(5)
 592
 593   grid1 = [1._k, 2._k, 3._k, 4._k, 5._k, 6._k, 7._k, 8._k]
 594   grid2 = [2._k, 3._k, 4._k, 5._k, 6._k]
 595   data1 = [101._k, 102._k, 103._k, 104._k, 105._k, 106._k, 107._k, 108._k]
 596   data2 = [102._k, 103._k, 104._k, 105._k, 106._k]
 597
 598   call interp_1d_1d_assert(grid1, grid2, data1, data2, "linear")
 599   call interp_1d_1d_assert(grid1, grid2, data1, data2, "cubic_spline")
 600 end subroutine
 601
 602 subroutine interp_1d_1d_assert(grid1, grid2, data1, data2_expected, method, yp1, yp2)
 603   real(k), intent(in), dimension(:) :: grid1, grid2, data1, data2_expected
 604   character(*), intent(in), optional :: method
 605   real(k), intent(in), optional :: yp1, yp2
 606   real(k) :: data2_test(size(data2_expected))
 607   character(:), allocatable :: test_name
 608
 609   test_name = "interp_1d_1d(" // &
 610               stringify(grid1) // ", " // &
 611               stringify(grid2) // ", " // &
 612               stringify(data1) // ", data2"
 613
 614   if (present(method)) then
 615     test_name = test_name // ", method=" // method
 616   endif
 617
 618   if (present(yp1)) then
 619     test_name = test_name // ", yp1=" // string(yp1)
 620   endif
 621
 622   if (present(yp2)) then
 623     test_name = test_name // ", yp2=" // string(yp2)
 624   endif
 625
 626   test_name = test_name // ")"
 627
 628   call interp_1d(grid1, grid2, data1, data2_test, method, yp1, yp2)
 629   call array_compare_1d(data2_test, data2_expected, test_name // " unit test failed")
 630 end subroutine
 631
 632 ! Status: SKELETAL
 633 ! TODO: More comprehensive interp_1d_2d test
 634 subroutine test_interp_1d_2d
 635   real(k) :: grid1(2,4), grid2(2,2), data1(2,4), data2(2,2)
 636
 637   grid1(1,:) = [1._k, 2._k, 3._k, 4._k]
 638   grid1(2,:) = [5._k, 6._k, 7._k, 8._k]
 639
 640   grid2(1,:) = [2._k, 3._k]
 641   grid2(2,:) = [6._k, 7._k]
 642
 643   data1(1,:) = [101._k, 102._k, 103._k, 104._k]
 644   data1(2,:) = [105._k, 106._k, 107._k, 108._k]
 645
 646   data2(1,:) = [102._k, 103._k]
 647   data2(2,:) = [106._k, 107._k]
 648
 649   call interp_1d_2d_assert(grid1, grid2, data1, data2)
 650 end subroutine
 651
 652 subroutine interp_1d_2d_assert(grid1, grid2, data1, data2_expected)
 653   real(k), intent(in), dimension(:,:) :: grid1, grid2, data1, data2_expected
 654   real(k) :: data2_test(size(data2_expected,1), size(data2_expected,2))
 655   character(:), allocatable :: test_name
 656
 657   test_name = "interp_1d_2d(" // &
 658               stringify(grid1) // ", " // &
 659               stringify(grid2) // ", " // &
 660               stringify(data1) // ", data2)"
 661
 662   call interp_1d(grid1, grid2, data1, data2_test)
 663   call array_compare_2d(data2_test, data2_expected, test_name // " unit test failed")
 664 end subroutine
 665
 666 ! Status: SKELETAL
 667 ! TODO: More comprehensive interp_1d_3d test
 668 subroutine test_interp_1d_3d
 669   real(k) :: grid1(2,2,4), grid2(2,2,2), data1(2,2,4), data2(2,2,2)
 670
 671   grid1(1,1,:) = [1._k, 2._k, 3._k, 4._k]
 672   grid1(1,2,:) = [5._k, 6._k, 7._k, 8._k]
 673   grid1(2,1,:) = [21._k, 22._k, 23._k, 24._k]
 674   grid1(2,2,:) = [25._k, 26._k, 27._k, 28._k]
 675
 676   grid2(1,1,:) = [2._k, 3._k]
 677   grid2(1,2,:) = [6._k, 7._k]
 678   grid2(2,1,:) = [22._k, 23._k]
 679   grid2(2,2,:) = [26._k, 27._k]
 680
 681   data1(1,1,:) = [101._k, 102._k, 103._k, 104._k]
 682   data1(1,2,:) = [105._k, 106._k, 107._k, 108._k]
 683   data1(2,1,:) = [201._k, 202._k, 203._k, 204._k]
 684   data1(2,2,:) = [205._k, 206._k, 207._k, 208._k]
 685
 686   data2(1,1,:) = [102._k, 103._k]
 687   data2(1,2,:) = [106._k, 107._k]
 688   data2(2,1,:) = [202._k, 203._k]
 689   data2(2,2,:) = [206._k, 207._k]
 690
 691   call interp_1d_3d_assert(grid1, grid2, data1, data2)
 692   call interp_1d_3d_assert(grid1, grid2, data1, data2, "linear")
 693   call interp_1d_3d_assert(grid1, grid2, data1, data2, "cubic_spline")
 694 end subroutine
 695
 696 subroutine interp_1d_3d_assert(grid1, grid2, data1, data2_expected, method, yp1, yp2)
 697   real(k), intent(in), dimension(:,:,:) :: grid1, grid2, data1, data2_expected
 698   character(*), intent(in), optional :: method
 699   real(k), intent(in), optional :: yp1, yp2
 700   real(k) :: data2_test(size(data2_expected,1), size(data2_expected,2), size(data2_expected,3))
 701   integer :: i,i2,i3
 702   character(:), allocatable :: test_name
 703
 704   test_name = "interp_1d_3d(" // &
 705               stringify(grid1) // ", " // &
 706               stringify(grid2) // ", " // &
 707               stringify(data1) // ", data2"
 708
 709   if (present(method)) then
 710     test_name = test_name // ", method=" // method
 711   endif
 712
 713   if (present(yp1)) then
 714     test_name = test_name // ", yp1=" // string(yp1)
 715   endif
 716
 717   if (present(yp2)) then
 718     test_name = test_name // ", yp2=" // string(yp2)
 719   endif
 720
 721   test_name = test_name // ")"
 722
 723   call interp_1d(grid1, grid2, data1, data2_test, method, yp1, yp2)
 724   call array_compare_3d(data2_test, data2_expected, test_name // " unit test failed")
 725 end subroutine
 726
 727 !
 728 ! Supporting utilities
 729 !
 730
 731 subroutine open_netcdf_w(fileobj)
 732   type(FmsNetcdfFile_t), intent(out) :: fileobj
 733
 734   if (.not.open_file(fileobj, "test_axis_utils.nc", "overwrite")) then
 735     call mpp_error(FATAL, "Error opening test_axis_utils.nc to write")
 736   endif
 737 end subroutine
 738
 739 subroutine open_netcdf_r(fileobj)
 740   type(FmsNetcdfFile_t), intent(out) :: fileobj
 741
 742   if (.not.open_file(fileobj, "test_axis_utils.nc", "read")) then
 743     call mpp_error(FATAL, "Error opening test_axis_utils.nc to read")
 744   endif
 745 end subroutine
 746
 747 subroutine array_compare_1d(arr1, arr2, msg)
 748   real(k), intent(in), dimension(:) :: arr1, arr2
 749   character(*), intent(in) :: msg
 750   integer :: i, m, n
 751
 752   m = size(arr1)
 753   n = size(arr2)
 754
 755   if (m.ne.n) then
 756     write(stderr(), "(A)") "1D array comparison failed due to incompatible array sizes"
 757     write(stderr(), "(A)") "Array 1 has size " // string(m) // " and array 2 has size " // string(n)
 758     call mpp_error(FATAL, msg)
 759   endif
 760
 761   do i=1,m
 762     if (arr1(i).ne.arr2(i)) then
 763       write(stderr(), "(A)") "1D array comparison failed due to element " // string(i)
 764       write(stderr(), "(A)") "Array 1 has value " // string(arr1(i)) // &
 765                            & " and array 2 has value " // string(arr2(i))
 766       call mpp_error(FATAL, msg)
 767     endif
 768   enddo
 769 end subroutine
 770
 771 subroutine array_compare_2d(arr1, arr2, msg)
 772   real(k), intent(in), dimension(:,:) :: arr1, arr2
 773   character(*), intent(in) :: msg
 774   integer :: i1, i2, m1, m2, n1, n2
 775
 776   m1 = size(arr1, 1)
 777   m2 = size(arr1, 2)
 778
 779   n1 = size(arr2, 1)
 780   n2 = size(arr2, 2)
 781
 782   if (m1.ne.n1 .or. m2.ne.n2) then
 783     write(stderr(), "(A)") "2D array comparison failed due to incompatible array sizes"
 784     write(stderr(), "(A)") "Array 1 has size " // string(m1) // "x" // string(m2) // &
 785                           & " and array 2 has size " // string(n1) // "x" // string(n2)
 786     call mpp_error(FATAL, msg)
 787   endif
 788
 789   do i2=1,m2
 790     do i1=1,m1
 791       if (arr1(i1,i2).ne.arr2(i1,i2)) then
 792         write(stderr(), "(A)") "2D array comparison failed due to element " // string(i1) // "," // string(i2)
 793         write(stderr(), "(A)") "Array 1 has value " // string(arr1(i1,i2)) // &
 794                              & " and array 2 has value " // string(arr2(i1,i2))
 795         call mpp_error(FATAL, msg)
 796       endif
 797     enddo
 798   enddo
 799 end subroutine
 800
 801 subroutine array_compare_3d(arr1, arr2, msg)
 802   real(k), intent(in), dimension(:,:,:) :: arr1, arr2
 803   character(*), intent(in) :: msg
 804   integer :: i1, i2, i3, m1, m2, m3, n1, n2, n3
 805
 806   m1 = size(arr1, 1)
 807   m2 = size(arr1, 2)
 808   m3 = size(arr1, 3)
 809
 810   n1 = size(arr2, 1)
 811   n2 = size(arr2, 2)
 812   n3 = size(arr2, 3)
 813
 814   if (m1.ne.n1 .or. m2.ne.n2 .or. m3.ne.n3) then
 815     write(stderr(), "(A)") "3D array comparison failed due to incompatible array sizes"
 816     write(stderr(), "(A)") "Array 1 has size " // string(m1) // "x" // string(m2) // "x" // string(m3) // &
 817                            & " and array 2 has size " // string(n1) // "x" // string(n2) // "x" // string(n3)
 818     call mpp_error(FATAL, msg)
 819   endif
 820
 821   do i3=1,m3
 822     do i2=1,m2
 823       do i1=1,m1
 824         if (arr1(i1,i2,i3).ne.arr2(i1,i2,i3)) then
 825           write(stderr(), "(A)") "3D array comparison failed due to element " // &
 826                                & string(i1) // "," // string(i2) // "," // string(i3)
 827           write(stderr(), "(A)") "Array 1 has value " // string(arr1(i1,i2,i3)) // &
 828                                & " and array 2 has value " // string(arr2(i1,i2,i3))
 829           call mpp_error(FATAL, msg)
 830         endif
 831       enddo
 832     enddo
 833   enddo
 834 end subroutine
 835
 836 end program test_axis_utils