libsim Versione 7.2.0
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
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
156
157! compute length of cumulation step in seconds
158CALL getval(step, asec=steps)
159
160lstart = datetime_miss
161IF (PRESENT(start)) lstart = start
162lfull_steps = optio_log(full_steps)
163
164! create a mask of suitable timeranges
165ALLOCATE(mask_timerange(SIZE(itimerange)))
166mask_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
171IF (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
176ENDIF
177! mask_timerange includes all candidate timeranges
178
179nitr = count(mask_timerange)
180ALLOCATE(f(nitr))
181j = 1
182DO i = 1, nitr
183 DO WHILE(.NOT.mask_timerange(j))
184 j = j + 1
185 ENDDO
186 f(i) = j
187 j = j + 1
188ENDDO
189
190! now we have to evaluate time/timerage pairs which do not need processing
191ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
192CALL compute_keep_tr()
193
194ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
195map_tr(:,:,:,:,:) = imiss
196
197! scan through all possible combinations of time and timerange
198DO 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
277ENDDO
278
279! we have to repeat the computation with sorted arrays
280CALL compute_keep_tr()
281
282otime => a_otime%array
283otimerange => a_otimerange%array
284! delete local objects keeping the contents
285CALL delete(a_otime, nodealloc=.true.)
286CALL delete(a_otimerange, nodealloc=.true.)
287
288#ifdef DEBUG
289CALL 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))))
294CALL l4f_log(l4f_debug, &
295 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr))/2)//', '// &
296 t2c(count(c_e(map_tr))/2))
297CALL l4f_log(l4f_debug, &
298 'recompute_stat_proc_diff, nitr: '//t2c(nitr))
299CALL l4f_log(l4f_debug, &
300 'recompute_stat_proc_diff, good timeranges: '//t2c(count(c_e(keep_tr))/2))
301CALL l4f_log(l4f_debug, &
302 'recompute_stat_proc_diff, output times: '//t2c(SIZE(otime)))
303CALL l4f_log(l4f_debug, &
304 'recompute_stat_proc_diff, output timeranges: '//t2c(SIZE(otimerange)))
305#endif
306
307CONTAINS
308
309SUBROUTINE compute_keep_tr()
310
311keep_tr(:,:,:) = imiss
312DO 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
347ENDDO
349END SUBROUTINE compute_keep_tr
350
351END SUBROUTINE recompute_stat_proc_diff_common
352
353
354! common operations for statistical processing by metamorphosis
355SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
356 otimerange, map_tr)
357INTEGER,INTENT(in) :: istat_proc
358TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
359INTEGER,INTENT(in) :: ostat_proc
360TYPE(vol7d_timerange),POINTER :: otimerange(:)
361INTEGER,POINTER :: map_tr(:)
362
363INTEGER :: i
364LOGICAL :: tr_mask(SIZE(itimerange))
366IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
367 ALLOCATE(otimerange(0), map_tr(0))
368 RETURN
369ENDIF
370
371! useful timeranges
372tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
373 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
374ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
375
376otimerange = pack(itimerange, mask=tr_mask)
377otimerange(:)%timerange = ostat_proc
378map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
379
380END SUBROUTINE compute_stat_proc_metamorph_common
381
382
383! common operations for statistical processing by aggregation
384SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
385 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
386TYPE(datetime),INTENT(in) :: itime(:)
387TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
388INTEGER,INTENT(in) :: stat_proc
389INTEGER,INTENT(in) :: tri
390TYPE(timedelta),INTENT(in) :: step
391INTEGER,INTENT(in) :: time_definition
392TYPE(datetime),POINTER :: otime(:)
393TYPE(vol7d_timerange),POINTER :: otimerange(:)
394TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
395INTEGER,POINTER,OPTIONAL :: dtratio(:)
396TYPE(datetime),INTENT(in),OPTIONAL :: start
397LOGICAL,INTENT(in),OPTIONAL :: full_steps
398
399INTEGER :: i, j, k, l, na, nf, n
400INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
401INTEGER(kind=int_ll) :: stepms, mstepms
402LOGICAL :: lforecast
403TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
404TYPE(arrayof_datetime) :: a_otime
405TYPE(arrayof_vol7d_timerange) :: a_otimerange
406TYPE(arrayof_integer) :: a_dtratio
407LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
408TYPE(ttr_mapper) :: lmapper
409CHARACTER(len=8) :: env_var
410LOGICAL :: climat_behavior
411
412
413! enable bad behavior for climat database (only for instantaneous data)
414env_var = ''
415CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
416climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
417
418! compute length of cumulation step in seconds
419CALL getval(timedelta_depop(step), asec=steps)
420
421! create a mask of suitable timeranges
422ALLOCATE(mask_timerange(SIZE(itimerange)))
423mask_timerange(:) = itimerange(:)%timerange == tri &
424 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
425
426IF (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
432ELSE ! instantaneous
433 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
434ENDIF
435
436#ifdef DEBUG
437CALL 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
444na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
445nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
446lforecast = nf >= na
447#ifdef DEBUG
448CALL 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?)
453IF (lforecast) THEN
454 CALL l4f_log(l4f_info, &
455 'recompute_stat_proc_agg, processing in forecast mode')
456ELSE
457 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
458 CALL l4f_log(l4f_info, &
459 'recompute_stat_proc_agg, processing in analysis mode')
460ENDIF
461
462#ifdef DEBUG
463CALL l4f_log(l4f_debug, &
464 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
465 t2c(count(mask_timerange)))
466#endif
467
468IF (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
472ENDIF
473
474! determine start and end of processing period, should work with p2==0
475lstart = datetime_miss
476IF (PRESENT(start)) lstart = start
477lend = itime(SIZE(itime))
478! compute some quantities
479maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
480maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
481minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
482IF (time_definition == 0) THEN ! reference time
483 lend = lend + timedelta_new(sec=maxp1)
484ENDIF
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
488lend = lend + step
489IF (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
504ENDIF
505
506#ifdef DEBUG
507CALL 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
513IF (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
565ELSE ! 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
573ENDIF
574
575CALL packarray(a_otime)
576CALL packarray(a_otimerange)
577otime => a_otime%array
578otimerange => a_otimerange%array
579CALL sort(otime)
580CALL sort(otimerange)
581! delete local objects keeping the contents
582CALL delete(a_otime, nodealloc=.true.)
583CALL delete(a_otimerange, nodealloc=.true.)
584
585#ifdef DEBUG
586CALL l4f_log(l4f_debug, &
587 'recompute_stat_proc_agg, output time and timerange: '//&
588 t2c(SIZE(otime))//', '//t2c(size(otimerange)))
589#endif
590
591IF (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.)
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
637ELSE
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
673ENDIF
674
675END SUBROUTINE recompute_stat_proc_agg_common
676
677
678SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
679 max_step, weights)
680TYPE(datetime),INTENT(in) :: vertime(:)
681TYPE(datetime),INTENT(in) :: pstart
682TYPE(datetime),INTENT(in) :: pend
683LOGICAL,INTENT(in) :: time_mask(:)
684TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
685DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
686
687INTEGER :: i, nt
688TYPE(datetime),ALLOCATABLE :: lvertime(:)
689TYPE(datetime) :: half, nexthalf
690INTEGER(kind=int_ll) :: dt, tdt
691
692nt = count(time_mask)
693ALLOCATE(lvertime(nt))
694lvertime = pack(vertime, mask=time_mask)
695
696IF (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
719ENDIF
720
721IF (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
738ENDIF
739
740END 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.
744SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
745 pstart, pend, reftime)
746TYPE(datetime),INTENT(in) :: time
747TYPE(vol7d_timerange),INTENT(in) :: timerange
748INTEGER,INTENT(in) :: time_definition
749TYPE(datetime),INTENT(out) :: reftime
750TYPE(datetime),INTENT(out) :: pstart
751TYPE(datetime),INTENT(out) :: pend
752
753TYPE(timedelta) :: p1, p2
754
755
756p1 = timedelta_new(sec=timerange%p1) ! end of period
757p2 = timedelta_new(sec=timerange%p2) ! length of period
758
759IF (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
766ENDIF
767
768IF (time_definition == 0) THEN ! time == reference time
769 reftime = time
770 pend = time + p1
771 pstart = pend - p2
772ELSE IF (time_definition == 1) THEN ! time == verification time
773 pend = time
774 pstart = time - p2
775 reftime = time - p1
776ELSE
777 pstart = datetime_miss
778 pend = datetime_miss
779 reftime = datetime_miss
780ENDIF
781
782END 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"
788SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
789 pstart, pend, reftime)
790TYPE(datetime),INTENT(in) :: time
791TYPE(vol7d_timerange),INTENT(in) :: timerange
792TYPE(timedelta),INTENT(in) :: step
793INTEGER,INTENT(in) :: time_definition
794TYPE(datetime),INTENT(out) :: reftime
795TYPE(datetime),INTENT(out) :: pstart
796TYPE(datetime),INTENT(out) :: pend
797
798TYPE(timedelta) :: p1
799
800
801p1 = timedelta_new(sec=timerange%p1) ! end of period
802
803IF (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
810ENDIF
811
812IF (time_definition == 0) THEN ! time == reference time
813 reftime = time
814 pend = time + p1
815 pstart = pend - step
816ELSE IF (time_definition == 1) THEN ! time == verification time
817 pend = time
818 pstart = time - step
819 reftime = time - p1
820ELSE
821 pstart = datetime_miss
822 pend = datetime_miss
823 reftime = datetime_miss
824ENDIF
825
826END SUBROUTINE time_timerange_get_period_pop
827
828
829! set time, timerange%p1, timerange%p2 according to pstart, pend,
830! reftime and time_definition.
831SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
832 pstart, pend, reftime)
833TYPE(datetime),INTENT(out) :: time
834TYPE(vol7d_timerange),INTENT(inout) :: timerange
835INTEGER,INTENT(in) :: time_definition
836TYPE(datetime),INTENT(in) :: reftime
837TYPE(datetime),INTENT(in) :: pstart
838TYPE(datetime),INTENT(in) :: pend
839
840TYPE(timedelta) :: p1, p2
841INTEGER(kind=int_ll) :: dmsec
842
843
844IF (time_definition == 0) THEN ! time == reference time
845 time = reftime
846 p1 = pend - reftime
847 p2 = pend - pstart
848ELSE IF (time_definition == 1) THEN ! time == verification time
849 time = pend
850 p1 = pend - reftime
851 p2 = pend - pstart
852ELSE
853 time = datetime_miss
854ENDIF
856IF (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)
861ELSE
862 timerange%p1 = imiss
863 timerange%p2 = imiss
864ENDIF
865
866END SUBROUTINE time_timerange_set_period
867
868
869END 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.