libsim Versione 7.2.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!> This module contains functions that are only for internal use of
20!! the library. It should not be used by user procedures because it is
21!! subject to change
22!! \ingroup vol7d
26IMPLICIT 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
30TYPE 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
35END TYPE ttr_mapper
36
37INTERFACE OPERATOR (==)
38 MODULE PROCEDURE ttr_mapper_eq
39END INTERFACE
40
41INTERFACE OPERATOR (/=)
42 MODULE PROCEDURE ttr_mapper_ne
43END INTERFACE
44
45INTERFACE OPERATOR (>)
46 MODULE PROCEDURE ttr_mapper_gt
47END INTERFACE
48
49INTERFACE OPERATOR (<)
50 MODULE PROCEDURE ttr_mapper_lt
51END INTERFACE
52
53INTERFACE OPERATOR (>=)
54 MODULE PROCEDURE ttr_mapper_ge
55END INTERFACE
56
57INTERFACE OPERATOR (<=)
58 MODULE PROCEDURE ttr_mapper_le
59END 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
77CONTAINS
78
79! simplified comparisons only on time
80ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res)
81TYPE(ttr_mapper),INTENT(IN) :: this, that
82LOGICAL :: res
83
84res = this%time == that%time
85
86END FUNCTION ttr_mapper_eq
87
88ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res)
89TYPE(ttr_mapper),INTENT(IN) :: this, that
90LOGICAL :: res
91
92res = this%time /= that%time
93
94END FUNCTION ttr_mapper_ne
95
96ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res)
97TYPE(ttr_mapper),INTENT(IN) :: this, that
98LOGICAL :: res
99
100res = this%time > that%time
101
102END FUNCTION ttr_mapper_gt
103
104ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res)
105TYPE(ttr_mapper),INTENT(IN) :: this, that
106LOGICAL :: res
107
108res = this%time < that%time
109
110END FUNCTION ttr_mapper_lt
111
112ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res)
113TYPE(ttr_mapper),INTENT(IN) :: this, that
114LOGICAL :: res
115
116res = this%time >= that%time
117
118END FUNCTION ttr_mapper_ge
119
120ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res)
121TYPE(ttr_mapper),INTENT(IN) :: this, that
122LOGICAL :: res
123
124res = this%time <= that%time
125
126END 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
133SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, &
134 otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, &
135 start)
136TYPE(datetime),INTENT(in) :: itime(:)
137TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
138INTEGER,INTENT(in) :: stat_proc
139TYPE(timedelta),INTENT(in) :: step
140TYPE(datetime),POINTER :: otime(:)
141TYPE(vol7d_timerange),POINTER :: otimerange(:)
142INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
143INTEGER :: nitr
144LOGICAL,ALLOCATABLE :: mask_timerange(:)
145INTEGER,INTENT(in) :: time_definition
146LOGICAL,INTENT(in),OPTIONAL :: full_steps
147TYPE(datetime),INTENT(in),OPTIONAL :: start
148
149INTEGER :: i, j, k, l, dirtyrep
150INTEGER :: steps
151LOGICAL :: lfull_steps, useful
152TYPE(datetime) :: lstart, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
153TYPE(vol7d_timerange) :: tmptimerange
154TYPE(arrayof_datetime) :: a_otime
155TYPE(arrayof_vol7d_timerange) :: a_otimerange
156TYPE(timedelta) :: start_delta
157
158! compute length of cumulation step in seconds
159CALL getval(step, asec=steps)
160
161lstart = datetime_miss
162IF (PRESENT(start)) lstart = start
163lfull_steps = optio_log(full_steps)
164
165! create a mask of suitable timeranges
166ALLOCATE(mask_timerange(SIZE(itimerange)))
167mask_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
172IF (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
177ENDIF
178! mask_timerange includes all candidate timeranges
179
180nitr = count(mask_timerange)
181ALLOCATE(f(nitr))
182j = 1
183DO i = 1, nitr
184 DO WHILE(.NOT.mask_timerange(j))
185 j = j + 1
186 ENDDO
187 f(i) = j
188 j = j + 1
189ENDDO
190
191! now we have to evaluate time/timerage pairs which do not need processing
192ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
193CALL compute_keep_tr()
194
195ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
196map_tr(:,:,:,:,:) = imiss
197
198! scan through all possible combinations of time and timerange
199DO 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
296ENDDO
297
298! we have to repeat the computation with sorted arrays
299CALL compute_keep_tr()
301otime => a_otime%array
302otimerange => a_otimerange%array
303! delete local objects keeping the contents
304CALL delete(a_otime, nodealloc=.true.)
305CALL delete(a_otimerange, nodealloc=.true.)
306
307#ifdef DEBUG
308CALL 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))))
313CALL l4f_log(l4f_debug, &
314 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr))/2)//', '// &
315 t2c(count(c_e(map_tr))/2))
316CALL l4f_log(l4f_debug, &
317 'recompute_stat_proc_diff, nitr: '//t2c(nitr))
318CALL l4f_log(l4f_debug, &
319 'recompute_stat_proc_diff, good timeranges: '//t2c(count(c_e(keep_tr))/2))
320CALL l4f_log(l4f_debug, &
321 'recompute_stat_proc_diff, output times: '//t2c(SIZE(otime)))
322CALL l4f_log(l4f_debug, &
323 'recompute_stat_proc_diff, output timeranges: '//t2c(SIZE(otimerange)))
324#endif
325
326CONTAINS
327
328SUBROUTINE compute_keep_tr()
329INTEGER :: start_deltas
330
331keep_tr(:,:,:) = imiss
332DO 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
377ENDDO
379END SUBROUTINE compute_keep_tr
380
381END SUBROUTINE recompute_stat_proc_diff_common
382
383
384! common operations for statistical processing by metamorphosis
385SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
386 otimerange, map_tr)
387INTEGER,INTENT(in) :: istat_proc
388TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
389INTEGER,INTENT(in) :: ostat_proc
390TYPE(vol7d_timerange),POINTER :: otimerange(:)
391INTEGER,POINTER :: map_tr(:)
392
393INTEGER :: i
394LOGICAL :: tr_mask(SIZE(itimerange))
395
396IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
397 ALLOCATE(otimerange(0), map_tr(0))
398 RETURN
399ENDIF
400
401! useful timeranges
402tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
403 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
404ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
405
406otimerange = pack(itimerange, mask=tr_mask)
407otimerange(:)%timerange = ostat_proc
408map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
409
410END SUBROUTINE compute_stat_proc_metamorph_common
411
412
413! common operations for statistical processing by aggregation
414SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
415 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
416TYPE(datetime),INTENT(in) :: itime(:)
417TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
418INTEGER,INTENT(in) :: stat_proc
419INTEGER,INTENT(in) :: tri
420TYPE(timedelta),INTENT(in) :: step
421INTEGER,INTENT(in) :: time_definition
422TYPE(datetime),POINTER :: otime(:)
423TYPE(vol7d_timerange),POINTER :: otimerange(:)
424TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
425INTEGER,POINTER,OPTIONAL :: dtratio(:)
426TYPE(datetime),INTENT(in),OPTIONAL :: start
427LOGICAL,INTENT(in),OPTIONAL :: full_steps
428
429INTEGER :: i, j, k, l, na, nf, n
430INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
431INTEGER(kind=int_ll) :: stepms, mstepms
432LOGICAL :: lforecast
433TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
434TYPE(arrayof_datetime) :: a_otime
435TYPE(arrayof_vol7d_timerange) :: a_otimerange
436TYPE(arrayof_integer) :: a_dtratio
437LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
438TYPE(ttr_mapper) :: lmapper
439CHARACTER(len=8) :: env_var
440LOGICAL :: climat_behavior
441
442
443! enable bad behavior for climat database (only for instantaneous data)
444env_var = ''
445CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
446climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
447
448! compute length of cumulation step in seconds
449CALL getval(timedelta_depop(step), asec=steps)
450
451! create a mask of suitable timeranges
452ALLOCATE(mask_timerange(SIZE(itimerange)))
453mask_timerange(:) = itimerange(:)%timerange == tri &
454 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
455
456IF (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
462ELSE ! instantaneous
463 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
464ENDIF
465
466#ifdef DEBUG
467CALL 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
474na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
475nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
476lforecast = nf >= na
477#ifdef DEBUG
478CALL 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?)
483IF (lforecast) THEN
484 CALL l4f_log(l4f_info, &
485 'recompute_stat_proc_agg, processing in forecast mode')
486ELSE
487 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
488 CALL l4f_log(l4f_info, &
489 'recompute_stat_proc_agg, processing in analysis mode')
490ENDIF
491
492#ifdef DEBUG
493CALL l4f_log(l4f_debug, &
494 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
495 t2c(count(mask_timerange)))
496#endif
497
498IF (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
502ENDIF
503
504! determine start and end of processing period, should work with p2==0
505lstart = datetime_miss
506IF (PRESENT(start)) lstart = start
507lend = itime(SIZE(itime))
508! compute some quantities
509maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
510maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
511minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
512IF (time_definition == 0) THEN ! reference time
513 lend = lend + timedelta_new(sec=maxp1)
514ENDIF
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
518lend = lend + step
519IF (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
534ENDIF
535
536#ifdef DEBUG
537CALL 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
543IF (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
595ELSE ! 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))
603ENDIF
604
605CALL packarray(a_otime)
606CALL packarray(a_otimerange)
607otime => a_otime%array
608otimerange => a_otimerange%array
609CALL sort(otime)
610CALL sort(otimerange)
611! delete local objects keeping the contents
612CALL delete(a_otime, nodealloc=.true.)
613CALL delete(a_otimerange, nodealloc=.true.)
614
615#ifdef DEBUG
616CALL l4f_log(l4f_debug, &
617 'recompute_stat_proc_agg, output time and timerange: '//&
618 t2c(SIZE(otime))//', '//t2c(size(otimerange)))
619#endif
620
621IF (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
667ELSE
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
703ENDIF
704
705END SUBROUTINE recompute_stat_proc_agg_common
706
707
708SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
709 max_step, weights)
710TYPE(datetime),INTENT(in) :: vertime(:)
711TYPE(datetime),INTENT(in) :: pstart
712TYPE(datetime),INTENT(in) :: pend
713LOGICAL,INTENT(in) :: time_mask(:)
714TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
715DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
716
717INTEGER :: i, nt
718TYPE(datetime),ALLOCATABLE :: lvertime(:)
719TYPE(datetime) :: half, nexthalf
720INTEGER(kind=int_ll) :: dt, tdt
721
722nt = count(time_mask)
723ALLOCATE(lvertime(nt))
724lvertime = pack(vertime, mask=time_mask)
725
726IF (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
749ENDIF
751IF (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
768ENDIF
769
770END 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.
774SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
775 pstart, pend, reftime)
776TYPE(datetime),INTENT(in) :: time
777TYPE(vol7d_timerange),INTENT(in) :: timerange
778INTEGER,INTENT(in) :: time_definition
779TYPE(datetime),INTENT(out) :: reftime
780TYPE(datetime),INTENT(out) :: pstart
781TYPE(datetime),INTENT(out) :: pend
782
783TYPE(timedelta) :: p1, p2
784
785
786p1 = timedelta_new(sec=timerange%p1) ! end of period
787p2 = timedelta_new(sec=timerange%p2) ! length of period
788
789IF (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
796ENDIF
797
798IF (time_definition == 0) THEN ! time == reference time
799 reftime = time
800 pend = time + p1
801 pstart = pend - p2
802ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
803 pend = time
804 pstart = time - p2
805 reftime = time - p1
806ELSE
807 pstart = datetime_miss
808 pend = datetime_miss
809 reftime = datetime_miss
810ENDIF
811
812END 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"
818SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
819 pstart, pend, reftime)
820TYPE(datetime),INTENT(in) :: time
821TYPE(vol7d_timerange),INTENT(in) :: timerange
822TYPE(timedelta),INTENT(in) :: step
823INTEGER,INTENT(in) :: time_definition
824TYPE(datetime),INTENT(out) :: reftime
825TYPE(datetime),INTENT(out) :: pstart
826TYPE(datetime),INTENT(out) :: pend
827
828TYPE(timedelta) :: p1
829
830
831p1 = timedelta_new(sec=timerange%p1) ! end of period
832
833IF (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
840ENDIF
841
842IF (time_definition == 0) THEN ! time == reference time
843 reftime = time
844 pend = time + p1
845 pstart = pend - step
846ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
847 pend = time
848 pstart = time - step
849 reftime = time - p1
850ELSE
851 pstart = datetime_miss
852 pend = datetime_miss
853 reftime = datetime_miss
854ENDIF
856END SUBROUTINE time_timerange_get_period_pop
857
858
859! set time, timerange%p1, timerange%p2 according to pstart, pend,
860! reftime and time_definition.
861SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
862 pstart, pend, reftime)
863TYPE(datetime),INTENT(out) :: time
864TYPE(vol7d_timerange),INTENT(inout) :: timerange
865INTEGER,INTENT(in) :: time_definition
866TYPE(datetime),INTENT(in) :: reftime
867TYPE(datetime),INTENT(in) :: pstart
868TYPE(datetime),INTENT(in) :: pend
869
870TYPE(timedelta) :: p1, p2
871INTEGER(kind=int_ll) :: dmsec
872
873
874IF (time_definition == 0) THEN ! time == reference time
875 time = reftime
876 p1 = pend - reftime
877 p2 = pend - pstart
878ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
879 time = pend
880 p1 = pend - reftime
881 p2 = pend - pstart
882ELSE
883 time = datetime_miss
884ENDIF
885
886IF (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)
891ELSE
892 timerange%p1 = imiss
893 timerange%p2 = imiss
894ENDIF
895
896END SUBROUTINE time_timerange_set_period
897
898
899END 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.
Index method.
Quick method to append an element to the array.
Destructor for finalizing an array object.
Index method with sorted array.
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...
Method for removing elements of the array at a desired position.
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.
Derived type defining a dynamically extensible array of TYPE(datetime) elements.
Class for expressing an absolute time value.
Class for expressing a relative time interval.
Derived type defining a dynamically extensible array of TYPE(ttr_mapper) elements.

Generated with Doxygen.