test_fms/time_interp/test_time_interp.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_time_interp
  21
  22  use          fms_mod, only: fms_init, fms_end, stdout, stdlog, FATAL, mpp_error
  23  use time_manager_mod, only: get_date, set_time, set_date, time_manager_init, set_calendar_type, operator(+)
  24  use time_manager_mod, only: JULIAN, time_type, increment_time, NOLEAP, print_date
  25  use  time_interp_mod, only: time_interp_init, time_interp, NONE, YEAR, MONTH, DAY
  26  use time_manager_mod, only: operator(<=), operator(>=), operator(==)
  27  use platform_mod
  28
  29  implicit none
  30
  31  integer, parameter :: num_Time=6, kindl = TI_TEST_KIND_
  32  type(time_type) :: Time_beg, Time_end, Time(num_Time)
  33  type(time_type), allocatable, dimension(:) :: Timelist
  34  integer :: index1, index2, mo, yr, outunit, ntest, nline
  35  real(TI_TEST_KIND_) :: weight
  36  real(TI_TEST_KIND_) :: ref_weights(num_Time), ref_weights_leap(num_Time)
  37  real(TI_TEST_KIND_), parameter :: SMALL = 1.0e-7_kindl ! r4 will fail with 8
  38  real(TI_TEST_KIND_), parameter :: midpoint = 0.483870967741935_kindl
  39  real(TI_TEST_KIND_), parameter :: day_before_leap_day =  0.964285714285714_kindl
  40  real(TI_TEST_KIND_), parameter :: day_before_leap_day_with_ly =  0.931034482758621_kindl
  41
  42  integer :: nmin, nmax
  43
  44  call fms_init
  45  outunit = stdout()
  46  call set_calendar_type(JULIAN)
  47  call time_interp_init
  48
  49  Time_beg = set_date(1, 1, 1)
  50  Time_end = set_date(2, 1, 1)
  51  Time(1) = Time_beg
  52  Time(2) = set_date(1, 1,16)
  53  Time(3) = set_date(1, 2, 1)
  54  Time(4) = set_date(1,12, 1)
  55  Time(5) = set_date(1,12,16)
  56  Time(6) = Time_end
  57
  58  ref_weights(1) = 0.0_kindl ! on 'edge' (timeList value)
  59  ref_weights(2) = midpoint ! rough midpoint of a month ie. jan 16
  60  ref_weights(3) = 0.0_kindl
  61  ref_weights(4) = 0.0_kindl
  62  ref_weights(5) = midpoint
  63  ref_weights(6) = 0.0_kindl
  64
  65  ref_weights_leap(1) = 0.0_kindl ! on 'edge' (timeList value)
  66  ref_weights_leap(2) = day_before_leap_day ! feb 28th
  67  ref_weights_leap(3) = midpoint
  68  ref_weights_leap(4) = 0.0_kindl
  69  ref_weights_leap(5) = day_before_leap_day
  70  ref_weights_leap(6) = day_before_leap_day ! checks that 29th gives same result
  71
  72 ! Tests with modulo time
  73  do nline=1,3
  74
  75    if(nline == 1) then
  76      allocate(Timelist(12))
  77      do mo=1,12
  78        Timelist(mo) = set_date(1, mo, 1)
  79      enddo
  80    else if(nline == 2) then
  81      allocate(Timelist(13))
  82      do mo=1,12
  83        Timelist(mo) = set_date(1, mo, 1)
  84      enddo
  85      Timelist(13) = set_date(2, 1, 1)
  86    else if(nline == 3) then
  87      allocate(Timelist(12))
  88      do mo=2,12
  89        Timelist(mo-1) = set_date(1, mo, 1)
  90      enddo
  91      Timelist(12) = set_date(2, 1, 1)
  92    endif
  93
  94    do ntest=1,num_Time
  95      print *, ntest
  96      call diagram(nline,ntest,modulo_time=.true.)
  97      call time_interp(Time(ntest), Time_beg, Time_end, Timelist, weight, index1, index2)
  98      write(outunit,*) 'time_interp_modulo:'
  99      write(outunit,'()')
 100      call print_date(Time(ntest),                'Time       =')
 101      call print_date(Time_beg,                   'Time_beg   =')
 102      call print_date(Time_end,                   'Time_end   =')
 103      call print_date(Timelist(1),                'Timelist(1)=')
 104      call print_date(Timelist(size(Timelist(:))),'Timelist(n)=')
 105      write(outunit,99) index1,index2,weight
 106      write(outunit,'()')
 107
 108      if(.not. is_valid_indices(index1, index2, Timelist, Time(ntest), weight, YEAR)) &
 109          call mpp_error(FATAL, "test_time_interp: invalid indices from time_interp_timelist")
 110      if(abs(weight - ref_weights(ntest)) .gt. SMALL) &
 111          call mpp_error(FATAL, "test_time_interp: incorrect weight value with reference")
 112
 113      call time_interp(Time(ntest), Timelist, weight, index1, index2, modtime=YEAR)
 114      write(outunit,*) 'time_interp_list with modtime=YEAR:'
 115      write(outunit,'()')
 116      call print_date(Time(ntest),                'Time       =')
 117      call print_date(Timelist(1),                'Timelist(1)=')
 118      call print_date(Timelist(size(Timelist(:))),'Timelist(n)=')
 119      write(outunit,99) index1,index2,weight
 120
 121      if(.not. is_valid_indices(index1, index2, Timelist, Time(ntest), weight, YEAR)) &
 122          call mpp_error(FATAL, "test_time_interp: invalid indices from time_interp_modulo")
 123      if(abs(weight - ref_weights(ntest)) .gt. SMALL) &
 124          call mpp_error(FATAL, "test_time_interp: incorrect weight value with reference")
 125
 126    enddo
 127    deallocate(Timelist)
 128  enddo
 129
 130
 131
 132 ! Tests without modulo time
 133  do nline=1,3
 134    if(nline == 1) then
 135      allocate(Timelist(12))
 136      do mo=1,12
 137        Timelist(mo) = set_date(1, mo, 1)
 138      enddo
 139    else if(nline == 2) then
 140      allocate(Timelist(13))
 141      do mo=1,12
 142        Timelist(mo) = set_date(1, mo, 1)
 143      enddo
 144      Timelist(13) = set_date(2, 1, 1)
 145    else if(nline == 3) then
 146      allocate(Timelist(12))
 147      do mo=2,12
 148        Timelist(mo-1) = set_date(1, mo, 1)
 149      enddo
 150      Timelist(12) = set_date(2, 1, 1)
 151    endif
 152
 153    if(nline == 1) then
 154      nmin = 1; nmax = 4
 155    else if(nline == 2) then
 156      nmin = 1; nmax = num_Time
 157    else if(nline == 3) then
 158      nmin = 3; nmax = num_Time
 159    endif
 160    do ntest=nmin,nmax
 161      call diagram(nline,ntest,modulo_time=.false.)
 162      call time_interp(Time(ntest), Timelist, weight, index1, index2, modtime=NONE)
 163      write(outunit,*) 'time_interp_list with modtime=NONE:'
 164      write(outunit,'()')
 165      call print_date(Time(ntest),                'Time       =')
 166      call print_date(Timelist(1),                'Timelist(1)=')
 167      call print_date(Timelist(size(Timelist(:))),'Timelist(n)=')
 168      write(outunit,99) index1,index2,weight
 169
 170      if( .not. is_valid_indices(index1, index2, TimeList, Time(ntest), weight, NONE)) &
 171         call mpp_error(FATAL, "invalid result without modtime")
 172      if(abs(weight - ref_weights(ntest)) .gt. SMALL) &
 173          call mpp_error(FATAL, "test_time_interp: incorrect weight value with reference")
 174
 175    enddo
 176    deallocate(Timelist)
 177  enddo
 178
 179 ! More tests with modulo time
 180  Time_beg = set_date(1999, 1, 1)
 181  Time_end = set_date(2000, 1, 1)
 182  Time(1)  = set_date(1998, 1, 1)
 183  Time(2)  = set_date(1998, 2,28)
 184  Time(3)  = set_date(1998,12,16)
 185  Time(4)  = set_date(2000, 1, 1)
 186  Time(5)  = set_date(2000, 2,28)
 187  Time(6)  = set_date(2000, 2,29)
 188
 189  allocate(Timelist(13))
 190  do mo=1,12
 191    Timelist(mo) = set_date(1999, mo, 1)
 192  enddo
 193  Timelist(13) = set_date(2000, 1, 1)
 194
 195  write(outunit,'("<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>",/)')
 196  write(outunit,'()')
 197  write(outunit,*) 'time_interp_modulo with correct_leap_year_inconsistency=.true.'
 198  write(outunit,'()')
 199  write(outunit,'(" Jan 1 1999                                     Jan 1 2000")')
 200  write(outunit,'("    |                                               |")')
 201  write(outunit,'("    v                                               v")')
 202  write(outunit,'("    x---x---x---x---x---x---x---x---x---x---x---x---x")')
 203  write(outunit,'("    ^                                               ^")')
 204  write(outunit,'("    |                                               |")')
 205  write(outunit,'(" Time_beg                                        Time_end ")')
 206  write(outunit,'()')
 207
 208  do ntest=1,num_Time
 209    call time_interp(Time(ntest), Time_beg, Time_end, Timelist, weight, index1, index2, &
 210                    &  correct_leap_year_inconsistency=.true.)
 211    call print_date(Time(ntest),' Time =')
 212    write(outunit,99) index1,index2,weight
 213    write(outunit,'()')
 214    if( .not. is_valid_indices(index1, index2, Timelist, Time(ntest), weight, YEAR)) &
 215        call mpp_error(FATAL, 'invalid results for indices with leap year correction')
 216    if(abs(weight - ref_weights_leap(ntest)) .gt. SMALL) &
 217        call mpp_error(FATAL, "test_time_interp: incorrect weight value with reference")
 218  enddo
 219  deallocate(Timelist)
 220
 221  ! swap around ref numbers for different data set
 222  ref_weights_leap(1) = day_before_leap_day
 223  ref_weights_leap(2) = day_before_leap_day ! feb 28th
 224  ref_weights_leap(3) = 0.0_kindl
 225  ref_weights_leap(4) = day_before_leap_day_with_ly
 226  ref_weights_leap(5) = 0.0_kindl
 227  ref_weights_leap(6) = 0.0_kindl
 228 ! Tests of modulo time and leap year inconsistency
 229  Time_beg = set_date(1978, 1, 1)
 230  Time_end = set_date(1981, 1, 1)
 231  Time(1)  = set_date(1976, 2,28)
 232  Time(2)  = set_date(1976, 2,29)
 233  Time(3)  = set_date(1976, 3, 1)
 234  Time(4)  = set_date(1983, 2,28)
 235  Time(5)  = set_date(1983, 3, 1)
 236  Time(6)  = set_date(1981, 1, 1)
 237  allocate(Timelist(37))
 238  do yr=1978,1980
 239    do mo=1,12
 240      Timelist(12*(yr-1978)+mo) = set_date(yr, mo, 1)
 241    enddo
 242  enddo
 243  Timelist(37) = set_date(1981, 1, 1)
 244
 245  write(outunit,'("<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>")')
 246  write(outunit,'()')
 247  write(outunit,*) 'time_interp_modulo with correct_leap_year_inconsistency=.true.'
 248  write(outunit,'()')
 249  write(outunit,'(" Jan 1 1978              Jan 1 1979              Jan 1 1980              Jan 1 1981")')
 250  write(outunit,'("     |                       |                       | <---- leap year ----> |")')
 251  write(outunit,'("     v                       v                       v                       v")')
 252  write(outunit,'("     x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x")')
 253  write(outunit,'("     ^                                                                       ^")')
 254  write(outunit,'("     |                                                                       |")')
 255  write(outunit,'("  Time_beg                                                               Time_end")')
 256  write(outunit,'()')
 257
 258  do ntest=1,num_Time
 259    call time_interp(Time(ntest), Time_beg, Time_end, Timelist, weight, index1, index2, &
 260                    &  correct_leap_year_inconsistency=.true.)
 261    call print_date(Time(ntest),' Time=')
 262    write(outunit,99) index1,index2,weight
 263    write(outunit,'()')
 264    if( .not. is_valid_indices(index1, index2, Timelist, Time(ntest), weight, YEAR)) &
 265        call mpp_error(FATAL, 'invalid results for indices with leap year correction')
 266    if(abs(weight - ref_weights_leap(ntest)) .gt. SMALL) &
 267        call mpp_error(FATAL, "test_time_interp: incorrect weight value with reference")
 268  enddo
 269  deallocate(Timelist)
 270
 271  allocate(Timelist(12))
 272  Timelist( 1) = set_date(1,  1, 16, hour=12) ! Jan midmonth
 273  Timelist( 2) = set_date(1,  2, 15, hour= 0) ! Feb midmonth (common year)
 274  Timelist( 3) = set_date(1,  3, 16, hour=12) ! Mar midmonth
 275  Timelist( 4) = set_date(1,  4, 16, hour= 0) ! Apr midmonth
 276  Timelist( 5) = set_date(1,  5, 16, hour=12) ! May midmonth
 277  Timelist( 6) = set_date(1,  6, 16, hour= 0) ! Jun midmonth
 278  Timelist( 7) = set_date(1,  7, 16, hour=12) ! Jul midmonth
 279  Timelist( 8) = set_date(1,  8, 16, hour=12) ! Aug midmonth
 280  Timelist( 9) = set_date(1,  9, 16, hour= 0) ! Sep midmonth
 281  Timelist(10) = set_date(1, 10, 16, hour=12) ! Oct midmonth
 282  Timelist(11) = set_date(1, 11, 16, hour= 0) ! Nov midmonth
 283  Timelist(12) = set_date(1, 12, 16, hour=12) ! Dec midmonth
 284  Time_beg = set_date(1, 1, 1)
 285  Time_end = set_date(2, 1, 1)
 286  call diagram(nline=4, ntest=0, modulo_time=.true.)
 287  do ntest=0,73
 288    Time(1) = set_date(1996, 1, 1) + set_time(seconds=0, days=5*ntest)
 289    call print_date(Time(1),' Time=')
 290    call time_interp(Time(1), Timelist, weight, index1, index2, modtime=YEAR)
 291    write(outunit,89) 'time_interp_list with modtime=YEAR:   ', index1,index2,weight
 292    call time_interp(Time(1), Time_beg, Time_end, Timelist, weight, index1, index2, &
 293                    &  correct_leap_year_inconsistency=.true.)
 294    write(outunit,89) 'time_interp_modulo: ', index1,index2,weight
 295    write(outunit,'()')
 296  enddo
 297
 298  99 format(' index1=',i3,'  index2=',i3,'  weight=',f18.15)
 299  89 format(a20,' index1=',i3,'  index2=',i3,'  weight=',f18.15)
 300  call fms_end
 301
 302  contains
 303
 304  subroutine diagram(nline,ntest,modulo_time)
 305  integer, intent(in) :: nline,ntest
 306  logical, intent(in) :: modulo_time
 307
 308  write(outunit,'("<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>")')
 309  write(outunit,'()')
 310  if(modulo_time) then
 311    write(outunit,'(" Time_beg                                      Time_end")')
 312    write(outunit,'("  |                                               |")')
 313    write(outunit,'("  v                                               v")')
 314  endif
 315
 316  if(nline == 1) then
 317    write(outunit,'("  x---x---x---x---x---x---x---x---x---x---x---x----")')
 318  else if(nline == 2) then
 319    write(outunit,'("  x---x---x---x---x---x---x---x---x---x---x---x---x")')
 320  else if(nline == 3) then
 321    write(outunit,'("  ----x---x---x---x---x---x---x---x---x---x---x---x")')
 322  else if(nline == 4) then
 323    write(outunit,'("  --x---x---x---x---x---x---x---x---x---x---x---x--")')
 324  endif
 325
 326  if(ntest == 1) then
 327    write(outunit,'("  ^")  ')
 328    write(outunit,'("  |")  ')
 329    write(outunit,'(" Time")')
 330  else if(ntest == 2) then
 331    write(outunit,'("    ^")  ')
 332    write(outunit,'("    |")  ')
 333    write(outunit,'("   Time")')
 334  else if(ntest == 3) then
 335    write(outunit,'("      ^")  ')
 336    write(outunit,'("      |")  ')
 337    write(outunit,'("     Time")')
 338  else if(ntest == 4) then
 339    write(outunit,'("                                              ^")  ')
 340    write(outunit,'("                                              |")  ')
 341    write(outunit,'("                                             Time")')
 342  else if(ntest == 5) then
 343    write(outunit,'("                                                ^")  ')
 344    write(outunit,'("                                                |")  ')
 345    write(outunit,'("                                               Time")')
 346  else if(ntest == 6) then
 347    write(outunit,'("                                                  ^")  ')
 348    write(outunit,'("                                                  |")  ')
 349    write(outunit,'("                                                 Time")')
 350  endif
 351  write(outunit,'()')
 352
 353  end subroutine diagram
 354
 355    !> helper function to check results
 356    !! true if invalid , false for valid
 357    logical function is_valid_indices(ind1, ind2, tList, tintv, res_weight, mtime)
 358        integer, intent(in) :: ind1, ind2
 359        type(time_type), intent(in) :: tList(:), tintv
 360        real(TI_TEST_KIND_), intent(in) :: res_weight
 361        integer, intent(in) :: mtime
 362        integer :: i
 363
 364         ! modulo_time determines wrap around
 365         if( mtime .eq. NONE) then
 366            if (ind1 .eq. SIZE(tList)) then
 367                 is_valid_indices = ind2 .eq. ind1
 368            else
 369                 is_valid_indices = ind2 .eq. ind1+1
 370            endif
 371         else ! YEAR, default
 372            if (ind1 .eq. 12 ) then
 373                 is_valid_indices = ind2 .eq. 1
 374            else
 375                 is_valid_indices = ind2 .eq. ind1+1
 376            endif
 377         endif
 378
 379     end function is_valid_indices
 380
 381  end program test_time_interp