Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! module where input data utility classes reside
3 : !
4 : ! classes:
5 : ! time_coordinate -- manages the time coordinate of input data sets
6 : !--------------------------------------------------------------------------------
7 : module input_data_utils
8 : use shr_kind_mod, only : r8 => shr_kind_r8, cs => shr_kind_cs, cl=> shr_kind_cl
9 : use cam_abortutils, only : endrun
10 : use cam_logfile, only : iulog
11 : use pio, only : file_desc_t, pio_inq_dimid, pio_inq_dimlen, pio_get_att
12 : use pio, only : pio_seterrorhandling, pio_get_var, pio_inq_varid
13 : use pio, only : PIO_NOWRITE, PIO_BCAST_ERROR, PIO_INTERNAL_ERROR, PIO_NOERR
14 : use time_manager, only : timemgr_get_calendar_cf, set_time_float_from_date, get_curr_date
15 : use spmd_utils, only : masterproc
16 :
17 : implicit none
18 :
19 : private
20 : public :: time_coordinate
21 :
22 : type :: time_coordinate
23 : integer :: ntimes
24 : real(r8) :: wghts(2)
25 : integer :: indxs(2)
26 : real(r8), allocatable :: times(:)
27 : real(r8), allocatable :: time_bnds(:,:)
28 : logical :: time_interp = .true.
29 : logical :: fixed = .false.
30 : integer :: fixed_ymd, fixed_tod
31 : character(len=cl) :: filename
32 : real(r8) :: dtime ! time shift in interpolation point (days)
33 : contains
34 : procedure :: initialize
35 : procedure :: advance
36 : procedure :: read_more
37 : procedure :: times_check
38 : procedure :: copy
39 : procedure :: destroy
40 : end type time_coordinate
41 :
42 : contains
43 :
44 : !-----------------------------------------------------------------------------
45 : ! initializer
46 : !-----------------------------------------------------------------------------
47 0 : subroutine initialize( this, filepath, fixed, fixed_ymd, fixed_tod, force_time_interp, set_weights, try_dates, delta_days )
48 : use ioFileMod, only : getfil
49 : use cam_pio_utils, only : cam_pio_openfile, cam_pio_closefile
50 : use string_utils, only : to_upper
51 : use shr_const_mod, only : SHR_CONST_CDAY
52 :
53 : class(time_coordinate), intent(inout) :: this
54 : character(len=*), intent(in) :: filepath
55 : logical, optional,intent(in) :: fixed
56 : integer, optional,intent(in) :: fixed_ymd
57 : integer, optional,intent(in) :: fixed_tod
58 : logical, optional,intent(in) :: force_time_interp
59 : logical, optional,intent(in) :: set_weights
60 : logical, optional,intent(in) :: try_dates
61 : real(r8), optional,intent(in) :: delta_days ! time shift in interpolation point (days) -- for previous day set this to -1.
62 :
63 : character(len=cl) :: filen
64 : character(len=cl) :: time_units, err_str
65 : character(len=cs) :: time_calendar, model_calendar
66 : character(len=4) :: yr_str
67 : character(len=2) :: mon_str, day_str, hr_str, min_str, sec_str
68 : integer :: ref_yr, ref_mon, ref_day, ref_hr, ref_min, ref_sec, tod
69 : integer :: varid, ierr
70 : real(r8) :: ref_time
71 :
72 0 : integer, allocatable :: dates(:)
73 0 : integer, allocatable :: datesecs(:)
74 0 : integer, allocatable :: year(:)
75 0 : integer, allocatable :: month(:)
76 0 : integer, allocatable :: day(:)
77 0 : real(r8), allocatable :: ut(:)
78 :
79 0 : real(r8), allocatable :: times_file(:)
80 0 : real(r8), allocatable :: times_modl(:)
81 0 : real(r8), allocatable :: time_bnds_file(:,:)
82 : type(file_desc_t) :: fileid
83 : logical :: force_interp, time_in_secs
84 : logical :: set_wghts
85 : logical :: use_time, adj_times, use_time_bnds
86 : integer :: i, yri, moni, dayi, hri, mini, seci
87 : integer :: err_handling
88 :
89 : character(len=*), parameter :: prefix = 'time_coordinate%initialize: '
90 :
91 0 : if (present(fixed)) this%fixed = fixed
92 0 : if (present(fixed_ymd)) this%fixed_ymd = fixed_ymd
93 0 : if (present(fixed_tod)) this%fixed_tod = fixed_tod
94 :
95 0 : this%dtime = 0._r8
96 0 : if (present(delta_days)) this%dtime = delta_days
97 :
98 0 : if (present(force_time_interp)) then
99 0 : force_interp = force_time_interp
100 : else
101 : force_interp = .false.
102 : endif
103 :
104 0 : if (present(set_weights)) then
105 0 : set_wghts = set_weights
106 : else
107 : set_wghts = .true.
108 : endif
109 :
110 0 : this%filename = trim(filepath)
111 :
112 0 : call getfil( filepath, filen, 0 )
113 0 : call cam_pio_openfile( fileid, filen, PIO_NOWRITE )
114 :
115 0 : call pio_seterrorhandling( fileid, PIO_BCAST_ERROR, oldmethod=err_handling )
116 :
117 0 : call get_dimension( fileid, 'time', this%ntimes )
118 0 : allocate ( times_file( this%ntimes ) )
119 0 : allocate ( times_modl( this%ntimes ) )
120 0 : allocate ( this%times( this%ntimes ) )
121 :
122 0 : ierr = pio_inq_varid( fileid, 'time', varid )
123 0 : use_time = ierr.eq.PIO_NOERR
124 0 : ierr = pio_get_att( fileid, varid, 'calendar', time_calendar)
125 0 : use_time = ierr.eq.PIO_NOERR .and. use_time
126 0 : ierr = pio_get_att( fileid, varid, 'units', time_units)
127 0 : use_time = ierr.eq.PIO_NOERR .and. use_time
128 0 : if (use_time) then
129 : use_time = time_units(1:10) == 'days since' .or. &
130 0 : time_units(1:13) == 'seconds since'
131 : endif
132 :
133 0 : if (present(try_dates)) then
134 0 : if (try_dates) then
135 0 : ierr = pio_inq_varid( fileid, 'date', varid )
136 0 : use_time = ierr .ne. PIO_NOERR
137 : endif
138 : endif
139 :
140 0 : adj_times = .false.
141 0 : use_time_bnds = .false.
142 :
143 0 : time_var_use: if (use_time) then
144 :
145 : ! check the calendar attribute - must match model calendar
146 0 : model_calendar = timemgr_get_calendar_cf()
147 :
148 0 : if (this%ntimes>2) then
149 : ! if only 2 time records then it is assumed that the input has 2 identical time records
150 : ! -- climatological or solar-cycle avaraged
151 :
152 0 : adj_times = (to_upper(time_calendar(1:6)) .ne. to_upper(model_calendar(1:6)))
153 :
154 0 : if (adj_times .and. masterproc) then
155 : write(iulog,*) prefix//'model calendar '//trim(model_calendar)// &
156 0 : ' does not match input data calendar '//trim(time_calendar)
157 0 : write(iulog,*) ' -- will try to use date and datesec in the input file to adjust the time coordinate.'
158 : end if
159 : end if
160 :
161 0 : time_in_secs = .false.
162 :
163 0 : if (index( time_units, 'days since') > 0) then
164 : ! time:units = "days since YYYY-MM-DD hh:mm:ss"
165 : ! 1234567890123456789012345678901234567890
166 : ! 1 2 3
167 : yri = 12
168 : moni = 17
169 : dayi = 20
170 : hri = 23
171 : mini = 26
172 : seci = 29
173 0 : else if(index( time_units, 'seconds since') > 0) then
174 : ! time:units = "seconds since YYYY-MM-DD hh:mm:ss"
175 : ! 1234567890123456789012345678901234567890
176 : ! 1 2 3
177 : yri = 15
178 : moni = 20
179 : dayi = 23
180 : hri = 26
181 : mini = 29
182 : seci = 32
183 : time_in_secs = .true.
184 : else
185 0 : call endrun('time units not recognized')
186 : end if
187 :
188 :
189 : ! parse out ref date and time
190 :
191 0 : yr_str = time_units(yri:yri+3)
192 0 : mon_str = time_units(moni:moni+1)
193 0 : day_str = time_units(dayi:dayi+1)
194 0 : hr_str = time_units(hri:hri+1)
195 0 : min_str = time_units(mini:mini+1)
196 :
197 0 : read( yr_str, * ) ref_yr
198 0 : read( mon_str, * ) ref_mon
199 0 : read( day_str, * ) ref_day
200 0 : read( hr_str, * ) ref_hr
201 0 : read( min_str, * ) ref_min
202 0 : if (len_trim(time_units)>seci) then
203 0 : sec_str = time_units(seci:seci+1)
204 0 : read( sec_str, * ) ref_sec
205 : else
206 0 : ref_sec = 0
207 : endif
208 :
209 0 : tod = ref_hr*3600 + ref_min*60 + ref_sec
210 0 : call set_time_float_from_date( ref_time, ref_yr, ref_mon, ref_day, tod )
211 :
212 0 : ierr = pio_get_var( fileid, varid, times_file )
213 0 : if (time_in_secs) then
214 0 : times_file = times_file/SHR_CONST_CDAY
215 : endif
216 :
217 0 : if (ierr.ne.PIO_NOERR) then
218 0 : call endrun(prefix//'not able to read times')
219 : endif
220 :
221 0 : times_file = times_file + ref_time
222 :
223 0 : ierr = pio_inq_varid( fileid, 'time_bnds', varid )
224 :
225 0 : use_time_bnds = (ierr==PIO_NOERR .and. .not.force_interp)
226 :
227 0 : if (use_time_bnds) then
228 0 : allocate ( this%time_bnds( 2, this%ntimes ) )
229 0 : allocate ( time_bnds_file( 2, this%ntimes ) )
230 0 : ierr = pio_get_var( fileid, varid, time_bnds_file )
231 0 : time_bnds_file = time_bnds_file + ref_time
232 0 : this%time_interp = .false.
233 0 : do i = 1,this%ntimes
234 0 : if (.not. (time_bnds_file(1,i)<times_file(i) &
235 0 : .and. time_bnds_file(2,i)>times_file(i)) ) then
236 0 : write(err_str,*) 'incorrect time_bnds -- time index: ',i,' file: '//trim(filepath)
237 0 : call endrun(err_str)
238 : endif
239 : enddo
240 : else
241 0 : this%time_interp = .true.
242 : endif
243 : else
244 0 : this%time_interp = .true.
245 : endif time_var_use
246 :
247 0 : read_dates: if (adj_times .or. .not.use_time) then
248 :
249 : ! try using date and datesec
250 0 : allocate(dates(this%ntimes), stat=ierr )
251 0 : if( ierr /= 0 ) then
252 0 : write(iulog,*) prefix//'failed to allocate dates; error = ',ierr
253 0 : call endrun(prefix//'failed to allocate dates')
254 : end if
255 :
256 0 : allocate(datesecs(this%ntimes), stat=ierr )
257 0 : if( ierr /= 0 ) then
258 0 : write(iulog,*) prefix//'failed to allocate datesecs; error = ',ierr
259 0 : call endrun(prefix//'failed to allocate datesecs')
260 : end if
261 :
262 0 : ierr = pio_inq_varid( fileid, 'date', varid )
263 0 : if (ierr==PIO_NOERR) then
264 0 : ierr = pio_get_var( fileid, varid, dates )
265 0 : if (ierr/=PIO_NOERR) then
266 0 : call endrun(prefix//' error reading date in '//trim(filepath))
267 : endif
268 : else
269 : ! try year, month, day
270 0 : allocate(year(this%ntimes), stat=ierr )
271 0 : if (ierr/=0) then
272 0 : call endrun(prefix//'issue with allocation of year array')
273 : endif
274 0 : allocate(month(this%ntimes), stat=ierr )
275 0 : if (ierr/=0) then
276 0 : call endrun(prefix//'issue with allocation of month array')
277 : endif
278 0 : allocate(day(this%ntimes), stat=ierr )
279 0 : if (ierr/=0) then
280 0 : call endrun(prefix//'issue with allocation of day array')
281 : endif
282 :
283 0 : ierr = pio_inq_varid( fileid, 'year', varid )
284 0 : if (ierr/=PIO_NOERR) then
285 0 : call endrun(prefix//' error inquiring year var in '//trim(filepath))
286 : endif
287 0 : ierr = pio_get_var( fileid, varid, year )
288 0 : if (ierr/=PIO_NOERR) then
289 0 : call endrun(prefix//' error reading year in '//trim(filepath))
290 : endif
291 :
292 0 : ierr = pio_inq_varid( fileid, 'month', varid )
293 0 : if (ierr/=PIO_NOERR) then
294 0 : call endrun(prefix//' error inquiring month var in '//trim(filepath))
295 : endif
296 0 : ierr = pio_get_var( fileid, varid, month )
297 0 : if (ierr/=PIO_NOERR) then
298 0 : call endrun(prefix//' error reading month in '//trim(filepath))
299 : endif
300 :
301 0 : ierr = pio_inq_varid( fileid, 'day', varid )
302 0 : if (ierr/=PIO_NOERR) then
303 0 : call endrun(prefix//' error inquiring day var in '//trim(filepath))
304 : endif
305 0 : ierr = pio_get_var( fileid, varid, day )
306 0 : if (ierr/=PIO_NOERR) then
307 0 : call endrun(prefix//' error reading day in '//trim(filepath))
308 : endif
309 :
310 0 : dates(:) = year(:)*10000 + month(:)*100 + day(:)
311 :
312 0 : deallocate(year,month,day)
313 : endif
314 0 : ierr = pio_inq_varid( fileid, 'datesec', varid )
315 0 : if (ierr==PIO_NOERR) then
316 0 : ierr = pio_get_var( fileid, varid, datesecs )
317 0 : if (ierr/=PIO_NOERR) then
318 0 : call endrun(prefix//' error reading datesec in '//trim(filepath))
319 : endif
320 : else
321 : ! try ut
322 :
323 0 : allocate(ut(this%ntimes), stat=ierr )
324 0 : ierr = pio_inq_varid( fileid, 'ut', varid ) ! fractional hours
325 0 : if (ierr==PIO_NOERR) then
326 0 : ierr = pio_get_var( fileid, varid, ut )
327 0 : if (ierr/=PIO_NOERR) then
328 0 : call endrun(prefix//' error reading ut in '//trim(filepath))
329 : endif
330 0 : datesecs = int(3600._r8*ut) ! hours -> secs
331 : else
332 0 : datesecs(:) = 0
333 : endif
334 0 : deallocate(ut)
335 : endif
336 :
337 0 : call convert_dates( dates, datesecs, times_modl )
338 :
339 0 : deallocate( dates, datesecs )
340 :
341 : endif read_dates
342 :
343 0 : if (adj_times) then
344 : ! time_bnds_modl - time_bnds_file = times_modl - times_file
345 : ! time_bnds_modl = time_bnds_file + times_modl - times_file
346 0 : this%times(:) = times_modl(:)
347 0 : if (use_time_bnds) then
348 0 : this%time_bnds(1,:) = time_bnds_file(1,:) + times_modl(:) - times_file(:)
349 0 : this%time_bnds(2,:) = time_bnds_file(2,:) + times_modl(:) - times_file(:)
350 : endif
351 0 : else if (use_time) then
352 0 : this%times(:) = times_file(:)
353 0 : if (use_time_bnds) then
354 0 : this%time_bnds(1,:) = time_bnds_file(1,:)
355 0 : this%time_bnds(2,:) = time_bnds_file(2,:)
356 : endif
357 : else
358 0 : this%times(:) = times_modl(:)
359 : endif
360 :
361 0 : deallocate( times_modl, times_file )
362 0 : if (use_time_bnds) deallocate(time_bnds_file)
363 :
364 0 : call pio_seterrorhandling( fileid, err_handling )
365 :
366 0 : call cam_pio_closefile(fileid)
367 :
368 0 : this%indxs(1)=1
369 0 : if (set_wghts) call set_wghts_indices(this)
370 :
371 0 : end subroutine initialize
372 :
373 : !-----------------------------------------------------------------------------
374 : ! advance the time coordinate
375 : !-----------------------------------------------------------------------------
376 0 : subroutine advance( this )
377 : class(time_coordinate) :: this
378 :
379 0 : if (.not.this%fixed) call set_wghts_indices(this)
380 :
381 0 : end subroutine advance
382 :
383 : !-----------------------------------------------------------------------------
384 : ! determine if need to read more data from input data set
385 : !-----------------------------------------------------------------------------
386 0 : function read_more(this) result(check)
387 : class(time_coordinate), intent(in) :: this
388 : logical :: check
389 :
390 : real(r8) :: model_time
391 :
392 0 : model_time = get_model_time() + this%dtime
393 :
394 0 : if (.not.this%fixed) then
395 0 : if (allocated(this%time_bnds)) then
396 0 : check = model_time > this%time_bnds(2,this%indxs(1))
397 : else
398 0 : check = model_time > this%times(this%indxs(2))
399 : endif
400 : else
401 : check = .false.
402 : endif
403 :
404 0 : end function read_more
405 :
406 : !-----------------------------------------------------------------------------
407 : ! times_check -- returns timing status indicator
408 : ! -1 : current model time is before the data times
409 : ! 0 : current model time is within the data times
410 : ! 1 : current model time is after the data times
411 : !-----------------------------------------------------------------------------
412 0 : integer function times_check(this)
413 : class(time_coordinate), intent(in) :: this
414 :
415 : real(r8) :: model_time
416 :
417 0 : model_time = get_model_time()
418 :
419 0 : times_check = 0
420 0 : if (model_time<this%times(1)) then
421 : times_check = -1
422 0 : else if (model_time>this%times(this%ntimes)) then
423 0 : times_check = 1
424 : end if
425 :
426 0 : end function times_check
427 :
428 : !-----------------------------------------------------------------------------
429 : ! destroy method -- deallocate memory and revert to default settings
430 : !-----------------------------------------------------------------------------
431 0 : subroutine destroy( this )
432 : class(time_coordinate), intent(inout) :: this
433 :
434 0 : if (allocated(this%times)) deallocate(this%times)
435 0 : if (allocated(this%time_bnds)) deallocate(this%time_bnds)
436 0 : this%ntimes = 0
437 0 : this%filename='NONE'
438 :
439 0 : end subroutine destroy
440 :
441 : !-----------------------------------------------------------------------------
442 : ! produce a duplicate time coordinate object
443 : !-----------------------------------------------------------------------------
444 0 : subroutine copy( this, obj )
445 : class(time_coordinate), intent(inout) :: this
446 : class(time_coordinate), intent(in) :: obj
447 :
448 0 : call this%destroy()
449 :
450 0 : this%ntimes = obj%ntimes
451 0 : this%fixed = obj%fixed
452 0 : this%fixed_ymd = obj%fixed_ymd
453 0 : this%fixed_tod = obj%fixed_tod
454 :
455 0 : allocate ( this%times( this%ntimes ) )
456 0 : this%times = obj%times
457 :
458 0 : if (allocated( obj%time_bnds )) then
459 0 : allocate ( this%time_bnds( 2, this%ntimes ) )
460 0 : this%time_bnds = obj%time_bnds
461 : endif
462 0 : this%filename = obj%filename
463 :
464 0 : end subroutine copy
465 :
466 : ! private methods
467 :
468 : !-----------------------------------------------------------------------
469 : ! set time interpolation weights
470 : !-----------------------------------------------------------------------
471 0 : subroutine set_wghts_indices(obj)
472 :
473 : class(time_coordinate), intent(inout) :: obj
474 :
475 : real(r8) :: model_time
476 : real(r8) :: datatm, datatp
477 : integer :: yr, mon, day
478 : integer :: index, i
479 : character(len=cl) :: errmsg
480 :
481 : ! set time indices and time-interpolation weights
482 0 : fixed_time: if (obj%fixed) then
483 0 : yr = obj%fixed_ymd/10000
484 0 : mon = (obj%fixed_ymd-yr*10000) / 100
485 0 : day = obj%fixed_ymd-yr*10000-mon*100
486 0 : call set_time_float_from_date( model_time, yr, mon, day, obj%fixed_tod )
487 0 : model_time = model_time + obj%dtime
488 : else
489 0 : model_time = get_model_time() + obj%dtime
490 : endif fixed_time
491 :
492 0 : index = -1
493 :
494 0 : findtimes: do i = obj%indxs(1), obj%ntimes
495 0 : if (allocated(obj%time_bnds)) then
496 0 : datatm = obj%time_bnds(1,i)
497 0 : datatp = obj%time_bnds(2,i)
498 : else
499 0 : if (i .ge. obj%ntimes) then
500 : errmsg = 'input_data_utils::set_wghts_indices cannot not find model time in: '&
501 0 : // trim(obj%filename)
502 0 : write(iulog,*) trim(errmsg)
503 0 : call endrun(trim(errmsg))
504 : endif
505 0 : datatm = obj%times(i)
506 0 : datatp = obj%times(i+1)
507 : endif
508 0 : if ( model_time .lt. datatm ) then
509 : errmsg = 'input_data_utils::set_wghts_indices cannot not find model time in: '&
510 0 : // trim(obj%filename)
511 0 : write(iulog,*) trim(errmsg)
512 0 : call endrun(trim(errmsg))
513 : endif
514 0 : if ( model_time .ge. datatm .and. model_time .le. datatp ) then
515 0 : index = i
516 0 : obj%indxs(1) = i
517 0 : obj%indxs(2) = i+1
518 0 : exit findtimes
519 : endif
520 : enddo findtimes
521 :
522 0 : if ((allocated(obj%time_bnds)) .and. (i<obj%ntimes)) then
523 0 : if (.not.(obj%time_bnds(1,i+1) > obj%time_bnds(1,i))) then
524 0 : obj%indxs = obj%indxs+1 ! skip 29 Feb when calendar is noleap
525 : endif
526 : endif
527 :
528 0 : if (.not.(index>0.and.index<=obj%ntimes)) then
529 : errmsg = 'input_data_utils::set_wghts_indices cannot not find time indices for input file: '&
530 0 : // trim(obj%filename)
531 0 : write(iulog,*) trim(errmsg)
532 0 : call endrun(trim(errmsg))
533 : endif
534 :
535 0 : if (obj%time_interp) then
536 0 : obj%wghts(2) = ( model_time - obj%times(index) ) / ( obj%times(index+1) - obj%times(index) )
537 0 : obj%wghts(1) = 1._r8 - obj%wghts(2)
538 : else
539 0 : obj%wghts(1) = 1._r8
540 0 : obj%wghts(2) = 0._r8
541 : endif
542 :
543 0 : end subroutine set_wghts_indices
544 :
545 : !-----------------------------------------------------------------------
546 : ! returns dimension size
547 : !-----------------------------------------------------------------------
548 0 : subroutine get_dimension( fid, dname, dsize )
549 : type(file_desc_t), intent(in) :: fid
550 : character(*), intent(in) :: dname
551 : integer, intent(out) :: dsize
552 :
553 : integer :: dimid, ierr
554 :
555 0 : ierr = pio_inq_dimid( fid, dname, dimid )
556 0 : ierr = pio_inq_dimlen( fid, dimid, dsize )
557 :
558 0 : end subroutine get_dimension
559 :
560 : !-----------------------------------------------------------------------
561 : ! returns a real which represents the current model time
562 : !-----------------------------------------------------------------------
563 0 : function get_model_time() result(time)
564 :
565 : real(r8) :: time
566 :
567 : integer yr, mon, day, ncsec ! components of a date
568 :
569 0 : call get_curr_date(yr, mon, day, ncsec)
570 :
571 0 : call set_time_float_from_date( time, yr, mon, day, ncsec )
572 :
573 0 : end function get_model_time
574 :
575 : !---------------------------------------------------------------------------
576 : ! convert a collection of dates and times to reals
577 : !---------------------------------------------------------------------------
578 0 : subroutine convert_dates( dates, secs, times )
579 :
580 : use time_manager, only: set_time_float_from_date
581 :
582 : integer, intent(in) :: dates(:)
583 : integer, intent(in) :: secs(:)
584 :
585 : real(r8), intent(out) :: times(:)
586 :
587 : integer :: year, month, day, sec,n ,i
588 :
589 0 : n = size( dates )
590 :
591 0 : do i=1,n
592 0 : year = dates(i)/10000
593 0 : month = (dates(i)-year*10000)/100
594 0 : day = dates(i)-year*10000-month*100
595 0 : sec = secs(i)
596 0 : call set_time_float_from_date( times(i), year, month, day, sec )
597 : enddo
598 :
599 0 : end subroutine convert_dates
600 :
601 0 : end module input_data_utils
|
