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