libsim Versione 7.1.11
volgrid6d_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
33IMPLICIT NONE
34
35CONTAINS
36
102SUBROUTINE volgrid6d_compute_stat_proc(this, that, stat_proc_input, stat_proc, &
103 step, start, full_steps, frac_valid, max_step, weighted, clone)
104TYPE(volgrid6d),INTENT(inout) :: this
105TYPE(volgrid6d),INTENT(out) :: that
106INTEGER,INTENT(in) :: stat_proc_input
107INTEGER,INTENT(in) :: stat_proc
108TYPE(timedelta),INTENT(in) :: step
109TYPE(datetime),INTENT(in),OPTIONAL :: start
110LOGICAL,INTENT(in),OPTIONAL :: full_steps
111REAL,INTENT(in),OPTIONAL :: frac_valid
112TYPE(timedelta),INTENT(in),OPTIONAL :: max_step ! maximum allowed distance in time between two single valid data within a dataset, for the dataset to be eligible for statistical processing
113LOGICAL,INTENT(in),OPTIONAL :: weighted
114LOGICAL , INTENT(in),OPTIONAL :: clone
115
116INTEGER :: dtmax, dtstep
117
118
119IF (stat_proc_input == 254) THEN
120 CALL l4f_category_log(this%category, l4f_info, &
121 'computing statistical processing by aggregation '//&
122 trim(to_char(stat_proc_input))//':'//trim(to_char(stat_proc)))
123
124 CALL volgrid6d_compute_stat_proc_agg(this, that, stat_proc, &
125 step, start, full_steps, max_step, clone)
126
127ELSE IF (stat_proc == 254) THEN
128 CALL l4f_category_log(this%category, l4f_error, &
129 'statistical processing to instantaneous data not implemented for gridded fields')
130 CALL raise_error()
131
132ELSE IF (stat_proc_input /= stat_proc) THEN
133 IF ((stat_proc_input == 0 .AND. stat_proc == 1) .OR. &
134 (stat_proc_input == 1 .AND. stat_proc == 0)) THEN
135 CALL l4f_category_log(this%category, l4f_info, &
136 'computing statistically processed data by integration/differentiation '// &
137 t2c(stat_proc_input)//':'//t2c(stat_proc))
138 CALL volgrid6d_compute_stat_proc_metamorph(this, that, stat_proc_input, &
139 stat_proc, clone)
140 ELSE
141 CALL l4f_category_log(this%category, l4f_error, &
142 'statistical processing '//t2c(stat_proc_input)//':'//t2c(stat_proc)// &
143 ' not implemented or does not make sense')
144 CALL raise_error()
145 ENDIF
146
147ELSE IF (count(this%timerange(:)%timerange == stat_proc) == 0) THEN
148 CALL l4f_category_log(this%category, l4f_warn, &
149 'no timeranges of the desired statistical processing type '//t2c(stat_proc)//' available')
150! return an empty volume, without signaling error
151 CALL init(that)
152 CALL volgrid6d_alloc_vol(that)
153
154ELSE
155! euristically determine whether aggregation or difference is more suitable
156 dtmax = maxval(this%timerange(:)%p2, &
157 mask=(this%timerange(:)%timerange == stat_proc))
158 CALL getval(step, asec=dtstep)
159
160#ifdef DEBUG
161 CALL l4f_category_log(this%category, l4f_debug, &
162 'stat_proc='//t2c(stat_proc)//' dtmax='//t2c(dtmax)//' dtstep='//t2c(dtstep))
163#endif
164
165 IF (dtstep < dtmax) THEN
166 CALL l4f_category_log(this%category, l4f_info, &
167 'recomputing statistically processed data by difference '// &
168 t2c(stat_proc_input)//':'//t2c(stat_proc))
169 CALL volgrid6d_recompute_stat_proc_diff(this, that, stat_proc, step, &
170 full_steps, start, clone)
171 ELSE
172 CALL l4f_category_log(this%category, l4f_info, &
173 'recomputing statistically processed data by aggregation '// &
174 t2c(stat_proc_input)//':'//t2c(stat_proc))
175 CALL volgrid6d_recompute_stat_proc_agg(this, that, stat_proc, step, start, &
176 full_steps, frac_valid, clone)
177 ENDIF
178
179ENDIF
180
181END SUBROUTINE volgrid6d_compute_stat_proc
182
183
226SUBROUTINE volgrid6d_recompute_stat_proc_agg(this, that, stat_proc, &
227 step, start, full_steps, frac_valid, clone)
228TYPE(volgrid6d),INTENT(inout) :: this
229TYPE(volgrid6d),INTENT(out) :: that
230INTEGER,INTENT(in) :: stat_proc
231TYPE(timedelta),INTENT(in) :: step
232TYPE(datetime),INTENT(in),OPTIONAL :: start
233LOGICAL,INTENT(in),OPTIONAL :: full_steps
234REAL,INTENT(in),OPTIONAL :: frac_valid
235LOGICAL, INTENT(in),OPTIONAL :: clone
236
237INTEGER :: tri
238INTEGER i, j, n, n1, ndtr, i3, i6
239TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
240INTEGER,POINTER :: dtratio(:)
241REAL :: lfrac_valid
242LOGICAL :: lclone
243REAL,POINTER :: voldatiin(:,:), voldatiout(:,:)
244
245
246NULLIFY(voldatiin, voldatiout)
247tri = stat_proc
248IF (PRESENT(frac_valid)) THEN
249 lfrac_valid = frac_valid
250ELSE
251 lfrac_valid = 1.0
252ENDIF
253
254CALL init(that)
255! be safe
256CALL volgrid6d_alloc_vol(this)
257
258! when volume is not decoded it is better to clone anyway to avoid
259! overwriting fields
260lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati)
261! initialise the output volume
262CALL init(that, griddim=this%griddim, time_definition=this%time_definition)
263CALL volgrid6d_alloc(that, dim=this%griddim%dim, ntimerange=1, &
264 nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.)
265that%level = this%level
266that%var = this%var
267
268CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
269 step, this%time_definition, that%time, that%timerange, map_ttr, &
270 dtratio=dtratio, start=start, full_steps=full_steps)
271
272CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati))
273
274do_otimerange: DO j = 1, SIZE(that%timerange)
275 do_otime: DO i = 1, SIZE(that%time)
276
277 DO n1 = 1, SIZE(dtratio)
278 IF (dtratio(n1) <= 0) cycle ! safety check
279
280 DO i6 = 1, SIZE(this%var)
281 DO i3 = 1, SIZE(this%level)
282 CALL volgrid_get_vol_2d(that, i3, i, j, i6, voldatiout)
283 ndtr = 0
284 DO n = 1, map_ttr(i,j)%arraysize
285 IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
286 ndtr = ndtr + 1
287 CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(n)%it, &
288 map_ttr(i,j)%array(n)%itr, i6, voldatiin)
289
290 IF (ndtr == 1) THEN
291 voldatiout = voldatiin
292 IF (lclone) THEN
293 CALL copy(this%gaid(i3, map_ttr(i,j)%array(n)%it,&
294 map_ttr(i,j)%array(n)%itr,i6), that%gaid(i3,i,j,i6))
295 ELSE
296 that%gaid(i3,i,j,i6) = this%gaid(i3, map_ttr(i,j)%array(n)%it, &
297 map_ttr(i,j)%array(n)%itr,i6)
298 ENDIF
299
300 ELSE ! second or more time
301 SELECT CASE(stat_proc)
302 CASE (0, 200, 1, 4) ! average, vectorial mean, accumulation, difference
303 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
304 voldatiout(:,:) = voldatiout(:,:) + voldatiin(:,:)
305 ELSEWHERE
306 voldatiout(:,:) = rmiss
307 END WHERE
308 CASE(2) ! maximum
309 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
310 voldatiout(:,:) = max(voldatiout(:,:), voldatiin(:,:))
311 ELSEWHERE
312 voldatiout(:,:) = rmiss
313 END WHERE
314 CASE(3) ! minimum
315 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
316 voldatiout(:,:) = min(voldatiout(:,:), voldatiin(:,:))
317 ELSEWHERE
318 voldatiout(:,:) = rmiss
319 END WHERE
320 END SELECT
321
322 ENDIF ! first time
323 ENDIF ! dtratio(n1)
324 ENDDO ! ttr
325
326#ifdef DEBUG
327 CALL l4f_log(l4f_debug, &
328 'compute_stat_proc_agg, ndtr/dtratio/frac_valid: '// &
329 t2c(ndtr)//'/'//t2c(dtratio(n1))//'/'//t2c(lfrac_valid))
330#endif
331 IF (ndtr > 0) THEN ! why this condition was not here before?
332 IF (real(ndtr)/real(dtratio(n1)) >= lfrac_valid) THEN ! success
333 IF (stat_proc == 0) THEN ! average
334 WHERE(c_e(voldatiout(:,:)))
335 voldatiout(:,:) = voldatiout(:,:)/ndtr
336 END WHERE
337 ENDIF
338 CALL volgrid_set_vol_2d(that, i3, i, j, i6, voldatiout)
339#ifdef DEBUG
340 CALL l4f_log(l4f_debug, &
341 'compute_stat_proc_agg, coding lev/t/tr/var: '// &
342 t2c(i3)//'/'//t2c(i)//'/'//t2c(j)//'/'//t2c(i6))
343#endif
344 ELSE
345! must nullify the output gaid here, otherwise an incomplete field will be output
346 IF (lclone) THEN
347 CALL delete(that%gaid(i3,i,j,i6))
348 ELSE
349 CALL init(that%gaid(i3,i,j,i6)) ! grid_id lacks a nullify method
350 ENDIF
351#ifdef DEBUG
352 CALL l4f_log(l4f_debug, &
353 'compute_stat_proc_agg, skipping lev/t/tr/var: '// &
354 t2c(i3)//'/'//t2c(i)//'/'//t2c(j)//'/'//t2c(i6))
355#endif
356 ENDIF
357 ENDIF ! ndtr > 0
358
359 ENDDO ! level
360 ENDDO ! var
361 ENDDO ! dtratio
362 CALL delete(map_ttr(i,j))
363 ENDDO do_otime
364ENDDO do_otimerange
365
366DEALLOCATE(dtratio, map_ttr)
367
368END SUBROUTINE volgrid6d_recompute_stat_proc_agg
369
370
394SUBROUTINE volgrid6d_compute_stat_proc_agg(this, that, stat_proc, &
395 step, start, full_steps, max_step, clone)
396TYPE(volgrid6d),INTENT(inout) :: this
397TYPE(volgrid6d),INTENT(out) :: that
398INTEGER,INTENT(in) :: stat_proc
399TYPE(timedelta),INTENT(in) :: step
400TYPE(datetime),INTENT(in),OPTIONAL :: start
401LOGICAL,INTENT(in),OPTIONAL :: full_steps
402TYPE(timedelta),INTENT(in),OPTIONAL :: max_step
403LOGICAL , INTENT(in),OPTIONAL :: clone
404
405INTEGER :: tri
406INTEGER i, j, n, ninp, i3, i6
407TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
408TYPE(timedelta) :: lmax_step
409LOGICAL :: lclone
410REAL,POINTER :: voldatiin(:,:), voldatiout(:,:)
411
412
413NULLIFY(voldatiin, voldatiout)
414tri = 254
415IF (PRESENT(max_step)) THEN
416 lmax_step = max_step
417ELSE
418 lmax_step = timedelta_max
419ENDIF
421CALL init(that)
422! be safe
423CALL volgrid6d_alloc_vol(this)
424
425! when volume is not decoded it is better to clone anyway to avoid
426! overwriting fields
427lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati)
428! initialise the output volume
429CALL init(that, griddim=this%griddim, time_definition=this%time_definition)
430CALL volgrid6d_alloc(that, dim=this%griddim%dim, ntimerange=1, &
431 nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.)
432that%level = this%level
433that%var = this%var
434
435CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
436 step, this%time_definition, that%time, that%timerange, map_ttr, &
437 start=start, full_steps=full_steps)
438
439CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati))
440
441do_otimerange: DO j = 1, SIZE(that%timerange)
442 do_otime: DO i = 1, SIZE(that%time)
443 ninp = map_ttr(i,j)%arraysize
444 IF (ninp <= 0) cycle do_otime
445
446 IF (stat_proc == 4) THEN ! check validity for difference
447 IF (map_ttr(i,j)%array(1)%extra_info /= 1 .OR. &
448 map_ttr(i,j)%array(ninp)%extra_info /= 2) THEN
449 CALL delete(map_ttr(i,j))
450 cycle do_otime
451 ENDIF
452 ELSE
453! check validity condition (missing values in volume are not accounted for)
454 DO n = 2, ninp
455 IF (map_ttr(i,j)%array(n)%time - map_ttr(i,j)%array(n-1)%time > &
456 lmax_step) THEN
457 CALL delete(map_ttr(i,j))
458 cycle do_otime
459 ENDIF
460 ENDDO
461 ENDIF
462
463 DO i6 = 1, SIZE(this%var)
464 DO i3 = 1, SIZE(this%level)
465 CALL volgrid_get_vol_2d(that, i3, i, j, i6, voldatiout)
466
467 IF (stat_proc == 4) THEN ! special treatment for difference
468 IF (lclone) THEN
469 CALL copy(this%gaid(i3, map_ttr(i,j)%array(1)%it,&
470 map_ttr(i,j)%array(1)%itr,i6), that%gaid(i3,i,j,i6))
471 ELSE
472 that%gaid(i3,i,j,i6) = this%gaid(i3, map_ttr(i,j)%array(1)%it, &
473 map_ttr(i,j)%array(1)%itr,i6)
474 ENDIF
475! improve the next workflow?
476 CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(ninp)%it, &
477 map_ttr(i,j)%array(ninp)%itr, i6, voldatiin)
478 voldatiout = voldatiin
479 CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(1)%it, &
480 map_ttr(i,j)%array(1)%itr, i6, voldatiin)
481
482 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
483 voldatiout(:,:) = voldatiout(:,:) - voldatiin(:,:)
484 ELSEWHERE
485 voldatiout(:,:) = rmiss
486 END WHERE
487
488 ELSE ! other stat_proc
489 DO n = 1, ninp
490 CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(n)%it, &
491 map_ttr(i,j)%array(n)%itr, i6, voldatiin)
492
493 IF (n == 1) THEN
494 voldatiout = voldatiin
495 IF (lclone) THEN
496 CALL copy(this%gaid(i3, map_ttr(i,j)%array(n)%it,&
497 map_ttr(i,j)%array(n)%itr,i6), that%gaid(i3,i,j,i6))
498 ELSE
499 that%gaid(i3,i,j,i6) = this%gaid(i3, map_ttr(i,j)%array(n)%it, &
500 map_ttr(i,j)%array(n)%itr,i6)
501 ENDIF
502
503 ELSE ! second or more time
504 SELECT CASE(stat_proc)
505 CASE (0, 1) ! average, accumulation
506 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
507 voldatiout(:,:) = voldatiout(:,:) + voldatiin(:,:)
508 ELSEWHERE
509 voldatiout(:,:) = rmiss
510 END WHERE
511 CASE(2) ! maximum
512 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
513 voldatiout(:,:) = max(voldatiout(:,:), voldatiin(:,:))
514 ELSEWHERE
515 voldatiout(:,:) = rmiss
516 END WHERE
517 CASE(3) ! minimum
518 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
519 voldatiout(:,:) = min(voldatiout(:,:), voldatiin(:,:))
520 ELSEWHERE
521 voldatiout(:,:) = rmiss
522 END WHERE
523 END SELECT
524
525 ENDIF ! first time
526 ENDDO
527 IF (stat_proc == 0) THEN ! average
528 WHERE(c_e(voldatiout(:,:)))
529 voldatiout(:,:) = voldatiout(:,:)/ninp
530 END WHERE
531 ENDIF
532 ENDIF
533 CALL volgrid_set_vol_2d(that, i3, i, j, i6, voldatiout)
534 ENDDO ! level
535 ENDDO ! var
536 CALL delete(map_ttr(i,j))
537 ENDDO do_otime
538ENDDO do_otimerange
539
540DEALLOCATE(map_ttr)
541
542
543END SUBROUTINE volgrid6d_compute_stat_proc_agg
544
545
570SUBROUTINE volgrid6d_recompute_stat_proc_diff(this, that, stat_proc, step, full_steps, start, clone)
571TYPE(volgrid6d),INTENT(inout) :: this
572TYPE(volgrid6d),INTENT(out) :: that
573INTEGER,INTENT(in) :: stat_proc
574TYPE(timedelta),INTENT(in) :: step
575LOGICAL,INTENT(in),OPTIONAL :: full_steps
576TYPE(datetime),INTENT(in),OPTIONAL :: start
577LOGICAL,INTENT(in),OPTIONAL :: clone
578INTEGER :: i3, i4, i6, i, j, k, l, nitr, steps
579INTEGER,ALLOCATABLE :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
580REAL,POINTER :: voldatiin1(:,:), voldatiin2(:,:), voldatiout(:,:)
581!LOGICAL,POINTER :: mask_timerange(:)
582LOGICAL :: lclone
583TYPE(vol7d_var),ALLOCATABLE :: varbufr(:)
584
585
586! be safe
587CALL volgrid6d_alloc_vol(this)
588! when volume is not decoded it is better to clone anyway to avoid
589! overwriting fields
590lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati)
591! initialise the output volume
592CALL init(that, griddim=this%griddim, time_definition=this%time_definition)
593CALL volgrid6d_alloc(that, dim=this%griddim%dim, &
594 nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.)
595that%level = this%level
596that%var = this%var
597
598! compute length of cumulation step in seconds
599CALL getval(step, asec=steps)
600
601! compute the statistical processing relations, output time and
602! timerange are defined here
603CALL recompute_stat_proc_diff_common(this%time, this%timerange, stat_proc, step, &
604 that%time, that%timerange, map_tr, f, keep_tr, &
605 this%time_definition, full_steps, start)
606nitr = SIZE(f)
607
608! complete the definition of the output volume
609CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati))
610! allocate workspace once
611IF (.NOT.ASSOCIATED(that%voldati)) THEN
612 ALLOCATE(voldatiin1(this%griddim%dim%nx, this%griddim%dim%ny), &
613 voldatiin2(this%griddim%dim%nx, this%griddim%dim%ny), &
614 voldatiout(this%griddim%dim%nx, this%griddim%dim%ny))
615ENDIF
616
617! copy the timeranges already satisfying the requested step, if any
618DO i4 = 1, SIZE(this%time)
619 DO i = 1, nitr
620 IF (c_e(keep_tr(i, i4, 2))) THEN
621 l = keep_tr(i, i4, 1)
622 k = keep_tr(i, i4, 2)
623#ifdef DEBUG
624 CALL l4f_category_log(this%category, l4f_debug, &
625 'volgrid6d_recompute_stat_proc_diff, good timerange: '//t2c(f(i))// &
626 '->'//t2c(k))
627#endif
628 DO i6 = 1, SIZE(this%var)
629 DO i3 = 1, SIZE(this%level)
630 IF (c_e(this%gaid(i3,i4,f(i),i6))) THEN
631 IF (lclone) THEN
632 CALL copy(this%gaid(i3,i4,f(i),i6), that%gaid(i3,l,k,i6))
633 ELSE
634 that%gaid(i3,l,k,i6) = this%gaid(i3,i4,f(i),i6)
635 ENDIF
636 IF (ASSOCIATED(that%voldati)) THEN
637 that%voldati(:,:,i3,l,k,i6) = this%voldati(:,:,i3,i4,f(i),i6)
638 ELSE
639 CALL volgrid_get_vol_2d(this, i3, i4, f(i), i6, voldatiout)
640 CALL volgrid_set_vol_2d(that, i3, l, k, i6, voldatiout)
641 ENDIF
642 ENDIF
643 ENDDO
644 ENDDO
645 ENDIF
646 ENDDO
647ENDDO
648
649! varbufr required for setting posdef, optimize with an array
650ALLOCATE(varbufr(SIZE(this%var)))
651DO i6 = 1, SIZE(this%var)
652 varbufr(i6) = convert(this%var(i6))
653ENDDO
654! compute statistical processing
655DO l = 1, SIZE(this%time)
656 DO k = 1, nitr
657 DO j = 1, SIZE(this%time)
658 DO i = 1, nitr
659 IF (c_e(map_tr(i,j,k,l,1))) THEN
660 DO i6 = 1, SIZE(this%var)
661 DO i3 = 1, SIZE(this%level)
662
663 IF (c_e(this%gaid(i3,j,f(i),i6)) .AND. &
664 c_e(this%gaid(i3,l,f(k),i6))) THEN
665! take the gaid from the second time/timerange contributing to the
666! result (l,f(k))
667 IF (lclone) THEN
668 CALL copy(this%gaid(i3,l,f(k),i6), &
669 that%gaid(i3,map_tr(i,j,k,l,1),map_tr(i,j,k,l,2),i6))
670 ELSE
671 that%gaid(i3,map_tr(i,j,k,l,1),map_tr(i,j,k,l,2),i6) = &
672 this%gaid(i3,l,f(k),i6)
673 ENDIF
674
675! get/set 2d sections API is used
676 CALL volgrid_get_vol_2d(this, i3, l, f(k), i6, voldatiin1)
677 CALL volgrid_get_vol_2d(this, i3, j, f(i), i6, voldatiin2)
678 IF (ASSOCIATED(that%voldati)) &
679 CALL volgrid_get_vol_2d(that, i3, &
680 map_tr(i,j,k,l,1), map_tr(i,j,k,l,2), i6, voldatiout)
681
682 IF (stat_proc == 0) THEN ! average
683 WHERE(c_e(voldatiin1(:,:)) .AND. c_e(voldatiin2(:,:)))
684 voldatiout(:,:) = &
685 (voldatiin1(:,:)*this%timerange(f(k))%p2 - &
686 voldatiin2(:,:)*this%timerange(f(i))%p2)/ &
687 steps
688 ELSEWHERE
689 voldatiout(:,:) = rmiss
690 END WHERE
691 ELSE IF (stat_proc == 1 .OR. stat_proc == 4) THEN ! acc, diff
692 WHERE(c_e(voldatiin1(:,:)) .AND. c_e(voldatiin2(:,:)))
693 voldatiout(:,:) = voldatiin1(:,:) - voldatiin2(:,:)
694 ELSEWHERE
695 voldatiout(:,:) = rmiss
696 END WHERE
697 IF (stat_proc == 1) THEN
698 CALL vol7d_var_features_posdef_apply(varbufr(i6), voldatiout)
699 ENDIF
700 ENDIF
701
702 CALL volgrid_set_vol_2d(that, i3, &
703 map_tr(i,j,k,l,1), map_tr(i,j,k,l,2), i6, voldatiout)
704
705 ENDIF
706 ENDDO
707 ENDDO
708 ENDIF
709 ENDDO
710 ENDDO
711 ENDDO
712ENDDO
713
714IF (.NOT.ASSOCIATED(that%voldati)) THEN
715 DEALLOCATE(voldatiin1, voldatiin2, voldatiout)
716ENDIF
717
718END SUBROUTINE volgrid6d_recompute_stat_proc_diff
719
720
748SUBROUTINE volgrid6d_compute_stat_proc_metamorph(this, that, stat_proc_input, stat_proc, clone)
749TYPE(volgrid6d),INTENT(inout) :: this
750TYPE(volgrid6d),INTENT(out) :: that
751INTEGER,INTENT(in) :: stat_proc_input
752INTEGER,INTENT(in) :: stat_proc
753LOGICAL , INTENT(in),OPTIONAL :: clone
754
755INTEGER :: j, i3, i4, i6
756INTEGER,POINTER :: map_tr(:)
757REAL,POINTER :: voldatiin(:,:), voldatiout(:,:)
758REAL,ALLOCATABLE :: int_ratio(:)
759LOGICAL :: lclone
760
761NULLIFY(voldatiin, voldatiout)
762
763! be safe
764CALL volgrid6d_alloc_vol(this)
765! when volume is not decoded it is better to clone anyway to avoid
766! overwriting fields
767lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati)
768
769IF (.NOT.((stat_proc_input == 0 .AND. stat_proc == 1) .OR. &
770 (stat_proc_input == 1 .AND. stat_proc == 0))) THEN
771
772 CALL l4f_category_log(this%category, l4f_warn, &
773 'compute_stat_proc_metamorph, can only be applied to average->accumulated timerange and viceversa')
774! return an empty volume, without signaling error
775 CALL init(that)
776 CALL volgrid6d_alloc_vol(that)
777 RETURN
778ENDIF
779
780! initialise the output volume
781CALL init(that, griddim=this%griddim, time_definition=this%time_definition)
782CALL volgrid6d_alloc(that, dim=this%griddim%dim, ntime=SIZE(this%time), &
783 nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.)
784that%time = this%time
785that%level = this%level
786that%var = this%var
787
788CALL compute_stat_proc_metamorph_common(stat_proc_input, this%timerange, stat_proc, &
789 that%timerange, map_tr)
790
791! complete the definition of the output volume
792CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati))
793
794IF (stat_proc == 0) THEN ! average -> integral
795 int_ratio = 1./real(that%timerange(:)%p2)
796ELSE ! cumulation
797 int_ratio = real(that%timerange(:)%p2)
798ENDIF
799
800DO i6 = 1, SIZE(this%var)
801 DO j = 1, SIZE(map_tr)
802 DO i4 = 1, SIZE(that%time)
803 DO i3 = 1, SIZE(this%level)
804
805 IF (lclone) THEN
806 CALL copy(this%gaid(i3,i4,map_tr(j),i6), that%gaid(i3,i4,j,i6))
807 ELSE
808 that%gaid(i3,i4,map_tr(j),i6) = this%gaid(i3,i4,j,i6)
809 ENDIF
810 CALL volgrid_get_vol_2d(this, i3, i4, map_tr(j), i6, voldatiin)
811 CALL volgrid_get_vol_2d(that, i3, i4, j, i6, voldatiout)
812 WHERE (c_e(voldatiin))
813 voldatiout = voldatiin*int_ratio(j)
814 ELSEWHERE
815 voldatiout = rmiss
816 END WHERE
817 CALL volgrid_set_vol_2d(that, i3, i4, j, i6, voldatiout)
818 ENDDO
819 ENDDO
820 ENDDO
821ENDDO
822
823
824END SUBROUTINE volgrid6d_compute_stat_proc_metamorph
825
840SUBROUTINE volgrid6d_compute_vert_coord_var(this, level, volgrid_lev)
841TYPE(volgrid6d),INTENT(in) :: this
842TYPE(vol7d_level),INTENT(in) :: level
843TYPE(volgrid6d),INTENT(out) :: volgrid_lev
844
845INTEGER :: nlev, i, ii, iii, iiii
846TYPE(grid_id) :: out_gaid
847LOGICAL,ALLOCATABLE :: levmask(:)
848TYPE(volgrid6d_var) :: lev_var
849
850CALL init(volgrid_lev) ! initialise to null
851IF (.NOT.ASSOCIATED(this%gaid)) THEN
852 CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: input volume not allocated')
853 RETURN
854ENDIF
855! if layer, both surfaces must be of the same type
856IF (c_e(level%level2) .AND. level%level1 /= level%level2) THEN
857 CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: requested (mixed) layer type not valid')
858 RETURN
859ENDIF
860
861! look for valid levels to be converted to vars
862ALLOCATE(levmask(SIZE(this%level)))
863levmask = this%level%level1 == level%level1 .AND. &
864 this%level%level2 == level%level2 .AND. c_e(this%level%l1)
865IF (c_e(level%level2)) levmask = levmask .AND. c_e(this%level%l2)
866nlev = count(levmask)
867IF (nlev == 0) THEN
868 CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: requested level type not available')
869 RETURN
870ENDIF
871
872out_gaid = grid_id_new()
873gaidloop: DO i=1 ,SIZE(this%gaid,1)
874 DO ii=1 ,SIZE(this%gaid,2)
875 DO iii=1 ,SIZE(this%gaid,3)
876 DO iiii=1 ,SIZE(this%gaid,4)
877 IF (c_e(this%gaid(i,ii,iii,iiii))) THEN ! conserve first valid gaid
878 CALL copy(this%gaid(i,ii,iii,iiii), out_gaid)
879 EXIT gaidloop
880 ENDIF
881 ENDDO
882 ENDDO
883 ENDDO
884ENDDO gaidloop
885
886! look for variable corresponding to level
887lev_var = convert(vol7d_var_new(btable=vol7d_level_to_var(level)), &
888 grid_id_template=out_gaid)
889IF (.NOT.c_e(lev_var)) THEN
890 CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: no variable corresponds to requested level type')
891 RETURN
892ENDIF
893
894! prepare output volume
895CALL init(volgrid_lev, griddim=this%griddim, &
896 time_definition=this%time_definition) !, categoryappend=categoryappend)
897CALL volgrid6d_alloc(volgrid_lev, ntime=SIZE(this%time), nlevel=nlev, &
898 ntimerange=SIZE(this%timerange), nvar=1)
899! fill metadata
900volgrid_lev%time = this%time
901volgrid_lev%level = pack(this%level, mask=levmask)
902volgrid_lev%timerange = this%timerange
903volgrid_lev%var(1) = lev_var
904
905CALL volgrid6d_alloc_vol(volgrid_lev, decode=.true.)
906! fill data
907DO i = 1, nlev
908 IF (c_e(level%level2)) THEN
909 volgrid_lev%voldati(:,:,i,:,:,:) = real(volgrid_lev%level(i)%l1 + &
910 volgrid_lev%level(i)%l2)* &
911 vol7d_level_to_var_factor(volgrid_lev%level(i))/2.
912 ELSE
913 volgrid_lev%voldati(:,:,i,:,:,:) = real(volgrid_lev%level(i)%l1)* &
914 vol7d_level_to_var_factor(volgrid_lev%level(i))
915 ENDIF
916ENDDO
917! fill gaid for subsequent export
918IF (c_e(out_gaid)) THEN
919 DO i=1 ,SIZE(volgrid_lev%gaid,1)
920 DO ii=1 ,SIZE(volgrid_lev%gaid,2)
921 DO iii=1 ,SIZE(volgrid_lev%gaid,3)
922 DO iiii=1 ,SIZE(volgrid_lev%gaid,4)
923 CALL copy(out_gaid, volgrid_lev%gaid(i,ii,iii,iiii))
924 ENDDO
925 ENDDO
926 ENDDO
927 ENDDO
928 CALL delete(out_gaid)
929ENDIF
930
931END SUBROUTINE volgrid6d_compute_vert_coord_var
932
Distruttori per le 2 classi.
Restituiscono il valore dell'oggetto nella forma desiderata.
Costruttori per le classi datetime e timedelta.
Functions that return a trimmed CHARACTER representation of the input variable.
Restituiscono il valore dell'oggetto in forma di stringa stampabile.
Make a deep copy, if possible, of the grid identifier.
Apply the conversion function this to values.
Classi per la gestione delle coordinate temporali.
This module defines an abstract interface to different drivers for access to files containing gridded...
Module for basic statistical computations taking into account missing data.
Definition: simple_stat.f90:25
This module contains functions that are only for internal use of the library.
Extension of volgrid6d_class with methods for performing simple statistical operations on entire volu...
This module defines objects and methods for managing data volumes on rectangular georeferenced grids.
Class for managing physical variables in a grib 1/2 fashion.
Class for expressing an absolute time value.
Class for expressing a relative time interval.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...
Derived type defining a dynamically extensible array of TYPE(ttr_mapper) elements.
Object describing a rectangular, homogeneous gridded dataset.
Definition of a physical variable in grib coding style.

Generated with Doxygen.