1module timeCoord_mod
2 ! MODULE timeCoord_mod (prefix='tim' category='7. Low-level data objects')
3 !
4 !:Purpose: To store public variables and procedures related to the time
5 ! coordinate.
6 !
7 use midasMpi_mod
8 use varNameList_mod
9 use utilities_mod
10 implicit none
11 save
12 private
13
14 ! public variables
15 public :: tim_dstepobs, tim_dstepobsinc, tim_windowsize
16 public :: tim_nstepobs, tim_nstepobsinc, tim_referencetime
17 public :: tim_fullyUseExtremeTimeBins
18 ! public procedures
19 public :: tim_setup, tim_initialized
20 public :: tim_getDateStamp, tim_setDateStamp, tim_getStampList, tim_getStepObsIndex
21 public :: tim_getDateStampFromFile, tim_dateStampToYYYYMMDDHH, tim_getValidDateTimeFromList
22
23 character(len=4) :: varNameForDate
24 character(len=6) :: tim_referencetime
25 real(8) :: tim_dstepobs
26 real(8) :: tim_dstepobsinc
27 real(8) :: tim_windowsize
28 integer :: tim_nstepobs
29 integer :: tim_nstepobsinc
30 logical :: tim_fullyUseExtremeTimeBins
31 integer :: datestamp = 0 ! datestamp is usually the centre of time window
32 logical :: initialized = .false.
33
34 integer, external :: get_max_rss
35
36contains
37
38 subroutine tim_readNml()
39 !
40 !:Purpose: Read the namelist block NAMTIME.
41 !
42 !:Namelist parameters:
43 ! :dstepobs: time step (hrs) between successive trial fields
44 ! for use in OmP determation. Set to dwindowsize for
45 ! single trial field, i.e. use of 3dvar instead of 3dvar-FGAT.
46 ! nstepobs = number of trial fields
47 ! :dstepobsinc: time step (hrs) between obs groupings in time. Set to
48 ! dwindowsize for use of a single obs group.
49 ! nstepobsinc = number of obs time intervals
50 ! :dwindowsize: Time window size (hrs).
51 !
52 !:Comment:
53 !
54 ! Provided dates and number of provided trial field files must be
55 ! consistent with nstepobs, dstepobs and dwindowsize with reference datestamp
56 ! corresponding to the date of the middle trial field file.
57 !
58 implicit none
59
60 ! Locals:
61 integer :: nulnam, ierr, fnom, fclos
62 logical, save :: firstCall = .true.
63
64 ! Namelist variables:
65 real(8) :: dstepobs ! time step length for background state (in hours)
66 real(8) :: dstepobsinc ! time step length for increment and/or B matrix (in hours)
67 real(8) :: dwindowsize ! length of assimilation window (in hours)
68 character(len=6) :: referencetime ! location of 'date' within the window: 'middle' or 'start'
69 logical :: fullyUseExtremeTimeBins ! choose to use full-size bins at both ends of window (usually only half size)
70
71 NAMELIST /NAMTIME/dstepobs, dstepobsinc, dwindowsize, referencetime, fullyUseExtremeTimeBins
72
73 if (.not.firstCall) then
74 write(*,*) 'tim_readNml: already initialized, just return'
75 return
76 else
77 firstCall = .false.
78 end if
79
80 ! Set default values for namelist variables
81 dstepobs = 6.0d0
82 dstepobsinc = 6.0d0
83 dwindowsize = 6.0d0
84 referenceTime = 'middle'
85 fullyUseExtremeTimeBins = .false.
86
87 ! Read the namelist
88 nulnam = 0
89 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
90 read(nulnam, nml=namtime, iostat=ierr)
91 if (ierr /= 0) call utl_abort('tim_readNml: Error reading namelist')
92 if (mmpi_myid == 0) write(*,nml=namtime)
93 ierr = fclos(nulnam)
94
95 ! Set the module variables for timestep length, number of timesteps and window length
96 tim_dstepobs = dstepobs
97 tim_dstepobsinc = dstepobsinc
98 tim_windowsize = dwindowsize
99 tim_referencetime = referenceTime
100 tim_fullyUseExtremeTimeBins = fullyUseExtremeTimeBins
101
102 if (tim_fullyUseExtremeTimeBins .and. tim_referencetime=="middle") then
103 write(*,*) "Warning: tim_fullyUseExtremeTimeBins==.true. and tim_referencetime=='middle' is a non-standard combination"
104 write(*,*) "Is it really what you want ?"
105 end if
106
107 if (dstepobs > dwindowsize) then
108 if (mmpi_myid == 0) write(*,*) 'tim_readNml: dstepobs>dwindowsize. Reset to dwindowsize value.'
109 tim_dstepobs = tim_windowsize
110 end if
111 if (dstepobsinc > dwindowsize) then
112 if (mmpi_myid == 0) write(*,*) 'tim_readNml: dstepobsinc>dwindowsize. Reset to dwindowsize value.'
113 tim_dstepobsinc = tim_windowsize
114 end if
115 if (tim_referencetime == "middle") then
116 tim_nstepobs = 2*nint(((tim_windowsize - tim_dstepobs)/2.d0)/tim_dstepobs) + 1
117 tim_nstepobsinc = 2*nint(((tim_windowsize - tim_dstepobsinc)/2.d0)/tim_dstepobsinc) + 1
118 end if
119
120 if (trim(tim_referencetime) == "start") then
121 tim_nstepobs = max(nint(tim_windowsize/tim_dstepobs), 1)
122 tim_nstepobsinc = max(nint(tim_windowsize/tim_dstepobsinc), 1)
123 end if
124
125 end subroutine tim_readNml
126
127
128 subroutine tim_getDateStampFromEnvVar(dateStamp)
129 !
130 !:Purpose: Determine the date from the environment variable MIDAS_DATE.
131 !
132 implicit none
133
134 ! Arguments:
135 integer, intent(inout) :: dateStamp
136
137 ! Locals:
138 integer :: lengthValidDateStr, status, datePrint, timePrint, imode, ierr
139 integer(8) :: dateTimePrint
140 character(len=256) :: validDateStr
141 integer :: newdate
142
143 status = 0
144 call get_environment_variable('MIDAS_DATE',validDateStr,lengthValidDateStr,status,.true.)
145
146 if (status > 1) then
147 call utl_abort('tim_getDateStampFromEnvVar: Problem when getting the environment variable MIDAS_DATE')
148 end if
149 if (status == 1) then
150 write(*,*) 'tim_getDateStampFromEnvVar: WARNING: The environment variable MIDAS_DATE has not been detected!'
151 !call utl_abort('tim_getDateStampFromEnvVar: The environment variable MIDAS_DATE has not been detected!')
152 return
153 end if
154
155 write(*,*)
156 write(*,*) 'tim_getDateStampFromEnvVar: The environment variable MIDAS_DATE has correctly been detected'
157
158 ! convert string to long integer
159 read (validDateStr,*) dateTimePrint
160
161 ! split dateTime long integer into separate date and time values
162 if (lengthValidDateStr == 10) then
163 datePrint = dateTimePrint/100
164 timePrint = (dateTimePrint - datePrint*100) * 1000000
165 else if (lengthValidDateStr == 12) then
166 datePrint = dateTimePrint/10000
167 timePrint = (dateTimePrint - datePrint*10000) * 10000
168 else if (lengthValidDateStr == 14) then
169 datePrint = dateTimePrint/1000000
170 timePrint = (dateTimePrint - datePrint*1000000) * 100
171 else if (lengthValidDateStr == 16) then
172 datePrint = dateTimePrint/100000000
173 timePrint = (dateTimePrint - datePrint*100000000)
174 else
175 write(*,*) 'length of MIDAS_DATE = ', lengthValidDateStr
176 call utl_abort('tim_getDateStampFromEnvVar: Unexpected length of variable MIDAS_DATE')
177 end if
178
179 ! convert to CMC dateStamp
180 imode = 3 ! printable to stamp
181 ierr = newdate(datestamp, datePrint, timePrint, imode)
182
183 write(*,*) 'tim_getDateStampFromEnvVar: envVar, validDate, dateStamp = ', trim(validDateStr), dateTimePrint, dateStamp
184
185 end subroutine tim_getDateStampFromEnvVar
186
187
188 subroutine tim_setup(fileNameForDate_opt)
189 !
190 !:Purpose: Setup of obs time window size and related trial field
191 ! time step for OmP determination.
192 !
193 implicit none
194
195 ! Arguments:
196 character(len=*), optional, intent(in) :: fileNameForDate_opt
197
198 ! Locals:
199 integer :: ierr, newdate, imode, prntdate, prnttime
200 integer :: dateStampEnvVar
201
202 call tim_readNml()
203
204 if (initialized) then
205 write(*,*) 'tim_setup: already initialized, just return'
206 return
207 end if
208
209 ! First try to set dateStamp from MIDAS_DATE
210 dateStampEnvVar = 0
211 call tim_getDateStampFromEnvVar(dateStampEnvVar)
212
213 ! Possibly set the datestamp (except when set later from burp files)
214 if (dateStampEnvVar /= 0) then
215 write(*,*) 'tim_setup: ====================================================='
216 write(*,*) 'tim_setup: DATESTAMP set by value in supplied MIDAS_DATE'
217 write(*,*) 'tim_setup: ====================================================='
218 dateStamp = dateStampEnvVar
219 imode = -3 ! stamp to printable
220 ierr = newdate(datestamp, prntdate, prnttime, imode)
221 write(*,*) 'tim_setup: printdate = ',prntdate
222 write(*,*) 'tim_setup: printtime = ',prnttime
223 write(*,*) 'tim_setup: datestamp = ',datestamp
224 else if (present(fileNameForDate_opt)) then
225 write(*,*) 'tim_setup: ====================================================='
226 write(*,*) 'tim_setup: DATESTAMP set by value in supplied file'
227 write(*,*) 'tim_setup: ====================================================='
228 datestamp = tim_getDatestampFromFile(fileNameForDate_opt)
229 imode = -3 ! stamp to printable
230 ierr = newdate(datestamp, prntdate, prnttime, imode)
231 write(*,*) 'tim_setup: printdate = ',prntdate
232 write(*,*) 'tim_setup: printtime = ',prnttime
233 write(*,*) 'tim_setup: datestamp = ',datestamp
234 else
235 write(*,*) 'tim_setup: =========================================================='
236 write(*,*) 'tim_setup: DATESTAMP not set in this subroutine, use tim_setDateStamp'
237 write(*,*) 'tim_setup: =========================================================='
238 end if
239
240 if (mmpi_myid == 0) write(*,*) 'tim_setup: dobs_windowsize=',tim_windowsize
241 if (mmpi_myid == 0) write(*,*) 'tim_setup: dstepobs =',tim_dstepobs
242 if (mmpi_myid == 0) write(*,*) 'tim_setup: nstepobs =',tim_nstepobs
243 if (mmpi_myid == 0) write(*,*) 'tim_setup: dstepobsinc =',tim_dstepobsinc
244 if (mmpi_myid == 0) write(*,*) 'tim_setup: nstepobsinc =',tim_nstepobsinc
245 if (mmpi_myid == 0) write(*,*) 'tim_setup: tim_referencetime =',tim_referencetime
246
247 initialized = .true.
248
249 end subroutine tim_setup
250
251
252 function tim_initialized() result(initialized_out)
253 implicit none
254
255 ! Result:
256 logical initialized_out
257
258 initialized_out = initialized
259
260 end function tim_initialized
261
262
263 function tim_getDatestampFromFile(fileName, varNameForDate_opt) result(dateStamp_out)
264 !
265 !:Purpose: to extract the dateStamp from the supplied file.
266 !
267 implicit none
268
269 ! Arguments:
270 character(len=*), intent(in) :: fileName
271 character(len=*), optional, intent(in) :: varNameForDate_opt
272 ! Result:
273 integer :: dateStamp_out
274
275 ! Locals:
276 integer :: nulFile, ierr
277 integer, parameter :: maxNumDates = 2000
278 integer :: numDates, ikeys(maxNumDates), varIndex
279 integer :: fnom, fstouv, fstinl, fstprm, fstfrm, fclos, newdate
280 integer :: prntdate, prnttime, imode, windowIndex, windowsPerDay, dateStamp_tmp
281 logical :: fileExists, foundWindow, foundVarNameInFile
282 real(8) :: leadTimeInHours, windowBegHour, windowEndHour, fileHour, middleHour
283 integer :: ideet, inpas, dateStamp_origin, ini, inj, ink, inbits, idatyp
284 integer :: ip1, ip2, ip3, ig1, ig2, ig3, ig4, iswa, ilng, idltf, iubc
285 integer :: iextra1, iextra2, iextra3
286 character(len=2) :: typvar
287 character(len=4) :: nomvar
288 character(len=12) :: etiket
289 character(len=1) :: grtyp
290
291 if (mmpi_myid == 0) then
292
293 ! Check if file for any date within the analysis window (except the last) exists
294 inquire(file=trim(fileName), exist=fileExists)
295 if (.not.fileExists) then
296 call utl_abort('tim_getDateStampFromFile: file not found '//trim(fileName))
297 end if
298
299 ! Determine variable to use for the date (default is P0)
300 varNameForDate = 'P0'
301 if (present(varNameForDate_opt)) then
302
303 varNameForDate = trim(varNameForDate_opt)
304 write(*,*) 'tim_getDateStampFromFile: defining dateStamp from the variable = ', varNameForDate
305
306 ! If P0 not present, look for another suitable variable in the file
307 else if (.not. utl_varNamePresentInFile(varNameForDate,fileName_opt=trim(fileName))) then
308
309 foundVarNameInFile = .false.
310 do varIndex = 1, vnl_numvarmax
311 varNameForDate = vnl_varNameList(varIndex)
312 ! check if variable is in the file
313 if (.not. utl_varNamePresentInFile(varNameForDate,fileName_opt=trim(fileName))) cycle
314 foundVarNameInFile = .true.
315 exit
316 end do
317
318 if (.not. foundVarNameInFile) then
319 call utl_abort('tim_getDateStampFromFile: NO variables found in the file!!!')
320 end if
321
322 end if
323
324 ! Extract the datestamp from the file
325 nulFile = 0
326 ierr = fnom(nulFile,trim(fileName),'RND+OLD+R/O',0)
327 ierr = fstouv(nulFile,'RND+OLD')
328 ierr = fstinl(nulFile,ini,inj,ink,-1,' ',-1,-1,-1,' ', &
329 trim(varNameForDate),ikeys,numdates,maxnumdates)
330 if (ikeys(1) <= 0) then
331 call utl_abort('tim_getDateStampFromFile: Could not find variable ' // &
332 trim(varNameForDate) // ' in the supplied file')
333 end if
334 write(*,*) 'tim_getDateStampFromFile: number of dates found = ', numDates
335 ierr = fstprm(ikeys(1), dateStamp_origin, ideet, inpas, ini, inj, &
336 ink, inbits, idatyp, ip1, ip2, ip3, &
337 typvar, nomvar, etiket, grtyp, ig1, ig2, ig3, ig4, &
338 iswa, ilng, idltf, iubc, iextra1, iextra2, iextra3)
339 leadTimeInHours = real(ideet*inpas,8)/3600.0d0
340 call incdatr(dateStamp_out, dateStamp_origin, leadTimeInHours)
341
342 ierr = fstfrm(nulFile)
343 ierr = fclos(nulFile)
344
345 end if
346
347 call rpn_comm_bcast(dateStamp_out, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
348
349 if (tim_referenceTime == 'middle') then
350 ! Modify date to ensure that it corresponds to the middle of the window
351 ! Note: For this, we have to assume that the date in the file
352 ! does NOT correspond to the final time of the window
353 imode = -3 ! stamp to printable, time is HHMMSShh
354 ierr = newdate(datestamp_out, prntdate, prnttime, imode)
355 fileHour = real(prnttime,8)/1000000.0d0
356 windowsPerDay = nint(24.0d0 / tim_windowsize)
357 foundWindow = .false.
358 window_loop: do windowIndex = 0, windowsPerDay
359 windowBegHour = (real(windowIndex,8) * tim_windowsize) - (tim_windowsize/2.0d0)
360 windowEndHour = (real(windowIndex,8) * tim_windowsize) + (tim_windowsize/2.0d0)
361 if (fileHour >= windowBegHour .and. fileHour < windowEndHour) then
362 foundWindow = .true.
363 middleHour = real(windowIndex,8) * tim_windowsize
364 exit window_loop
365 end if
366 end do window_loop
367
368 if (.not. foundWindow) then
369 write(*,*) 'windowsPerDay, fileHour=', windowsPerDay, fileHour
370 call utl_abort('tim_getDateStampFromFile: could not determine assimilation window position')
371 end if
372
373 ! handle special case when window centered on hour 24
374 if (nint(middleHour) == 24) then
375 ! add 24h to dateStamp and recompute prntdate
376 dateStamp_tmp = dateStamp_out
377 call incdatr(dateStamp_out, dateStamp_tmp, 24.0d0)
378 imode = -3 ! stamp to printable, only prntdate will be used
379 ierr = newdate(datestamp_out, prntdate, prnttime, imode)
380
381 ! subtract 24h from middleHour
382 middleHour = 0.0d0
383 end if
384
385 prnttime = nint(middleHour) * 1000000
386 imode = 3 ! printable to stamp
387 ierr = newdate(datestamp_out, prntdate, prnttime, imode)
388
389 end if
390
391 end function tim_getDateStampFromFile
392
393
394 subroutine tim_setDatestamp(datestamp_in)
395 !
396 !:Purpose: to control access to the minimization object. Sets the date
397 ! of the window centre of analysis validity to the indicated value.
398 !
399 implicit none
400
401 ! Arguments:
402 integer, intent(in) :: datestamp_in
403
404 if (.not.initialized) call utl_abort('tim_setDateStamp: module not initialized')
405
406 datestamp = datestamp_in
407
408 end subroutine tim_setDatestamp
409
410
411 function tim_getDatestamp() result(datestamp_out)
412 !
413 !:Purpose: to control access to the minimization object. Returns the date
414 ! of the window centre of analysis validity.
415 !
416 implicit none
417
418 ! Result:
419 integer :: datestamp_out
420
421 if (.not.initialized) call utl_abort('tim_getDateStamp: module not initialized')
422
423 datestamp_out = datestamp
424
425 end function tim_getDatestamp
426
427
428 subroutine tim_getStampList(datestamplist, numStep, referenceDateStamp)
429 !
430 !:Purpose: Compute a list of STAMPS corresponding to stepobs time
431 !
432 implicit none
433
434 ! Arguments:
435 integer, intent(in) :: numStep ! number of step obs
436 integer, intent(in) :: referenceDateStamp ! Synoptic time
437 integer, intent(out) :: datestamplist(numStep)
438
439 ! Locals:
440 integer :: stepIndex
441 real(8) :: dldelt ! delta time in hours between middle time and each step
442 real(8) :: dtstep ! delta time in hours between step obs
443
444 if (.not. initialized) call utl_abort('tim_getStampList: module not initialized')
445
446 if (referenceDateStamp == -1) then
447
448 if (mmpi_myid == 0) write(*,*) 'tim_getStampList: datestamp is not specified, keep as -1'
449 datestamplist(:) = -1
450
451 else
452
453 if (tim_referencetime == "middle") then
454 if (numStep > 1) then
455 dtstep = tim_windowsize/(real(numStep-1,8))
456 else
457 dtstep = tim_windowsize
458 end if
459
460 do stepIndex = 1, numStep
461 dldelt = (stepIndex - ((numStep-1)/2 + 1)) * dtstep
462 call incdatr(datestamplist(stepIndex), referenceDateStamp, dldelt)
463 end do
464 end if
465
466 if (trim(tim_referencetime) == "start") then
467
468 dtstep = tim_windowsize/(real(numStep,8))
469
470 do stepIndex = 1, numStep
471 dldelt = (stepIndex - 1) * dtstep
472 call incdatr(datestamplist(stepIndex), referenceDateStamp, dldelt)
473 end do
474
475 end if
476
477 end if
478
479 end subroutine tim_getStampList
480
481
482 subroutine tim_getStepObsIndex(dnstepobs, referenceDateStamp, obsYYYMMDD, obsHHMM, numStep)
483 !
484 !:Purpose: Return the stepobs index as a real number (-1.0 if out of range)
485 !
486 implicit none
487
488 ! Arguments:
489 real(8), intent(out) :: dnstepobs ! number of stepobs from reference time
490 integer, intent(in) :: referenceDateStamp ! Synop CMC date-time stamp
491 integer, intent(in) :: obsYYYMMDD ! Obs date YYYYMMDD
492 integer, intent(in) :: obsHHMM ! Obs time HHMM
493 integer, intent(in) :: numStep ! number of stepobs in assimilation window
494
495 ! Locals:
496 integer :: newdate, istat, imode
497 real(8) :: dddt ! delta time in hours
498 integer :: istobs ! obs CMC date-time stamp
499 integer :: itobs ! obs time HHMMSShh
500 real(8) :: dlhours ! delta time from synop time
501
502 if (.not. initialized) call utl_abort('tim_getStepObsIndex: module not initialized')
503
504 ! Building observation stamp
505 imode = 3 ! printable to stamp
506 itobs = obsHHMM * 10000
507 istat = newdate(istobs, obsYYYMMDD, itobs, imode)
508
509 ! Difference (in hours) between obs time and reference time
510 call difdatr(istobs, referenceDateStamp, dlhours)
511
512 if (numStep > 1) then
513 ! FGAT: more than 1 trial field in time window
514 if (tim_referencetime == "middle") then
515 dddt = tim_windowsize / (real(numStep-1,8))
516 else
517 dddt = tim_windowsize / (real(numStep,8))
518 end if
519 dnstepobs = dlhours / dddt ! number of step obs from reference (e.g. synoptic)
520 if (tim_referencetime == "middle") then
521 dnstepobs = dnstepobs + real((numStep+1)/2,8)
522 end if
523 if (trim(tim_referencetime) == "start") then
524 dnstepobs = dnstepobs + 1.d0
525 end if
526 if (dnstepobs < 0.5d0.or.dnstepobs > (0.5d0+real(numStep,8))) dnstepobs = -1.0d0
527
528 else
529 ! 3D: only 1 trial field in time window
530 if (tim_referencetime == "middle") then
531 if (dlhours < -tim_windowsize/2.0D0 .or. dlhours > tim_windowsize/2.0D0) then
532 ! outside time window
533 dnstepobs = -1.0d0
534 else
535 ! inside time window
536 dnstepobs = 1.0d0
537 end if
538 else
539 dddt = tim_windowsize
540 if (dlhours < -dddt/2.0d0 .or. dlhours > tim_windowsize + dddt/2.d0) then
541 ! outside time window
542 dnstepobs = -1.0d0
543 else
544 ! inside time window
545 dnstepobs = 1.0d0
546 end if
547 endif
548
549 end if
550
551 end subroutine tim_getStepObsIndex
552
553 !----------------------------------------------------------------------------------------
554 ! tim_dateStampToDDMMYYYY
555 !----------------------------------------------------------------------------------------
556 subroutine tim_dateStampToYYYYMMDDHH(dateStamp, prnttime, dd, mm, ndays, yyyy, verbose_opt)
557 !
558 !:Purpose: to get day (DD), month (MM), number of days in this month
559 ! and year (YYYY) from dateStamp
560 !
561 implicit none
562
563 ! Arguments:
564 integer, intent(in) :: dateStamp
565 integer, intent(inout) :: prnttime, dd, mm, ndays, yyyy
566 logical, optional, intent(in) :: verbose_opt
567
568 ! Locals:
569 character(len=8) :: yyyymmdd
570 character(len=3), parameter :: months(12) = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
571 integer :: ndaysM(12)
572 integer :: imode, ierr, newdate, prntdate
573 logical :: verbose = .True.
574
575 ndaysM(:) = [ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
576
577 imode = -3 ! stamp to printable
578 ierr = newdate(dateStamp, prntdate, prnttime, imode)
579 write(yyyymmdd,'(i8)') prntdate
580 read (yyyymmdd(7:8), '(i)') dd
581 read (yyyymmdd(5:6), '(i)') mm
582 read (yyyymmdd(1:4), '(i)') yyyy
583
584 ! leap year
585 if (mm == 2 .and. mod(yyyy,4)==0) ndaysM(mm) = 29
586 ndays = ndaysM(mm)
587
588 if (present(verbose_opt)) verbose = verbose_opt
589
590 if(verbose) then
591 write(*,*) 'tim_dateStampToYYYYMMDDHH: date = ', prntdate
592 write(*,*) 'tim_dateStampToYYYYMMDDHH: year = ', yyyy
593 write(*,'(a,i5,a,i5,a)') 'tim_dateStampToYYYYMMDDHH: month = ', mm, ' ( '// months(mm)//' where there are ', ndays, ' days)'
594 write(*,*) 'tim_dateStampToYYYYMMDDHH: day = ', dd
595 write(*,*) 'tim_dateStampToYYYYMMDDHH: time = ', prnttime
596 end if
597
598 end subroutine tim_dateStampToYYYYMMDDHH
599
600 !----------------------------------------------------------------------------------------
601 ! tim_getValidDateTimeFromList
602 !----------------------------------------------------------------------------------------
603 subroutine tim_getValidDateTimeFromList(headDateValues, headTimeValues, validDate, validtime)
604 implicit none
605
606 ! Arguments:
607 integer, intent(in) :: headDateValues(:)
608 integer, intent(in) :: headTimeValues(:)
609 integer, intent(out) :: validDate
610 integer, intent(out) :: validTime
611
612 ! Locals:
613 integer :: numDates, numWindowsPerDay, windowIndex, timeMin, timeMax, dateMin, dateMax
614 integer :: windowBoundaryMin, windowBoundaryMax, validTimeMin, validTimeMax, validDateMin, validDateMax
615 integer(8) :: dateTimeMin, dateTimeMax
616 integer(8), allocatable :: dateTimeValues(:), windowBoundaries(:)
617 integer :: ier, imode
618 integer :: newdate, dateStampIn, dateStampOut
619
620 call tim_readNml()
621
622 numDates = size(headDateValues)
623 checkNumDates: if (numDates > 0) then
624 write(*,*) 'tim_getValidDateTimeFromList: check inputs: time min/max = ', minval(headTimeValues), maxval(headTimeValues)
625 write(*,*) 'tim_getValidDateTimeFromList: check inputs: date min/max = ', minval(headDateValues), maxval(headDateValues)
626 allocate(dateTimeValues(size(headDateValues)))
627 dateTimeValues(:) = headDateValues(:)
628 dateTimeValues(:) = 10000*dateTimeValues(:) + headTimeValues(:)
629
630 dateTimeMin = minval(dateTimeValues(:))
631 dateTimeMax = maxval(dateTimeValues(:))
632 deallocate(dateTimeValues)
633 dateMin = dateTimeMin/10000
634 dateMax = dateTimeMax/10000
635 timeMin = dateTimeMin - 10000*(dateTimeMin/10000)
636 timeMax = dateTimeMax - 10000*(dateTimeMax/10000)
637 ! convert from hhmm to just minutes: hhmm - 100*(hhmm/100) + 60*(hhmm/100)
638 timeMin = timeMin - 100*(timeMin/100) + 60*(timeMin/100)
639 timeMax = timeMax - 100*(timeMax/100) + 60*(timeMax/100)
640 write(*,*) 'tim_getValidDateTimeFromList: min/max DateTime = ', dateTimeMin, dateTimeMax
641 write(*,*) 'tim_getValidDateTimeFromList: min/max time (in minutes) = ', timeMin, timeMax
642 if (tim_windowSize < 24.0d0) then
643 numWindowsPerDay = nint(24.0/tim_windowSize)
644
645 ! define boundaries between assimilation windows in minutes relative to 0UTC
646 allocate(windowBoundaries(0:numWindowsPerDay))
647 do windowIndex = 0, numWindowsPerDay
648 ! example for windowSize=6h, boundaries = -3h,+3h,+9h,+15h,+21h
649 windowBoundaries(windowIndex) = nint(windowIndex*60.0*tim_windowSize - 60.0*tim_windowSize/2.0)
650 end do
651 write(*,*) 'tim_getValidDateTimeFromList: boundaries (in minutes) = ', windowBoundaries(:)
652
653 ! find left boundary of window where timeMin/Max are located
654 windowBoundaryMin = -1
655 windowBoundaryMax = -1
656 do windowIndex = 0, numWindowsPerDay
657 if (timeMin >= windowBoundaries(windowIndex)) then
658 windowBoundaryMin = windowIndex
659 end if
660 if (timeMax >= windowBoundaries(windowIndex)) then
661 windowBoundaryMax = windowIndex
662 end if
663 end do
664
665 ! find validTimeMin/Max from left boundary
666 validTimeMin = nint((windowBoundaries(windowBoundaryMin) + 60.0*tim_windowSize/2.0)/60.0)
667 if (validTimeMin >= 24) then
668 validTimeMin = 0
669 imode = 3
670 ier = newdate(dateStampIn, dateMin, validTimeMin, imode)
671 call incdat(dateStampOut, dateStampIn, 24) ! add 1 day to get validDate
672 imode = -3
673 ier = newdate(dateStampOut, validDateMin, validTimeMin, imode)
674 validTimeMin = 0
675 else
676 validDateMin = dateMin
677 end if
678 validTimeMax = nint((windowBoundaries(windowBoundaryMax) + 60.0*tim_windowSize/2.0)/60.0)
679 if (validTimeMax >= 24) then
680 validTimeMax = 0
681 imode = 3
682 ier = newdate(dateStampIn, dateMax, validTimeMax, imode)
683 call incdat(dateStampOut, dateStampIn, 24) ! add 1 day to get validDate
684 imode = -3
685 ier = newdate(dateStampOut, validDateMax, validTimeMax, imode)
686 validTimeMax = 0
687 else
688 validDateMax = dateMax
689 end if
690 write(*,*) 'tim_getValidDateTimeFromList: date from Min/Max = ', validDateMin, validDateMax
691 write(*,*) 'tim_getValidDateTimeFromList: hour from Min/Max = ', validTimeMin, validTimeMax
692 if (validTimeMin /= validTimeMax) call utl_abort('validTimeMin/Max not equal')
693 validTime = validTimeMin
694 validDate = validDateMin
695 deallocate(windowBoundaries)
696 else
697 write(*,*) 'tim_getValidDateTimeFromList: WARNING: window size equal or greater than 1 day, cannot get dateStamp'
698 validTime = 0
699 validDate = 0
700 end if
701 end if checkNumDates
702
703 end subroutine tim_getValidDateTimeFromList
704
705end module timeCoord_mod
