libsim  Versione 7.2.4
stat_proc_engine.F90
1 ! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2 ! authors:
3 ! Davide Cesari <dcesari@arpa.emr.it>
4 ! Paolo Patruno <ppatruno@arpa.emr.it>
5 
6 ! This program is free software; you can redistribute it and/or
7 ! modify it under the terms of the GNU General Public License as
8 ! published by the Free Software Foundation; either version 2 of
9 ! the License, or (at your option) any later version.
10 
11 ! This program is distributed in the hope that it will be useful,
12 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
13 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 ! GNU General Public License for more details.
15 
16 ! You should have received a copy of the GNU General Public License
17 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
18 #include "config.h"
19 
23 MODULE stat_proc_engine
25 USE vol7d_class
26 IMPLICIT NONE
27 
28 ! per ogni elemento i,j dell'output, definire n elementi di input ad
29 ! esso contribuenti (arrayof_ttr_mapper) con le seguenti informazioni
30 TYPE ttr_mapper
31  INTEGER :: it=imiss ! k
32  INTEGER :: itr=imiss ! l
33  INTEGER :: extra_info=imiss ! dtratio for intervals, extreme for point in time
34  TYPE(datetime) :: time=datetime_miss ! time for point in time
35 END TYPE ttr_mapper
36 
37 INTERFACE OPERATOR (==)
38  MODULE PROCEDURE ttr_mapper_eq
39 END INTERFACE
40 
41 INTERFACE OPERATOR (/=)
42  MODULE PROCEDURE ttr_mapper_ne
43 END INTERFACE
44 
45 INTERFACE OPERATOR (>)
46  MODULE PROCEDURE ttr_mapper_gt
47 END INTERFACE
48 
49 INTERFACE OPERATOR (<)
50  MODULE PROCEDURE ttr_mapper_lt
51 END INTERFACE
52 
53 INTERFACE OPERATOR (>=)
54  MODULE PROCEDURE ttr_mapper_ge
55 END INTERFACE
56 
57 INTERFACE OPERATOR (<=)
58  MODULE PROCEDURE ttr_mapper_le
59 END INTERFACE
60 
61 #undef VOL7D_POLY_TYPE
62 #undef VOL7D_POLY_TYPES
63 #undef ENABLE_SORT
64 #define VOL7D_POLY_TYPE TYPE(ttr_mapper)
65 #define VOL7D_POLY_TYPES _ttr_mapper
66 #define ENABLE_SORT
67 #include "array_utilities_pre.F90"
68 
69 #define ARRAYOF_ORIGTYPE TYPE(ttr_mapper)
70 #define ARRAYOF_TYPE arrayof_ttr_mapper
71 #define ARRAYOF_ORIGEQ 1
72 #define ARRAYOF_ORIGGT 1
73 #include "arrayof_pre.F90"
74 ! from arrayof
75 
76 
77 CONTAINS
78 
79 ! simplified comparisons only on time
80 ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res)
81 TYPE(ttr_mapper),INTENT(IN) :: this, that
82 LOGICAL :: res
83 
84 res = this%time == that%time
85 
86 END FUNCTION ttr_mapper_eq
87 
88 ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res)
89 TYPE(ttr_mapper),INTENT(IN) :: this, that
90 LOGICAL :: res
91 
92 res = this%time /= that%time
93 
94 END FUNCTION ttr_mapper_ne
95 
96 ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res)
97 TYPE(ttr_mapper),INTENT(IN) :: this, that
98 LOGICAL :: res
99 
100 res = this%time > that%time
101 
102 END FUNCTION ttr_mapper_gt
103 
104 ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res)
105 TYPE(ttr_mapper),INTENT(IN) :: this, that
106 LOGICAL :: res
107 
108 res = this%time < that%time
109 
110 END FUNCTION ttr_mapper_lt
111 
112 ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res)
113 TYPE(ttr_mapper),INTENT(IN) :: this, that
114 LOGICAL :: res
115 
116 res = this%time >= that%time
117 
118 END FUNCTION ttr_mapper_ge
119 
120 ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res)
121 TYPE(ttr_mapper),INTENT(IN) :: this, that
122 LOGICAL :: res
123 
124 res = this%time <= that%time
125 
126 END FUNCTION ttr_mapper_le
127 
128 #include "arrayof_post.F90"
129 #include "array_utilities_inc.F90"
130 
131 
132 ! common operations for statistical processing by differences
133 SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, &
134  otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, &
135  start)
136 TYPE(datetime),INTENT(in) :: itime(:)
137 TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
138 INTEGER,INTENT(in) :: stat_proc
139 TYPE(timedelta),INTENT(in) :: step
140 TYPE(datetime),POINTER :: otime(:)
141 TYPE(vol7d_timerange),POINTER :: otimerange(:)
142 INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
143 INTEGER :: nitr
144 LOGICAL,ALLOCATABLE :: mask_timerange(:)
145 INTEGER,INTENT(in) :: time_definition
146 LOGICAL,INTENT(in),OPTIONAL :: full_steps
147 TYPE(datetime),INTENT(in),OPTIONAL :: start
148 
149 INTEGER :: i, j, k, l, dirtyrep
150 INTEGER :: steps
151 LOGICAL :: lfull_steps, useful
152 TYPE(datetime) :: lstart, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
153 TYPE(vol7d_timerange) :: tmptimerange
154 TYPE(arrayof_datetime) :: a_otime
155 TYPE(arrayof_vol7d_timerange) :: a_otimerange
156 TYPE(timedelta) :: start_delta
157 
158 ! compute length of cumulation step in seconds
159 CALL getval(step, asec=steps)
160 
161 lstart = datetime_miss
162 IF (PRESENT(start)) lstart = start
163 lfull_steps = optio_log(full_steps)
164 
165 ! create a mask of suitable timeranges
166 ALLOCATE(mask_timerange(SIZE(itimerange)))
167 mask_timerange(:) = itimerange(:)%timerange == stat_proc &
168  .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
169  .AND. itimerange(:)%p1 >= 0 &
170  .AND. itimerange(:)%p2 > 0
171 
172 IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
173  mask_timerange(:) = mask_timerange(:) .AND. &
174  (itimerange(:)%p1 == 0 .OR. & ! all analyses
175  mod(itimerange(:)%p1, steps) == 0 .OR. & ! end time is mod step
176  mod(itimerange(:)%p1 - itimerange(:)%p2, steps) == 0) ! start time is mod step
177 ENDIF
178 ! mask_timerange includes all candidate timeranges
179 
180 nitr = count(mask_timerange)
181 ALLOCATE(f(nitr))
182 j = 1
183 DO i = 1, nitr
184  DO WHILE(.NOT.mask_timerange(j))
185  j = j + 1
186  ENDDO
187  f(i) = j
188  j = j + 1
189 ENDDO
190 
191 ! now we have to evaluate time/timerage pairs which do not need processing
192 ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
193 CALL compute_keep_tr()
194 
195 ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
196 map_tr(:,:,:,:,:) = imiss
197 
198 ! scan through all possible combinations of time and timerange
199 DO dirtyrep = 1, 2
200  IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
201  CALL packarray(a_otime)
202  CALL packarray(a_otimerange)
203  CALL sort(a_otime%array)
204  CALL sort(a_otimerange%array)
205  ENDIF
206  DO l = 1, SIZE(itime)
207  DO k = 1, nitr
208  CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
209  time_definition, pstart2, pend2, reftime2)
210 
211  DO j = 1, SIZE(itime)
212  DO i = 1, nitr
213  useful = .false.
214  CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
215  time_definition, pstart1, pend1, reftime1)
216  tmptimerange = vol7d_timerange_new(timerange=stat_proc)
217 
218  IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
219  IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
220  CALL time_timerange_set_period(tmptime, tmptimerange, &
221  time_definition, pend1, pend2, reftime2)
222  IF (lfull_steps) THEN
223  IF (mod(reftime2, step) == timedelta_0) THEN
224  useful = .true.
225  ENDIF
226  ELSE
227  useful = .true.
228  ENDIF
229 
230  ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
231  CALL time_timerange_set_period(tmptime, tmptimerange, &
232  time_definition, pstart2, pstart1, pstart1)
233  IF (lfull_steps) THEN
234  IF (mod(pstart1, step) == timedelta_0) THEN
235  useful = .true.
236  ENDIF
237  ELSE
238  useful = .true.
239  ENDIF
240  ENDIF
241 
242  ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
243  IF (lfull_steps) THEN
244  IF (c_e(lstart)) THEN
245 ! lstart shifts the interval for computing modulo step, this does not
246 ! remove data before lstart but just shifts the phase
247  start_delta = lstart-reftime2
248  ELSE
249  start_delta = timedelta_0
250  ENDIF
251  ENDIF
252 
253  IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
254  CALL time_timerange_set_period(tmptime, tmptimerange, &
255  time_definition, pend1, pend2, reftime2)
256  IF (lfull_steps) THEN
257  IF (mod(pend2-reftime2-start_delta, step) == timedelta_0) THEN
258  useful = .true.
259  ENDIF
260  ELSE
261  useful = .true.
262  ENDIF
263 ! keep only data after lstart
264  IF (c_e(lstart)) THEN
265  IF (lstart > pend1) useful = .false.
266  ENDIF
267 
268  ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
269  CALL time_timerange_set_period(tmptime, tmptimerange, &
270  time_definition, pstart2, pstart1, reftime2)
271  IF (lfull_steps) THEN
272  IF (mod(pstart1-reftime2-start_delta, step) == timedelta_0) THEN
273  useful = .true.
274  ENDIF
275  ELSE
276  useful = .true.
277  ENDIF
278 ! keep only data after lstart
279  IF (c_e(lstart)) THEN
280  IF (lstart > pstart2) useful = .false.
281  ENDIF
282 
283  ENDIF
284  ENDIF
285  useful = useful .AND. tmptime /= datetime_miss .AND. &
286  tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
287 
288  IF (useful) THEN ! add a_otime, a_otimerange
289  map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
290  map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
291  ENDIF
292  ENDDO
293  ENDDO
294  ENDDO
295  ENDDO
296 ENDDO
297 
298 ! we have to repeat the computation with sorted arrays
299 CALL compute_keep_tr()
300 
301 otime => a_otime%array
302 otimerange => a_otimerange%array
303 ! delete local objects keeping the contents
304 CALL delete(a_otime, nodealloc=.true.)
305 CALL delete(a_otimerange, nodealloc=.true.)
306 
307 #ifdef DEBUG
308 CALL l4f_log(l4f_debug, &
309  'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr,1)))//', '// &
310  t2c((SIZE(map_tr,2)))//', '// &
311  t2c((SIZE(map_tr,3)))//', '// &
312  t2c((SIZE(map_tr,4))))
313 CALL l4f_log(l4f_debug, &
314  'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr))/2)//', '// &
315  t2c(count(c_e(map_tr))/2))
316 CALL l4f_log(l4f_debug, &
317  'recompute_stat_proc_diff, nitr: '//t2c(nitr))
318 CALL l4f_log(l4f_debug, &
319  'recompute_stat_proc_diff, good timeranges: '//t2c(count(c_e(keep_tr))/2))
320 CALL l4f_log(l4f_debug, &
321  'recompute_stat_proc_diff, output times: '//t2c(SIZE(otime)))
322 CALL l4f_log(l4f_debug, &
323  'recompute_stat_proc_diff, output timeranges: '//t2c(SIZE(otimerange)))
324 #endif
325 
326 CONTAINS
327 
328 SUBROUTINE compute_keep_tr()
329 INTEGER :: start_deltas
330 
331 keep_tr(:,:,:) = imiss
332 DO l = 1, SIZE(itime)
333  itrloop: DO k = 1, nitr
334  IF (itimerange(f(k))%p2 == steps) THEN
335  CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
336  time_definition, pstart2, pend2, reftime2)
337  useful = .false.
338 ! keep only data after lstart
339  IF (c_e(lstart)) THEN
340  IF (lstart > pstart2) cycle itrloop
341  ENDIF
342  IF (reftime2 == pend2) THEN ! analysis
343  IF (c_e(lstart)) THEN ! in analysis mode start wins over full_steps
344  IF (mod(reftime2-lstart, step) == timedelta_0) THEN
345  useful = .true.
346  ENDIF
347  ELSE IF (lfull_steps) THEN
348  IF (mod(reftime2, step) == timedelta_0) THEN
349  useful = .true.
350  ENDIF
351  ELSE
352  useful = .true.
353  ENDIF
354  ELSE ! forecast
355  IF (lfull_steps) THEN
356 ! same as above for start_delta, but in seconds and not in timerange/timedelta units
357  IF (c_e(lstart)) THEN
358  start_deltas = timedelta_getamsec(lstart-reftime2)/1000_int_ll
359  ELSE
360  start_deltas = 0
361  ENDIF
362  IF (mod(itimerange(f(k))%p1 - start_deltas, steps) == 0) THEN
363  useful = .true.
364  ENDIF
365  ELSE
366  useful = .true.
367  ENDIF
368  ENDIF
369  IF (useful) THEN
370 ! CALL time_timerange_set_period(tmptime, tmptimerange, &
371 ! time_definition, pstart2, pend2, reftime2)
372  keep_tr(k,l,1) = append_unique(a_otime, itime(l))
373  keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
374  ENDIF
375  ENDIF
376  ENDDO itrloop
377 ENDDO
378 
379 END SUBROUTINE compute_keep_tr
380 
381 END SUBROUTINE recompute_stat_proc_diff_common
382 
383 
384 ! common operations for statistical processing by metamorphosis
385 SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
386  otimerange, map_tr)
387 INTEGER,INTENT(in) :: istat_proc
388 TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
389 INTEGER,INTENT(in) :: ostat_proc
390 TYPE(vol7d_timerange),POINTER :: otimerange(:)
391 INTEGER,POINTER :: map_tr(:)
392 
393 INTEGER :: i
394 LOGICAL :: tr_mask(SIZE(itimerange))
395 
396 IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
397  ALLOCATE(otimerange(0), map_tr(0))
398  RETURN
399 ENDIF
400 
401 ! useful timeranges
402 tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
403  .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
404 ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
405 
406 otimerange = pack(itimerange, mask=tr_mask)
407 otimerange(:)%timerange = ostat_proc
408 map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
409 
410 END SUBROUTINE compute_stat_proc_metamorph_common
411 
412 
413 ! common operations for statistical processing by aggregation
414 SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
415  step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
416 TYPE(datetime),INTENT(in) :: itime(:)
417 TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
418 INTEGER,INTENT(in) :: stat_proc
419 INTEGER,INTENT(in) :: tri
420 TYPE(timedelta),INTENT(in) :: step
421 INTEGER,INTENT(in) :: time_definition
422 TYPE(datetime),POINTER :: otime(:)
423 TYPE(vol7d_timerange),POINTER :: otimerange(:)
424 TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
425 INTEGER,POINTER,OPTIONAL :: dtratio(:)
426 TYPE(datetime),INTENT(in),OPTIONAL :: start
427 LOGICAL,INTENT(in),OPTIONAL :: full_steps
428 
429 INTEGER :: i, j, k, l, na, nf, n
430 INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
431 INTEGER(kind=int_ll) :: stepms, mstepms
432 LOGICAL :: lforecast
433 TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
434 TYPE(arrayof_datetime) :: a_otime
435 TYPE(arrayof_vol7d_timerange) :: a_otimerange
436 TYPE(arrayof_integer) :: a_dtratio
437 LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
438 TYPE(ttr_mapper) :: lmapper
439 CHARACTER(len=8) :: env_var
440 LOGICAL :: climat_behavior
441 
442 
443 ! enable bad behavior for climat database (only for instantaneous data)
444 env_var = ''
445 CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
446 climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
447 
448 ! compute length of cumulation step in seconds
449 CALL getval(timedelta_depop(step), asec=steps)
450 
451 ! create a mask of suitable timeranges
452 ALLOCATE(mask_timerange(SIZE(itimerange)))
453 mask_timerange(:) = itimerange(:)%timerange == tri &
454  .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
455 
456 IF (PRESENT(dtratio)) THEN
457  WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
458  mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
459  ELSEWHERE
460  mask_timerange(:) = .false.
461  END WHERE
462 ELSE ! instantaneous
463  mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
464 ENDIF
465 
466 #ifdef DEBUG
467 CALL l4f_log(l4f_debug, &
468  '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
469  t2c(count(mask_timerange)))
470 #endif
471 
472 ! euristically determine whether we are dealing with an
473 ! analysis/observation or a forecast dataset
474 na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
475 nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
476 lforecast = nf >= na
477 #ifdef DEBUG
478 CALL l4f_log(l4f_debug, &
479  'recompute_stat_proc_agg, na: '//t2c(na)//', nf: '//t2c(nf))
480 #endif
481 
482 ! keep only timeranges of one type (really necessary?)
483 IF (lforecast) THEN
484  CALL l4f_log(l4f_info, &
485  'recompute_stat_proc_agg, processing in forecast mode')
486 ELSE
487  mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
488  CALL l4f_log(l4f_info, &
489  'recompute_stat_proc_agg, processing in analysis mode')
490 ENDIF
491 
492 #ifdef DEBUG
493 CALL l4f_log(l4f_debug, &
494  '(re)compute_stat_proc_agg, number of useful timeranges: '// &
495  t2c(count(mask_timerange)))
496 #endif
497 
498 IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
499  ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
500  IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
501  RETURN
502 ENDIF
503 
504 ! determine start and end of processing period, should work with p2==0
505 lstart = datetime_miss
506 IF (PRESENT(start)) lstart = start
507 lend = itime(SIZE(itime))
508 ! compute some quantities
509 maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
510 maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
511 minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
512 IF (time_definition == 0) THEN ! reference time
513  lend = lend + timedelta_new(sec=maxp1)
514 ENDIF
515 ! extend interval at the end in order to include all the data in case
516 ! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
517 ! below in order to exclude the last full interval which would be empty
518 lend = lend + step
519 IF (lstart == datetime_miss) THEN ! autodetect
520  lstart = itime(1)
521 ! if autodetected, adjust to obtain real absolute start of data
522  IF (time_definition == 0) THEN ! reference time
523  lstart = lstart + timedelta_new(sec=minp1mp2)
524  ELSE ! verification time
525 ! go back to start of longest processing interval
526  lstart = lstart - timedelta_new(sec=maxp2)
527  ENDIF
528 ! apply full_steps in analysis mode and when start is not specified
529 ! (start by itself in analysis mode implies full_steps with respect to
530 ! start instead of absolute full steps), todo in forecast mode
531  IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
532  lstart = lstart - (mod(lstart, step)) ! round to step, (should be MODULO, not MOD)
533  ENDIF
534 ENDIF
535 
536 #ifdef DEBUG
537 CALL l4f_log(l4f_debug, &
538  'recompute_stat_proc_agg, processing period: '//t2c(lstart)//' - '//t2c(lend))
539 #endif
540 
541 ! create output time and timerange lists
542 
543 IF (lforecast) THEN ! forecast mode
544  IF (time_definition == 0) THEN ! reference time
545  CALL insert(a_otime, itime) ! should I limit to elements itime >= lstart?
546 
547 ! apply start shift to timerange, not to reftime
548 ! why did we use itime(SIZE(itime)) (last ref time)?
549 ! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
550  CALL getval(lstart-itime(1), asec=dstart)
551 ! allow starting before first reftime but restrict dtstart to -steps
552 ! dstart = MAX(0, dstart)
553  IF (dstart < 0) dstart = mod(dstart, steps)
554  DO p1 = steps + dstart, maxp1, steps
555  CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
556  ENDDO
557 
558  ELSE ! verification time
559 
560 ! verification time in forecast mode is the ugliest case, because we
561 ! don't know a priori how many different (thus incompatible) reference
562 ! times we have, so some assumption of regularity has to be made. For
563 ! this reason msteps, the minimum step between two times, is
564 ! computed. We choose to compute it as a difference between itime
565 ! elements, it could be also computed as difference of itimerange%p1
566 ! elements. But what if it is not modulo steps?
567  mstepms = steps*1000_int_ll
568  DO i = 2, SIZE(itime)
569  CALL getval(itime(i)-itime(i-1), amsec=stepms)
570  IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
571  msteps = stepms/1000_int_ll
572  IF (mod(steps, msteps) == 0) mstepms = stepms
573  ENDIF
574  ENDDO
575  msteps = mstepms/1000_int_ll
576 
577  tmptime = lstart + step
578  DO WHILE(tmptime < lend) ! < since lend has been extended
579  CALL insert_unique(a_otime, tmptime)
580  tmptime = tmptime + step
581  ENDDO
582 
583 ! in next loop, we used initially steps instead of msteps (see comment
584 ! above), this gave much less combinations of time/timeranges with
585 ! possible empty output; we start from msteps instead of steps in
586 ! order to include partial initial intervals in case frac_valid<1;
587 ! however, as a gemeral rule, for aggregation of forecasts it's better
588 ! to use reference time
589  DO p1 = msteps, maxp1, msteps
590  CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
591  ENDDO
592 
593  ENDIF
594 
595 ELSE ! analysis mode
596  tmptime = lstart + step
597  DO WHILE(tmptime < lend) ! < since lend has been extended
598  CALL insert_unique(a_otime, tmptime)
599  tmptime = tmptime + step
600  ENDDO
601  CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
602 
603 ENDIF
604 
605 CALL packarray(a_otime)
606 CALL packarray(a_otimerange)
607 otime => a_otime%array
608 otimerange => a_otimerange%array
609 CALL sort(otime)
610 CALL sort(otimerange)
611 ! delete local objects keeping the contents
612 CALL delete(a_otime, nodealloc=.true.)
613 CALL delete(a_otimerange, nodealloc=.true.)
614 
615 #ifdef DEBUG
616 CALL l4f_log(l4f_debug, &
617  'recompute_stat_proc_agg, output time and timerange: '//&
618  t2c(SIZE(otime))//', '//t2c(size(otimerange)))
619 #endif
620 
621 IF (PRESENT(dtratio)) THEN
622 ! count the possible i/o interval ratios
623  DO k = 1, SIZE(itimerange)
624  IF (itimerange(k)%p2 /= 0) &
625  CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
626  ENDDO
627  CALL packarray(a_dtratio)
628  dtratio => a_dtratio%array
629  CALL sort(dtratio)
630 ! delete local object keeping the contents
631  CALL delete(a_dtratio, nodealloc=.true.)
632 
633 #ifdef DEBUG
634  CALL l4f_log(l4f_debug, &
635  'recompute_stat_proc_agg, found '//t2c(size(dtratio))// &
636  ' possible aggregation ratios, from '// &
637  t2c(dtratio(1))//' to '//t2c(dtratio(SIZE(dtratio))))
638 #endif
639 
640  ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
641  do_itimerange1: DO l = 1, SIZE(itimerange)
642  IF (.NOT.mask_timerange(l)) cycle do_itimerange1
643  do_itime1: DO k = 1, SIZE(itime)
644  CALL time_timerange_get_period(itime(k), itimerange(l), &
645  time_definition, pstart1, pend1, reftime1)
646  do_otimerange1: DO j = 1, SIZE(otimerange)
647  do_otime1: DO i = 1, SIZE(otime)
648  CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
649  time_definition, pstart2, pend2, reftime2)
650  IF (lforecast) THEN
651  IF (reftime1 /= reftime2) cycle do_otime1
652  ENDIF
653 
654  IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
655  mod(pstart1-pstart2, pend1-pstart1) == timedelta_0) THEN ! useful
656  lmapper%it = k
657  lmapper%itr = l
658  lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
659  n = append(map_ttr(i,j), lmapper)
660  cycle do_itime1 ! can contribute only to a single interval
661  ENDIF
662  ENDDO do_otime1
663  ENDDO do_otimerange1
664  ENDDO do_itime1
665  ENDDO do_itimerange1
666 
667 ELSE
668 
669  ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
670  do_itimerange2: DO l = 1, SIZE(itimerange)
671  IF (.NOT.mask_timerange(l)) cycle do_itimerange2
672  do_itime2: DO k = 1, SIZE(itime)
673  CALL time_timerange_get_period(itime(k), itimerange(l), &
674  time_definition, pstart1, pend1, reftime1)
675  do_otimerange2: DO j = 1, SIZE(otimerange)
676  do_otime2: DO i = 1, SIZE(otime)
677  CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
678  time_definition, pstart2, pend2, reftime2)
679  IF (lforecast) THEN
680  IF (reftime1 /= reftime2) cycle do_otime2
681  ENDIF
682 
683  IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
684  IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
685  lmapper%it = k
686  lmapper%itr = l
687  IF (pstart1 == pstart2) THEN
688  lmapper%extra_info = 1 ! start of interval
689  ELSE IF (pend1 == pend2) THEN
690  lmapper%extra_info = 2 ! end of interval
691  ELSE
692  lmapper%extra_info = imiss
693  ENDIF
694  lmapper%time = pstart1 ! = pend1, order by time?
695  n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
696 ! no CYCLE, a single input can contribute to multiple output intervals
697  ENDIF
698  ENDDO do_otime2
699  ENDDO do_otimerange2
700  ENDDO do_itime2
701  ENDDO do_itimerange2
702 
703 ENDIF
704 
705 END SUBROUTINE recompute_stat_proc_agg_common
706 
707 
708 SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
709  max_step, weights)
710 TYPE(datetime),INTENT(in) :: vertime(:)
711 TYPE(datetime),INTENT(in) :: pstart
712 TYPE(datetime),INTENT(in) :: pend
713 LOGICAL,INTENT(in) :: time_mask(:)
714 TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
715 DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
716 
717 INTEGER :: i, nt
718 TYPE(datetime),ALLOCATABLE :: lvertime(:)
719 TYPE(datetime) :: half, nexthalf
720 INTEGER(kind=int_ll) :: dt, tdt
721 
722 nt = count(time_mask)
723 ALLOCATE(lvertime(nt))
724 lvertime = pack(vertime, mask=time_mask)
725 
726 IF (PRESENT(max_step)) THEN
727 ! new way
728 ! max_step = timedelta_0
729 ! DO i = 1, nt - 1
730 ! IF (lvertime(i+1) - lvertime(i) > max_step) &
731 ! max_step = lvertime(i+1) - lvertime(i)
732 ! ENDDO
733 ! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
734 ! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
735 ! old way
736  IF (nt == 1) THEN
737  max_step = pend - pstart
738  ELSE
739  half = lvertime(1) + (lvertime(2) - lvertime(1))/2
740  max_step = half - pstart
741  DO i = 2, nt - 1
742  nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
743  IF (nexthalf - half > max_step) max_step = nexthalf - half
744  half = nexthalf
745  ENDDO
746  IF (pend - half > max_step) max_step = pend - half
747  ENDIF
748 
749 ENDIF
750 
751 IF (PRESENT(weights)) THEN
752  IF (nt == 1) THEN
753  weights(1) = 1.0
754  ELSE
755  CALL getval(pend - pstart, amsec=tdt)
756  half = lvertime(1) + (lvertime(2) - lvertime(1))/2
757  CALL getval(half - pstart, amsec=dt)
758  weights(1) = dble(dt)/dble(tdt)
759  DO i = 2, nt - 1
760  nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
761  CALL getval(nexthalf - half, amsec=dt)
762  weights(i) = dble(dt)/dble(tdt)
763  half = nexthalf
764  ENDDO
765  CALL getval(pend - half, amsec=dt)
766  weights(nt) = dble(dt)/dble(tdt)
767  ENDIF
768 ENDIF
769 
770 END SUBROUTINE compute_stat_proc_agg_sw
771 
772 ! get start of period, end of period and reference time from time,
773 ! timerange, according to time_definition.
774 SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
775  pstart, pend, reftime)
776 TYPE(datetime),INTENT(in) :: time
777 TYPE(vol7d_timerange),INTENT(in) :: timerange
778 INTEGER,INTENT(in) :: time_definition
779 TYPE(datetime),INTENT(out) :: reftime
780 TYPE(datetime),INTENT(out) :: pstart
781 TYPE(datetime),INTENT(out) :: pend
782 
783 TYPE(timedelta) :: p1, p2
784 
785 
786 p1 = timedelta_new(sec=timerange%p1) ! end of period
787 p2 = timedelta_new(sec=timerange%p2) ! length of period
788 
789 IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
790 ! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
791  timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
792  pstart = datetime_miss
793  pend = datetime_miss
794  reftime = datetime_miss
795  RETURN
796 ENDIF
797 
798 IF (time_definition == 0) THEN ! time == reference time
799  reftime = time
800  pend = time + p1
801  pstart = pend - p2
802 ELSE IF (time_definition == 1) THEN ! time == verification time
803  pend = time
804  pstart = time - p2
805  reftime = time - p1
806 ELSE
807  pstart = datetime_miss
808  pend = datetime_miss
809  reftime = datetime_miss
810 ENDIF
811 
812 END SUBROUTINE time_timerange_get_period
813 
814 
815 ! get start of period, end of period and reference time from time,
816 ! timerange, according to time_definition. step is taken separately
817 ! from timerange and can be "popular"
818 SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
819  pstart, pend, reftime)
820 TYPE(datetime),INTENT(in) :: time
821 TYPE(vol7d_timerange),INTENT(in) :: timerange
822 TYPE(timedelta),INTENT(in) :: step
823 INTEGER,INTENT(in) :: time_definition
824 TYPE(datetime),INTENT(out) :: reftime
825 TYPE(datetime),INTENT(out) :: pstart
826 TYPE(datetime),INTENT(out) :: pend
827 
828 TYPE(timedelta) :: p1
829 
830 
831 p1 = timedelta_new(sec=timerange%p1) ! end of period
832 
833 IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
834 ! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
835  timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
836  pstart = datetime_miss
837  pend = datetime_miss
838  reftime = datetime_miss
839  RETURN
840 ENDIF
841 
842 IF (time_definition == 0) THEN ! time == reference time
843  reftime = time
844  pend = time + p1
845  pstart = pend - step
846 ELSE IF (time_definition == 1) THEN ! time == verification time
847  pend = time
848  pstart = time - step
849  reftime = time - p1
850 ELSE
851  pstart = datetime_miss
852  pend = datetime_miss
853  reftime = datetime_miss
854 ENDIF
855 
856 END SUBROUTINE time_timerange_get_period_pop
857 
858 
859 ! set time, timerange%p1, timerange%p2 according to pstart, pend,
860 ! reftime and time_definition.
861 SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
862  pstart, pend, reftime)
863 TYPE(datetime),INTENT(out) :: time
864 TYPE(vol7d_timerange),INTENT(inout) :: timerange
865 INTEGER,INTENT(in) :: time_definition
866 TYPE(datetime),INTENT(in) :: reftime
867 TYPE(datetime),INTENT(in) :: pstart
868 TYPE(datetime),INTENT(in) :: pend
869 
870 TYPE(timedelta) :: p1, p2
871 INTEGER(kind=int_ll) :: dmsec
872 
873 
874 IF (time_definition == 0) THEN ! time == reference time
875  time = reftime
876  p1 = pend - reftime
877  p2 = pend - pstart
878 ELSE IF (time_definition == 1) THEN ! time == verification time
879  time = pend
880  p1 = pend - reftime
881  p2 = pend - pstart
882 ELSE
883  time = datetime_miss
884 ENDIF
885 
886 IF (time /= datetime_miss) THEN
887  CALL getval(p1, amsec=dmsec) ! end of period
888  timerange%p1 = int(dmsec/1000_int_ll)
889  CALL getval(p2, amsec=dmsec) ! length of period
890  timerange%p2 = int(dmsec/1000_int_ll)
891 ELSE
892  timerange%p1 = imiss
893  timerange%p2 = imiss
894 ENDIF
895 
896 END SUBROUTINE time_timerange_set_period
897 
898 
899 END MODULE stat_proc_engine
Restituiscono il valore dell'oggetto nella forma desiderata.
Operatore di resto della divisione.
Functions that return a trimmed CHARACTER representation of the input variable.
Quick method to append an element to the array.
Destructor for finalizing an array object.
Method for inserting elements of the array at a desired position.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Classi per la gestione delle coordinate temporali.
This module contains functions that are only for internal use of the library.
Classe per la gestione di un volume completo di dati osservati.
Class for expressing an absolute time value.
Class for expressing a relative time interval.

Generated with Doxygen.