libsim  Versione6.3.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 
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  nitr, otime, otimerange, map_tr, f, mask_timerange, 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 INTEGER,INTENT(out) :: nitr
141 TYPE(datetime),POINTER :: otime(:)
142 TYPE(vol7d_timerange),POINTER :: otimerange(:)
143 INTEGER,POINTER :: map_tr(:,:,:,:,:), f(:)
144 LOGICAL,POINTER :: 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 :: 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) > 1 .AND. c_e(start)) THEN ! security check
163  CALL getval(start-itime(1), asec=deltas)
164  ENDIF
165 ENDIF
166 
167 ! create a mask of suitable timeranges
168 ALLOCATE(mask_timerange(SIZE(itimerange)))
169 mask_timerange(:) = itimerange(:)%timerange == stat_proc &
170  .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
171  .AND. itimerange(:)%p1 >= 0 &
172  .AND. itimerange(:)%p2 > 0
173 
174 IF (optio_log(full_steps) .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer steps, check better steps /= 0
175  mask_timerange(:) = mask_timerange(:) .AND. (mod(itimerange(:)%p2-deltas, steps) == 0)
176 ENDIF
177 nitr = count(mask_timerange)
178 ALLOCATE(f(nitr))
179 j = 1
180 DO i = 1, nitr
181  DO WHILE(.NOT.mask_timerange(j))
182  j = j + 1
183  ENDDO
184  f(i) = j
185  j = j + 1
186 ENDDO
187 
188 ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
189 map_tr(:,:,:,:,:) = imiss
190 
191 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == steps
192 DO i = 1, SIZE(mask_timerange)
193  IF (mask_timerange(i)) THEN
194  j = append_unique(a_otimerange, itimerange(i))
195  ENDIF
196 ENDDO
197 
198 ! scan through all possible combinations of time and timerange
199 DO dirtyrep = 1, 2
200  IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
201  CALL packarray(a_otime)
202  CALL packarray(a_otimerange)
203  CALL sort(a_otime%array)
204  CALL sort(a_otimerange%array)
205  ENDIF
206  DO l = 1, SIZE(itime)
207  DO k = 1, nitr
208  CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
209  time_definition, pstart2, pend2, reftime2)
210 
211  DO j = 1, SIZE(itime)
212  DO i = 1, nitr
213  useful = .false.
214  CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
215  time_definition, pstart1, pend1, reftime1)
216  tmptimerange = vol7d_timerange_new(timerange=stat_proc)
217 
218  IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
219  IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
220  CALL time_timerange_set_period(tmptime, tmptimerange, &
221  time_definition, pend1, pend2, reftime2)
222  useful = .true.
223 
224  ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
225  CALL time_timerange_set_period(tmptime, tmptimerange, &
226  time_definition, pstart2, pstart1, pstart1)
227  useful = .true.
228  ENDIF
229 
230  ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
231  IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
232  CALL time_timerange_set_period(tmptime, tmptimerange, &
233  time_definition, pend1, pend2, reftime2)
234  useful = .true.
235 
236  ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
237  CALL time_timerange_set_period(tmptime, tmptimerange, &
238  time_definition, pstart2, pstart1, reftime2)
239  useful = .true.
240  ENDIF
241 
242  ENDIF
243  useful = useful .AND. tmptime /= datetime_miss .AND. &
244  tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
245 
246  IF (useful) THEN ! add a_otime, a_otimerange
247  map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
248  map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
249  ENDIF
250  ENDDO
251  ENDDO
252  ENDDO
253  ENDDO
254 ENDDO
255 
256 otime => a_otime%array
257 otimerange => a_otimerange%array
258 ! delete local objects keeping the contents
259 CALL delete(a_otime, nodealloc=.true.)
260 CALL delete(a_otimerange, nodealloc=.true.)
261 
262 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == steps
263 
264 #ifdef DEBUG
265 CALL l4f_log(l4f_debug, &
266  'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr,1)))//', '// &
267  t2c((SIZE(map_tr,2)))//', '// &
268  t2c((SIZE(map_tr,3)))//', '// &
269  t2c((SIZE(map_tr,4))))
270 CALL l4f_log(l4f_debug, &
271  'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr)))//', '// &
272  t2c(count(c_e(map_tr))))
273 CALL l4f_log(l4f_debug, &
274  'recompute_stat_proc_diff, nitr: '//t2c(nitr))
275 CALL l4f_log(l4f_debug, &
276  'recompute_stat_proc_diff, good timeranges: '//t2c(count(mask_timerange)))
277 CALL l4f_log(l4f_debug, &
278  'recompute_stat_proc_diff, output times: '//t2c(SIZE(otime)))
279 CALL l4f_log(l4f_debug, &
280  'recompute_stat_proc_diff, output timeranges: '//t2c(SIZE(otimerange)))
281 #endif
282 
283 END SUBROUTINE recompute_stat_proc_diff_common
284 
285 
286 ! common operations for statistical processing by metamorphosis
287 SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
288  otimerange, map_tr)
289 INTEGER,INTENT(in) :: istat_proc
290 TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
291 INTEGER,INTENT(in) :: ostat_proc
292 !TYPE(timedelta),INTENT(in) :: step
293 TYPE(vol7d_timerange),POINTER :: otimerange(:)
294 INTEGER,POINTER :: map_tr(:)
295 
296 INTEGER :: i
297 LOGICAL :: tr_mask(SIZE(itimerange))
298 
299 IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
300  ALLOCATE(otimerange(0), map_tr(0))
301  RETURN
302 ENDIF
303 
304 ! useful timeranges
305 tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
306  .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
307 ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
308 
309 otimerange = pack(itimerange, mask=tr_mask)
310 otimerange(:)%timerange = ostat_proc
311 map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
312 
313 END SUBROUTINE compute_stat_proc_metamorph_common
314 
315 
316 ! common operations for statistical processing by aggregation
317 SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
318  step, time_definition, otime, otimerange, map_ttr, dtratio, start)
319 TYPE(datetime),INTENT(in) :: itime(:)
320 TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
321 INTEGER,INTENT(in) :: stat_proc
322 INTEGER,INTENT(in) :: tri
323 TYPE(timedelta),INTENT(in) :: step
324 INTEGER,INTENT(in) :: time_definition
325 TYPE(datetime),POINTER :: otime(:)
326 TYPE(vol7d_timerange),POINTER :: otimerange(:)
327 TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
328 INTEGER,POINTER,OPTIONAL :: dtratio(:)
329 TYPE(datetime),INTENT(in),OPTIONAL :: start
331 INTEGER :: i, j, k, l, na, nf, n
332 INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart
333 LOGICAL :: lforecast
334 TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
335 TYPE(arrayof_datetime) :: a_otime
336 TYPE(arrayof_vol7d_timerange) :: a_otimerange
337 TYPE(arrayof_integer) :: a_dtratio
338 LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
339 TYPE(ttr_mapper) :: lmapper
340 CHARACTER(len=8) :: env_var
341 LOGICAL :: climat_behavior
342 
343 
344 ! enable bad behavior for climat database (only for instantaneous data)
345 env_var = ''
346 CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
347 climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
348 
349 ! compute length of cumulation step in seconds
350 CALL getval(timedelta_depop(step), asec=steps)
351 
352 ! create a mask of suitable timeranges
353 ALLOCATE(mask_timerange(SIZE(itimerange)))
354 mask_timerange(:) = itimerange(:)%timerange == tri &
355  .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
356 
357 IF (PRESENT(dtratio)) THEN
358  mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 /= imiss &
359  .AND. itimerange(:)%p2 > 0 .AND. mod(steps, itimerange(:)%p2) == 0
360 ELSE ! instantaneous
361  mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
362 ENDIF
363 
364 #ifdef DEBUG
365 CALL l4f_log(l4f_debug, &
366  '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
367  t2c(count(mask_timerange)))
368 #endif
369 
370 ! euristically determine whether we are dealing with an
371 ! analysis/observation or a forecast dataset
372 na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
373 nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
374 lforecast = nf >= na
375 #ifdef DEBUG
376 CALL l4f_log(l4f_debug, &
377  'recompute_stat_proc_agg, na: '//t2c(na)//', nf: '//t2c(nf))
378 #endif
379 
380 ! keep only timeranges of one type (really necessary?)
381 IF (lforecast) THEN
382  CALL l4f_log(l4f_info, &
383  'recompute_stat_proc_agg, processing in forecast mode')
384 ELSE
385  mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
386  CALL l4f_log(l4f_info, &
387  'recompute_stat_proc_agg, processing in analysis mode')
388 ENDIF
389 
390 #ifdef DEBUG
391 CALL l4f_log(l4f_debug, &
392  '(re)compute_stat_proc_agg, number of useful timeranges: '// &
393  t2c(count(mask_timerange)))
394 #endif
395 
396 IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
397  ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
398  IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
399  RETURN
400 ENDIF
401 
402 ! determine start and end of processing period, should work with p2==0
403 lstart = datetime_miss
404 IF (PRESENT(start)) lstart = start
405 lend = itime(SIZE(itime))
406 ! compute some quantities
407 maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
408 maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
409 minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
410 IF (time_definition == 0) THEN ! reference time
411  lend = lend + timedelta_new(sec=maxp1)
412 ENDIF
413 ! extend interval at the end in order to include al the data, must use
414 ! < and not <= in the DO WHILE loops some lines below
415 lend = lend + step
416 IF (lstart == datetime_miss) THEN ! autodetect
417  lstart = itime(1)
418 ! if autodetected, adjust to obtain real absolute start of data
419  IF (time_definition == 0) THEN ! reference time
420  lstart = lstart + timedelta_new(sec=minp1mp2)
421  ELSE ! verification time
422 ! go back to start of longest processing interval
423  lstart = lstart - timedelta_new(sec=maxp2)
424 ! lstart = lstart - (MOD(lstart, step)) ! round to step, + or - ?
425  ENDIF
426 ENDIF
427 
428 #ifdef DEBUG
429 CALL l4f_log(l4f_debug, &
430  'recompute_stat_proc_agg, processing period: '//t2c(lstart)//' - '//t2c(lend))
431 #endif
432 
433 ! create output time and timerange lists
434 
435 IF (lforecast) THEN ! forecast mode
436  IF (time_definition == 0) THEN ! reference time
437  CALL insert(a_otime, itime) ! should I limit to elements itime >= lstart?
438 
439 ! apply start shift to timerange, not to reftime
440 ! why did we use itime(SIZE(itime)) (last ref time)?
441 ! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
442  CALL getval(lstart-itime(1), asec=dstart)
443 ! allow starting before first reftime but restrict dtstart to -steps
444 ! dstart = MAX(0, dstart)
445  IF (dstart < 0) dstart = mod(dstart, steps)
446  DO p1 = steps + dstart, maxp1, steps
447  CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
448  ENDDO
449 
450  ELSE ! verification time
451  tmptime = lstart + step
452  DO WHILE(tmptime < lend) ! < since lend has been extended
453  CALL insert_unique(a_otime, tmptime)
454  tmptime = tmptime + step
455  ENDDO
456  DO p1 = steps, maxp1, steps
457  CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
458  ENDDO
459 
460  ENDIF
461 
462 ELSE ! analysis mode
463  tmptime = lstart + step
464  DO WHILE(tmptime < lend) ! < since lend has been extended
465  CALL insert_unique(a_otime, tmptime)
466  tmptime = tmptime + step
467  ENDDO
468  CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
469 
470 ENDIF
471 
472 CALL packarray(a_otime)
473 CALL packarray(a_otimerange)
474 otime => a_otime%array
475 otimerange => a_otimerange%array
476 CALL sort(otime)
477 CALL sort(otimerange)
478 ! delete local objects keeping the contents
479 CALL delete(a_otime, nodealloc=.true.)
480 CALL delete(a_otimerange, nodealloc=.true.)
481 
482 #ifdef DEBUG
483 CALL l4f_log(l4f_debug, &
484  'recompute_stat_proc_agg, output time and timerange: '//&
485  t2c(SIZE(otime))//', '//t2c(size(otimerange)))
486 #endif
487 
488 IF (PRESENT(dtratio)) THEN
489 ! count the possible i/o interval ratios
490  DO k = 1, SIZE(itimerange)
491  IF (itimerange(k)%p2 /= 0) &
492  CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
493  ENDDO
494  CALL packarray(a_dtratio)
495  dtratio => a_dtratio%array
496  CALL sort(dtratio)
497 ! delete local object keeping the contents
498  CALL delete(a_dtratio, nodealloc=.true.)
499 
500 #ifdef DEBUG
501  CALL l4f_log(l4f_debug, &
502  'recompute_stat_proc_agg, found '//t2c(size(dtratio))// &
503  ' possible aggregation ratios, from '// &
504  t2c(dtratio(1))//' to '//t2c(dtratio(SIZE(dtratio))))
505 #endif
506 
507  ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
508  do_itimerange1: DO l = 1, SIZE(itimerange)
509  IF (.NOT.mask_timerange(l)) cycle do_itimerange1
510  do_itime1: DO k = 1, SIZE(itime)
511  CALL time_timerange_get_period(itime(k), itimerange(l), &
512  time_definition, pstart1, pend1, reftime1)
513  do_otimerange1: DO j = 1, SIZE(otimerange)
514  do_otime1: DO i = 1, SIZE(otime)
515  CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
516  time_definition, pstart2, pend2, reftime2)
517  IF (lforecast) THEN
518  IF (reftime1 /= reftime2) cycle do_otime1
519  ENDIF
520 
521  IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
522  mod(pstart1-pstart2, pend1-pstart1) == timedelta_0) THEN ! useful
523  lmapper%it = k
524  lmapper%itr = l
525  lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
526  n = append(map_ttr(i,j), lmapper)
527  cycle do_itime1 ! can contribute only to a single interval
528  ENDIF
529  ENDDO do_otime1
530  ENDDO do_otimerange1
531  ENDDO do_itime1
532  ENDDO do_itimerange1
533 
534 ELSE
535 
536  ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
537  do_itimerange2: DO l = 1, SIZE(itimerange)
538  IF (.NOT.mask_timerange(l)) cycle do_itimerange2
539  do_itime2: DO k = 1, SIZE(itime)
540  CALL time_timerange_get_period(itime(k), itimerange(l), &
541  time_definition, pstart1, pend1, reftime1)
542  do_otimerange2: DO j = 1, SIZE(otimerange)
543  do_otime2: DO i = 1, SIZE(otime)
544  CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
545  time_definition, pstart2, pend2, reftime2)
546  IF (lforecast) THEN
547  IF (reftime1 /= reftime2) cycle do_otime2
548  ENDIF
549 
550  IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
551  IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
552  lmapper%it = k
553  lmapper%itr = l
554  IF (pstart1 == pstart2) THEN
555  lmapper%extra_info = 1 ! start of interval
556  ELSE IF (pend1 == pend2) THEN
557  lmapper%extra_info = 2 ! end of interval
558  ELSE
559  lmapper%extra_info = imiss
560  ENDIF
561  lmapper%time = pstart1 ! = pend1, order by time?
562  n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
563 ! no CYCLE, a single input can contribute to multiple output intervals
564  ENDIF
565  ENDDO do_otime2
566  ENDDO do_otimerange2
567  ENDDO do_itime2
568  ENDDO do_itimerange2
570 ENDIF
571 
572 END SUBROUTINE recompute_stat_proc_agg_common
573 
574 
575 SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
576  max_step, weights)
577 TYPE(datetime),INTENT(in) :: vertime(:)
578 TYPE(datetime),INTENT(in) :: pstart
579 TYPE(datetime),INTENT(in) :: pend
580 LOGICAL,INTENT(in) :: time_mask(:)
581 TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
582 DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
583 
584 INTEGER :: i, nt
585 TYPE(datetime),ALLOCATABLE :: lvertime(:)
586 TYPE(datetime) :: half, nexthalf
587 INTEGER(kind=int_ll) :: dt, tdt
588 
589 nt = count(time_mask)
590 ALLOCATE(lvertime(nt))
591 lvertime = pack(vertime, mask=time_mask)
592 
593 IF (PRESENT(max_step)) THEN
594 ! new way
595 ! max_step = timedelta_0
596 ! DO i = 1, nt - 1
597 ! IF (lvertime(i+1) - lvertime(i) > max_step) &
598 ! max_step = lvertime(i+1) - lvertime(i)
599 ! ENDDO
600 ! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
601 ! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
602 ! old way
603  IF (nt == 1) THEN
604  max_step = pend - pstart
605  ELSE
606  half = lvertime(1) + (lvertime(2) - lvertime(1))/2
607  max_step = half - pstart
608  DO i = 2, nt - 1
609  nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
610  IF (nexthalf - half > max_step) max_step = nexthalf - half
611  half = nexthalf
612  ENDDO
613  IF (pend - half > max_step) max_step = pend - half
614  ENDIF
615 
616 ENDIF
617 
618 IF (PRESENT(weights)) THEN
619  IF (nt == 1) THEN
620  weights(1) = 1.0
621  ELSE
622  CALL getval(pend - pstart, amsec=tdt)
623  half = lvertime(1) + (lvertime(2) - lvertime(1))/2
624  CALL getval(half - pstart, amsec=dt)
625  weights(1) = dble(dt)/dble(tdt)
626  DO i = 2, nt - 1
627  nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
628  CALL getval(nexthalf - half, amsec=dt)
629  weights(i) = dble(dt)/dble(tdt)
630  half = nexthalf
631  ENDDO
632  CALL getval(pend - half, amsec=dt)
633  weights(nt) = dble(dt)/dble(tdt)
634  ENDIF
635 ENDIF
636 
637 END SUBROUTINE compute_stat_proc_agg_sw
638 
639 ! get start of period, end of period and reference time from time,
640 ! timerange, according to time_definition.
641 SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
642  pstart, pend, reftime)
643 TYPE(datetime),INTENT(in) :: time
644 TYPE(vol7d_timerange),INTENT(in) :: timerange
645 INTEGER,INTENT(in) :: time_definition
646 TYPE(datetime),INTENT(out) :: reftime
647 TYPE(datetime),INTENT(out) :: pstart
648 TYPE(datetime),INTENT(out) :: pend
649 
650 TYPE(timedelta) :: p1, p2
651 
652 
653 p1 = timedelta_new(sec=timerange%p1) ! end of period
654 p2 = timedelta_new(sec=timerange%p2) ! length of period
655 
656 IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
657 ! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
658  timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
659  pstart = datetime_miss
660  pend = datetime_miss
661  reftime = datetime_miss
662  RETURN
663 ENDIF
664 
665 IF (time_definition == 0) THEN ! time == reference time
666  reftime = time
667  pend = time + p1
668  pstart = pend - p2
669 ELSE IF (time_definition == 1) THEN ! time == verification time
670  pend = time
671  pstart = time - p2
672  reftime = time - p1
673 ELSE
674  pstart = datetime_miss
675  pend = datetime_miss
676  reftime = datetime_miss
677 ENDIF
678 
679 END SUBROUTINE time_timerange_get_period
680 
681 
682 ! get start of period, end of period and reference time from time,
683 ! timerange, according to time_definition. step is taken separately
684 ! from timerange and can be "popular"
685 SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
686  pstart, pend, reftime)
687 TYPE(datetime),INTENT(in) :: time
688 TYPE(vol7d_timerange),INTENT(in) :: timerange
689 TYPE(timedelta),INTENT(in) :: step
690 INTEGER,INTENT(in) :: time_definition
691 TYPE(datetime),INTENT(out) :: reftime
692 TYPE(datetime),INTENT(out) :: pstart
693 TYPE(datetime),INTENT(out) :: pend
694 
695 TYPE(timedelta) :: p1
696 
697 
698 p1 = timedelta_new(sec=timerange%p1) ! end of period
699 
700 IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
701 ! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
702  timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
703  pstart = datetime_miss
704  pend = datetime_miss
705  reftime = datetime_miss
706  RETURN
707 ENDIF
708 
709 IF (time_definition == 0) THEN ! time == reference time
710  reftime = time
711  pend = time + p1
712  pstart = pend - step
713 ELSE IF (time_definition == 1) THEN ! time == verification time
714  pend = time
715  pstart = time - step
716  reftime = time - p1
717 ELSE
718  pstart = datetime_miss
719  pend = datetime_miss
720  reftime = datetime_miss
721 ENDIF
722 
723 END SUBROUTINE time_timerange_get_period_pop
725 
726 ! set time, timerange%p1, timerange%p2 according to pstart, pend,
727 ! reftime and time_definition.
728 SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
729  pstart, pend, reftime)
730 TYPE(datetime),INTENT(out) :: time
731 TYPE(vol7d_timerange),INTENT(inout) :: timerange
732 INTEGER,INTENT(in) :: time_definition
733 TYPE(datetime),INTENT(in) :: reftime
734 TYPE(datetime),INTENT(in) :: pstart
735 TYPE(datetime),INTENT(in) :: pend
736 
737 TYPE(timedelta) :: p1, p2
738 INTEGER(kind=int_ll) :: dmsec
739 
740 
741 IF (time_definition == 0) THEN ! time == reference time
742  time = reftime
743  p1 = pend - reftime
744  p2 = pend - pstart
745 ELSE IF (time_definition == 1) THEN ! time == verification time
746  time = pend
747  p1 = pend - reftime
748  p2 = pend - pstart
749 ELSE
750  time = datetime_miss
751 ENDIF
752 
753 IF (time /= datetime_miss) THEN
754  CALL getval(p1, amsec=dmsec) ! end of period
755  timerange%p1 = int(dmsec/1000_int_ll)
756  CALL getval(p2, amsec=dmsec) ! length of period
757  timerange%p2 = int(dmsec/1000_int_ll)
758 ELSE
759  timerange%p1 = imiss
760  timerange%p2 = imiss
761 ENDIF
762 
763 END SUBROUTINE time_timerange_set_period
764 
765 
766 END MODULE stat_proc_engine
Functions that return a trimmed CHARACTER representation of the input variable.
Class for expressing an absolute time value.
This module contains functions that are only for internal use of the library.
Classi per la gestione delle coordinate temporali.
Quick method to append an element to the array.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Restituiscono il valore dell&#39;oggetto nella forma desiderata.
Operatore di resto della divisione.
Classe per la gestione di un volume completo di dati osservati.
Class for expressing a relative time interval.
Destructor for finalizing an array object.
Method for inserting elements of the array at a desired position.

Generated with Doxygen.