libsim  Versione 7.2.4
vol7d_class_compute.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 
28 USE vol7d_class
30 USE simple_stat
31 IMPLICIT NONE
32 
33 CONTAINS
34 
35 
106 SUBROUTINE vol7d_compute_stat_proc(this, that, stat_proc_input, stat_proc, &
107  step, start, full_steps, frac_valid, max_step, weighted, other)
108 TYPE(vol7d),INTENT(inout) :: this
109 TYPE(vol7d),INTENT(out) :: that
110 INTEGER,INTENT(in) :: stat_proc_input
111 INTEGER,INTENT(in) :: stat_proc
112 TYPE(timedelta),INTENT(in) :: step
113 TYPE(datetime),INTENT(in),OPTIONAL :: start
114 LOGICAL,INTENT(in),OPTIONAL :: full_steps
115 REAL,INTENT(in),OPTIONAL :: frac_valid
116 TYPE(timedelta),INTENT(in),OPTIONAL :: max_step
117 LOGICAL,INTENT(in),OPTIONAL :: weighted
118 TYPE(vol7d),INTENT(inout),OPTIONAL :: other
119 
120 TYPE(vol7d) :: that1, that2, other1
121 INTEGER :: steps
122 
123 IF (stat_proc_input == 254) THEN
124  CALL l4f_log(l4f_info, 'computing statistical processing by aggregation '//&
125  trim(to_char(stat_proc_input))//':'//trim(to_char(stat_proc)))
126 
127  CALL vol7d_compute_stat_proc_agg(this, that, stat_proc, &
128  step, start, full_steps, max_step, weighted, other)
129 
130 ELSE IF (stat_proc == 254) THEN
131  CALL l4f_log(l4f_info, &
132  'computing instantaneous data from statistically processed '//&
133  trim(to_char(stat_proc_input))//':'//trim(to_char(stat_proc)))
134 
135 ! compute length of cumulation step in seconds
136  CALL getval(step, asec=steps)
137 
138  IF (any(this%timerange(:)%p2 == steps)) THEN ! data is ready
139  CALL vol7d_decompute_stat_proc(this, that, step, other, stat_proc_input)
140  ELSE
141  IF (any(this%timerange(:)%p2 == steps/2)) THEN ! need to average
142 ! average twice on step interval, with a shift of step/2, check full_steps
143  CALL vol7d_recompute_stat_proc_agg(this, that1, stat_proc_input, &
144  step, full_steps=.false., frac_valid=1.0)
145  CALL vol7d_recompute_stat_proc_agg(this, that2, stat_proc_input, &
146  step, start=that1%time(1)+step/2, full_steps=.false., frac_valid=1.0)
147 ! merge the result
148  CALL vol7d_append(that1, that2, sort=.true., lanasimple=.true.)
149 ! and process it
150  CALL vol7d_decompute_stat_proc(that1, that, step, other, stat_proc_input)
151  CALL delete(that1)
152  CALL delete(that2)
153  ELSE
154 ! do nothing
155  ENDIF
156  ENDIF
157 
158 ELSE IF (stat_proc_input == stat_proc .OR. &
159  (stat_proc == 0 .OR. stat_proc == 2 .OR. stat_proc == 3)) THEN
160 ! avg, min and max can be computed from any input, with care
161  CALL l4f_log(l4f_info, &
162  'recomputing statistically processed data by aggregation and difference '//&
163  trim(to_char(stat_proc_input))//':'//trim(to_char(stat_proc)))
164 
165  IF (PRESENT(other)) THEN
166  CALL vol7d_recompute_stat_proc_agg(this, that1, stat_proc, &
167  step, start, full_steps, frac_valid, &
168  other=other, stat_proc_input=stat_proc_input)
169  CALL vol7d_recompute_stat_proc_diff(this, that2, stat_proc, &
170  step, full_steps, start, other=other1)
171  CALL vol7d_merge(other, other1, sort=.true.)
172  ELSE
173  CALL vol7d_recompute_stat_proc_agg(this, that1, stat_proc, &
174  step, start, full_steps, frac_valid, stat_proc_input=stat_proc_input)
175  CALL vol7d_recompute_stat_proc_diff(this, that2, stat_proc, step, full_steps, &
176  start)
177  ENDIF
178 
179  CALL vol7d_merge(that1, that2, sort=.true., bestdata=.true.)
180  CALL delete(that2)
181  that = that1
182 ELSE ! IF (stat_proc_input /= stat_proc) THEN
183  IF ((stat_proc_input == 0 .AND. stat_proc == 1) .OR. &
184  (stat_proc_input == 1 .AND. stat_proc == 0)) THEN
185  CALL l4f_log(l4f_info, &
186  'computing statistically processed data by integration/differentiation '// &
187  t2c(stat_proc_input)//':'//t2c(stat_proc))
188  CALL vol7d_compute_stat_proc_metamorph(this, that, stat_proc_input, &
189  stat_proc)
190  ELSE
191  CALL l4f_log(l4f_error, &
192  'statistical processing '//t2c(stat_proc_input)//':'//t2c(stat_proc)// &
193  ' not implemented or does not make sense')
194  CALL raise_error()
195  ENDIF
196 ENDIF
197 
198 END SUBROUTINE vol7d_compute_stat_proc
199 
200 
246 SUBROUTINE vol7d_recompute_stat_proc_agg(this, that, stat_proc, &
247  step, start, full_steps, frac_valid, other, stat_proc_input)
248 TYPE(vol7d),INTENT(inout) :: this
249 TYPE(vol7d),INTENT(out) :: that
250 INTEGER,INTENT(in) :: stat_proc
251 TYPE(timedelta),INTENT(in) :: step
252 TYPE(datetime),INTENT(in),OPTIONAL :: start
253 LOGICAL,INTENT(in),OPTIONAL :: full_steps
254 REAL,INTENT(in),OPTIONAL :: frac_valid
255 TYPE(vol7d),INTENT(inout),OPTIONAL :: other
256 INTEGER,INTENT(in),OPTIONAL :: stat_proc_input
257 
258 INTEGER :: tri
259 INTEGER :: i, j, n, n1, ndtr, i1, i3, i5, i6
260 INTEGER :: linshape(1)
261 REAL :: lfrac_valid, frac_c, frac_m
262 LOGICAL,ALLOCATABLE :: ttr_mask(:,:)
263 TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
264 INTEGER,POINTER :: dtratio(:)
265 
266 
267 IF (PRESENT(stat_proc_input)) THEN
268  tri = stat_proc_input
269 ELSE
270  tri = stat_proc
271 ENDIF
272 IF (PRESENT(frac_valid)) THEN
273  lfrac_valid = frac_valid
274 ELSE
275  lfrac_valid = 1.0
276 ENDIF
277 
278 ! be safe
279 CALL vol7d_alloc_vol(this)
280 ! initial check
281 
282 ! cleanup the original volume
283 CALL vol7d_smart_sort(this, lsort_time=.true.) ! time-ordered volume needed
284 CALL vol7d_reform(this, miss=.false., sort=.false., unique=.true.)
285 
286 CALL init(that, time_definition=this%time_definition)
287 CALL vol7d_alloc(that, nana=SIZE(this%ana), nlevel=SIZE(this%level), &
288  nnetwork=SIZE(this%network))
289 IF (ASSOCIATED(this%dativar%r)) THEN
290  CALL vol7d_alloc(that, ndativarr=SIZE(this%dativar%r))
291  that%dativar%r = this%dativar%r
292 ENDIF
293 IF (ASSOCIATED(this%dativar%d)) THEN
294  CALL vol7d_alloc(that, ndativard=SIZE(this%dativar%d))
295  that%dativar%d = this%dativar%d
296 ENDIF
297 that%ana = this%ana
298 that%level = this%level
299 that%network = this%network
300 
301 ! compute the output time and timerange and all the required mappings
302 CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
303  step, this%time_definition, that%time, that%timerange, map_ttr, dtratio, &
304  start, full_steps)
305 CALL vol7d_alloc_vol(that)
306 
307 ALLOCATE(ttr_mask(SIZE(this%time), SIZE(this%timerange)))
308 linshape = (/SIZE(ttr_mask)/)
309 ! finally perform computations
310 IF (ASSOCIATED(this%voldatir)) THEN
311  DO j = 1, SIZE(that%timerange)
312  DO i = 1, SIZE(that%time)
313 
314  DO i1 = 1, SIZE(this%ana)
315  DO i3 = 1, SIZE(this%level)
316  DO i6 = 1, SIZE(this%network)
317  DO i5 = 1, SIZE(this%dativar%r)
318 
319  frac_m = 0.
320  DO n1 = SIZE(dtratio), 1, -1 ! precedence to longer periods
321  IF (dtratio(n1) <= 0) cycle ! safety check
322  ttr_mask = .false.
323  DO n = 1, map_ttr(i,j)%arraysize
324  IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
325  IF (c_e(this%voldatir(i1,map_ttr(i,j)%array(n)%it,i3, &
326  map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
327  ttr_mask(map_ttr(i,j)%array(n)%it, &
328  map_ttr(i,j)%array(n)%itr) = .true.
329  ENDIF
330  ENDIF
331  ENDDO
332 
333  ndtr = count(ttr_mask)
334  frac_c = real(ndtr)/real(dtratio(n1))
335 
336  IF (ndtr > 0 .AND. frac_c >= max(lfrac_valid, frac_m)) THEN
337  frac_m = frac_c
338  SELECT CASE(stat_proc)
339  CASE (0, 200) ! average, vectorial mean
340  that%voldatir(i1,i,i3,j,i5,i6) = &
341  sum(this%voldatir(i1,:,i3,:,i5,i6), &
342  mask=ttr_mask)/ndtr
343  CASE (1, 4) ! accumulation, difference
344  that%voldatir(i1,i,i3,j,i5,i6) = &
345  sum(this%voldatir(i1,:,i3,:,i5,i6), &
346  mask=ttr_mask)
347  CASE (2) ! maximum
348  that%voldatir(i1,i,i3,j,i5,i6) = &
349  maxval(this%voldatir(i1,:,i3,:,i5,i6), &
350  mask=ttr_mask)
351  CASE (3) ! minimum
352  that%voldatir(i1,i,i3,j,i5,i6) = &
353  minval(this%voldatir(i1,:,i3,:,i5,i6), &
354  mask=ttr_mask)
355  CASE (6) ! stddev
356  that%voldatir(i1,i,i3,j,i5,i6) = &
357  stat_stddev( &
358  reshape(this%voldatir(i1,:,i3,:,i5,i6), shape=linshape), &
359  mask=reshape(ttr_mask, shape=linshape))
360  END SELECT
361  ENDIF
362 
363  ENDDO ! dtratio
364  ENDDO ! var
365  ENDDO ! network
366  ENDDO ! level
367  ENDDO ! ana
368  CALL delete(map_ttr(i,j))
369  ENDDO ! otime
370  ENDDO ! otimerange
371 ENDIF
372 
373 IF (ASSOCIATED(this%voldatid)) THEN
374  DO j = 1, SIZE(that%timerange)
375  DO i = 1, SIZE(that%time)
376 
377  DO i1 = 1, SIZE(this%ana)
378  DO i3 = 1, SIZE(this%level)
379  DO i6 = 1, SIZE(this%network)
380  DO i5 = 1, SIZE(this%dativar%d)
381 
382  frac_m = 0.
383  DO n1 = SIZE(dtratio), 1, -1 ! precedence to longer periods
384  IF (dtratio(n1) <= 0) cycle ! safety check
385  ttr_mask = .false.
386  DO n = 1, map_ttr(i,j)%arraysize
387  IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
388  IF (c_e(this%voldatid(i1,map_ttr(i,j)%array(n)%it,i3, &
389  map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
390  ttr_mask(map_ttr(i,j)%array(n)%it, &
391  map_ttr(i,j)%array(n)%itr) = .true.
392  ENDIF
393  ENDIF
394  ENDDO
395 
396  ndtr = count(ttr_mask)
397  frac_c = real(ndtr)/real(dtratio(n1))
398 
399  IF (ndtr > 0 .AND. frac_c >= max(lfrac_valid, frac_m)) THEN
400  frac_m = frac_c
401  SELECT CASE(stat_proc)
402  CASE (0) ! average
403  that%voldatid(i1,i,i3,j,i5,i6) = &
404  sum(this%voldatid(i1,:,i3,:,i5,i6), &
405  mask=ttr_mask)/ndtr
406  CASE (1, 4) ! accumulation, difference
407  that%voldatid(i1,i,i3,j,i5,i6) = &
408  sum(this%voldatid(i1,:,i3,:,i5,i6), &
409  mask=ttr_mask)
410  CASE (2) ! maximum
411  that%voldatid(i1,i,i3,j,i5,i6) = &
412  maxval(this%voldatid(i1,:,i3,:,i5,i6), &
413  mask=ttr_mask)
414  CASE (3) ! minimum
415  that%voldatid(i1,i,i3,j,i5,i6) = &
416  minval(this%voldatid(i1,:,i3,:,i5,i6), &
417  mask=ttr_mask)
418  CASE (6) ! stddev
419  that%voldatid(i1,i,i3,j,i5,i6) = &
420  stat_stddev( &
421  reshape(this%voldatid(i1,:,i3,:,i5,i6), shape=linshape), &
422  mask=reshape(ttr_mask, shape=linshape))
423  END SELECT
424  ENDIF
425 
426  ENDDO ! dtratio
427  ENDDO ! var
428  ENDDO ! network
429  ENDDO ! level
430  ENDDO ! ana
431  CALL delete(map_ttr(i,j))
432  ENDDO ! otime
433  ENDDO ! otimerange
434 ENDIF
435 
436 DEALLOCATE(ttr_mask)
437 
438 CALL makeother()
439 
440 CONTAINS
441 
442 SUBROUTINE makeother()
443 IF (PRESENT(other)) THEN ! create volume with the remaining data for further processing
444  CALL vol7d_copy(this, other, miss=.false., sort=.false., unique=.false., &
445  ltimerange=(this%timerange(:)%timerange /= tri .OR. this%timerange(:)%p2 == imiss &
446  .OR. this%timerange(:)%p2 == 0)) ! or MOD(steps, this%timerange(:)%p2) == 0
447 ENDIF
448 END SUBROUTINE makeother
449 
450 END SUBROUTINE vol7d_recompute_stat_proc_agg
451 
452 
484 SUBROUTINE vol7d_compute_stat_proc_agg(this, that, stat_proc, &
485  step, start, full_steps, max_step, weighted, other)
486 TYPE(vol7d),INTENT(inout) :: this
487 TYPE(vol7d),INTENT(out) :: that
488 INTEGER,INTENT(in) :: stat_proc
489 TYPE(timedelta),INTENT(in) :: step
490 TYPE(datetime),INTENT(in),OPTIONAL :: start
491 LOGICAL,INTENT(in),OPTIONAL :: full_steps
492 TYPE(timedelta),INTENT(in),OPTIONAL :: max_step
493 LOGICAL,INTENT(in),OPTIONAL :: weighted
494 TYPE(vol7d),INTENT(inout),OPTIONAL :: other
495 !INTEGER,INTENT(in),OPTIONAL :: stat_proc_input !< to be used with care, type of statistical processing of data that has to be processed (from grib2 table), only data having timerange of this type will be recomputed, the actual statistical processing performed and which will appear in the output volume, is however determined by \a stat_proc argument
496 
497 TYPE(vol7d) :: v7dtmp
498 INTEGER :: tri
499 INTEGER :: i, j, n, ninp, ndtr, i1, i3, i5, i6, vartype, maxsize
500 TYPE(timedelta) :: lmax_step, act_max_step
501 TYPE(datetime) :: pstart, pend, reftime
502 TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
503 REAL,ALLOCATABLE :: tmpvolr(:)
504 DOUBLE PRECISION,ALLOCATABLE :: tmpvold(:), weights(:)
505 LOGICAL,ALLOCATABLE :: lin_mask(:)
506 LOGICAL :: lweighted
507 CHARACTER(len=8) :: env_var
508 
509 IF (PRESENT(max_step)) THEN
510  lmax_step = max_step
511 ELSE
512  lmax_step = timedelta_max
513 ENDIF
514 lweighted = optio_log(weighted)
515 tri = 254
516 ! enable bad behavior for climat database
517 env_var = ''
518 CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
519 lweighted = lweighted .AND. len_trim(env_var) == 0
520 ! only average can be weighted
521 lweighted = lweighted .AND. stat_proc == 0
522 
523 ! be safe
524 CALL vol7d_alloc_vol(this)
525 ! initial check
526 
527 ! cleanup the original volume
528 CALL vol7d_smart_sort(this, lsort_time=.true.) ! time-ordered volume needed
529 CALL vol7d_reform(this, miss=.false., sort=.false., unique=.true.)
530 ! copy everything except time and timerange
531 CALL vol7d_copy(this, v7dtmp, ltime=(/.false./), ltimerange=(/.false./))
532 
533 ! create new volume
534 CALL init(that, time_definition=this%time_definition)
535 ! compute the output time and timerange and all the required mappings
536 CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
537  step, this%time_definition, that%time, that%timerange, map_ttr, start=start, &
538  full_steps=full_steps)
539 ! merge with information from original volume
540 CALL vol7d_merge(that, v7dtmp)
541 
542 maxsize = maxval(map_ttr(:,:)%arraysize)
543 ALLOCATE(tmpvolr(maxsize), tmpvold(maxsize), lin_mask(maxsize), weights(maxsize))
544 do_otimerange: DO j = 1, SIZE(that%timerange)
545  do_otime: DO i = 1, SIZE(that%time)
546  ninp = map_ttr(i,j)%arraysize
547  IF (ninp <= 0) cycle do_otime
548 ! required for some computations
549  CALL time_timerange_get_period(that%time(i), that%timerange(j), &
550  that%time_definition, pstart, pend, reftime)
551 
552  IF (ASSOCIATED(this%voldatir)) THEN
553  DO i1 = 1, SIZE(this%ana)
554  DO i3 = 1, SIZE(this%level)
555  DO i6 = 1, SIZE(this%network)
556  DO i5 = 1, SIZE(this%dativar%r)
557 ! stat_proc difference treated separately here
558  IF (stat_proc == 4) THEN
559  IF (ninp >= 2) THEN
560  IF (map_ttr(i,j)%array(1)%extra_info == 1 .AND. &
561  map_ttr(i,j)%array(ninp)%extra_info == 2) THEN
562  IF (c_e(this%voldatir(i1,map_ttr(i,j)%array(1)%it,i3, &
563  map_ttr(i,j)%array(1)%itr,i5,i6)) .AND. &
564  c_e(this%voldatir(i1,map_ttr(i,j)%array(ninp)%it,i3, &
565  map_ttr(i,j)%array(ninp)%itr,i5,i6))) THEN
566  that%voldatir(i1,i,i3,j,i5,i6) = &
567  this%voldatir(i1,map_ttr(i,j)%array(ninp)%it,i3, &
568  map_ttr(i,j)%array(ninp)%itr,i5,i6) - &
569  this%voldatir(i1,map_ttr(i,j)%array(1)%it,i3, &
570  map_ttr(i,j)%array(1)%itr,i5,i6)
571  ENDIF
572  ENDIF
573  ENDIF
574  cycle
575  ENDIF
576 ! other stat_proc
577  vartype = vol7d_vartype(this%dativar%r(i5))
578  lin_mask = .false.
579  ndtr = 0
580  DO n = 1, ninp
581  IF (c_e(this%voldatir(i1,map_ttr(i,j)%array(n)%it,i3, &
582  map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
583  ndtr = ndtr + 1
584  tmpvolr(ndtr) = this%voldatir(i1,map_ttr(i,j)%array(n)%it,i3, &
585  map_ttr(i,j)%array(n)%itr,i5,i6)
586  lin_mask(n) = .true.
587  ENDIF
588  ENDDO
589  IF (ndtr == 0) cycle
590  IF (lweighted) THEN
591  CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
592  pstart, pend, lin_mask(1:ninp), act_max_step, weights)
593  ELSE
594  CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
595  pstart, pend, lin_mask(1:ninp), act_max_step)
596  ENDIF
597  IF (act_max_step > lmax_step) cycle
598 
599  SELECT CASE(stat_proc)
600  CASE (0) ! average
601  IF (lweighted) THEN
602  that%voldatir(i1,i,i3,j,i5,i6) = &
603  sum(real(weights(1:ndtr))*tmpvolr(1:ndtr))
604  ELSE
605  that%voldatir(i1,i,i3,j,i5,i6) = &
606  sum(tmpvolr(1:ndtr))/ndtr
607  ENDIF
608  CASE (2) ! maximum
609  that%voldatir(i1,i,i3,j,i5,i6) = &
610  maxval(tmpvolr(1:ndtr))
611  CASE (3) ! minimum
612  that%voldatir(i1,i,i3,j,i5,i6) = &
613  minval(tmpvolr(1:ndtr))
614  CASE (6) ! stddev
615  that%voldatir(i1,i,i3,j,i5,i6) = &
616  stat_stddev(tmpvolr(1:ndtr))
617  CASE (201) ! mode
618 ! mode only for angles at the moment, with predefined histogram
619  IF (vartype == var_dir360) THEN
620 ! remove undefined wind direction (=0), improve check?
621 ! and reduce to interval [22.5,382.5[
622  WHERE (tmpvolr(1:ndtr) == 0.0)
623  tmpvolr(1:ndtr) = rmiss
624  ELSE WHERE (tmpvolr(1:ndtr) < 22.5 .AND. tmpvolr(1:ndtr) > 0.0)
625  tmpvolr(1:ndtr) = tmpvolr(1:ndtr) + 360.
626  END WHERE
627  that%voldatir(i1,i,i3,j,i5,i6) = &
628  stat_mode_histogram(tmpvolr(1:ndtr), &
629  8, 22.5, 382.5)
630  ENDIF
631  END SELECT
632  ENDDO
633  ENDDO
634  ENDDO
635  ENDDO
636  ENDIF
637 
638  IF (ASSOCIATED(this%voldatid)) THEN
639  DO i1 = 1, SIZE(this%ana)
640  DO i3 = 1, SIZE(this%level)
641  DO i6 = 1, SIZE(this%network)
642  DO i5 = 1, SIZE(this%dativar%d)
643 ! stat_proc difference treated separately here
644  IF (stat_proc == 4) THEN
645  IF (n >= 2) THEN
646  IF (map_ttr(i,j)%array(1)%extra_info == 1 .AND. &
647  map_ttr(i,j)%array(ninp)%extra_info == 2) THEN
648  IF (c_e(this%voldatid(i1,map_ttr(i,j)%array(1)%it,i3, &
649  map_ttr(i,j)%array(1)%itr,i5,i6)) .AND. &
650  c_e(this%voldatid(i1,map_ttr(i,j)%array(ninp)%it,i3, &
651  map_ttr(i,j)%array(ninp)%itr,i5,i6))) THEN
652  that%voldatid(i1,i,i3,j,i5,i6) = &
653  this%voldatid(i1,map_ttr(i,j)%array(ninp)%it,i3, &
654  map_ttr(i,j)%array(ninp)%itr,i5,i6) - &
655  this%voldatid(i1,map_ttr(i,j)%array(1)%it,i3, &
656  map_ttr(i,j)%array(1)%itr,i5,i6)
657  ENDIF
658  ENDIF
659  ENDIF
660  cycle
661  ENDIF
662 ! other stat_proc
663  vartype = vol7d_vartype(this%dativar%d(i5))
664  lin_mask = .false.
665  ndtr = 0
666  DO n = 1, ninp
667  IF (c_e(this%voldatid(i1,map_ttr(i,j)%array(n)%it,i3, &
668  map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
669  ndtr = ndtr + 1
670  tmpvold(ndtr) = this%voldatid(i1,map_ttr(i,j)%array(n)%it,i3, &
671  map_ttr(i,j)%array(n)%itr,i5,i6)
672  lin_mask(n) = .true.
673  ENDIF
674  ENDDO
675  IF (ndtr == 0) cycle
676  IF (lweighted) THEN
677  CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
678  pstart, pend, lin_mask(1:ninp), act_max_step, weights)
679  ELSE
680  CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
681  pstart, pend, lin_mask(1:ninp), act_max_step)
682  ENDIF
683  IF (act_max_step > lmax_step) cycle
684 
685  SELECT CASE(stat_proc)
686  CASE (0) ! average
687  IF (lweighted) THEN
688  that%voldatid(i1,i,i3,j,i5,i6) = &
689  sum(real(weights(1:ndtr))*tmpvold(1:ndtr))
690  ELSE
691  that%voldatid(i1,i,i3,j,i5,i6) = &
692  sum(tmpvold(1:ndtr))/ndtr
693  ENDIF
694  CASE (2) ! maximum
695  that%voldatid(i1,i,i3,j,i5,i6) = &
696  maxval(tmpvold(1:ndtr))
697  CASE (3) ! minimum
698  that%voldatid(i1,i,i3,j,i5,i6) = &
699  minval(tmpvold(1:ndtr))
700  CASE (6) ! stddev
701  that%voldatid(i1,i,i3,j,i5,i6) = &
702  stat_stddev(tmpvold(1:ndtr))
703  CASE (201) ! mode
704 ! mode only for angles at the moment, with predefined histogram
705  IF (vartype == var_dir360) THEN
706 ! remove undefined wind direction (=0), improve check?
707 ! and reduce to interval [22.5,382.5[
708  WHERE (tmpvold(1:ndtr) == 0.0d0)
709  tmpvold(1:ndtr) = dmiss
710  ELSE WHERE (tmpvold(1:ndtr) < 22.5d0 .AND. tmpvold(1:ndtr) > 0.0d0)
711  tmpvold(1:ndtr) = tmpvold(1:ndtr) + 360.0d0
712  END WHERE
713  that%voldatid(i1,i,i3,j,i5,i6) = &
714  stat_mode_histogram(tmpvold(1:ndtr), &
715  8, 22.5d0, 382.5d0)
716  ENDIF
717  END SELECT
718  ENDDO
719  ENDDO
720  ENDDO
721  ENDDO
722  ENDIF
723 
724  CALL delete(map_ttr(i,j))
725  ENDDO do_otime
726 ENDDO do_otimerange
727 
728 DEALLOCATE(map_ttr)
729 DEALLOCATE(tmpvolr, tmpvold, lin_mask, weights)
730 
731 IF (PRESENT(other)) THEN ! create volume with the remaining data for further processing
732  CALL vol7d_copy(this, other, miss=.false., sort=.false., unique=.false., &
733  ltimerange=(this%timerange(:)%timerange /= tri))
734 ENDIF
735 
736 END SUBROUTINE vol7d_compute_stat_proc_agg
737 
738 
754 SUBROUTINE vol7d_decompute_stat_proc(this, that, step, other, stat_proc_input)
755 TYPE(vol7d),INTENT(inout) :: this
756 TYPE(vol7d),INTENT(out) :: that
757 TYPE(timedelta),INTENT(in) :: step
758 TYPE(vol7d),INTENT(inout),OPTIONAL :: other
759 INTEGER,INTENT(in),OPTIONAL :: stat_proc_input
760 
761 INTEGER :: i, tri, steps
762 
763 
764 IF (PRESENT(stat_proc_input)) THEN
765  tri = stat_proc_input
766 ELSE
767  tri = 0
768 ENDIF
769 ! be safe
770 CALL vol7d_alloc_vol(this)
771 
772 ! compute length of cumulation step in seconds
773 CALL getval(step, asec=steps)
774 
775 ! filter requested data
776 CALL vol7d_copy(this, that, miss=.false., sort=.false., unique=.false., &
777  ltimerange=(this%timerange(:)%timerange == tri .AND. &
778  this%timerange(:)%p1 == 0 .AND. this%timerange(:)%p2 == steps))
779 
780 ! convert timerange to instantaneous and go back half step in time
781 that%timerange(:)%timerange = 254
782 that%timerange(:)%p2 = 0
783 DO i = 1, SIZE(that%time(:))
784  that%time(i) = that%time(i) - step/2
785 ENDDO
786 
787 IF (PRESENT(other)) THEN ! create volume with the remaining data for further processing
788  CALL vol7d_copy(this, other, miss=.false., sort=.false., unique=.false., &
789  ltimerange=(this%timerange(:)%timerange /= tri .OR. &
790  this%timerange(:)%p1 /= 0 .OR. this%timerange(:)%p2 /= steps))
791 ENDIF
792 
793 END SUBROUTINE vol7d_decompute_stat_proc
794 
795 
822 SUBROUTINE vol7d_recompute_stat_proc_diff(this, that, stat_proc, step, full_steps, start, other)
823 TYPE(vol7d),INTENT(inout) :: this
824 TYPE(vol7d),INTENT(out) :: that
825 INTEGER,INTENT(in) :: stat_proc
826 TYPE(timedelta),INTENT(in) :: step
827 LOGICAL,INTENT(in),OPTIONAL :: full_steps
828 TYPE(datetime),INTENT(in),OPTIONAL :: start
829 TYPE(vol7d),INTENT(out),OPTIONAL :: other
830 
831 INTEGER :: i1, i3, i5, i6, i, j, k, l, nitr, steps
832 INTEGER,ALLOCATABLE :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
833 LOGICAL,ALLOCATABLE :: mask_timerange(:)
834 LOGICAL,ALLOCATABLE :: mask_time(:)
835 TYPE(vol7d) :: v7dtmp
836 
837 
838 ! be safe
839 CALL vol7d_alloc_vol(this)
840 ! initialise the template of the output volume
841 CALL init(that, time_definition=this%time_definition)
842 
843 ! compute length of cumulation step in seconds
844 CALL getval(step, asec=steps)
845 
846 ! compute the statistical processing relations, output time and
847 ! timerange are defined here
848 CALL recompute_stat_proc_diff_common(this%time, this%timerange, stat_proc, step, &
849  that%time, that%timerange, map_tr, f, keep_tr, &
850  this%time_definition, full_steps, start)
851 nitr = SIZE(f)
852 
853 ! complete the definition of the empty output template
854 CALL vol7d_alloc(that, nana=0, nlevel=0, nnetwork=0)
855 CALL vol7d_alloc_vol(that)
856 ! shouldn't we exit here with an empty volume if stat_proc/=0,1 ?
857 ALLOCATE(mask_time(SIZE(this%time)), mask_timerange(SIZE(this%timerange)))
858 DO l = 1, SIZE(this%time)
859  mask_time(l) = any(this%time(l) == that%time(:))
860 ENDDO
861 DO l = 1, SIZE(this%timerange)
862  mask_timerange(l) = any(this%timerange(l) == that%timerange(:))
863 ENDDO
864 ! create template for the output volume, keep all ana, level, network
865 ! and variables; copy only the timeranges already satisfying the
866 ! requested step, if any and only the times already existing in the
867 ! output
868 CALL vol7d_copy(this, v7dtmp, miss=.false., sort=.false., unique=.false., &
869  ltimerange=mask_timerange(:), ltime=mask_time(:))
870 ! merge output created so far with template
871 CALL vol7d_merge(that, v7dtmp, lanasimple=.true., llevelsimple=.true.)
872 
873 ! compute statistical processing
874 IF (ASSOCIATED(this%voldatir)) THEN
875  DO l = 1, SIZE(this%time)
876  DO k = 1, nitr
877  DO j = 1, SIZE(this%time)
878  DO i = 1, nitr
879  IF (c_e(map_tr(i,j,k,l,1))) THEN
880  DO i6 = 1, SIZE(this%network)
881  DO i5 = 1, SIZE(this%dativar%r)
882  DO i3 = 1, SIZE(this%level)
883  DO i1 = 1, SIZE(this%ana)
884  IF (c_e(this%voldatir(i1,l,i3,f(k),i5,i6)) .AND. &
885  c_e(this%voldatir(i1,j,i3,f(i),i5,i6))) THEN
886 
887  IF (stat_proc == 0) THEN ! average
888  that%voldatir( &
889  i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
890  (this%voldatir(i1,l,i3,f(k),i5,i6)*this%timerange(f(k))%p2 - &
891  this%voldatir(i1,j,i3,f(i),i5,i6)*this%timerange(f(i))%p2)/ &
892  steps ! optimize avoiding conversions
893  ELSE IF (stat_proc == 1 .OR. stat_proc == 4) THEN ! acc, diff
894  that%voldatir( &
895  i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
896  this%voldatir(i1,l,i3,f(k),i5,i6) - &
897  this%voldatir(i1,j,i3,f(i),i5,i6)
898  ENDIF
899 
900  ENDIF
901  ENDDO
902  ENDDO
903  ENDDO
904  ENDDO
905  ENDIF
906  ENDDO
907  ENDDO
908  ENDDO
909  ENDDO
910 ENDIF
911 
912 IF (ASSOCIATED(this%voldatid)) THEN
913  DO l = 1, SIZE(this%time)
914  DO k = 1, nitr
915  DO j = 1, SIZE(this%time)
916  DO i = 1, nitr
917  IF (c_e(map_tr(i,j,k,l,1))) THEN
918  DO i6 = 1, SIZE(this%network)
919  DO i5 = 1, SIZE(this%dativar%d)
920  DO i3 = 1, SIZE(this%level)
921  DO i1 = 1, SIZE(this%ana)
922  IF (c_e(this%voldatid(i1,l,i3,f(k),i5,i6)) .AND. &
923  c_e(this%voldatid(i1,j,i3,f(i),i5,i6))) THEN
924 ! IF (.NOT.c_e(that%voldatid( &
925 ! i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6))) THEN
926 
927  IF (stat_proc == 0) THEN ! average
928  that%voldatid( &
929  i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
930  (this%voldatid(i1,l,i3,f(k),i5,i6)*this%timerange(f(k))%p2 - &
931  this%voldatid(i1,j,i3,f(i),i5,i6)*this%timerange(f(i))%p2)/ &
932  steps ! optimize avoiding conversions
933  ELSE IF (stat_proc == 1 .OR. stat_proc == 4) THEN ! acc, diff
934  that%voldatid( &
935  i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
936  this%voldatid(i1,l,i3,f(k),i5,i6) - &
937  this%voldatid(i1,j,i3,f(i),i5,i6)
938  ENDIF
939 
940 ! ENDIF
941  ENDIF
942  ENDDO
943  ENDDO
944  ENDDO
945  ENDDO
946  ENDIF
947  ENDDO
948  ENDDO
949  ENDDO
950  ENDDO
951 ENDIF
952 
953 ! this should be avoided by sorting descriptors upstream
954 ! descriptors now are sorted upstream with a dirty and expensive trick
955 ! but the order may be scrambled in the call to vol7d_merge above
956 CALL vol7d_smart_sort(that, lsort_time=.true., lsort_timerange=.true.)
957 
958 CALL makeother(.true.)
959 
960 CONTAINS
961 
962 SUBROUTINE makeother(filter)
963 LOGICAL,INTENT(in) :: filter
964 IF (PRESENT(other)) THEN
965  IF (filter) THEN ! create volume with the remaining data for further processing
966  CALL vol7d_copy(this, other, miss=.false., sort=.false., unique=.false., &
967  ltimerange=(this%timerange(:)%timerange /= stat_proc))
968  ELSE
969  CALL vol7d_copy(this, other, miss=.false., sort=.false., unique=.false.)
970  ENDIF
971 ENDIF
972 END SUBROUTINE makeother
973 
974 END SUBROUTINE vol7d_recompute_stat_proc_diff
975 
976 
1004 SUBROUTINE vol7d_compute_stat_proc_metamorph(this, that, stat_proc_input, stat_proc)
1005 TYPE(vol7d),INTENT(inout) :: this
1006 TYPE(vol7d),INTENT(out) :: that
1007 INTEGER,INTENT(in) :: stat_proc_input
1008 INTEGER,INTENT(in) :: stat_proc
1009 
1010 INTEGER :: j
1011 LOGICAL,ALLOCATABLE :: tr_mask(:)
1012 REAL,ALLOCATABLE :: int_ratio(:)
1013 DOUBLE PRECISION,ALLOCATABLE :: int_ratiod(:)
1014 
1015 IF (.NOT.((stat_proc_input == 0 .AND. stat_proc == 1) .OR. &
1016  (stat_proc_input == 1 .AND. stat_proc == 0))) THEN
1017 
1018  CALL l4f_log(l4f_warn, &
1019  'compute_stat_proc_metamorph, can only be applied to average->accumulated timerange and viceversa')
1020 ! return an empty volume, without signaling error
1021  CALL init(that)
1022  CALL vol7d_alloc_vol(that)
1023  RETURN
1024 ENDIF
1025 
1026 ! be safe
1027 CALL vol7d_alloc_vol(this)
1028 
1029 ! useful timeranges
1030 tr_mask = this%timerange(:)%timerange == stat_proc_input .AND. this%timerange(:)%p2 /= imiss &
1031  .AND. this%timerange(:)%p2 /= 0
1032 
1033 ! initial check (necessary?)
1034 IF (count(tr_mask) == 0) THEN
1035  CALL l4f_log(l4f_warn, &
1036  'vol7d_compute, no timeranges suitable for statistical processing by metamorphosis')
1037  CALL init(that)
1038 ! CALL makeother()
1039  RETURN
1040 ENDIF
1041 
1042 ! copy required data and reset timerange
1043 CALL vol7d_copy(this, that, ltimerange=tr_mask)
1044 that%timerange(:)%timerange = stat_proc
1045 ! why next automatic f2003 allocation does not always work?
1046 ALLOCATE(int_ratio(SIZE(that%timerange)), int_ratiod(SIZE(that%timerange)))
1047 
1048 IF (stat_proc == 0) THEN ! average -> integral
1049  int_ratio = 1./real(that%timerange(:)%p2)
1050  int_ratiod = 1./dble(that%timerange(:)%p2)
1051 ELSE ! cumulation
1052  int_ratio = real(that%timerange(:)%p2)
1053  int_ratiod = dble(that%timerange(:)%p2)
1054 ENDIF
1055 
1056 IF (ASSOCIATED(that%voldatir)) THEN
1057  DO j = 1, SIZE(that%timerange)
1058  WHERE(c_e(that%voldatir(:,:,:,j,:,:)))
1059  that%voldatir(:,:,:,j,:,:) = that%voldatir(:,:,:,j,:,:)*int_ratio(j)
1060  ELSEWHERE
1061  that%voldatir(:,:,:,j,:,:) = rmiss
1062  END WHERE
1063  ENDDO
1064 ENDIF
1065 
1066 IF (ASSOCIATED(that%voldatid)) THEN
1067  DO j = 1, SIZE(that%timerange)
1068  WHERE(c_e(that%voldatid(:,:,:,j,:,:)))
1069  that%voldatid(:,:,:,j,:,:) = that%voldatid(:,:,:,j,:,:)*int_ratiod(j)
1070  ELSEWHERE
1071  that%voldatid(:,:,:,j,:,:) = rmiss
1072  END WHERE
1073  ENDDO
1074 ENDIF
1075 
1076 
1077 END SUBROUTINE vol7d_compute_stat_proc_metamorph
1078 
1079 
1080 SUBROUTINE vol7d_recompute_stat_proc_agg_multiv(this, that, &
1081  step, start, frac_valid, multiv_proc)
1082 TYPE(vol7d),INTENT(inout) :: this
1083 TYPE(vol7d),INTENT(out) :: that
1084 !INTEGER,INTENT(in) :: stat_proc !< type of statistical processing to be recomputed (from grib2 table), only data having timerange of this type will be recomputed and will appear in the output volume
1085 TYPE(timedelta),INTENT(in) :: step
1086 TYPE(datetime),INTENT(in),OPTIONAL :: start
1087 REAL,INTENT(in),OPTIONAL :: frac_valid
1088 !TYPE(vol7d),INTENT(inout),OPTIONAL :: other !< optional volume that, on exit, is going to contain the data that did not contribute to the statistical processing
1089 !INTEGER,INTENT(in),OPTIONAL :: stat_proc_input !< to be used with care, type of statistical processing of data that has to be processed (from grib2 table), only data having timerange of this type will be recomputed, the actual statistical processing performed and which will appear in the output volume, is however determined by \a stat_proc argument
1090 INTEGER,INTENT(in) :: multiv_proc
1091 
1092 INTEGER :: tri
1093 INTEGER :: i, j, n, n1, ndtr, i1, i3, i5, i6
1094 INTEGER :: linshape(1)
1095 REAL :: lfrac_valid, frac_c, frac_m
1096 LOGICAL,ALLOCATABLE :: ttr_mask(:,:)
1097 TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
1098 INTEGER,POINTER :: dtratio(:)
1099 INTEGER :: stat_proc_input, stat_proc
1100 
1101 SELECT CASE(multiv_proc)
1102 CASE (1) ! direction of maximum
1103  stat_proc_input = 205
1104  stat_proc=205
1105 END SELECT
1106 
1107 tri = stat_proc_input
1108 IF (PRESENT(frac_valid)) THEN
1109  lfrac_valid = frac_valid
1110 ELSE
1111  lfrac_valid = 1.0
1112 ENDIF
1113 
1114 ! be safe
1115 CALL vol7d_alloc_vol(this)
1116 ! initial check
1117 
1118 ! cleanup the original volume
1119 CALL vol7d_smart_sort(this, lsort_time=.true.) ! time-ordered volume needed
1120 CALL vol7d_reform(this, miss=.false., sort=.false., unique=.true.)
1121 
1122 CALL init(that, time_definition=this%time_definition)
1123 CALL vol7d_alloc(that, nana=SIZE(this%ana), nlevel=SIZE(this%level), &
1124  nnetwork=SIZE(this%network))
1125 IF (ASSOCIATED(this%dativar%r)) THEN
1126  CALL vol7d_alloc(that, ndativarr=SIZE(this%dativar%r))
1127  that%dativar%r = this%dativar%r
1128 ENDIF
1129 IF (ASSOCIATED(this%dativar%d)) THEN
1130  CALL vol7d_alloc(that, ndativard=SIZE(this%dativar%d))
1131  that%dativar%d = this%dativar%d
1132 ENDIF
1133 that%ana = this%ana
1134 that%level = this%level
1135 that%network = this%network
1136 
1137 ! compute the output time and timerange and all the required mappings
1138 CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
1139  step, this%time_definition, that%time, that%timerange, map_ttr, &
1140  dtratio=dtratio, start=start)
1141 CALL vol7d_alloc_vol(that)
1142 
1143 ALLOCATE(ttr_mask(SIZE(this%time), SIZE(this%timerange)))
1144 linshape = (/SIZE(ttr_mask)/)
1145 ! finally perform computations
1146 IF (ASSOCIATED(this%voldatir)) THEN
1147  DO j = 1, SIZE(that%timerange)
1148  DO i = 1, SIZE(that%time)
1149 
1150  DO i1 = 1, SIZE(this%ana)
1151  DO i3 = 1, SIZE(this%level)
1152  DO i6 = 1, SIZE(this%network)
1153  DO i5 = 1, SIZE(this%dativar%r)
1154 
1155  frac_m = 0.
1156  DO n1 = SIZE(dtratio), 1, -1 ! precedence to longer periods
1157  IF (dtratio(n1) <= 0) cycle ! safety check
1158  ttr_mask = .false.
1159  DO n = 1, map_ttr(i,j)%arraysize
1160  IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
1161  IF (c_e(this%voldatir(i1,map_ttr(i,j)%array(n)%it,i3, &
1162  map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
1163  ttr_mask(map_ttr(i,j)%array(n)%it, &
1164  map_ttr(i,j)%array(n)%itr) = .true.
1165  ENDIF
1166  ENDIF
1167  ENDDO
1168 
1169  ndtr = count(ttr_mask)
1170  frac_c = real(ndtr)/real(dtratio(n1))
1171 
1172  IF (ndtr > 0 .AND. frac_c >= max(lfrac_valid, frac_m)) THEN
1173  frac_m = frac_c
1174  SELECT CASE(multiv_proc)
1175  CASE (1) ! average, vectorial mean
1176  that%voldatir(i1,i,i3,j,i5,i6) = &
1177  sum(this%voldatir(i1,:,i3,:,i5,i6), &
1178  mask=ttr_mask)/ndtr
1179  END SELECT
1180  ENDIF
1181 
1182  ENDDO ! dtratio
1183  ENDDO ! var
1184  ENDDO ! network
1185  ENDDO ! level
1186  ENDDO ! ana
1187  CALL delete(map_ttr(i,j))
1188  ENDDO ! otime
1189  ENDDO ! otimerange
1190 ENDIF
1191 
1192 IF (ASSOCIATED(this%voldatid)) THEN
1193  DO j = 1, SIZE(that%timerange)
1194  DO i = 1, SIZE(that%time)
1195 
1196  DO i1 = 1, SIZE(this%ana)
1197  DO i3 = 1, SIZE(this%level)
1198  DO i6 = 1, SIZE(this%network)
1199  DO i5 = 1, SIZE(this%dativar%d)
1200 
1201  frac_m = 0.
1202  DO n1 = SIZE(dtratio), 1, -1 ! precedence to longer periods
1203  IF (dtratio(n1) <= 0) cycle ! safety check
1204  ttr_mask = .false.
1205  DO n = 1, map_ttr(i,j)%arraysize
1206  IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
1207  IF (c_e(this%voldatid(i1,map_ttr(i,j)%array(n)%it,i3, &
1208  map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
1209  ttr_mask(map_ttr(i,j)%array(n)%it, &
1210  map_ttr(i,j)%array(n)%itr) = .true.
1211  ENDIF
1212  ENDIF
1213  ENDDO
1214 
1215  ndtr = count(ttr_mask)
1216  frac_c = real(ndtr)/real(dtratio(n1))
1217 
1218  IF (ndtr > 0 .AND. frac_c >= max(lfrac_valid, frac_m)) THEN
1219  frac_m = frac_c
1220  SELECT CASE(stat_proc)
1221  CASE (0) ! average
1222  that%voldatid(i1,i,i3,j,i5,i6) = &
1223  sum(this%voldatid(i1,:,i3,:,i5,i6), &
1224  mask=ttr_mask)/ndtr
1225  CASE (1, 4) ! accumulation, difference
1226  that%voldatid(i1,i,i3,j,i5,i6) = &
1227  sum(this%voldatid(i1,:,i3,:,i5,i6), &
1228  mask=ttr_mask)
1229  CASE (2) ! maximum
1230  that%voldatid(i1,i,i3,j,i5,i6) = &
1231  maxval(this%voldatid(i1,:,i3,:,i5,i6), &
1232  mask=ttr_mask)
1233  CASE (3) ! minimum
1234  that%voldatid(i1,i,i3,j,i5,i6) = &
1235  minval(this%voldatid(i1,:,i3,:,i5,i6), &
1236  mask=ttr_mask)
1237  CASE (6) ! stddev
1238  that%voldatid(i1,i,i3,j,i5,i6) = &
1239  stat_stddev( &
1240  reshape(this%voldatid(i1,:,i3,:,i5,i6), shape=linshape), &
1241  mask=reshape(ttr_mask, shape=linshape))
1242  END SELECT
1243  ENDIF
1244 
1245  ENDDO ! dtratio
1246  ENDDO ! var
1247  ENDDO ! network
1248  ENDDO ! level
1249  ENDDO ! ana
1250  CALL delete(map_ttr(i,j))
1251  ENDDO ! otime
1252  ENDDO ! otimerange
1253 ENDIF
1254 
1255 DEALLOCATE(ttr_mask)
1256 
1257 END SUBROUTINE vol7d_recompute_stat_proc_agg_multiv
1258 
1275 SUBROUTINE vol7d_fill_time(this, that, step, start, stopp, cyclicdt)
1276 TYPE(vol7d),INTENT(inout) :: this
1277 TYPE(vol7d),INTENT(inout) :: that
1278 TYPE(timedelta),INTENT(in) :: step
1279 TYPE(datetime),INTENT(in),OPTIONAL :: start
1280 TYPE(datetime),INTENT(in),OPTIONAL :: stopp
1281 TYPE(cyclicdatetime),INTENT(in),OPTIONAL :: cyclicdt
1282 
1283 TYPE(cyclicdatetime) :: lcyclicdt
1284 TYPE(datetime) :: counter, lstart, lstop
1285 INTEGER :: i, naddtime
1286 
1287 CALL safe_start_stop(this, lstart, lstop, start, stopp)
1288 IF (.NOT. c_e(lstart) .OR. .NOT. c_e(lstop) .OR. .NOT. c_e(step)) RETURN
1289 
1290 lcyclicdt=cyclicdatetime_miss
1291 if (present(cyclicdt)) then
1292  if(c_e(cyclicdt)) lcyclicdt=cyclicdt
1293 end if
1294 
1295 CALL l4f_log(l4f_info, 'vol7d_fill_time: time interval '//trim(to_char(lstart))// &
1296  ' '//trim(to_char(lstop)))
1297 
1298 ! Count the number of time levels required for completing the series
1299 ! valid also in the case (SIZE(this%time) == 0)
1300 naddtime = 0
1301 counter = lstart
1302 i = 1
1303 naddcount: DO WHILE(counter <= lstop)
1304  DO WHILE(i <= SIZE(this%time)) ! this%time(i) chases counter
1305  IF (counter < this%time(i)) THEN ! this%time(i) overtook counter
1306  i = max(i-1,1) ! go back if possible
1307  EXIT
1308  ELSE IF (counter == this%time(i) .OR. .NOT. counter == lcyclicdt) THEN ! found matching time
1309  counter = counter + step
1310  cycle naddcount
1311  ENDIF
1312  i = i + 1
1313  ENDDO
1314  naddtime = naddtime + 1
1315  counter = counter + step
1316 ENDDO naddcount
1317 
1318 ! old universal algorithm, not optimized, check that the new one is equivalent
1319 !naddtime = 0
1320 !counter = lstart
1321 !DO WHILE(counter <= lstop)
1322 ! IF (.NOT.ANY(counter == this%time(:))) THEN
1323 ! naddtime = naddtime + 1
1324 ! ENDIF
1325 ! counter = counter + step
1326 !ENDDO
1327 
1328 IF (naddtime > 0) THEN
1329 
1330  CALL init(that)
1331  CALL vol7d_alloc(that, ntime=naddtime)
1332  CALL vol7d_alloc_vol(that)
1333 
1334  ! Repeat the count loop setting the time levels to be added
1335  naddtime = 0
1336  counter = lstart
1337  i = 1
1338  naddadd: DO WHILE(counter <= lstop)
1339  DO WHILE(i <= SIZE(this%time)) ! this%time(i) chases counter
1340  IF (counter < this%time(i)) THEN ! this%time(i) overtook counter
1341  i = max(i-1,1) ! go back if possible
1342  EXIT
1343  ELSE IF (counter == this%time(i) .OR. .NOT. counter == lcyclicdt) THEN ! found matching time
1344  counter = counter + step
1345  cycle naddadd
1346  ENDIF
1347  i = i + 1
1348  ENDDO
1349  naddtime = naddtime + 1
1350  that%time(naddtime) = counter ! only difference
1351  counter = counter + step
1352  ENDDO naddadd
1353 
1354  CALL vol7d_append(that, this, sort=.true.)
1355 
1356 ELSE
1357 !! ? why sort all dimension ?
1358 !! CALL vol7d_copy(this, that, lsort_time=.TRUE.)
1359  CALL vol7d_copy(this, that, sort=.true.)
1360 ENDIF
1361 
1362 
1363 END SUBROUTINE vol7d_fill_time
1364 
1365 
1377 SUBROUTINE vol7d_filter_time(this, that, step, start, stopp, cyclicdt)
1378 TYPE(vol7d),INTENT(inout) :: this
1379 TYPE(vol7d),INTENT(inout) :: that
1380 TYPE(timedelta),INTENT(in),optional :: step
1381 TYPE(datetime),INTENT(in),OPTIONAL :: start
1382 TYPE(datetime),INTENT(in),OPTIONAL :: stopp
1383 TYPE(cyclicdatetime),INTENT(in),OPTIONAL :: cyclicdt
1384 
1385 TYPE(datetime) :: lstart, lstop
1386 LOGICAL, ALLOCATABLE :: time_mask(:)
1387 
1388 CALL safe_start_stop(this, lstart, lstop, start, stopp)
1389 IF (.NOT. c_e(lstart) .OR. .NOT. c_e(lstop)) RETURN
1390 
1391 CALL l4f_log(l4f_info, 'vol7d_filter_time: time interval '//trim(to_char(lstart))// &
1392  ' '//trim(to_char(lstop)))
1393 
1394 ALLOCATE(time_mask(SIZE(this%time)))
1395 
1396 time_mask = this%time >= lstart .AND. this%time <= lstop
1397 
1398 IF (PRESENT(cyclicdt)) THEN
1399  IF (c_e(cyclicdt)) THEN
1400  time_mask = time_mask .AND. this%time == cyclicdt
1401  ENDIF
1402 ENDIF
1403 
1404 IF (PRESENT(step)) THEN
1405  IF (c_e(step)) THEN
1406  time_mask = time_mask .AND. mod(this%time - lstart, step) == timedelta_0
1407  ENDIF
1408 ENDIF
1409 
1410 CALL vol7d_copy(this,that, ltime=time_mask)
1411 
1412 DEALLOCATE(time_mask)
1413 
1414 END SUBROUTINE vol7d_filter_time
1415 
1416 
1420 SUBROUTINE vol7d_fill_data(this, step, start, stopp, tolerance)
1421 TYPE(vol7d),INTENT(inout) :: this
1422 TYPE(timedelta),INTENT(in) :: step
1423 TYPE(datetime),INTENT(in),OPTIONAL :: start
1424 TYPE(datetime),INTENT(in),OPTIONAL :: stopp
1425 TYPE(timedelta),INTENT(in),optional :: tolerance
1426 
1427 TYPE(datetime) :: lstart, lstop
1428 integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork, iindtime
1429 type(timedelta) :: deltato,deltat, ltolerance
1430 
1431 CALL safe_start_stop(this, lstart, lstop, start, stopp)
1432 IF (.NOT. c_e(lstart) .OR. .NOT. c_e(lstop)) RETURN
1433 
1434 CALL l4f_log(l4f_info, 'vol7d_fill_data: time interval '//trim(to_char(lstart))// &
1435  ' '//trim(to_char(lstop)))
1436 
1437 
1438 ltolerance=step/2
1439 
1440 if (present(tolerance))then
1441  if (c_e(tolerance)) ltolerance=tolerance
1442 end if
1443 
1444 
1445 do indtime=1,size(this%time)
1446 
1447  IF (this%time(indtime) < lstart .OR. this%time(indtime) > lstop .OR. &
1448  mod(this%time(indtime) - lstart, step) /= timedelta_0) cycle
1449  do indtimerange=1,size(this%timerange)
1450  if (this%timerange(indtimerange)%timerange /= 254) cycle
1451  do indnetwork=1,size(this%network)
1452  do inddativarr=1,size(this%dativar%r)
1453  do indlevel=1,size(this%level)
1454  do indana=1,size(this%ana)
1455 
1456  !find the nearest data in time if data is missing
1457  if (.not. c_e(this%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)))then
1458  deltato=timedelta_miss
1459 
1460  !do iindtime=max(indtime-20,1),min(indtime+20,size(this%time)) !check on a chunk: 20 should be enought
1461 
1462  do iindtime=indtime+1,size(this%time) !check forward
1463 
1464  if (c_e(this%voldatir (indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork )))then
1465  deltat=this%time(iindtime)-this%time(indtime)
1466 
1467  if (deltat >= ltolerance) exit
1468 
1469  if (deltat < deltato) then
1470  this%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork) = &
1471  this%voldatir(indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork)
1472  deltato=deltat
1473  end if
1474  end if
1475  end do
1476 
1477  do iindtime=indtime-1,1,-1 !check backward
1478 
1479  if (c_e(this%voldatir (indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork )))then
1480  if (iindtime < indtime) then
1481  deltat=this%time(indtime)-this%time(iindtime)
1482  else if (iindtime > indtime) then
1483  deltat=this%time(iindtime)-this%time(indtime)
1484  else
1485  cycle
1486  end if
1487 
1488  if (deltat >= ltolerance) exit
1489 
1490  if (deltat < deltato) then
1491  this%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork) = &
1492  this%voldatir(indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork)
1493  deltato=deltat
1494  end if
1495  end if
1496  end do
1497 
1498  end if
1499  end do
1500  end do
1501  end do
1502  end do
1503  end do
1504 end do
1505 
1506 END SUBROUTINE vol7d_fill_data
1507 
1508 
1509 ! private utility routine for checking interval and start-stop times
1510 ! in input missing start-stop values are treated as not present
1511 ! in output missing start-stop values mean "do nothing"
1512 SUBROUTINE safe_start_stop(this, lstart, lstop, start, stopp)
1513 TYPE(vol7d),INTENT(inout) :: this
1514 TYPE(datetime),INTENT(out) :: lstart
1515 TYPE(datetime),INTENT(out) :: lstop
1516 TYPE(datetime),INTENT(in),OPTIONAL :: start
1517 TYPE(datetime),INTENT(in),OPTIONAL :: stopp
1518 
1519 lstart = datetime_miss
1520 lstop = datetime_miss
1521 ! initial safety operation
1522 CALL vol7d_alloc_vol(this)
1523 IF (SIZE(this%time) == 0) RETURN ! avoid segmentation fault in case of empty volume
1524 CALL vol7d_smart_sort(this, lsort_time=.true.)
1525 
1526 IF (PRESENT(start)) THEN
1527  IF (c_e(start)) THEN
1528  lstart = start
1529  ELSE
1530  lstart = this%time(1)
1531  ENDIF
1532 ELSE
1533  lstart = this%time(1)
1534 ENDIF
1535 IF (PRESENT(stopp)) THEN
1536  IF (c_e(stopp)) THEN
1537  lstop = stopp
1538  ELSE
1539  lstop = this%time(SIZE(this%time))
1540  ENDIF
1541 ELSE
1542  lstop = this%time(SIZE(this%time))
1543 ENDIF
1544 
1545 END SUBROUTINE safe_start_stop
1546 
1547 
1554 SUBROUTINE vol7d_normalize_vcoord(this,that,ana,time,timerange,network)
1555 TYPE(vol7d),INTENT(INOUT) :: this
1556 TYPE(vol7d),INTENT(OUT) :: that
1557 integer,intent(in) :: time,ana,timerange,network
1558 
1559 character(len=1) :: type
1560 integer :: ind
1561 TYPE(vol7d_var) :: var
1562 LOGICAL,allocatable :: ltime(:),ltimerange(:),lana(:),lnetwork(:)
1563 logical,allocatable :: maschera(:)
1564 
1565 
1566 allocate(ltime(size(this%time)))
1567 allocate(ltimerange(size(this%timerange)))
1568 allocate(lana(size(this%ana)))
1569 allocate(lnetwork(size(this%network)))
1570 
1571 ltime=.false.
1572 ltimerange=.false.
1573 lana=.false.
1574 lnetwork=.false.
1576 ltime(time)=.true.
1577 ltimerange(timerange)=.true.
1578 lana(ana)=.true.
1579 lnetwork(network)=.true.
1580 
1581 call vol7d_copy(this, that,unique=.true.,&
1582  ltime=ltime,ltimerange=ltimerange,lana=lana,lnetwork=lnetwork )
1583 
1584 call init(var, btable="B10004") ! Pressure
1585 type=cmiss
1586 !type="i"
1587 ind = index(that%dativar, var, type=type)
1588 
1589 allocate(maschera(size(that%level)))
1590 
1591 maschera = (&
1592  (that%level%level1 == 105.and.that%level%level2 == 105) .or. &
1593  (that%level%level1 == 103 .and. that%level%level2 == imiss ) .or. &
1594  (that%level%level1 == 102 .and. that%level%level2 == imiss )) &
1595  .and. c_e(that%voldatic(1,1,:,1,ind,1))
1596 
1597 
1598 select case (type)
1599 
1600 case("d")
1601 
1602  where (maschera)
1603  that%level%level1 = 100
1604  that%level%l1 = int(realdat(that%voldatid(1,1,:,1,ind,1),that%dativar%d(ind)))
1605  that%level%l1 = int(that%voldatid(1,1,:,1,ind,1))
1606  that%level%level2 = imiss
1607  that%level%l2 = imiss
1608  end where
1609 
1610 case("r")
1611 
1612  where (maschera)
1613  that%level%level1 = 100
1614  that%level%l1 = int(realdat(that%voldatir(1,1,:,1,ind,1),that%dativar%r(ind)))
1615  that%level%level2 = imiss
1616  that%level%l2 = imiss
1617  end where
1619 case("i")
1620 
1621  where (maschera)
1622  that%level%level1 = 100
1623  that%level%l1 = int(realdat(that%voldatii(1,1,:,1,ind,1),that%dativar%i(ind)))
1624  that%level%level2 = imiss
1625  that%level%l2 = imiss
1626  end where
1627 
1628 case("b")
1629 
1630  where (maschera)
1631  that%level%level1 = 100
1632  that%level%l1 = int(realdat(that%voldatib(1,1,:,1,ind,1),that%dativar%b(ind)))
1633  that%level%level2 = imiss
1634  that%level%l2 = imiss
1635  end where
1636 
1637 case("c")
1638 
1639  where (maschera)
1640  that%level%level1 = 100
1641  that%level%l1 = int(realdat(that%voldatic(1,1,:,1,ind,1),that%dativar%c(ind)))
1642  that%level%level2 = imiss
1643  that%level%l2 = imiss
1644  end where
1645 
1646 end select
1647 
1648 deallocate(ltime)
1649 deallocate(ltimerange)
1650 deallocate(lana)
1651 deallocate(lnetwork)
1652 
1653 END SUBROUTINE vol7d_normalize_vcoord
1654 
1655 
1656 !!$!> Metodo per calcolare variabili derivate.
1657 !!$!! TO DO !!
1658 !!$SUBROUTINE vol7d_compute_var(this,that,var)
1659 !!$TYPE(vol7d),INTENT(INOUT) :: this !< oggetto da normalizzare
1660 !!$TYPE(vol7d),INTENT(OUT) :: that !< oggetto normalizzato
1661 !!$
1662 !!$character(len=1) :: type
1663 !!$TYPE(vol7d_var),intent(in) :: var
1664 
1665 
1666 !!$call init(var, btable="B10004") ! Pressure
1667 !!$type=cmiss
1668 !!$call vol7d_varvect_index(that%dativar,var , type=type,index_v=ind)
1669 !!$
1670 !!$select case (type)
1671 !!$
1672 !!$case("d")
1673 !!$
1674 !!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
1675 !!$ that%level%level1 = 100
1676 !!$ that%level%l1 = realdat(that%voldatid(1,1,:,1,ind,1),that%dativar%d(ind))
1677 !!$ that%level%level2 = imiss
1678 !!$ that%level%l2 = imiss
1679 !!$ end where
1680 !!$
1681 !!$case("r")
1682 !!$
1683 !!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
1684 !!$ that%level%level1 = 100
1685 !!$ that%level%l1 = realdat(that%voldatir(1,1,:,1,ind,1),that%dativar%r(ind))
1686 !!$ that%level%level2 = imiss
1687 !!$ that%level%l2 = imiss
1688 !!$ end where
1689 !!$
1690 !!$case("i")
1691 !!$
1692 !!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
1693 !!$ that%level%level1 = 100
1694 !!$ that%level%l1 = realdat(that%voldatii(1,1,:,1,ind,1),that%dativar%i(ind))
1695 !!$ that%level%level2 = imiss
1696 !!$ that%level%l2 = imiss
1697 !!$ end where
1698 !!$
1699 !!$case("b")
1700 !!$
1701 !!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
1702 !!$ that%level%level1 = 100
1703 !!$ that%level%l1 = realdat(that%voldatib(1,1,:,1,ind,1),that%dativar%b(ind))
1704 !!$ that%level%level2 = imiss
1705 !!$ that%level%l2 = imiss
1706 !!$ end where
1707 !!$
1708 !!$case("c")
1709 !!$
1710 !!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
1711 !!$ that%level%level1 = 100
1712 !!$ that%level%l1 = realdat(that%voldatic(1,1,:,1,ind,1),that%dativar%c(ind))
1713 !!$ that%level%level2 = imiss
1714 !!$ that%level%l2 = imiss
1715 !!$ end where
1716 !!$
1717 !!$end select
1718 
1719 !!$
1720 !!$END SUBROUTINE vol7d_compute_var
1721 !!$
1722 
1723 
1724 
1725 END MODULE vol7d_class_compute
Distruttori per le 2 classi.
Restituiscono il valore dell'oggetto nella forma desiderata.
Costruttori per le classi datetime e timedelta.
Operatore di resto della divisione.
Functions that return a trimmed CHARACTER representation of the input variable.
Restituiscono il valore dell'oggetto in forma di stringa stampabile.
Index method.
Compute the mode of the random variable provided taking into account missing data.
Compute the standard deviation of the random variable provided, taking into account missing data.
Definition: simple_stat.f90:69
real data conversion
Classi per la gestione delle coordinate temporali.
Module for basic statistical computations taking into account missing data.
Definition: simple_stat.f90:25
This module contains functions that are only for internal use of the library.
Extension of vol7d_class with methods for performing simple statistical operations on entire volumes ...
Classe per la gestione di un volume completo di dati osservati.
Class for expressing a cyclic datetime.
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.
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...

Generated with Doxygen.