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