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