test_fms/interpolator/test_interpolator2.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 !> @file
  20 !! @brief unit tests for interpolator_mod
  21 !! @author MiKyung Lee
  22 !! @email gfdl.climate.model.info@noaa.gov
  23 !! @description This program tests the various subroutines in interpolator_mod.  The code in
  24 !! in test_interpolator.F90 shows how to use interpolator_mod.  This program contains the
  25 !! actual testing.
  26 !! TODO:
  27
  28 program test_interpolator2
  29
  30   use fms2_io_mod, only: FmsNetcdfFile_t, UNLIMITED,                  &
  31                          register_field, register_variable_attribute, &
  32                          register_axis,                               &
  33                          write_data, close_file, open_file
  34   use mpp_mod,          only: mpp_error, FATAL, WARNING
  35   use time_manager_mod, only: time_type, set_calendar_type, time_manager_init, &
  36                               get_date_no_leap, get_date_julian, set_date, set_time, &
  37                               set_date_no_leap, set_date_julian, operator(/), &
  38                               operator(+), operator(-), time_type_to_real, increment_time, &
  39                               leap_year, days_in_month, print_date, print_time
  40   use fms_mod,          only: fms_init
  41   use constants_mod,    only: PI
  42   use platform_mod,     only: r4_kind, r8_kind
  43
  44   use interpolator_mod
  45
  46   implicit none
  47
  48   character(100), parameter :: ncfile='immadeup.o3.climatology.nc' !< fake climatology file
  49   integer, parameter :: lkind=TEST_INTP_KIND_
  50   !> the interpolation methods are not perfect.Will not get perfectly agreeing answers
  51   real(r8_kind), parameter :: tol=0.1_lkind
  52   integer :: calendar_type
  53
  54   !> climatology related variables and arrays (made up data)
  55   integer :: nlonlat       !< number of latitude and longitudinal center coordinates in file
  56   integer :: nlonlatb      !< number of latitude and longitudinal boundary coordinates in file
  57   integer :: ntime         !< number of time slices
  58   integer :: npfull        !< number of p levels
  59   integer :: nphalf        !< number of half p levels
  60   real(TEST_INTP_KIND_), allocatable :: lat(:)   !< climatology coordinates
  61   real(TEST_INTP_KIND_), allocatable :: lon(:)   !< climatology coordinates
  62   real(TEST_INTP_KIND_), allocatable :: latb(:)  !< climatology coordinates
  63   real(TEST_INTP_KIND_), allocatable :: lonb(:)  !< climatology coordinates
  64   real(r8_kind), allocatable :: clim_time (:) !< climatology time
  65   real(TEST_INTP_KIND_), allocatable :: pfull(:) !< climatology p level
  66   real(TEST_INTP_KIND_), allocatable :: phalf(:) !< climatology p half level
  67   real(TEST_INTP_KIND_), allocatable :: ozone(:,:,:,:) !< climatology ozone data
  68
  69   !> model related variables and arrays
  70   integer :: nlonlat_mod, nlonlatb_mod !< number of latitude and longitude coordinates in the model
  71   real(TEST_INTP_KIND_), allocatable :: lat_mod(:,:)  !< model coordinates
  72   real(TEST_INTP_KIND_), allocatable :: lon_mod(:,:)  !< model coordinates
  73   real(TEST_INTP_KIND_), allocatable :: latb_mod(:,:) !< model coordinates
  74   real(TEST_INTP_KIND_), allocatable :: lonb_mod(:,:) !< model coordinates
  75
  76   !> array holding model times
  77   type(time_type), allocatable :: model_time_julian(:), model_time_noleap(:)
  78
  79   type(interpolate_type) :: o3 !< recyclable interpolate_type
  80
  81   !> whether the file input is yearly, daily data
  82   logical :: yearly, daily
  83   logical :: noleap
  84
  85   logical :: test_file_daily_julian=.true.,   test_file_daily_noleap=.false.
  86   logical :: test_file_yearly_noleap=.false., test_file_yearly_julian=.false.
  87   logical :: test_file_no_time=.false.
  88   integer :: nml_unit_var=99
  89   character(*), parameter :: nml_file='test_interpolator.nml'
  90   NAMELIST / test_interpolator_nml / test_file_daily_noleap, test_file_daily_julian, &
  91                                      test_file_yearly_noleap, test_file_yearly_julian, test_file_no_time
  92
  93   open(unit=nml_unit_var, file=nml_file)
  94   read(unit=nml_unit_var, nml=test_interpolator_nml)
  95   close(nml_unit_var)
  96
  97   call fms_init
  98   call time_manager_init
  99   call write_header
 100
 101   !> set data
 102   call set_parameters_wrapper
 103   call set_and_write_data
 104
 105   if(.not.test_file_no_time) then
 106      !> test interpolator when  model calendar is JULIAN
 107      calendar_type=2
 108      call set_calendar_type(calendar_type)
 109      call run_test_set
 110      !---------------------------------------------------------
 111      !> test interpolator when model calendar is NOLEAP
 112      calendar_type=4
 113      call set_calendar_type(calendar_type)
 114      call run_test_set
 115   end if
 116
 117   !> test interpolator_no_time_axis
 118   if(test_file_no_time) then
 119      write(*,*) 'test_intepolator_no_time_axis'
 120      calendar_type=2  !< still need to set model calendar
 121      call set_calendar_type(calendar_type)
 122      call test_interpolator_init(o3)
 123      call test_interpolator_no_time_axis(o3)
 124      call test_interpolator_end(o3)
 125   end if
 126
 127 contains
 128
 129 #include "test_interpolator_write_climatology.inc"
 130
 131   !===============================================!
 132   subroutine test_interpolator_init(clim_type)
 133
 134     !> interopolator_init initializes the interpolate_type after reading in the
 135     !! climatology data.  The required arguments are clim_type, file_name, lonb_mod, latb_mod
 136     !! where lonb_mod and latb_mod contain the model longitude and latitude values on the grid.
 137
 138     implicit none
 139     type(interpolate_type), intent(inout) :: clim_type
 140     integer, dimension(1) :: data_out_of_bounds
 141
 142     data_out_of_bounds(1)=CONSTANT
 143     call interpolator_init(clim_type,trim(ncfile), lonb_mod, latb_mod, data_out_of_bounds=data_out_of_bounds)
 144
 145   end subroutine test_interpolator_init
 146   !===============================================!
 147   subroutine test_interpolator(clim_type, model_time)
 148
 149     !> call the variants of interpolator (4D-2d) that interpolates data at a given time-point
 150     !! The tests here do not test the "no_axis" interpolator routines
 151     !! This subroutine also tests obtain_interpolator_time_slices for the 2D case.
 152
 153     implicit none
 154
 155     type(interpolate_type), intent(inout) :: clim_type
 156     type(time_type), dimension(ntime), intent(in) :: model_time
 157     type(time_type) :: test_time
 158     real(TEST_INTP_KIND_), dimension(nlonlat_mod,nlonlat_mod,npfull,1) :: interp_data !<only 1 field
 159     real(TEST_INTP_KIND_), dimension(nlonlat_mod,nlonlat_mod,nphalf) :: phalf_in
 160     integer :: itime, i, j, k, l
 161
 162     real(TEST_INTP_KIND_) :: answer
 163
 164
 165     do i=1, nphalf
 166        phalf_in(:,:,i)=phalf(i)
 167     end do
 168
 169     do itime=2, ntime-1
 170
 171        answer=0.5_lkind*ozone(1,1,1,itime-1)+0.5_lkind*ozone(1,1,1,itime)
 172        test_time=model_time(itime-1) + (model_time(itime)-model_time(itime-1))/2
 173
 174        !> test interpolator_4D_r4/8
 175        call interpolator(clim_type, test_time, phalf_in, interp_data, 'ozone')
 176        do i=1, npfull
 177           do j=1, nlonlat_mod
 178              do k=1, nlonlat_mod
 179                 call check_answers(interp_data(k,j,i,1), answer, tol, 'test interpolator_4D')
 180              end do
 181           end do
 182        end do
 183
 184        !> test interpolator_3_r4/8
 185        call interpolator(clim_type, test_time, phalf_in, interp_data(:,:,:,1), 'ozone')
 186        do i=1, npfull
 187           do j=1, nlonlat_mod
 188              do k=1, nlonlat_mod
 189                 call check_answers(interp_data(k,j,i,1), answer, tol, 'test interpolator_3D')
 190              end do
 191           end do
 192        end do
 193
 194        !> test interpolator_2D_r4/8
 195        call interpolator(clim_type, test_time, interp_data(:,:,1,1), 'ozone')
 196        do j=1, nlonlat_mod
 197           do k=1, nlonlat_mod
 198              call check_answers(interp_data(k,j,1,1), answer, tol, 'test interpolator_2D')
 199           end do
 200        end do
 201
 202        !> Test obtain_interpolator_time_slices
 203        call obtain_interpolator_time_slices(clim_type,test_time)
 204        call interpolator(clim_type, test_time, interp_data(:,:,1,1), 'ozone')
 205        call unset_interpolator_time_flag(clim_type)
 206        do j=1, nlonlat_mod
 207           do k=1, nlonlat_mod
 208              call check_answers(interp_data(k,j,1,1), answer, tol, 'test interpolator_2D')
 209           end do
 210        end do
 211
 212     end do
 213
 214   end subroutine test_interpolator
 215   !===============================================!
 216   subroutine test_interpolator_end(clim_type)
 217
 218     !> This subroutine tests interpolator_end
 219
 220     implicit none
 221
 222     type(interpolate_type) :: clim_type
 223
 224     call interpolator_end(clim_type)
 225
 226   end subroutine test_interpolator_end
 227   !===============================================!
 228   subroutine test_interpolator_no_time_axis(clim_type)
 229
 230     !> This subroutine tests the variants (42-2D) of interpolator_no_time_axis
 231
 232     implicit none
 233
 234     type(interpolate_type) :: clim_type
 235
 236     real(TEST_INTP_KIND_), dimension(nlonlat,nlonlat,npfull,1) :: interp_data !< last column, there is only one field
 237     real(TEST_INTP_KIND_), dimension(nlonlat,nlonlat,nphalf) :: phalf_in
 238     integer :: i, j, k
 239
 240     do i=1, nphalf
 241        phalf_in(:,:,i)=phalf(i)
 242     end do
 243
 244     !> test interpolator_4D_no_time_axis_r4/8
 245     call interpolator(clim_type, phalf_in, interp_data, 'ozone')
 246     do i=1, npfull
 247        do j=1, nlonlat
 248           do k=1, nlonlat
 249              call check_answers(interp_data(k,j,i,1), ozone(k,j,i,1), tol, 'test interpolator_4D_no_time_axis')
 250           end do
 251        end do
 252     end do
 253
 254     !> test interpolator_3D_no_time_axis_r4/8
 255     call interpolator(clim_type, phalf_in, interp_data(:,:,:,1), 'ozone')
 256     do i=1, npfull
 257        do j=1, nlonlat
 258           do k=1, nlonlat
 259              call check_answers(interp_data(k,j,i,1), ozone(k,j,i,1), tol, 'test interpolator_3D_no_time_axis')
 260           end do
 261        end do
 262     end do
 263
 264     !> test interpolator_2D_no_time_axis_r4/8
 265     call interpolator(clim_type, interp_data(:,:,1,1), 'ozone')
 266     do j=1, nlonlat
 267        do k=1, nlonlat
 268           call check_answers(interp_data(k,j,1,1), ozone(k,j,1,1), tol, 'test interpolator_2D_no_time_axis')
 269        end do
 270     end do
 271
 272   end subroutine test_interpolator_no_time_axis
 273   !===============================================!
 274   subroutine test_interpolate_type_eq
 275
 276     !> This subroutine tests interpolaote_type_eq (assignment = operator)
 277     !! The success of "=" is insured by checking to see if interpolation with o3_copy succeeds.
 278
 279     implicit none
 280
 281     type(interpolate_type) :: o3_copy
 282
 283     o3_copy = o3
 284     if(calendar_type==2) call test_interpolator(o3_copy, model_time_julian)
 285     if(calendar_type==4) call test_interpolator(o3_copy, model_time_noleap)
 286
 287   end subroutine test_interpolate_type_eq
 288   !===============================================!
 289   subroutine test_query_interpolator
 290
 291     !> This subroutne tests query_interpolator
 292
 293     implicit none
 294
 295     character(100), parameter :: answer_field_name='ozone'
 296     integer, parameter :: answer_nfields=1
 297
 298     character(100), allocatable :: field_names(:)
 299     integer :: nfields
 300
 301     call query_interpolator(o3,nfields)
 302     if( nfields .ne. answer_nfields) call mpp_error(FATAL, '')
 303
 304     allocate(field_names(nfields))
 305     call query_interpolator(o3,nfields,field_names)
 306     if(trim(answer_field_name).ne.trim(field_names(1))) call mpp_error(FATAL,'')
 307
 308     deallocate(field_names)
 309
 310   end subroutine test_query_interpolator
 311   !===============================================!
 312   subroutine run_test_set
 313
 314     implicit none
 315     integer :: i
 316
 317     if(calendar_type==2) write(*,*) "** MODEL CALENDAR JULIAN ** MODEL CALENDAR JULIAN ** MODEL CALENDAR JULIAN **"
 318     if(calendar_type==4) write(*,*) "** MODEL CALENDAR NOLEAP ** MODEL CALENDAR NOLEAP ** MODEL CALENDAR NOLEAP **"
 319
 320     write(*,*) '1.  interpolator_init'
 321     call test_interpolator_init(o3)
 322
 323     !> test interpolator 2D-4D
 324     write(*,*) '2.  interpolator4D-2D'
 325     if(calendar_type==2) call test_interpolator(o3, model_time_julian)
 326     if(calendar_type==4) call test_interpolator(o3, model_time_noleap)
 327
 328     !> test interpolate_type_eq
 329     !! This test has been commented out and will be included in the testing suite once fileobj cp is added into
 330     write(*,*) '3.  skipping interpolate_type_eq until bug in interpolator is fixed'
 331     !call test_interpolate_type_eq()
 332
 333     !> test query_interpolator
 334     write(*,*) '4.  query_interpolator'
 335     call test_query_interpolator()
 336
 337     !> test interpolator end
 338     write(*,*) '5.  interpolator_end'
 339     call test_interpolator_end(o3)
 340
 341   end subroutine run_test_set
 342   !===============================================!
 343   subroutine check_answers(results, answers, tol, whoami)
 344
 345     implicit none
 346     real(TEST_INTP_KIND_), intent(in) :: results, answers
 347     real(r8_kind), intent(in) :: tol
 348     character(*) :: whoami
 349
 350     if (real(abs(results-answers),r8_kind).gt.tol) then
 351     !if (results.ne.answers) then
 352        write(*,*) '      EXPECTED ', answers, ' but computed ', results
 353        call mpp_error(FATAL, trim(whoami))
 354     end if
 355
 356   end subroutine check_answers
 357   !===============================================!
 358   subroutine write_header
 359
 360
 361     if(test_file_daily_noleap) &
 362          write(*,"(////10x,a,i0////)") &
 363          " ** DAILY FILE CAL NOLEAP ** DAILY FILE CAL NOLEAP ** DAILY FILE CAL NOLEAP ** ", lkind
 364     if(test_file_daily_julian)  &
 365          write(*,"(////10x,a,i0/////)") &
 366          ' ** DAILY FILE CAL JULIAN ** DAILY FILE CAL JULIAN ** DAILY FILE CAL JULIAN ** ', lkind
 367     if(test_file_yearly_noleap) &
 368          write(*,"(////10x,a,i0/////)") &
 369          ' ** YEARLY FILE CAL NOLEAP ** YEARLY FILE CAL NOLEAP ** YEARLY FILE CAL NOLEAP ** ', lkind
 370     if(test_file_yearly_julian) &
 371          write(*,"(////10x,a,i0/////)") &
 372          ' ** YEARLY FILE CAL JULIAN ** YEARLY FILE CAL JULIAN ** YEARLY FILE CAL JULIAN ** ', lkind
 373     if(test_file_no_time) &
 374          write(*,"(////10x,a/////)") " ** NO TIME AXIS ** NO TIME AXIS ** NO TIME AXIS **"
 375
 376   end subroutine write_header
 377   !===============================================!
 378   subroutine set_parameters_wrapper
 379
 380     if(test_file_daily_noleap)  then
 381        call set_parameters(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=240, npfull_in=3, &
 382             daily_in=.true., yearly_in=.false., noleap_in=.true.)
 383     else if(test_file_daily_julian) then
 384        call set_parameters(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=240, npfull_in=3, &
 385             daily_in=.true., yearly_in=.false., noleap_in=.false.)
 386     else if(test_file_yearly_noleap) then
 387        call set_parameters(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=240, npfull_in=3, &
 388             daily_in=.false., yearly_in=.true., noleap_in=.true.)
 389     else if(test_file_yearly_julian) then
 390        call set_parameters(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=240, npfull_in=3, &
 391             daily_in=.false., yearly_in=.true., noleap_in=.false.)
 392     else if(test_file_no_time) then
 393        call set_parameters(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=0, npfull_in=3, &
 394             daily_in=.true., yearly_in=.false., noleap_in=.false.)
 395     end if
 396
 397   end subroutine set_parameters_wrapper
 398   !===============================================!
 399 end program test_interpolator2