FMS  2025.04
Flexible Modeling System
time_interp.inc
1 !***********************************************************************
2 !* Apache License 2.0
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* Licensed under the Apache License, Version 2.0 (the "License");
7 !* you may not use this file except in compliance with the License.
8 !* You may obtain a copy of the License at
9 !*
10 !* http://www.apache.org/licenses/LICENSE-2.0
11 !*
12 !* FMS is distributed in the hope that it will be useful, but WITHOUT
13 !* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
14 !* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 !* PARTICULAR PURPOSE. See the License for the specific language
16 !* governing permissions and limitations under the License.
17 !***********************************************************************
18 
19 !> @addtogroup time_interp_mod
20 !> @{
21 
22  !> @brief Calculates the fractional time into the current year
23  subroutine time_interp_frac_ ( Time, weight )
24 
25  type(time_type), intent(in) :: Time
26  real(FMS_TI_KIND_), intent(out) :: weight !< fractional time
27 
28  integer :: yr, mo, dy, hour, minute, second
29  type(time_type) :: Year_beg, Year_end
30 
31 
32  if ( .not. module_is_initialized ) call time_interp_init
33 
34 ! ---- compute fractional time of year -----
35 
36  call get_date (time, yr, mo, dy, hour, minute, second)
37 
38  year_beg = set_date(yr , 1, 1)
39  year_end = set_date(yr+1, 1, 1)
40 
41  weight = real( (time - year_beg) // (year_end - year_beg) , kind=fms_ti_kind_)
42 
43  end subroutine time_interp_frac_
44 
45 
46  !> @brief Calculates fractional time between mid points of consecutive years
47  subroutine time_interp_year_ ( Time, weight, year1, year2 )
48 
49  type(time_type), intent(in) :: Time
50  real(FMS_TI_KIND_), intent(out) :: weight !< fractional time between midpoints of year1 and year2
51  integer , intent(out) :: year1, year2
52 
53  integer :: yr, mo, dy, hour, minute, second
54  type (time_type) :: Mid_year, Mid_year1, Mid_year2
55 
56 
57  if ( .not. module_is_initialized ) call time_interp_init()
58 
59  call get_date (time, yr, mo, dy, hour, minute, second)
60 
61  ! mid point of current year
62  mid_year = year_midpt(yr)
63 
64  if ( time >= mid_year ) then
65  ! current time is after mid point of current year
66  year1 = yr
67  year2 = yr+1
68  mid_year2 = year_midpt(year2)
69  weight = real( (time - mid_year) // (mid_year2 - mid_year) , kind=fms_ti_kind_ )
70  else
71  ! current time is before mid point of current year
72  year2 = yr
73  year1 = yr-1
74  mid_year1 = year_midpt(year1)
75  weight = real( (time - mid_year1) // (mid_year - mid_year1), kind=fms_ti_kind_ )
76  endif
77 
78  end subroutine time_interp_year_
79 
80  !> @brief Calculates fractional time between mid points of consecutive months
81  subroutine time_interp_month_ ( Time, weight, year1, year2, month1, month2 )
82 
83  type(time_type), intent(in) :: Time
84  real(FMS_TI_KIND_) , intent(out) :: weight
85  integer , intent(out) :: year1, year2, month1, month2
86 
87  integer :: yr, mo, dy, hour, minute, second, &
88  mid_month, cur_month, mid1, mid2
89 
90  if ( .not. module_is_initialized ) call time_interp_init()
91 
92  call get_date (time, yr, mo, dy, hour, minute, second)
93 
94  ! mid point of current month in seconds
95  mid_month = days_in_month(time) * halfday
96  ! time into current month in seconds
97  cur_month = second + secmin*minute + sechour*hour + secday*(dy-1)
98 
99  if ( cur_month >= mid_month ) then
100  ! current time is after mid point of current month
101  year1 = yr; month1 = mo
102  year2 = yr; month2 = mo+1
103  if (month2 > monyear) then
104  year2 = year2+1; month2 = 1
105  endif
106  mid1 = mid_month
107  mid2 = days_in_month(set_date(year2,month2,2)) * halfday
108  weight = real(cur_month - mid1, fms_ti_kind_) / real(mid1+mid2, fms_ti_kind_)
109  else
110  ! current time is before mid point of current month
111  year2 = yr; month2 = mo
112  year1 = yr; month1 = mo-1
113  if (month1 < 1) then
114  year1 = year1-1; month1 = monyear
115  endif
116  if (year1>0) then
117  mid1 = days_in_month(set_date(year1,month1,2)) * halfday
118  else
119  ! this can happen if we are at the beginning of year 1. In this case
120  ! use December 0001 to calculate the duration of December 0000.
121  ! This should work for all calendars
122  mid1 = days_in_month(set_date(1,month1,2)) * halfday
123  endif
124  mid2 = mid_month
125  weight = real(cur_month + mid1, fms_ti_kind_) / real(mid1+mid2, fms_ti_kind_)
126  endif
127 
128  end subroutine time_interp_month_
129 
130  !> @brief Calculates fractional time between mid points of consecutive days
131  subroutine time_interp_day_ ( Time, weight, year1, year2, month1, month2, day1, day2 )
132 
133  type(time_type), intent(in) :: Time
134  real(FMS_TI_KIND_), intent(out) :: weight
135  integer , intent(out) :: year1, year2, month1, month2, day1, day2
136 
137  integer :: yr, mo, dy, hour, minute, second, sday
138 
139  if ( .not. module_is_initialized ) call time_interp_init()
140 
141  call get_date (time, yr, mo, dy, hour, minute, second)
142 
143  ! time into current day in seconds
144  sday = second + secmin*minute + sechour*hour
145 
146  if ( sday >= halfday ) then
147  ! current time is after mid point of day
148  year1 = yr; month1 = mo; day1 = dy
149  year2 = yr; month2 = mo; day2 = dy + 1
150  weight = real(sday - halfday, fms_ti_kind_) / real(secday, fms_ti_kind_)
151 
152  if (day2 > days_in_month(time)) then
153  month2 = month2 + 1
154  day2 = 1
155  if (month2 > monyear) then
156  month2 = 1; year2 = year2+1
157  endif
158  endif
159  else
160  ! current time is before mid point of day
161  year2 = yr; month2 = mo ; day2 = dy
162  year1 = yr; month1 = mo; day1 = dy - 1
163  weight = real(sday + halfday,fms_ti_kind_) / real(secday,fms_ti_kind_)
164 
165  if (day1 < 1) then
166  month1 = month1 - 1
167  if (month1 < 1) then
168  month1 = monyear; year1 = year1-1
169  endif
170  day1 = days_in_month(set_date(year1,month1,2))
171  endif
172  endif
173 
174  end subroutine time_interp_day_
175 
176  !> Part of the time_interp interface, calculates for cyclical data
177  !! Time_beg and Time_end mark a repeating period
178  !!
179  !! Finds mid points and fractional weight for a time perioid
180 subroutine time_interp_modulo_(Time, Time_beg, Time_end, Timelist, weight, index1, index2, &
181  correct_leap_year_inconsistency, err_msg)
182 type(time_type), intent(in) :: Time !< a specific time value
183 type(time_type), intent(in) :: Time_beg !< begining of period to search with
184 type(time_type), intent(in) :: Time_end !< end of period to search with
185 type(time_type), intent(in) :: Timelist(:) !< ascending time values to search between
186 real(FMS_TI_KIND_) , intent(out) :: weight
187 integer , intent(out) :: index1, index2 !< indices of bounding time values within Timelist
188 logical, intent(in), optional :: correct_leap_year_inconsistency!< When true turns on a kluge for an
189  !! inconsistency which may occur in a special case.
190  !! When the modulo time period (i.e. Time_end - Time_beg) is a
191  !! whole number of years and is not a multiple of 4, and the calendar
192  !! in use has leap years, then it is likely that the interpolation
193  !! will involve mapping a common year onto a leap year. In this case
194  !! it is often desirable, but not absolutely necessary, to use data
195  !! for Feb 28 of the leap year when it is mapped onto a common year.
196 character(len=*), intent(out), optional :: err_msg
197 
198  type(time_type) :: Period, T
199  integer :: is, ie,i1,i2
200  integer :: ys,ms,ds,hs,mins,ss ! components of the starting date
201  integer :: ye,me,de,he,mine,se ! components of the ending date
202  integer :: yt,mt,dt,ht,mint,st ! components of the current date
203  integer :: dt1 ! temporary value for day
204  integer :: n ! size of Timelist
205  integer :: stdoutunit
206  logical :: correct_lyr, calendar_has_leap_years, do_the_lyr_correction
207  integer, parameter :: kindl = fms_ti_kind_
208 
209  if ( .not. module_is_initialized ) call time_interp_init
210  if( present(err_msg) ) err_msg = ''
211 
212  stdoutunit = stdout()
213  n = size(timelist)
214 
215  if (time_beg>=time_end) then
216  if(fms_error_handler('time_interp_modulo', &
217  'end of the specified time loop interval must be later than its beginning',err_msg)) return
218  endif
219 
220  calendar_has_leap_years = (get_calendar_type() == julian .or. get_calendar_type() == gregorian)
221 
222  period = time_end-time_beg ! period of the time axis
223 
224  if(present(correct_leap_year_inconsistency)) then
225  correct_lyr = correct_leap_year_inconsistency
226  else
227  correct_lyr = .false.
228  endif
229 
230  ! bring the requested time inside the specified time period
231  t = time
232 
233  do_the_lyr_correction = .false.
234 
235  ! Determine if the leap year correction needs to be done.
236  ! It never needs to be done unless 3 conditions are met:
237  ! 1) We are using a calendar with leap years
238  ! 2) optional argument correct_leap_year_inconsistency is present and equals .true.
239  ! 3) The modulo time period is an integer number of years
240  ! If all of these are true then set do_the_lyr_correction to .true.
241 
242  if(calendar_has_leap_years .and. correct_lyr) then
243  call get_date(time_beg,ys,ms,ds,hs,mins,ss)
244  call get_date(time_end,ye,me,de,he,mine,se)
245  if(ms==me.and.ds==de.and.hs==he.and.mins==mine.and.ss==se) then
246  ! whole number of years
247  do_the_lyr_correction = .true.
248  endif
249  endif
250 
251  if(do_the_lyr_correction) then
252  call get_date(t,yt,mt,dt,ht,mint,st)
253  yt = ys+modulo(yt-ys,ye-ys)
254  dt1 = dt
255  ! If it is Feb 29, but we map into a common year, use Feb 28
256  if(mt==2.and.dt==29.and..not.leap_year(set_date(yt,1,1))) dt1=28
257  t = set_date(yt,mt,dt1,ht,mint,st)
258  if (t < time_beg) then
259  ! the requested time is within the first year,
260  ! but before the starting date. So we shift it to the last year.
261  if(mt==2.and.dt==29.and..not.leap_year(set_date(ye,1,1))) dt=28
262  t = set_date(ye,mt,dt,ht,mint,st)
263  endif
264  else
265  do while ( t >= time_end )
266  t = t-period
267  enddo
268  do while ( t < time_beg )
269  t = t+period
270  enddo
271  endif
272 
273  ! find indices of the first and last records in the Timelist that are within
274  ! the requested time period.
275  if (time_end<=timelist(1).or.time_beg>=timelist(n)) then
276  if(get_calendar_type() == no_calendar) then
277  call print_time(time_beg, 'Time_beg' )
278  call print_time(time_end, 'Time_end' )
279  call print_time(timelist(1), 'Timelist(1)' )
280  call print_time(timelist(n), 'Timelist(n)' )
281  else
282  call print_date(time_beg, 'Time_beg' )
283  call print_date(time_end, 'Time_end' )
284  call print_date(timelist(1), 'Timelist(1)' )
285  call print_date(timelist(n), 'Timelist(n)' )
286  endif
287  write(stdoutunit,*)'where n = size(Timelist) =',n
288  if(fms_error_handler('time_interp_modulo', &
289  'the entire time list is outside the specified time loop interval',err_msg)) return
290  endif
291 
292  call bisect(timelist,time_beg,index1=i1,index2=i2)
293  if (i1 < 1) then
294  is = 1 ! Time_beg before lower boundary
295  else if (time_beg == timelist(i1)) then
296  is = i1 ! Time_beg right on the lower boundary
297  else
298  is = i2 ! Time_beg inside the interval or on upper boundary
299  endif
300  call bisect(timelist,time_end,index1=i1,index2=i2)
301  if (time_end > timelist(i1)) then
302  ie = i1
303  else if (time_end == timelist(i1)) then
304  if(time_beg == timelist(is)) then
305  ! Timelist includes time levels at both the lower and upper ends of the period.
306  ! The endpoints of Timelist specify the same point in the cycle.
307  ! This ambiguity is resolved by ignoring the last time level.
308  ie = i1-1
309  else
310  ie = i1
311  endif
312  else
313 ! This should never happen because bisect does not return i1 such that Time_end < Timelist(i1)
314  endif
315  if (is>=ie) then
316  if(get_calendar_type() == no_calendar) then
317  call print_time(time_beg, 'Time_beg =')
318  call print_time(time_end, 'Time_end =')
319  call print_time(timelist(1), 'Timelist(1)=')
320  call print_time(timelist(n), 'Timelist(n)=')
321  else
322  call print_date(time_beg, 'Time_beg =')
323  call print_date(time_end, 'Time_end =')
324  call print_date(timelist(1), 'Timelist(1)=')
325  call print_date(timelist(n), 'Timelist(n)=')
326  endif
327  write(stdoutunit,*)'where n = size(Timelist) =',n
328  write(stdoutunit,*)'is =',is,'ie =',ie
329  if(fms_error_handler('time_interp_modulo', &
330  'error in calculation of time list bounds within the specified time loop interval',err_msg)) return
331  endif
332 
333  ! handle special cases:
334  if( t>=timelist(ie) ) then
335  ! time is after the end of the portion of the time list within the requested period
336  index1 = ie; index2 = is
337  weight = real((t-timelist(ie))//(period-(timelist(ie)-timelist(is))), fms_ti_kind_ )
338  else if (t<timelist(is)) then
339  ! time is before the beginning of the portion of the time list within the requested period
340  index1 = ie; index2 = is
341  weight = 1.0_kindl - real(((timelist(is)-t)//(period-(timelist(ie)-timelist(is)))), fms_ti_kind_ )
342  else
343  call bisect(timelist,t,index1,index2)
344  weight = real((t-timelist(index1)) // (timelist(index2)-timelist(index1)), fms_ti_kind_ )
345  endif
346 
347 end subroutine time_interp_modulo_
348 
349 subroutine time_interp_list_ ( Time, Timelist, weight, index1, index2, modtime, err_msg )
350 type(time_type) , intent(in) :: Time, Timelist(:)
351 real(FMS_TI_KIND_) , intent(out) :: weight
352 integer , intent(out) :: index1, index2
353 integer, optional, intent(in) :: modtime
354 character(len=*), intent(out), optional :: err_msg
355 
356 integer :: n, hr, mn, se, mtime
357 type(time_type) :: T, Ts, Te, Td, Period, Time_mod
358 character(len=:),allocatable :: terr, tserr, teerr
359 integer, parameter :: kindl = fms_ti_kind_
360 
361  if ( .not. module_is_initialized ) call time_interp_init
362 
363  if( present(err_msg) ) err_msg = ''
364 
365  weight = 0.0_kindl; index1 = 0; index2 = 0
366  n = size(timelist(:))
367 
368 ! setup modular time axis?
369  mtime = none
370  if (present(modtime)) then
371  mtime = modtime
372  time_mod = (timelist(1)+timelist(n))/2
373  call get_date (time_mod, yrmod, momod, dymod, hr, mn, se)
374  mod_leapyear = leap_year(time_mod)
375  endif
376 
377 ! set period for modulo axis
378  select case (mtime)
379  case (none)
380  ! do nothing
381  case (year)
382  period = set_time(0,days_in_year(time_mod))
383  case (month)
384  ! month length must be equal
385  if (days_in_month(time_mod) /= days_in_month(time)) then
386  if(fms_error_handler('time_interp_list','modulo months must have same length',err_msg)) return
387  endif
388  period = set_time(0,days_in_month(time_mod))
389  case (day)
390  period = set_time(0,1)
391  case default
392  if(fms_error_handler('time_interp_list','invalid value for argument modtime',err_msg)) return
393  end select
394 
395 ! If modulo time is in effect and Timelist spans a time interval exactly equal to
396 ! the modulo period, then the endpoints of Timelist specify the same point in the cycle.
397 ! This ambiguity is resolved by ignoring the last time level.
398  if (mtime /= none .and. timelist(size(timelist))-timelist(1) == period) then
399  n = size(timelist) - 1
400  else
401  n = size(timelist)
402  endif
403 
404 ! starting and ending times from list
405  ts = timelist(1)
406  te = timelist(n)
407  td = te-ts
408  t = set_modtime(time,mtime)
409 
410 ! Check that Timelist does not span a time interval greater than the modulo period
411  if (mtime /= none) then
412  if (td > period) then
413  if(fms_error_handler('time_interp_list','period of list exceeds modulo period',err_msg)) return
414  endif
415  endif
416 
417 ! time falls on start or between start and end list values
418  if ( t >= ts .and. t < te ) then
419  call bisect(timelist(1:n),t,index1,index2)
420  weight = real( (t-timelist(index1)) // (timelist(index2)-timelist(index1)), fms_ti_kind_)
421 
422 ! time falls before starting list value
423  else if ( t < ts ) then
424  if (mtime == none) then
425  call time_list_error(t,terr)
426  call time_list_error(ts,tserr)
427  call time_list_error(te,teerr)
428  if(fms_error_handler('time_interp_list',&
429  'time '//trim(terr)//' ('//date_to_string(t)//' is before range of list '//trim(tserr)//'-'//trim(teerr)//&
430  '('//date_to_string(ts)//' - '//date_to_string(te)//')',&
431  err_msg)) return
432  deallocate(terr,tserr,teerr)
433  endif
434  td = te-ts
435  weight = 1.0_kindl - real(((ts-t) // (period-td)), fms_ti_kind_ )
436  index1 = n
437  index2 = 1
438 
439 ! time falls on ending list value
440  else if ( t == te ) then
441  if(perthlike_behavior) then
442  weight = 1.0_kindl
443  index1 = n-1
444  index2 = n
445  else
446  weight = 0.0_kindl
447  index1 = n
448  if (mtime == none) then
449  index2 = n
450  else
451  index2 = 1
452  endif
453  endif
454 
455 ! time falls after ending list value
456  else if ( t > te ) then
457  if (mtime == none) then
458  call time_list_error(t,terr)
459  call time_list_error(ts,tserr)
460  call time_list_error(te,teerr)
461  if(fms_error_handler('time_interp_list',&
462  'time '//trim(terr)//' ('//date_to_string(t)//' is after range of list '//trim(tserr)//'-'//trim(teerr)//&
463  '('//date_to_string(ts)//' - '//date_to_string(te)//')',&
464  err_msg)) return
465  deallocate(terr,tserr,teerr)
466  endif
467  td = te-ts
468  weight = real( (t-te) // (period-td), fms_ti_kind_)
469  index1 = n
470  index2 = 1
471  endif
472 
473 end subroutine time_interp_list_
474 !> }
integer function stdout()
This function returns the current standard fortran unit numbers for output.
Definition: mpp_util.inc:42
subroutine time_interp_year_(Time, weight, year1, year2)
Calculates fractional time between mid points of consecutive years.
Definition: time_interp.inc:48
subroutine time_interp_month_(Time, weight, year1, year2, month1, month2)
Calculates fractional time between mid points of consecutive months.
Definition: time_interp.inc:82
subroutine time_interp_frac_(Time, weight)
Calculates the fractional time into the current year.
Definition: time_interp.inc:24
subroutine time_interp_modulo_(Time, Time_beg, Time_end, Timelist, weight, index1, index2, correct_leap_year_inconsistency, err_msg)
Part of the time_interp interface, calculates for cyclical data Time_beg and Time_end mark a repeatin...
subroutine time_interp_day_(Time, weight, year1, year2, month1, month2, day1, day2)
Calculates fractional time between mid points of consecutive days.