libsim  Versione 7.1.6
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, deltas
151 LOGICAL :: lfull_steps, useful
152 TYPE(datetime) :: 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 
157 ! compute length of cumulation step in seconds
158 CALL getval(step, asec=steps)
159 
160 deltas = 0
161 IF (PRESENT(start)) THEN
162  IF (SIZE(itime) > 0 .AND. c_e(start)) THEN ! security check
163  CALL getval(start-itime(1), asec=deltas)
164  ENDIF
165 ENDIF
166 
167 lfull_steps = optio_log(full_steps)
168 
169 ! create a mask of suitable timeranges
170 ALLOCATE(mask_timerange(SIZE(itimerange)))
171 mask_timerange(:) = itimerange(:)%timerange == stat_proc &
172  .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
173  .AND. itimerange(:)%p1 >= 0 &
174  .AND. itimerange(:)%p2 > 0
175 
176 IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
177  mask_timerange(:) = mask_timerange(:) .AND. &
178  (itimerange(:)%p1 == 0 .OR. & ! all analyses
179  mod(itimerange(:)%p1, steps) == 0 .OR. & ! end time is mod step
180  mod(itimerange(:)%p1 - itimerange(:)%p2, steps) == 0) ! start time is mod step
181 ENDIF
182 ! mask_timerange includes all candidate timeranges
183 
184 nitr = count(mask_timerange)
185 ALLOCATE(f(nitr))
186 j = 1
187 DO i = 1, nitr
188  DO WHILE(.NOT.mask_timerange(j))
189  j = j + 1
190  ENDDO
191  f(i) = j
192  j = j + 1
193 ENDDO
194 
195 ! now we have to evaluate time/timerage pairs which do not need processing
196 ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
197 CALL compute_keep_tr()
198 
199 ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
200 map_tr(:,:,:,:,:) = imiss
201 
202 ! scan through all possible combinations of time and timerange
203 DO dirtyrep = 1, 2
204  IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
205  CALL packarray(a_otime)
206  CALL packarray(a_otimerange)
207  CALL sort(a_otime%array)
208  CALL sort(a_otimerange%array)
209  ENDIF
210  DO l = 1, SIZE(itime)
211  DO k = 1, nitr
212  CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
213  time_definition, pstart2, pend2, reftime2)
214 
215  DO j = 1, SIZE(itime)
216  DO i = 1, nitr
217  useful = .false.
218  CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
219  time_definition, pstart1, pend1, reftime1)
220  tmptimerange = vol7d_timerange_new(timerange=stat_proc)
221 
222  IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
223  IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
224  CALL time_timerange_set_period(tmptime, tmptimerange, &
225  time_definition, pend1, pend2, reftime2)
226  IF (lfull_steps) THEN
227  IF (mod(reftime2, step) == timedelta_0) THEN
228  useful = .true.
229  ENDIF
230  ELSE
231  useful = .true.
232  ENDIF
233 
234  ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
235  CALL time_timerange_set_period(tmptime, tmptimerange, &
236  time_definition, pstart2, pstart1, pstart1)
237  IF (lfull_steps) THEN
238  IF (mod(pstart1, step) == timedelta_0) THEN
239  useful = .true.
240  ENDIF
241  ELSE
242  useful = .true.
243  ENDIF
244  ENDIF
245 
246  ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
247  IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
248  CALL time_timerange_set_period(tmptime, tmptimerange, &
249  time_definition, pend1, pend2, reftime2)
250  IF (lfull_steps) THEN
251  IF (mod(pend2-reftime2, step) == timedelta_0) THEN
252  useful = .true.
253  ENDIF
254  ELSE
255  useful = .true.
256  ENDIF
257 
258  ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
259  CALL time_timerange_set_period(tmptime, tmptimerange, &
260  time_definition, pstart2, pstart1, reftime2)
261  IF (lfull_steps) THEN
262  IF (mod(pstart1-reftime2, step) == timedelta_0) THEN
263  useful = .true.
264  ENDIF
265  ELSE
266  useful = .true.
267  ENDIF
268 
269  ENDIF
270  ENDIF
271  useful = useful .AND. tmptime /= datetime_miss .AND. &
272  tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
273 
274  IF (useful) THEN ! add a_otime, a_otimerange
275  map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
276  map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
277  ENDIF
278  ENDDO
279  ENDDO
280  ENDDO
281  ENDDO
282 ENDDO
283 
284 ! we have to repeat the computation with sorted arrays
285 CALL compute_keep_tr()
286 
287 otime => a_otime%array
288 otimerange => a_otimerange%array
289 ! delete local objects keeping the contents
290 CALL delete(a_otime, nodealloc=.true.)
291 CALL delete(a_otimerange, nodealloc=.true.)
292 
293 #ifdef DEBUG
294 CALL l4f_log(l4f_debug, &
295  'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr,1)))//', '// &
296  t2c((SIZE(map_tr,2)))//', '// &
297  t2c((SIZE(map_tr,3)))//', '// &
298  t2c((SIZE(map_tr,4))))
299 CALL l4f_log(l4f_debug, &
300  'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr))/2)//', '// &
301  t2c(count(c_e(map_tr))/2))
302 CALL l4f_log(l4f_debug, &
303  'recompute_stat_proc_diff, nitr: '//t2c(nitr))
304 CALL l4f_log(l4f_debug, &
305  'recompute_stat_proc_diff, good timeranges: '//t2c(count(c_e(keep_tr))/2))
306 CALL l4f_log(l4f_debug, &
307  'recompute_stat_proc_diff, output times: '//t2c(SIZE(otime)))
308 CALL l4f_log(l4f_debug, &
309  'recompute_stat_proc_diff, output timeranges: '//t2c(SIZE(otimerange)))
310 #endif
311 
312 CONTAINS
313 
314 SUBROUTINE compute_keep_tr()
315 
316 keep_tr(:,:,:) = imiss
317 DO l = 1, SIZE(itime)
318  DO k = 1, nitr
319  IF (itimerange(f(k))%p2 == steps) THEN
320  CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
321  time_definition, pstart2, pend2, reftime2)
322  useful = .false.
323  IF (reftime2 == pend2) THEN ! analysis
324  IF (lfull_steps) THEN
325  IF (mod(reftime2, step) == timedelta_0) THEN
326  useful = .true.
327  ENDIF
328  ELSE
329  useful = .true.
330  ENDIF
331  ELSE ! forecast
332  IF (lfull_steps) THEN
333  IF (mod(itimerange(f(k))%p1, steps) == 0) THEN
334  useful = .true.
335  ENDIF
336  ELSE
337  useful = .true.
338  ENDIF
339  ENDIF
340  IF (useful) THEN
341 ! CALL time_timerange_set_period(tmptime, tmptimerange, &
342 ! time_definition, pstart2, pend2, reftime2)
343  keep_tr(k,l,1) = append_unique(a_otime, itime(l))
344  keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
345  ENDIF
346  ENDIF
347  ENDDO
348 ENDDO
349 
350 END SUBROUTINE compute_keep_tr
351 
352 END SUBROUTINE recompute_stat_proc_diff_common
353 
354 
355 ! common operations for statistical processing by metamorphosis
356 SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
357  otimerange, map_tr)
358 INTEGER,INTENT(in) :: istat_proc
359 TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
360 INTEGER,INTENT(in) :: ostat_proc
361 TYPE(vol7d_timerange),POINTER :: otimerange(:)
362 INTEGER,POINTER :: map_tr(:)
363 
364 INTEGER :: i
365 LOGICAL :: tr_mask(SIZE(itimerange))
366 
367 IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
368  ALLOCATE(otimerange(0), map_tr(0))
369  RETURN
370 ENDIF
371 
372 ! useful timeranges
373 tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
374  .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
375 ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
376 
377 otimerange = pack(itimerange, mask=tr_mask)
378 otimerange(:)%timerange = ostat_proc
379 map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
380 
381 END SUBROUTINE compute_stat_proc_metamorph_common
382 
383 
384 ! common operations for statistical processing by aggregation
385 SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
386  step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
387 TYPE(datetime),INTENT(in) :: itime(:)
388 TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
389 INTEGER,INTENT(in) :: stat_proc
390 INTEGER,INTENT(in) :: tri
391 TYPE(timedelta),INTENT(in) :: step
392 INTEGER,INTENT(in) :: time_definition
393 TYPE(datetime),POINTER :: otime(:)
394 TYPE(vol7d_timerange),POINTER :: otimerange(:)
395 TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
396 INTEGER,POINTER,OPTIONAL :: dtratio(:)
397 TYPE(datetime),INTENT(in),OPTIONAL :: start
398 LOGICAL,INTENT(in),OPTIONAL :: full_steps
399 
400 INTEGER :: i, j, k, l, na, nf, n
401 INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
402 INTEGER(kind=int_ll) :: stepms, mstepms
403 LOGICAL :: lforecast
404 TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
405 TYPE(arrayof_datetime) :: a_otime
406 TYPE(arrayof_vol7d_timerange) :: a_otimerange
407 TYPE(arrayof_integer) :: a_dtratio
408 LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
409 TYPE(ttr_mapper) :: lmapper
410 CHARACTER(len=8) :: env_var
411 LOGICAL :: climat_behavior
412 
413 
414 ! enable bad behavior for climat database (only for instantaneous data)
415 env_var = ''
416 CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
417 climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
418 
419 ! compute length of cumulation step in seconds
420 CALL getval(timedelta_depop(step), asec=steps)
421 
422 ! create a mask of suitable timeranges
423 ALLOCATE(mask_timerange(SIZE(itimerange)))
424 mask_timerange(:) = itimerange(:)%timerange == tri &
425  .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
426 
427 IF (PRESENT(dtratio)) THEN
428  WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
429  mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
430  ELSEWHERE
431  mask_timerange(:) = .false.
432  END WHERE
433 ELSE ! instantaneous
434  mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
435 ENDIF
436 
437 #ifdef DEBUG
438 CALL l4f_log(l4f_debug, &
439  '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
440  t2c(count(mask_timerange)))
441 #endif
442 
443 ! euristically determine whether we are dealing with an
444 ! analysis/observation or a forecast dataset
445 na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
446 nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
447 lforecast = nf >= na
448 #ifdef DEBUG
449 CALL l4f_log(l4f_debug, &
450  'recompute_stat_proc_agg, na: '//t2c(na)//', nf: '//t2c(nf))
451 #endif
452 
453 ! keep only timeranges of one type (really necessary?)
454 IF (lforecast) THEN
455  CALL l4f_log(l4f_info, &
456  'recompute_stat_proc_agg, processing in forecast mode')
457 ELSE
458  mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
459  CALL l4f_log(l4f_info, &
460  'recompute_stat_proc_agg, processing in analysis mode')
461 ENDIF
462 
463 #ifdef DEBUG
464 CALL l4f_log(l4f_debug, &
465  '(re)compute_stat_proc_agg, number of useful timeranges: '// &
466  t2c(count(mask_timerange)))
467 #endif
468 
469 IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
470  ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
471  IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
472  RETURN
473 ENDIF
474 
475 ! determine start and end of processing period, should work with p2==0
476 lstart = datetime_miss
477 IF (PRESENT(start)) lstart = start
478 lend = itime(SIZE(itime))
479 ! compute some quantities
480 maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
481 maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
482 minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
483 IF (time_definition == 0) THEN ! reference time
484  lend = lend + timedelta_new(sec=maxp1)
485 ENDIF
486 ! extend interval at the end in order to include all the data in case
487 ! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
488 ! below in order to exclude the last full interval which would be empty
489 lend = lend + step
490 IF (lstart == datetime_miss) THEN ! autodetect
491  lstart = itime(1)
492 ! if autodetected, adjust to obtain real absolute start of data
493  IF (time_definition == 0) THEN ! reference time
494  lstart = lstart + timedelta_new(sec=minp1mp2)
495  ELSE ! verification time
496 ! go back to start of longest processing interval
497  lstart = lstart - timedelta_new(sec=maxp2)
498  ENDIF
499  IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
500  lstart = lstart - (mod(lstart, step)) ! round to step, (should be MODULO, not MOD)
501  ENDIF
502 ENDIF
503 
504 #ifdef DEBUG
505 CALL l4f_log(l4f_debug, &
506  'recompute_stat_proc_agg, processing period: '//t2c(lstart)//' - '//t2c(lend))
507 #endif
508 
509 ! create output time and timerange lists
510 
511 IF (lforecast) THEN ! forecast mode
512  IF (time_definition == 0) THEN ! reference time
513  CALL insert(a_otime, itime) ! should I limit to elements itime >= lstart?
514 
515 ! apply start shift to timerange, not to reftime
516 ! why did we use itime(SIZE(itime)) (last ref time)?
517 ! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
518  CALL getval(lstart-itime(1), asec=dstart)
519 ! allow starting before first reftime but restrict dtstart to -steps
520 ! dstart = MAX(0, dstart)
521  IF (dstart < 0) dstart = mod(dstart, steps)
522  DO p1 = steps + dstart, maxp1, steps
523  CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
524  ENDDO
525 
526  ELSE ! verification time
527 
528 ! verification time in forecast mode is the ugliest case, because we
529 ! don't know a priori how many different (thus incompatible) reference
530 ! times we have, so some assumption of regularity has to be made. For
531 ! this reason msteps, the minimum step between two times, is
532 ! computed. We choose to compute it as a difference between itime
533 ! elements, it could be also computed as difference of itimerange%p1
534 ! elements. But what if it is not modulo steps?
535  mstepms = steps*1000_int_ll
536  DO i = 2, SIZE(itime)
537  CALL getval(itime(i)-itime(i-1), amsec=stepms)
538  IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
539  msteps = stepms/1000_int_ll
540  IF (mod(steps, msteps) == 0) mstepms = stepms
541  ENDIF
542  ENDDO
543  msteps = mstepms/1000_int_ll
544 
545  tmptime = lstart + step
546  DO WHILE(tmptime < lend) ! < since lend has been extended
547  CALL insert_unique(a_otime, tmptime)
548  tmptime = tmptime + step
549  ENDDO
550 
551 ! in next loop, we used initially steps instead of msteps (see comment
552 ! above), this gave much less combinations of time/timeranges with
553 ! possible empty output; we start from msteps instead of steps in
554 ! order to include partial initial intervals in case frac_valid<1;
555 ! however, as a gemeral rule, for aggregation of forecasts it's better
556 ! to use reference time
557  DO p1 = msteps, maxp1, msteps
558  CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
559  ENDDO
560 
561  ENDIF
562 
563 ELSE ! analysis mode
564  tmptime = lstart + step
565  DO WHILE(tmptime < lend) ! < since lend has been extended
566  CALL insert_unique(a_otime, tmptime)
567  tmptime = tmptime + step
568  ENDDO
569  CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
570 
571 ENDIF
572 
573 CALL packarray(a_otime)
574 CALL packarray(a_otimerange)
575 otime => a_otime%array
576 otimerange => a_otimerange%array
577 CALL sort(otime)
578 CALL sort(otimerange)
579 ! delete local objects keeping the contents
580 CALL delete(a_otime, nodealloc=.true.)
581 CALL delete(a_otimerange, nodealloc=.true.)
582 
583 #ifdef DEBUG
584 CALL l4f_log(l4f_debug, &
585  'recompute_stat_proc_agg, output time and timerange: '//&
586  t2c(SIZE(otime))//', '//t2c(size(otimerange)))
587 #endif
588 
589 IF (PRESENT(dtratio)) THEN
590 ! count the possible i/o interval ratios
591  DO k = 1, SIZE(itimerange)
592  IF (itimerange(k)%p2 /= 0) &
593  CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
594  ENDDO
595  CALL packarray(a_dtratio)
596  dtratio => a_dtratio%array
597  CALL sort(dtratio)
598 ! delete local object keeping the contents
599  CALL delete(a_dtratio, nodealloc=.true.)
600 
601 #ifdef DEBUG
602  CALL l4f_log(l4f_debug, &
603  'recompute_stat_proc_agg, found '//t2c(size(dtratio))// &
604  ' possible aggregation ratios, from '// &
605  t2c(dtratio(1))//' to '//t2c(dtratio(SIZE(dtratio))))
606 #endif
607 
608  ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
609  do_itimerange1: DO l = 1, SIZE(itimerange)
610  IF (.NOT.mask_timerange(l)) cycle do_itimerange1
611  do_itime1: DO k = 1, SIZE(itime)
612  CALL time_timerange_get_period(itime(k), itimerange(l), &
613  time_definition, pstart1, pend1, reftime1)
614  do_otimerange1: DO j = 1, SIZE(otimerange)
615  do_otime1: DO i = 1, SIZE(otime)
616  CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
617  time_definition, pstart2, pend2, reftime2)
618  IF (lforecast) THEN
619  IF (reftime1 /= reftime2) cycle do_otime1
620  ENDIF
621 
622  IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
623  mod(pstart1-pstart2, pend1-pstart1) == timedelta_0) THEN ! useful
624  lmapper%it = k
625  lmapper%itr = l
626  lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
627  n = append(map_ttr(i,j), lmapper)
628  cycle do_itime1 ! can contribute only to a single interval
629  ENDIF
630  ENDDO do_otime1
631  ENDDO do_otimerange1
632  ENDDO do_itime1
633  ENDDO do_itimerange1
634 
635 ELSE
636 
637  ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
638  do_itimerange2: DO l = 1, SIZE(itimerange)
639  IF (.NOT.mask_timerange(l)) cycle do_itimerange2
640  do_itime2: DO k = 1, SIZE(itime)
641  CALL time_timerange_get_period(itime(k), itimerange(l), &
642  time_definition, pstart1, pend1, reftime1)
643  do_otimerange2: DO j = 1, SIZE(otimerange)
644  do_otime2: DO i = 1, SIZE(otime)
645  CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
646  time_definition, pstart2, pend2, reftime2)
647  IF (lforecast) THEN
648  IF (reftime1 /= reftime2) cycle do_otime2
649  ENDIF
650 
651  IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
652  IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
653  lmapper%it = k
654  lmapper%itr = l
655  IF (pstart1 == pstart2) THEN
656  lmapper%extra_info = 1 ! start of interval
657  ELSE IF (pend1 == pend2) THEN
658  lmapper%extra_info = 2 ! end of interval
659  ELSE
660  lmapper%extra_info = imiss
661  ENDIF
662  lmapper%time = pstart1 ! = pend1, order by time?
663  n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
664 ! no CYCLE, a single input can contribute to multiple output intervals
665  ENDIF
666  ENDDO do_otime2
667  ENDDO do_otimerange2
668  ENDDO do_itime2
669  ENDDO do_itimerange2
670 
671 ENDIF
672 
673 END SUBROUTINE recompute_stat_proc_agg_common
674 
675 
676 SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
677  max_step, weights)
678 TYPE(datetime),INTENT(in) :: vertime(:)
679 TYPE(datetime),INTENT(in) :: pstart
680 TYPE(datetime),INTENT(in) :: pend
681 LOGICAL,INTENT(in) :: time_mask(:)
682 TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
683 DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
684 
685 INTEGER :: i, nt
686 TYPE(datetime),ALLOCATABLE :: lvertime(:)
687 TYPE(datetime) :: half, nexthalf
688 INTEGER(kind=int_ll) :: dt, tdt
689 
690 nt = count(time_mask)
691 ALLOCATE(lvertime(nt))
692 lvertime = pack(vertime, mask=time_mask)
693 
694 IF (PRESENT(max_step)) THEN
695 ! new way
696 ! max_step = timedelta_0
697 ! DO i = 1, nt - 1
698 ! IF (lvertime(i+1) - lvertime(i) > max_step) &
699 ! max_step = lvertime(i+1) - lvertime(i)
700 ! ENDDO
701 ! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
702 ! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
703 ! old way
704  IF (nt == 1) THEN
705  max_step = pend - pstart
706  ELSE
707  half = lvertime(1) + (lvertime(2) - lvertime(1))/2
708  max_step = half - pstart
709  DO i = 2, nt - 1
710  nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
711  IF (nexthalf - half > max_step) max_step = nexthalf - half
712  half = nexthalf
713  ENDDO
714  IF (pend - half > max_step) max_step = pend - half
715  ENDIF
716 
717 ENDIF
718 
719 IF (PRESENT(weights)) THEN
720  IF (nt == 1) THEN
721  weights(1) = 1.0
722  ELSE
723  CALL getval(pend - pstart, amsec=tdt)
724  half = lvertime(1) + (lvertime(2) - lvertime(1))/2
725  CALL getval(half - pstart, amsec=dt)
726  weights(1) = dble(dt)/dble(tdt)
727  DO i = 2, nt - 1
728  nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
729  CALL getval(nexthalf - half, amsec=dt)
730  weights(i) = dble(dt)/dble(tdt)
731  half = nexthalf
732  ENDDO
733  CALL getval(pend - half, amsec=dt)
734  weights(nt) = dble(dt)/dble(tdt)
735  ENDIF
736 ENDIF
737 
738 END SUBROUTINE compute_stat_proc_agg_sw
739 
740 ! get start of period, end of period and reference time from time,
741 ! timerange, according to time_definition.
742 SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
743  pstart, pend, reftime)
744 TYPE(datetime),INTENT(in) :: time
745 TYPE(vol7d_timerange),INTENT(in) :: timerange
746 INTEGER,INTENT(in) :: time_definition
747 TYPE(datetime),INTENT(out) :: reftime
748 TYPE(datetime),INTENT(out) :: pstart
749 TYPE(datetime),INTENT(out) :: pend
750 
751 TYPE(timedelta) :: p1, p2
752 
753 
754 p1 = timedelta_new(sec=timerange%p1) ! end of period
755 p2 = timedelta_new(sec=timerange%p2) ! length of period
756 
757 IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
758 ! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
759  timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
760  pstart = datetime_miss
761  pend = datetime_miss
762  reftime = datetime_miss
763  RETURN
764 ENDIF
765 
766 IF (time_definition == 0) THEN ! time == reference time
767  reftime = time
768  pend = time + p1
769  pstart = pend - p2
770 ELSE IF (time_definition == 1) THEN ! time == verification time
771  pend = time
772  pstart = time - p2
773  reftime = time - p1
774 ELSE
775  pstart = datetime_miss
776  pend = datetime_miss
777  reftime = datetime_miss
778 ENDIF
779 
780 END SUBROUTINE time_timerange_get_period
781 
782 
783 ! get start of period, end of period and reference time from time,
784 ! timerange, according to time_definition. step is taken separately
785 ! from timerange and can be "popular"
786 SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
787  pstart, pend, reftime)
788 TYPE(datetime),INTENT(in) :: time
789 TYPE(vol7d_timerange),INTENT(in) :: timerange
790 TYPE(timedelta),INTENT(in) :: step
791 INTEGER,INTENT(in) :: time_definition
792 TYPE(datetime),INTENT(out) :: reftime
793 TYPE(datetime),INTENT(out) :: pstart
794 TYPE(datetime),INTENT(out) :: pend
795 
796 TYPE(timedelta) :: p1
797 
798 
799 p1 = timedelta_new(sec=timerange%p1) ! end of period
800 
801 IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
802 ! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
803  timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
804  pstart = datetime_miss
805  pend = datetime_miss
806  reftime = datetime_miss
807  RETURN
808 ENDIF
809 
810 IF (time_definition == 0) THEN ! time == reference time
811  reftime = time
812  pend = time + p1
813  pstart = pend - step
814 ELSE IF (time_definition == 1) THEN ! time == verification time
815  pend = time
816  pstart = time - step
817  reftime = time - p1
818 ELSE
819  pstart = datetime_miss
820  pend = datetime_miss
821  reftime = datetime_miss
822 ENDIF
823 
824 END SUBROUTINE time_timerange_get_period_pop
825 
826 
827 ! set time, timerange%p1, timerange%p2 according to pstart, pend,
828 ! reftime and time_definition.
829 SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
830  pstart, pend, reftime)
831 TYPE(datetime),INTENT(out) :: time
832 TYPE(vol7d_timerange),INTENT(inout) :: timerange
833 INTEGER,INTENT(in) :: time_definition
834 TYPE(datetime),INTENT(in) :: reftime
835 TYPE(datetime),INTENT(in) :: pstart
836 TYPE(datetime),INTENT(in) :: pend
837 
838 TYPE(timedelta) :: p1, p2
839 INTEGER(kind=int_ll) :: dmsec
840 
841 
842 IF (time_definition == 0) THEN ! time == reference time
843  time = reftime
844  p1 = pend - reftime
845  p2 = pend - pstart
846 ELSE IF (time_definition == 1) THEN ! time == verification time
847  time = pend
848  p1 = pend - reftime
849  p2 = pend - pstart
850 ELSE
851  time = datetime_miss
852 ENDIF
853 
854 IF (time /= datetime_miss) THEN
855  CALL getval(p1, amsec=dmsec) ! end of period
856  timerange%p1 = int(dmsec/1000_int_ll)
857  CALL getval(p2, amsec=dmsec) ! length of period
858  timerange%p2 = int(dmsec/1000_int_ll)
859 ELSE
860  timerange%p1 = imiss
861  timerange%p2 = imiss
862 ENDIF
863 
864 END SUBROUTINE time_timerange_set_period
865 
866 
867 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.