libsim  Versione6.3.0
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, 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
143  CALL vol7d_recompute_stat_proc_agg(this, that1, stat_proc_input, &
144  step, frac_valid=1.0)
145  CALL vol7d_recompute_stat_proc_agg(this, that2, stat_proc_input, &
146  step, start=that1%time(1)+step/2, 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 == 2 .OR. stat_proc == 3)) THEN ! min and max can be computed from any input
160  CALL l4f_log(l4f_info, &
161  'recomputing statistically processed data by aggregation and difference '//&
162  trim(to_char(stat_proc_input))//':'//trim(to_char(stat_proc)))
163 
164  IF (PRESENT(other)) THEN
165  CALL vol7d_recompute_stat_proc_agg(this, that1, stat_proc, &
166  step, start, frac_valid, other=other, stat_proc_input=stat_proc_input)
167  CALL vol7d_recompute_stat_proc_diff(this, that2, stat_proc, &
168  step, full_steps, other=other1)
169  CALL vol7d_merge(other, other1, sort=.true.)
170  ELSE
171  CALL vol7d_recompute_stat_proc_agg(this, that1, stat_proc, &
172  step, start, frac_valid, stat_proc_input=stat_proc_input)
173  CALL vol7d_recompute_stat_proc_diff(this, that2, stat_proc, step, full_steps)
174  ENDIF
175 
176  CALL vol7d_merge(that1, that2, sort=.true.)
177  CALL delete(that2)
178  that = that1
179 ELSE ! IF (stat_proc_input /= stat_proc) THEN
180  IF ((stat_proc_input == 0 .AND. stat_proc == 1) .OR. &
181  (stat_proc_input == 1 .AND. stat_proc == 0)) THEN
182  CALL l4f_log(l4f_info, &
183  'computing statistically processed data by integration/differentiation '// &
184  t2c(stat_proc_input)//':'//t2c(stat_proc))
185  CALL vol7d_compute_stat_proc_metamorph(this, that, stat_proc_input, &
186  stat_proc)
187  ELSE
188  CALL l4f_log(l4f_error, &
189  'statistical processing '//t2c(stat_proc_input)//':'//t2c(stat_proc)// &
190  ' not implemented or does not make sense')
191  CALL raise_error()
192  ENDIF
193 ENDIF
194 
195 END SUBROUTINE vol7d_compute_stat_proc
196 
197 
242 SUBROUTINE vol7d_recompute_stat_proc_agg(this, that, stat_proc, &
243  step, start, frac_valid, other, stat_proc_input)
244 TYPE(vol7d),INTENT(inout) :: this
245 TYPE(vol7d),INTENT(out) :: that
246 INTEGER,INTENT(in) :: stat_proc
247 TYPE(timedelta),INTENT(in) :: step
248 TYPE(datetime),INTENT(in),OPTIONAL :: start
249 REAL,INTENT(in),OPTIONAL :: frac_valid
250 TYPE(vol7d),INTENT(inout),OPTIONAL :: other
251 INTEGER,INTENT(in),OPTIONAL :: stat_proc_input
252 
253 INTEGER :: tri
254 INTEGER :: i, j, n, n1, ndtr, i1, i3, i5, i6
255 INTEGER :: linshape(1)
256 REAL :: lfrac_valid, frac_c, frac_m
257 LOGICAL,ALLOCATABLE :: ttr_mask(:,:)
258 TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
259 INTEGER,POINTER :: dtratio(:)
260 
261 
262 IF (PRESENT(stat_proc_input)) THEN
263  tri = stat_proc_input
264 ELSE
265  tri = stat_proc
266 ENDIF
267 IF (PRESENT(frac_valid)) THEN
268  lfrac_valid = frac_valid
269 ELSE
270  lfrac_valid = 1.0
271 ENDIF
272 
273 ! be safe
274 CALL vol7d_alloc_vol(this)
275 ! initial check
276 
277 ! cleanup the original volume
278 CALL vol7d_smart_sort(this, lsort_time=.true.) ! time-ordered volume needed
279 CALL vol7d_reform(this, miss=.false., sort=.false., unique=.true.)
280 
281 CALL init(that, time_definition=this%time_definition)
282 CALL vol7d_alloc(that, nana=SIZE(this%ana), nlevel=SIZE(this%level), &
283  nnetwork=SIZE(this%network))
284 IF (ASSOCIATED(this%dativar%r)) THEN
285  CALL vol7d_alloc(that, ndativarr=SIZE(this%dativar%r))
286  that%dativar%r = this%dativar%r
287 ENDIF
288 IF (ASSOCIATED(this%dativar%d)) THEN
289  CALL vol7d_alloc(that, ndativard=SIZE(this%dativar%d))
290  that%dativar%d = this%dativar%d
291 ENDIF
292 that%ana = this%ana
293 that%level = this%level
294 that%network = this%network
295 
296 ! compute the output time and timerange and all the required mappings
297 CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
298  step, this%time_definition, that%time, that%timerange, map_ttr, dtratio, start)
299 CALL vol7d_alloc_vol(that)
300 
301 ALLOCATE(ttr_mask(SIZE(this%time), SIZE(this%timerange)))
302 linshape = (/SIZE(ttr_mask)/)
303 ! finally perform computations
304 IF (ASSOCIATED(this%voldatir)) THEN
305  DO j = 1, SIZE(that%timerange)
306  DO i = 1, SIZE(that%time)
307 
308  DO i1 = 1, SIZE(this%ana)
309  DO i3 = 1, SIZE(this%level)
310  DO i6 = 1, SIZE(this%network)
311  DO i5 = 1, SIZE(this%dativar%r)
312 
313  frac_m = 0.
314  DO n1 = SIZE(dtratio), 1, -1 ! precedence to longer periods
315  IF (dtratio(n1) <= 0) cycle ! safety check
316  ttr_mask = .false.
317  DO n = 1, map_ttr(i,j)%arraysize
318  IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
319  IF (c_e(this%voldatir(i1,map_ttr(i,j)%array(n)%it,i3, &
320  map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
321  ttr_mask(map_ttr(i,j)%array(n)%it, &
322  map_ttr(i,j)%array(n)%itr) = .true.
323  ENDIF
324  ENDIF
325  ENDDO
326 
327  ndtr = count(ttr_mask)
328  frac_c = REAL(ndtr)/REAL(dtratio(n1))
329 
330  IF (ndtr > 0 .AND. frac_c >= max(lfrac_valid, frac_m)) THEN
331  frac_m = frac_c
332  SELECT CASE(stat_proc)
333  CASE (0) ! average
334  that%voldatir(i1,i,i3,j,i5,i6) = &
335  sum(this%voldatir(i1,:,i3,:,i5,i6), &
336  mask=ttr_mask)/ndtr
337  CASE (1, 4) ! accumulation, difference
338  that%voldatir(i1,i,i3,j,i5,i6) = &
339  sum(this%voldatir(i1,:,i3,:,i5,i6), &
340  mask=ttr_mask)
341  CASE (2) ! maximum
342  that%voldatir(i1,i,i3,j,i5,i6) = &
343  maxval(this%voldatir(i1,:,i3,:,i5,i6), &
344  mask=ttr_mask)
345  CASE (3) ! minimum
346  that%voldatir(i1,i,i3,j,i5,i6) = &
347  minval(this%voldatir(i1,:,i3,:,i5,i6), &
348  mask=ttr_mask)
349  CASE (6) ! stddev
350  that%voldatir(i1,i,i3,j,i5,i6) = &
351  stat_stddev( &
352  reshape(this%voldatir(i1,:,i3,:,i5,i6), shape=linshape), &
353  mask=reshape(ttr_mask, shape=linshape))
354  END SELECT
355  ENDIF
356 
357  ENDDO ! dtratio
358  ENDDO ! var
359  ENDDO ! network
360  ENDDO ! level
361  ENDDO ! ana
362  CALL delete(map_ttr(i,j))
363  ENDDO ! otime
364  ENDDO ! otimerange
365 ENDIF
366 
367 IF (ASSOCIATED(this%voldatid)) THEN
368  DO j = 1, SIZE(that%timerange)
369  DO i = 1, SIZE(that%time)
370 
371  DO i1 = 1, SIZE(this%ana)
372  DO i3 = 1, SIZE(this%level)
373  DO i6 = 1, SIZE(this%network)
374  DO i5 = 1, SIZE(this%dativar%d)
375 
376  frac_m = 0.
377  DO n1 = SIZE(dtratio), 1, -1 ! precedence to longer periods
378  IF (dtratio(n1) <= 0) cycle ! safety check
379  ttr_mask = .false.
380  DO n = 1, map_ttr(i,j)%arraysize
381  IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
382  IF (c_e(this%voldatid(i1,map_ttr(i,j)%array(n)%it,i3, &
383  map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
384  ttr_mask(map_ttr(i,j)%array(n)%it, &
385  map_ttr(i,j)%array(n)%itr) = .true.
386  ENDIF
387  ENDIF
388  ENDDO
389 
390  ndtr = count(ttr_mask)
391  frac_c = REAL(ndtr)/REAL(dtratio(n1))
392 
393  IF (ndtr > 0 .AND. frac_c >= max(lfrac_valid, frac_m)) THEN
394  frac_m = frac_c
395  SELECT CASE(stat_proc)
396  CASE (0) ! average
397  that%voldatid(i1,i,i3,j,i5,i6) = &
398  sum(this%voldatid(i1,:,i3,:,i5,i6), &
399  mask=ttr_mask)/ndtr
400  CASE (1, 4) ! accumulation, difference
401  that%voldatid(i1,i,i3,j,i5,i6) = &
402  sum(this%voldatid(i1,:,i3,:,i5,i6), &
403  mask=ttr_mask)
404  CASE (2) ! maximum
405  that%voldatid(i1,i,i3,j,i5,i6) = &
406  maxval(this%voldatid(i1,:,i3,:,i5,i6), &
407  mask=ttr_mask)
408  CASE (3) ! minimum
409  that%voldatid(i1,i,i3,j,i5,i6) = &
410  minval(this%voldatid(i1,:,i3,:,i5,i6), &
411  mask=ttr_mask)
412  CASE (6) ! stddev
413  that%voldatid(i1,i,i3,j,i5,i6) = &
414  stat_stddev( &
415  reshape(this%voldatid(i1,:,i3,:,i5,i6), shape=linshape), &
416  mask=reshape(ttr_mask, shape=linshape))
417  END SELECT
418  ENDIF
419 
420  ENDDO ! dtratio
421  ENDDO ! var
422  ENDDO ! network
423  ENDDO ! level
424  ENDDO ! ana
425  CALL delete(map_ttr(i,j))
426  ENDDO ! otime
427  ENDDO ! otimerange
428 ENDIF
429 
430 DEALLOCATE(ttr_mask)
431 
432 CALL makeother()
433 
434 CONTAINS
435 
436 SUBROUTINE makeother()
437 IF (PRESENT(other)) THEN ! create volume with the remaining data for further processing
438  CALL vol7d_copy(this, other, miss=.false., sort=.false., unique=.false., &
439  ltimerange=(this%timerange(:)%timerange /= tri .OR. this%timerange(:)%p2 == imiss &
440  .OR. this%timerange(:)%p2 == 0)) ! or MOD(steps, this%timerange(:)%p2) == 0
441 ENDIF
442 END SUBROUTINE makeother
443 
444 END SUBROUTINE vol7d_recompute_stat_proc_agg
445 
446 
478 SUBROUTINE vol7d_compute_stat_proc_agg(this, that, stat_proc, &
479  step, start, max_step, weighted, other)
480 TYPE(vol7d),INTENT(inout) :: this
481 TYPE(vol7d),INTENT(out) :: that
482 INTEGER,INTENT(in) :: stat_proc
483 TYPE(timedelta),INTENT(in) :: step
484 TYPE(datetime),INTENT(in),OPTIONAL :: start
485 TYPE(timedelta),INTENT(in),OPTIONAL :: max_step
486 LOGICAL,INTENT(in),OPTIONAL :: weighted
487 TYPE(vol7d),INTENT(inout),OPTIONAL :: other
488 !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
489 
490 TYPE(vol7d) :: v7dtmp
491 INTEGER :: tri
492 INTEGER :: i, j, n, ninp, ndtr, i1, i3, i5, i6, vartype, maxsize
493 TYPE(timedelta) :: lmax_step, act_max_step
494 TYPE(datetime) :: pstart, pend, reftime
495 TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
496 REAL,ALLOCATABLE :: tmpvolr(:)
497 DOUBLE PRECISION,ALLOCATABLE :: tmpvold(:), weights(:)
498 LOGICAL,ALLOCATABLE :: lin_mask(:)
499 LOGICAL :: lweighted
500 CHARACTER(len=8) :: env_var
501 
502 IF (PRESENT(max_step)) THEN
503  lmax_step = max_step
504 ELSE
505  lmax_step = timedelta_max
506 ENDIF
507 lweighted = optio_log(weighted)
508 tri = 254
509 ! enable bad behavior for climat database
510 env_var = ''
511 CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
512 lweighted = lweighted .AND. len_trim(env_var) == 0
513 ! only average can be weighted
514 lweighted = lweighted .AND. stat_proc == 0
515 
516 ! be safe
517 CALL vol7d_alloc_vol(this)
518 ! initial check
519 
520 ! cleanup the original volume
521 CALL vol7d_smart_sort(this, lsort_time=.true.) ! time-ordered volume needed
522 CALL vol7d_reform(this, miss=.false., sort=.false., unique=.true.)
523 ! copy everything except time and timerange
524 CALL vol7d_copy(this, v7dtmp, ltime=(/.false./), ltimerange=(/.false./))
525 
526 ! create new volume
527 CALL init(that, time_definition=this%time_definition)
528 ! compute the output time and timerange and all the required mappings
529 CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
530  step, this%time_definition, that%time, that%timerange, map_ttr, start=start)
531 ! merge with information from original volume
532 CALL vol7d_merge(that, v7dtmp)
533 
534 maxsize = maxval(map_ttr(:,:)%arraysize)
535 ALLOCATE(tmpvolr(maxsize), tmpvold(maxsize), lin_mask(maxsize), weights(maxsize))
536 do_otimerange: DO j = 1, SIZE(that%timerange)
537  do_otime: DO i = 1, SIZE(that%time)
538  ninp = map_ttr(i,j)%arraysize
539  IF (ninp <= 0) cycle do_otime
540 ! required for some computations
541  CALL time_timerange_get_period(that%time(i), that%timerange(j), &
542  that%time_definition, pstart, pend, reftime)
543 
544  IF (ASSOCIATED(this%voldatir)) THEN
545  DO i1 = 1, SIZE(this%ana)
546  DO i3 = 1, SIZE(this%level)
547  DO i6 = 1, SIZE(this%network)
548  DO i5 = 1, SIZE(this%dativar%r)
549 ! stat_proc difference treated separately here
550  IF (stat_proc == 4) THEN
551  IF (ninp >= 2) THEN
552  IF (map_ttr(i,j)%array(1)%extra_info == 1 .AND. &
553  map_ttr(i,j)%array(n)%extra_info == 2) THEN
554  IF (c_e(this%voldatir(i1,map_ttr(i,j)%array(1)%it,i3, &
555  map_ttr(i,j)%array(1)%itr,i5,i6)) .AND. &
556  c_e(this%voldatir(i1,map_ttr(i,j)%array(ninp)%it,i3, &
557  map_ttr(i,j)%array(ninp)%itr,i5,i6))) THEN
558  that%voldatir(i1,i,i3,j,i5,i6) = &
559  this%voldatir(i1,map_ttr(i,j)%array(ninp)%it,i3, &
560  map_ttr(i,j)%array(ninp)%itr,i5,i6) - &
561  this%voldatir(i1,map_ttr(i,j)%array(1)%it,i3, &
562  map_ttr(i,j)%array(1)%itr,i5,i6)
563  ENDIF
564  ENDIF
565  ENDIF
566  cycle
567  ENDIF
568 ! other stat_proc
569  vartype = vol7d_vartype(this%dativar%r(i5))
570  lin_mask = .false.
571  ndtr = 0
572  DO n = 1, ninp
573  IF (c_e(this%voldatir(i1,map_ttr(i,j)%array(n)%it,i3, &
574  map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
575  ndtr = ndtr + 1
576  tmpvolr(ndtr) = this%voldatir(i1,map_ttr(i,j)%array(n)%it,i3, &
577  map_ttr(i,j)%array(n)%itr,i5,i6)
578  lin_mask(n) = .true.
579  ENDIF
580  ENDDO
581  IF (ndtr == 0) cycle
582  IF (lweighted) THEN
583  CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
584  pstart, pend, lin_mask(1:ninp), act_max_step, weights)
585  ELSE
586  CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
587  pstart, pend, lin_mask(1:ninp), act_max_step)
588  ENDIF
589  IF (act_max_step > lmax_step) cycle
590 
591  SELECT CASE(stat_proc)
592  CASE (0) ! average
593  IF (lweighted) THEN
594  that%voldatir(i1,i,i3,j,i5,i6) = &
595  sum(REAL(weights(1:ndtr))*tmpvolr(1:ndtr))
596  ELSE
597  that%voldatir(i1,i,i3,j,i5,i6) = &
598  sum(tmpvolr(1:ndtr))/ndtr
599  ENDIF
600  CASE (2) ! maximum
601  that%voldatir(i1,i,i3,j,i5,i6) = &
602  maxval(tmpvolr(1:ndtr))
603  CASE (3) ! minimum
604  that%voldatir(i1,i,i3,j,i5,i6) = &
605  minval(tmpvolr(1:ndtr))
606  CASE (6) ! stddev
607  that%voldatir(i1,i,i3,j,i5,i6) = &
608  stat_stddev(tmpvolr(1:ndtr))
609  CASE (201) ! mode
610 ! mode only for angles at the moment, with predefined histogram
611  IF (vartype == var_dir360) THEN
612 ! reduce to interval [-22.5,337.5]
613  WHERE (tmpvolr(1:ndtr) > 337.5)
614  tmpvolr(1:ndtr) = tmpvolr(1:ndtr) - 360.
615  END WHERE
616  that%voldatir(i1,i,i3,j,i5,i6) = &
617  stat_mode_histogram(tmpvolr(1:ndtr), &
618  8, -22.5, 337.5)
619  ENDIF
620  END SELECT
621  ENDDO
622  ENDDO
623  ENDDO
624  ENDDO
625  ENDIF
626 
627  IF (ASSOCIATED(this%voldatid)) THEN
628  DO i1 = 1, SIZE(this%ana)
629  DO i3 = 1, SIZE(this%level)
630  DO i6 = 1, SIZE(this%network)
631  DO i5 = 1, SIZE(this%dativar%d)
632 ! stat_proc difference treated separately here
633  IF (stat_proc == 4) THEN
634  IF (n >= 2) THEN
635  IF (map_ttr(i,j)%array(1)%extra_info == 1 .AND. &
636  map_ttr(i,j)%array(ninp)%extra_info == 2) THEN
637  IF (c_e(this%voldatid(i1,map_ttr(i,j)%array(1)%it,i3, &
638  map_ttr(i,j)%array(1)%itr,i5,i6)) .AND. &
639  c_e(this%voldatid(i1,map_ttr(i,j)%array(ninp)%it,i3, &
640  map_ttr(i,j)%array(ninp)%itr,i5,i6))) THEN
641  that%voldatid(i1,i,i3,j,i5,i6) = &
642  this%voldatid(i1,map_ttr(i,j)%array(ninp)%it,i3, &
643  map_ttr(i,j)%array(ninp)%itr,i5,i6) - &
644  this%voldatid(i1,map_ttr(i,j)%array(1)%it,i3, &
645  map_ttr(i,j)%array(1)%itr,i5,i6)
646  ENDIF
647  ENDIF
648  ENDIF
649  cycle
650  ENDIF
651 ! other stat_proc
652  vartype = vol7d_vartype(this%dativar%d(i5))
653  lin_mask = .false.
654  ndtr = 0
655  DO n = 1, ninp
656  IF (c_e(this%voldatid(i1,map_ttr(i,j)%array(n)%it,i3, &
657  map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
658  ndtr = ndtr + 1
659  tmpvold(ndtr) = this%voldatid(i1,map_ttr(i,j)%array(n)%it,i3, &
660  map_ttr(i,j)%array(n)%itr,i5,i6)
661  lin_mask(n) = .true.
662  ENDIF
663  ENDDO
664  IF (ndtr == 0) cycle
665  IF (lweighted) THEN
666  CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
667  pstart, pend, lin_mask(1:ninp), act_max_step, weights)
668  ELSE
669  CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
670  pstart, pend, lin_mask(1:ninp), act_max_step)
671  ENDIF
672  IF (act_max_step > lmax_step) cycle
673 
674  SELECT CASE(stat_proc)
675  CASE (0) ! average
676  IF (lweighted) THEN
677  that%voldatid(i1,i,i3,j,i5,i6) = &
678  sum(REAL(weights(1:ndtr))*tmpvold(1:ndtr))
679  ELSE
680  that%voldatid(i1,i,i3,j,i5,i6) = &
681  sum(tmpvold(1:ndtr))/ndtr
682  ENDIF
683  CASE (2) ! maximum
684  that%voldatid(i1,i,i3,j,i5,i6) = &
685  maxval(tmpvold(1:ndtr))
686  CASE (3) ! minimum
687  that%voldatid(i1,i,i3,j,i5,i6) = &
688  minval(tmpvold(1:ndtr))
689  CASE (6) ! stddev
690  that%voldatid(i1,i,i3,j,i5,i6) = &
691  stat_stddev(tmpvold(1:ndtr))
692  CASE (201) ! mode
693 ! mode only for angles at the moment, with predefined histogram
694  IF (vartype == var_dir360) THEN
695 ! reduce to interval [-22.5,337.5]
696  WHERE (tmpvold(1:ndtr) > 337.5)
697  tmpvold(1:ndtr) = tmpvold(1:ndtr) - 360.
698  END WHERE
699  that%voldatid(i1,i,i3,j,i5,i6) = &
700  stat_mode_histogram(tmpvold(1:ndtr), &
701  8, -22.5d0, 337.5d0)
702  ENDIF
703  END SELECT
704  ENDDO
705  ENDDO
706  ENDDO
707  ENDDO
708  ENDIF
709 
710  CALL delete(map_ttr(i,j))
711  ENDDO do_otime
712 ENDDO do_otimerange
713 
714 DEALLOCATE(map_ttr)
715 DEALLOCATE(tmpvolr, tmpvold, lin_mask, weights)
716 
717 IF (PRESENT(other)) THEN ! create volume with the remaining data for further processing
718  CALL vol7d_copy(this, other, miss=.false., sort=.false., unique=.false., &
719  ltimerange=(this%timerange(:)%timerange /= tri))
720 ENDIF
721 
722 END SUBROUTINE vol7d_compute_stat_proc_agg
723 
724 
740 SUBROUTINE vol7d_decompute_stat_proc(this, that, step, other, stat_proc_input)
741 TYPE(vol7d),INTENT(inout) :: this
742 TYPE(vol7d),INTENT(out) :: that
743 TYPE(timedelta),INTENT(in) :: step
744 TYPE(vol7d),INTENT(inout),OPTIONAL :: other
745 INTEGER,INTENT(in),OPTIONAL :: stat_proc_input
746 
747 INTEGER :: i, tri, steps
748 
749 
750 IF (PRESENT(stat_proc_input)) THEN
751  tri = stat_proc_input
752 ELSE
753  tri = 0
754 ENDIF
755 ! be safe
756 CALL vol7d_alloc_vol(this)
757 
758 ! compute length of cumulation step in seconds
759 CALL getval(step, asec=steps)
760 
761 ! filter requested data
762 CALL vol7d_copy(this, that, miss=.false., sort=.false., unique=.false., &
763  ltimerange=(this%timerange(:)%timerange == tri .AND. &
764  this%timerange(:)%p1 == 0 .AND. this%timerange(:)%p2 == steps))
765 
766 ! convert timerange to instantaneous and go back half step in time
767 that%timerange(:)%timerange = 254
768 that%timerange(:)%p2 = 0
769 DO i = 1, SIZE(that%time(:))
770  that%time(i) = that%time(i) - step/2
771 ENDDO
772 
773 IF (PRESENT(other)) THEN ! create volume with the remaining data for further processing
774  CALL vol7d_copy(this, other, miss=.false., sort=.false., unique=.false., &
775  ltimerange=(this%timerange(:)%timerange /= tri .OR. &
776  this%timerange(:)%p1 /= 0 .OR. this%timerange(:)%p2 /= steps))
777 ENDIF
778 
779 END SUBROUTINE vol7d_decompute_stat_proc
780 
781 
808 SUBROUTINE vol7d_recompute_stat_proc_diff(this, that, stat_proc, step, full_steps, other)
809 TYPE(vol7d),INTENT(inout) :: this
810 TYPE(vol7d),INTENT(out) :: that
811 INTEGER,INTENT(in) :: stat_proc
812 TYPE(timedelta),INTENT(in) :: step
813 LOGICAL,INTENT(in),OPTIONAL :: full_steps
814 TYPE(vol7d),INTENT(out),OPTIONAL :: other
815 
816 INTEGER :: i1, i3, i5, i6, i, j, k, l, nitr, steps
817 INTEGER,POINTER :: map_tr(:,:,:,:,:), f(:)
818 LOGICAL,POINTER :: mask_timerange(:)
819 LOGICAL,ALLOCATABLE :: mask_time(:)
820 TYPE(vol7d) :: v7dtmp
821 
822 
823 ! be safe
824 CALL vol7d_alloc_vol(this)
825 ! initialise the template of the output volume
826 CALL init(that, time_definition=this%time_definition)
827 
828 ! compute length of cumulation step in seconds
829 CALL getval(step, asec=steps)
830 
831 ! compute the statistical processing relations, output time and
832 ! timerange are defined here
833 CALL recompute_stat_proc_diff_common(this%time, this%timerange, stat_proc, step, &
834  nitr, that%time, that%timerange, map_tr, f, mask_timerange, &
835  this%time_definition, full_steps)
836 
837 ! complete the definition of the empty output template
838 CALL vol7d_alloc(that, nana=0, nlevel=0, nnetwork=0)
839 CALL vol7d_alloc_vol(that)
840 ! shouldn't we exit here with an empty volume if stat_proc/=0,1 ?
841 ALLOCATE(mask_time(SIZE(this%time)))
842 DO l = 1, SIZE(this%time)
843  mask_time(l) = any(this%time(l) == that%time(:))
844 ENDDO
845 ! create template for the output volume, keep all ana, level, network
846 ! and variables; copy only the timeranges already satisfying the
847 ! requested step, if any and only the times already existing in the
848 ! output
849 CALL vol7d_copy(this, v7dtmp, miss=.false., sort=.false., unique=.false., &
850  ltimerange=mask_timerange(:), ltime=mask_time(:))
851 ! merge output created so far with template
852 CALL vol7d_merge(that, v7dtmp, lanasimple=.true., llevelsimple=.true.)
853 DEALLOCATE(mask_time)
854 
855 ! compute statistical processing
856 IF (ASSOCIATED(this%voldatir)) THEN
857  DO l = 1, SIZE(this%time)
858  DO k = 1, nitr
859  DO j = 1, SIZE(this%time)
860  DO i = 1, nitr
861  IF (c_e(map_tr(i,j,k,l,1))) THEN
862  DO i6 = 1, SIZE(this%network)
863  DO i5 = 1, SIZE(this%dativar%r)
864  DO i3 = 1, SIZE(this%level)
865  DO i1 = 1, SIZE(this%ana)
866  IF (c_e(this%voldatir(i1,l,i3,f(k),i5,i6)) .AND. &
867  c_e(this%voldatir(i1,j,i3,f(i),i5,i6))) THEN
868 
869  IF (stat_proc == 0) THEN ! average
870  that%voldatir( &
871  i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
872  (this%voldatir(i1,l,i3,f(k),i5,i6)*this%timerange(f(k))%p2 - &
873  this%voldatir(i1,j,i3,f(i),i5,i6)*this%timerange(f(i))%p2)/ &
874  steps ! optimize avoiding conversions
875  ELSE IF (stat_proc == 1 .OR. stat_proc == 4) THEN ! acc, diff
876  that%voldatir( &
877  i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
878  this%voldatir(i1,l,i3,f(k),i5,i6) - &
879  this%voldatir(i1,j,i3,f(i),i5,i6)
880  ENDIF
881 
882  ENDIF
883  ENDDO
884  ENDDO
885  ENDDO
886  ENDDO
887  ENDIF
888  ENDDO
889  ENDDO
890  ENDDO
891  ENDDO
892 ENDIF
893 
894 IF (ASSOCIATED(this%voldatid)) THEN
895  DO l = 1, SIZE(this%time)
896  DO k = 1, nitr
897  DO j = 1, SIZE(this%time)
898  DO i = 1, nitr
899  IF (c_e(map_tr(i,j,k,l,1))) THEN
900  DO i6 = 1, SIZE(this%network)
901  DO i5 = 1, SIZE(this%dativar%d)
902  DO i3 = 1, SIZE(this%level)
903  DO i1 = 1, SIZE(this%ana)
904  IF (c_e(this%voldatid(i1,l,i3,f(k),i5,i6)) .AND. &
905  c_e(this%voldatid(i1,j,i3,f(i),i5,i6))) THEN
906 ! IF (.NOT.c_e(that%voldatid( &
907 ! i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6))) THEN
908 
909  IF (stat_proc == 0) THEN ! average
910  that%voldatid( &
911  i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
912  (this%voldatid(i1,l,i3,f(k),i5,i6)*this%timerange(f(k))%p2 - &
913  this%voldatid(i1,j,i3,f(i),i5,i6)*this%timerange(f(i))%p2)/ &
914  steps ! optimize avoiding conversions
915  ELSE IF (stat_proc == 1 .OR. stat_proc == 4) THEN ! acc, diff
916  that%voldatid( &
917  i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
918  this%voldatid(i1,l,i3,f(k),i5,i6) - &
919  this%voldatid(i1,j,i3,f(i),i5,i6)
920  ENDIF
921 
922 ! ENDIF
923  ENDIF
924  ENDDO
925  ENDDO
926  ENDDO
927  ENDDO
928  ENDIF
929  ENDDO
930  ENDDO
931  ENDDO
932  ENDDO
933 ENDIF
934 
935 ! this should be avoided by sorting descriptors upstream
936 ! descriptors now are sorted upstream with a dirty and expensive trick
937 ! but the order may be scrambled in the call to vol7d_merge above
938 CALL vol7d_smart_sort(that, lsort_time=.true., lsort_timerange=.true.)
939 
940 DEALLOCATE(map_tr, f, mask_timerange)
941 CALL makeother(.true.)
942 
943 CONTAINS
944 
945 SUBROUTINE makeother(filter)
946 LOGICAL,INTENT(in) :: filter
947 IF (PRESENT(other)) THEN
948  IF (filter) THEN ! create volume with the remaining data for further processing
949  CALL vol7d_copy(this, other, miss=.false., sort=.false., unique=.false., &
950  ltimerange=(this%timerange(:)%timerange /= stat_proc))
951  ELSE
952  CALL vol7d_copy(this, other, miss=.false., sort=.false., unique=.false.)
953  ENDIF
954 ENDIF
955 END SUBROUTINE makeother
956 
957 END SUBROUTINE vol7d_recompute_stat_proc_diff
958 
959 
987 SUBROUTINE vol7d_compute_stat_proc_metamorph(this, that, stat_proc_input, stat_proc)
988 TYPE(vol7d),INTENT(inout) :: this
989 TYPE(vol7d),INTENT(out) :: that
990 INTEGER,INTENT(in) :: stat_proc_input
991 INTEGER,INTENT(in) :: stat_proc
992 
993 INTEGER :: j
994 LOGICAL,ALLOCATABLE :: tr_mask(:)
995 REAL,ALLOCATABLE :: int_ratio(:)
996 DOUBLE PRECISION,ALLOCATABLE :: int_ratiod(:)
997 
998 IF (.NOT.((stat_proc_input == 0 .AND. stat_proc == 1) .OR. &
999  (stat_proc_input == 1 .AND. stat_proc == 0))) THEN
1000 
1001  CALL l4f_log(l4f_warn, &
1002  'compute_stat_proc_metamorph, can only be applied to average->accumulated timerange and viceversa')
1003 ! return an empty volume, without signaling error
1004  CALL init(that)
1005  CALL vol7d_alloc_vol(that)
1006  RETURN
1007 ENDIF
1008 
1009 ! be safe
1010 CALL vol7d_alloc_vol(this)
1011 
1012 ! useful timeranges
1013 tr_mask = this%timerange(:)%timerange == stat_proc_input .AND. this%timerange(:)%p2 /= imiss &
1014  .AND. this%timerange(:)%p2 /= 0
1015 
1016 ! initial check (necessary?)
1017 IF (count(tr_mask) == 0) THEN
1018  CALL l4f_log(l4f_warn, &
1019  'vol7d_compute, no timeranges suitable for statistical processing by metamorphosis')
1020  CALL init(that)
1021 ! CALL makeother()
1022  RETURN
1023 ENDIF
1024 
1025 ! copy required data and reset timerange
1026 CALL vol7d_copy(this, that, ltimerange=tr_mask)
1027 that%timerange(:)%timerange = stat_proc
1028 ! why next automatic f2003 allocation does not always work?
1029 ALLOCATE(int_ratio(SIZE(that%timerange)), int_ratiod(SIZE(that%timerange)))
1030 
1031 IF (stat_proc == 0) THEN ! average -> integral
1032  int_ratio = 1./REAL(that%timerange(:)%p2)
1033  int_ratiod = 1./dble(that%timerange(:)%p2)
1034 ELSE ! cumulation
1035  int_ratio = REAL(that%timerange(:)%p2)
1036  int_ratiod = dble(that%timerange(:)%p2)
1037 ENDIF
1038 
1039 IF (ASSOCIATED(that%voldatir)) THEN
1040  DO j = 1, SIZE(that%timerange)
1041  WHERE(c_e(that%voldatir(:,:,:,j,:,:)))
1042  that%voldatir(:,:,:,j,:,:) = that%voldatir(:,:,:,j,:,:)*int_ratio(j)
1043  ELSEWHERE
1044  that%voldatir(:,:,:,j,:,:) = rmiss
1045  END WHERE
1046  ENDDO
1047 ENDIF
1048 
1049 IF (ASSOCIATED(that%voldatid)) THEN
1050  DO j = 1, SIZE(that%timerange)
1051  WHERE(c_e(that%voldatid(:,:,:,j,:,:)))
1052  that%voldatid(:,:,:,j,:,:) = that%voldatid(:,:,:,j,:,:)*int_ratiod(j)
1053  ELSEWHERE
1054  that%voldatid(:,:,:,j,:,:) = rmiss
1055  END WHERE
1056  ENDDO
1057 ENDIF
1058 
1059 
1060 END SUBROUTINE vol7d_compute_stat_proc_metamorph
1061 
1062 
1079 SUBROUTINE vol7d_fill_time(this, that, step, start, stopp, cyclicdt)
1080 TYPE(vol7d),INTENT(inout) :: this
1081 TYPE(vol7d),INTENT(inout) :: that
1082 TYPE(timedelta),INTENT(in) :: step
1083 TYPE(datetime),INTENT(in),OPTIONAL :: start
1084 TYPE(datetime),INTENT(in),OPTIONAL :: stopp
1085 TYPE(cyclicdatetime),INTENT(in),OPTIONAL :: cyclicdt
1086 
1087 TYPE(cyclicdatetime) :: lcyclicdt
1088 TYPE(datetime) :: counter, lstart, lstop
1089 INTEGER :: i, naddtime
1090 
1091 CALL safe_start_stop(this, lstart, lstop, start, stopp)
1092 IF (.NOT. c_e(lstart) .OR. .NOT. c_e(lstop) .OR. .NOT. c_e(step)) RETURN
1093 
1094 lcyclicdt=cyclicdatetime_miss
1095 if (present(cyclicdt)) then
1096  if(c_e(cyclicdt)) lcyclicdt=cyclicdt
1097 end if
1098 
1099 CALL l4f_log(l4f_info, 'vol7d_fill_time: time interval '//trim(to_char(lstart))// &
1100  ' '//trim(to_char(lstop)))
1101 
1102 ! Count the number of time levels required for completing the series
1103 ! valid also in the case (SIZE(this%time) == 0)
1104 naddtime = 0
1105 counter = lstart
1106 i = 1
1107 naddcount: DO WHILE(counter <= lstop)
1108  DO WHILE(i <= SIZE(this%time)) ! this%time(i) chases counter
1109  IF (counter < this%time(i)) THEN ! this%time(i) overtook counter
1110  i = max(i-1,1) ! go back if possible
1111  EXIT
1112  ELSE IF (counter == this%time(i) .OR. .NOT. counter == lcyclicdt) THEN ! found matching time
1113  counter = counter + step
1114  cycle naddcount
1115  ENDIF
1116  i = i + 1
1117  ENDDO
1118  naddtime = naddtime + 1
1119  counter = counter + step
1120 ENDDO naddcount
1121 
1122 ! old universal algorithm, not optimized, check that the new one is equivalent
1123 !naddtime = 0
1124 !counter = lstart
1125 !DO WHILE(counter <= lstop)
1126 ! IF (.NOT.ANY(counter == this%time(:))) THEN
1127 ! naddtime = naddtime + 1
1128 ! ENDIF
1129 ! counter = counter + step
1130 !ENDDO
1131 
1132 IF (naddtime > 0) THEN
1133 
1134  CALL init(that)
1135  CALL vol7d_alloc(that, ntime=naddtime)
1136  CALL vol7d_alloc_vol(that)
1137 
1138  ! Repeat the count loop setting the time levels to be added
1139  naddtime = 0
1140  counter = lstart
1141  i = 1
1142  naddadd: DO WHILE(counter <= lstop)
1143  DO WHILE(i <= SIZE(this%time)) ! this%time(i) chases counter
1144  IF (counter < this%time(i)) THEN ! this%time(i) overtook counter
1145  i = max(i-1,1) ! go back if possible
1146  EXIT
1147  ELSE IF (counter == this%time(i) .OR. .NOT. counter == lcyclicdt) THEN ! found matching time
1148  counter = counter + step
1149  cycle naddadd
1150  ENDIF
1151  i = i + 1
1152  ENDDO
1153  naddtime = naddtime + 1
1154  that%time(naddtime) = counter ! only difference
1155  counter = counter + step
1156  ENDDO naddadd
1157 
1158  CALL vol7d_append(that, this, sort=.true.)
1159 
1160 ELSE
1161 !! ? why sort all dimension ?
1162 !! CALL vol7d_copy(this, that, lsort_time=.TRUE.)
1163  CALL vol7d_copy(this, that, sort=.true.)
1164 ENDIF
1165 
1166 
1167 END SUBROUTINE vol7d_fill_time
1168 
1169 
1181 SUBROUTINE vol7d_filter_time(this, that, step, start, stopp, cyclicdt)
1182 TYPE(vol7d),INTENT(inout) :: this
1183 TYPE(vol7d),INTENT(inout) :: that
1184 TYPE(timedelta),INTENT(in),optional :: step
1185 TYPE(datetime),INTENT(in),OPTIONAL :: start
1186 TYPE(datetime),INTENT(in),OPTIONAL :: stopp
1187 TYPE(cyclicdatetime),INTENT(in),OPTIONAL :: cyclicdt
1188 
1189 TYPE(datetime) :: lstart, lstop
1190 LOGICAL, ALLOCATABLE :: time_mask(:)
1191 
1192 CALL safe_start_stop(this, lstart, lstop, start, stopp)
1193 IF (.NOT. c_e(lstart) .OR. .NOT. c_e(lstop)) RETURN
1194 
1195 CALL l4f_log(l4f_info, 'vol7d_filter_time: time interval '//trim(to_char(lstart))// &
1196  ' '//trim(to_char(lstop)))
1197 
1198 ALLOCATE(time_mask(SIZE(this%time)))
1199 
1200 time_mask = this%time >= lstart .AND. this%time <= lstop
1202 IF (PRESENT(cyclicdt)) THEN
1203  IF (c_e(cyclicdt)) THEN
1204  time_mask = time_mask .AND. this%time == cyclicdt
1205  ENDIF
1206 ENDIF
1207 
1208 IF (PRESENT(step)) THEN
1209  IF (c_e(step)) THEN
1210  time_mask = time_mask .AND. mod(this%time - lstart, step) == timedelta_0
1211  ENDIF
1212 ENDIF
1213 
1214 CALL vol7d_copy(this,that, ltime=time_mask)
1215 
1216 DEALLOCATE(time_mask)
1217 
1218 END SUBROUTINE vol7d_filter_time
1219 
1220 
1224 SUBROUTINE vol7d_fill_data(this, step, start, stopp, tolerance)
1225 TYPE(vol7d),INTENT(inout) :: this
1226 TYPE(timedelta),INTENT(in) :: step
1227 TYPE(datetime),INTENT(in),OPTIONAL :: start
1228 TYPE(datetime),INTENT(in),OPTIONAL :: stopp
1229 TYPE(timedelta),INTENT(in),optional :: tolerance
1230 
1231 TYPE(datetime) :: lstart, lstop
1232 integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork, iindtime
1233 type(timedelta) :: deltato,deltat, ltolerance
1234 
1235 CALL safe_start_stop(this, lstart, lstop, start, stopp)
1236 IF (.NOT. c_e(lstart) .OR. .NOT. c_e(lstop)) RETURN
1237 
1238 CALL l4f_log(l4f_info, 'vol7d_fill_data: time interval '//trim(to_char(lstart))// &
1239  ' '//trim(to_char(lstop)))
1240 
1241 
1242 ltolerance=step/2
1243 
1244 if (present(tolerance))then
1245  if (c_e(tolerance)) ltolerance=tolerance
1246 end if
1247 
1248 
1249 do indtime=1,size(this%time)
1250 
1251  IF (this%time(indtime) < lstart .OR. this%time(indtime) > lstop .OR. &
1252  mod(this%time(indtime) - lstart, step) /= timedelta_0) cycle
1253  do indtimerange=1,size(this%timerange)
1254  if (this%timerange(indtimerange)%timerange /= 254) cycle
1255  do indnetwork=1,size(this%network)
1256  do inddativarr=1,size(this%dativar%r)
1257  do indlevel=1,size(this%level)
1258  do indana=1,size(this%ana)
1259 
1260  !find the nearest data in time if data is missing
1261  if (.not. c_e(this%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)))then
1262  deltato=timedelta_miss
1263 
1264  !do iindtime=max(indtime-20,1),min(indtime+20,size(this%time)) !check on a chunk: 20 should be enought
1265 
1266  do iindtime=indtime+1,size(this%time) !check forward
1267 
1268  if (c_e(this%voldatir (indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork )))then
1269  deltat=this%time(iindtime)-this%time(indtime)
1270 
1271  if (deltat >= ltolerance) exit
1272 
1273  if (deltat < deltato) then
1274  this%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork) = &
1275  this%voldatir(indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork)
1276  deltato=deltat
1277  end if
1278  end if
1279  end do
1280 
1281  do iindtime=indtime-1,1,-1 !check backward
1282 
1283  if (c_e(this%voldatir (indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork )))then
1284  if (iindtime < indtime) then
1285  deltat=this%time(indtime)-this%time(iindtime)
1286  else if (iindtime > indtime) then
1287  deltat=this%time(iindtime)-this%time(indtime)
1288  else
1289  cycle
1290  end if
1291 
1292  if (deltat >= ltolerance) exit
1294  if (deltat < deltato) then
1295  this%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork) = &
1296  this%voldatir(indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork)
1297  deltato=deltat
1298  end if
1299  end if
1300  end do
1301 
1302  end if
1303  end do
1304  end do
1305  end do
1306  end do
1307  end do
1308 end do
1309 
1310 END SUBROUTINE vol7d_fill_data
1311 
1312 
1313 ! private utility routine for checking interval and start-stop times
1314 ! in input missing start-stop values are treated as not present
1315 ! in output missing start-stop values mean "do nothing"
1316 SUBROUTINE safe_start_stop(this, lstart, lstop, start, stopp)
1317 TYPE(vol7d),INTENT(inout) :: this
1318 TYPE(datetime),INTENT(out) :: lstart
1319 TYPE(datetime),INTENT(out) :: lstop
1320 TYPE(datetime),INTENT(in),OPTIONAL :: start
1321 TYPE(datetime),INTENT(in),OPTIONAL :: stopp
1322 
1323 lstart = datetime_miss
1324 lstop = datetime_miss
1325 ! initial safety operation
1326 CALL vol7d_alloc_vol(this)
1327 IF (SIZE(this%time) == 0) RETURN ! avoid segmentation fault in case of empty volume
1328 CALL vol7d_smart_sort(this, lsort_time=.true.)
1329 
1330 IF (PRESENT(start)) THEN
1331  IF (c_e(start)) THEN
1332  lstart = start
1333  ELSE
1334  lstart = this%time(1)
1335  ENDIF
1336 ELSE
1337  lstart = this%time(1)
1338 ENDIF
1339 IF (PRESENT(stopp)) THEN
1340  IF (c_e(stopp)) THEN
1341  lstop = stopp
1342  ELSE
1343  lstop = this%time(SIZE(this%time))
1344  ENDIF
1345 ELSE
1346  lstop = this%time(SIZE(this%time))
1347 ENDIF
1348 
1349 END SUBROUTINE safe_start_stop
1350 
1351 
1358 SUBROUTINE vol7d_normalize_vcoord(this,that,ana,time,timerange,network)
1359 TYPE(vol7d),INTENT(INOUT) :: this
1360 TYPE(vol7d),INTENT(OUT) :: that
1361 integer,intent(in) :: time,ana,timerange,network
1362 
1363 character(len=1) :: type
1364 integer :: ind
1365 TYPE(vol7d_var) :: var
1366 LOGICAL,allocatable :: ltime(:),ltimerange(:),lana(:),lnetwork(:)
1367 logical,allocatable :: maschera(:)
1368 
1369 
1370 allocate(ltime(size(this%time)))
1371 allocate(ltimerange(size(this%timerange)))
1372 allocate(lana(size(this%ana)))
1373 allocate(lnetwork(size(this%network)))
1374 
1375 ltime=.false.
1376 ltimerange=.false.
1377 lana=.false.
1378 lnetwork=.false.
1379 
1380 ltime(time)=.true.
1381 ltimerange(timerange)=.true.
1382 lana(ana)=.true.
1383 lnetwork(network)=.true.
1384 
1385 call vol7d_copy(this, that,unique=.true.,&
1386  ltime=ltime,ltimerange=ltimerange,lana=lana,lnetwork=lnetwork )
1387 
1388 call init(var, btable="B10004") ! Pressure
1389 type=cmiss
1390 !type="i"
1391 ind = index(that%dativar, var, type=type)
1392 
1393 allocate(maschera(size(that%level)))
1394 
1395 maschera = (&
1396  (that%level%level1 == 105.and.that%level%level2 == 105) .or. &
1397  (that%level%level1 == 103 .and. that%level%level2 == imiss ) .or. &
1398  (that%level%level1 == 102 .and. that%level%level2 == imiss )) &
1399  .and. c_e(that%voldatic(1,1,:,1,ind,1))
1400 
1401 
1402 select case (type)
1403 
1404 case("d")
1405 
1406  where (maschera)
1407  that%level%level1 = 100
1408  that%level%l1 = int(realdat(that%voldatid(1,1,:,1,ind,1),that%dativar%d(ind)))
1409  that%level%l1 = int(that%voldatid(1,1,:,1,ind,1))
1410  that%level%level2 = imiss
1411  that%level%l2 = imiss
1412  end where
1413 
1414 case("r")
1415 
1416  where (maschera)
1417  that%level%level1 = 100
1418  that%level%l1 = int(realdat(that%voldatir(1,1,:,1,ind,1),that%dativar%r(ind)))
1419  that%level%level2 = imiss
1420  that%level%l2 = imiss
1421  end where
1422 
1423 case("i")
1424 
1425  where (maschera)
1426  that%level%level1 = 100
1427  that%level%l1 = int(realdat(that%voldatii(1,1,:,1,ind,1),that%dativar%i(ind)))
1428  that%level%level2 = imiss
1429  that%level%l2 = imiss
1430  end where
1431 
1432 case("b")
1433 
1434  where (maschera)
1435  that%level%level1 = 100
1436  that%level%l1 = int(realdat(that%voldatib(1,1,:,1,ind,1),that%dativar%b(ind)))
1437  that%level%level2 = imiss
1438  that%level%l2 = imiss
1439  end where
1440 
1441 case("c")
1442 
1443  where (maschera)
1444  that%level%level1 = 100
1445  that%level%l1 = int(realdat(that%voldatic(1,1,:,1,ind,1),that%dativar%c(ind)))
1446  that%level%level2 = imiss
1447  that%level%l2 = imiss
1448  end where
1449 
1450 end select
1451 
1452 deallocate(ltime)
1453 deallocate(ltimerange)
1454 deallocate(lana)
1455 deallocate(lnetwork)
1456 
1457 END SUBROUTINE vol7d_normalize_vcoord
1458 
1459 
1460 !!$!> Metodo per calcolare variabili derivate.
1461 !!$!! TO DO !!
1462 !!$SUBROUTINE vol7d_compute_var(this,that,var)
1463 !!$TYPE(vol7d),INTENT(INOUT) :: this !< oggetto da normalizzare
1464 !!$TYPE(vol7d),INTENT(OUT) :: that !< oggetto normalizzato
1465 !!$
1466 !!$character(len=1) :: type
1467 !!$TYPE(vol7d_var),intent(in) :: var
1468 
1469 
1470 !!$call init(var, btable="B10004") ! Pressure
1471 !!$type=cmiss
1472 !!$call vol7d_varvect_index(that%dativar,var , type=type,index_v=ind)
1473 !!$
1474 !!$select case (type)
1475 !!$
1476 !!$case("d")
1477 !!$
1478 !!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
1479 !!$ that%level%level1 = 100
1480 !!$ that%level%l1 = realdat(that%voldatid(1,1,:,1,ind,1),that%dativar%d(ind))
1481 !!$ that%level%level2 = imiss
1482 !!$ that%level%l2 = imiss
1483 !!$ end where
1484 !!$
1485 !!$case("r")
1486 !!$
1487 !!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
1488 !!$ that%level%level1 = 100
1489 !!$ that%level%l1 = realdat(that%voldatir(1,1,:,1,ind,1),that%dativar%r(ind))
1490 !!$ that%level%level2 = imiss
1491 !!$ that%level%l2 = imiss
1492 !!$ end where
1493 !!$
1494 !!$case("i")
1495 !!$
1496 !!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
1497 !!$ that%level%level1 = 100
1498 !!$ that%level%l1 = realdat(that%voldatii(1,1,:,1,ind,1),that%dativar%i(ind))
1499 !!$ that%level%level2 = imiss
1500 !!$ that%level%l2 = imiss
1501 !!$ end where
1502 !!$
1503 !!$case("b")
1504 !!$
1505 !!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
1506 !!$ that%level%level1 = 100
1507 !!$ that%level%l1 = realdat(that%voldatib(1,1,:,1,ind,1),that%dativar%b(ind))
1508 !!$ that%level%level2 = imiss
1509 !!$ that%level%l2 = imiss
1510 !!$ end where
1511 !!$
1512 !!$case("c")
1513 !!$
1514 !!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
1515 !!$ that%level%level1 = 100
1516 !!$ that%level%l1 = realdat(that%voldatic(1,1,:,1,ind,1),that%dativar%c(ind))
1517 !!$ that%level%level2 = imiss
1518 !!$ that%level%l2 = imiss
1519 !!$ end where
1520 !!$
1521 !!$end select
1522 
1523 !!$
1524 !!$END SUBROUTINE vol7d_compute_var
1525 !!$
1526 
1527 
1528 
1529 END MODULE vol7d_class_compute
Functions that return a trimmed CHARACTER representation of the input variable.
Compute the mode of the random variable provided taking into account missing data.
This module contains functions that are only for internal use of the library.
Classi per la gestione delle coordinate temporali.
Compute the standard deviation of the random variable provided, taking into account missing data...
Definition: simple_stat.f90:69
Distruttori per le 2 classi.
Index method.
Restituiscono il valore dell&#39;oggetto in forma di stringa stampabile.
Restituiscono il valore dell&#39;oggetto nella forma desiderata.
Operatore di resto della divisione.
Classe per la gestione di un volume completo di dati osservati.
Costruttori per le classi datetime e timedelta.
Extension of vol7d_class with methods for performing simple statistical operations on entire volumes ...
Module for basic statistical computations taking into account missing data.
Definition: simple_stat.f90:25
real data conversion

Generated with Doxygen.