libsim Versione 7.2.6
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!> Extension of volgrid6d_class with methods for performing simple
20!! statistical operations on entire volumes of data. This module
21!! includes methods for performing operations such as (re)averaging
22!! and (re)cumulation in time of entire volgrid6d_class::volgrid6d
23!! objects.
24!!
25!! \ingroup volgrid6d
33IMPLICIT NONE
34
35CONTAINS
36
37!> General-purpose method for computing a statistical processing on
38!! data in a volgrid6d object already processed with the same statistical
39!! processing, on a different time interval specified by \a step and
40!! \a start. This method tries to apply all the suitable specialized
41!! statistical processing methods according to the input and output
42!! statistical processing requested. The argument \a stat_proc_input
43!! determines which data will be processed, while the \a stat_proc
44!! argument determines the type of statistical process to be applied
45!! and which will be owned by output data.
46!!
47!! The possible combinations are:
48!!
49!! - \a stat_proc_input = 254
50!! - \a stat_proc = 0 average instantaneous observations
51!! - \a stat_proc = 2 compute maximum of instantaneous observations
52!! - \a stat_proc = 3 compute minimum of instantaneous observations
53!! - \a stat_proc = 4 compute difference of instantaneous observations
54!!
55!! processing is computed on longer time intervals by aggregation,
56!! see the description of volgrid6d_compute_stat_proc_agg()
57!!
58!! - \a stat_proc_input = *
59!! - \a stat_proc = 254 consider statistically processed values as
60!! instantaneous without any extra processing
61!!
62!! see the description of volgrid6d_decompute_stat_proc()
63!!
64!! - \a stat_proc_input = 0, 1, 2, 3, 4, 200
65!! - \a stat_proc = \a stat_proc_input recompute input data on
66!! different intervals
67!!
68!! the same statistical processing is applied to obtain data
69!! processed on a different interval, either longer, by
70!! aggregation, or shorter, by differences, see the description of
71!! volgrid6d_recompute_stat_proc_agg() and
72!! volgrid6d_recompute_stat_proc_diff() respectively; it is also
73!! possible to provide \a stat_proc_input \a /= \a stat_proc, but
74!! it has to be used with care.
75!!
76!! - \a stat_proc_input = 0
77!! - \a stat_proc = 1
78!!
79!! a time-averaged rate or flux is transformed into a
80!! time-integrated value (sometimes called accumulated) on the same
81!! interval by multiplying the values by the length of the time
82!! interval in seconds, keeping constant all the rest, including
83!! the variable; the unit of the variable implicitly changes
84!! accordingly, this is supported officially in grib2 standard, in
85!! the other cases it is a forcing of the standards.
86!!
87!! - \a stat_proc_input = 1
88!! - \a stat_proc = 0
89!!
90!! a time-integrated value (sometimes called accumulated) is
91!! transformed into a time-averaged rate or flux on the same
92!! interval by dividing the values by the length of the time
93!! interval in seconds, see also the previous description of the
94!! opposite computation.
95!!
96!! If a particular statistical processing cannot be performed on the
97!! input data, the program continues with a warning and, if requested,
98!! the input data is passed over to the volume specified by the \a
99!! other argument, in order to allow continuation of processing. All
100!! the other parameters are passed over to the specifical statistical
101!! processing methods and are documented there.
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 !< volume providing data to be recomputed, it is not modified by the method, apart from performing a \a volgrid6d_alloc_vol on it
105TYPE(volgrid6d),INTENT(out) :: that !< output volume which will contain the recomputed data
106INTEGER,INTENT(in) :: stat_proc_input !< type of statistical processing of data that has to be processed (from grib2 table), only data having timerange of this type will be processed, the actual statistical processing performed and which will appear in the output volume, is however determined by \a stat_proc argument
107INTEGER,INTENT(in) :: stat_proc !< type of statistical processing to be recomputed (from grib2 table), data in output volume \a that will have a timerange of this type
108TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed
109TYPE(datetime),INTENT(in),OPTIONAL :: start !< start of statistical processing interval
110LOGICAL,INTENT(in),OPTIONAL :: full_steps !< if \a .TRUE. apply processing only on intervals starting at a forecast time or a reference time modulo \a step
111REAL,INTENT(in),OPTIONAL :: frac_valid !< minimum fraction of valid data required for considering acceptable a recomputed value, default=1.
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 !< if provided and \c .TRUE., the statistical process is computed, if possible, by weighting every value with a weight proportional to its validity interval
114LOGICAL , INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a this to \a that
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 .OR. &
133 (stat_proc == 0 .OR. stat_proc == 2 .OR. stat_proc == 3)) THEN
134! avg, min and max can be computed from any input, with care
135
136 IF (count(this%timerange(:)%timerange == stat_proc_input) == 0) THEN
137 CALL l4f_category_log(this%category, l4f_warn, &
138 'no timeranges of the desired statistical processing type '//t2c(stat_proc)//' available')
139! return an empty volume, without signaling error
140 CALL init(that)
141 CALL volgrid6d_alloc_vol(that)
142
143 ELSE
144! euristically determine whether aggregation or difference is more suitable
145 dtmax = maxval(this%timerange(:)%p2, &
146 mask=(this%timerange(:)%timerange == stat_proc))
147 CALL getval(step, asec=dtstep)
148
149#ifdef DEBUG
150 CALL l4f_category_log(this%category, l4f_debug, &
151 'stat_proc='//t2c(stat_proc)//' dtmax='//t2c(dtmax)//' dtstep='//t2c(dtstep))
152#endif
153
154 IF (dtstep <= dtmax) THEN
155 CALL l4f_category_log(this%category, l4f_info, &
156 'recomputing statistically processed data by difference '// &
157 t2c(stat_proc_input)//':'//t2c(stat_proc))
158 CALL volgrid6d_recompute_stat_proc_diff(this, that, stat_proc, step, &
159 full_steps, start, clone)
160 ELSE
161 CALL l4f_category_log(this%category, l4f_info, &
162 'recomputing statistically processed data by aggregation '// &
163 t2c(stat_proc_input)//':'//t2c(stat_proc))
164 CALL volgrid6d_recompute_stat_proc_agg(this, that, stat_proc, step, start, &
165 full_steps, frac_valid, clone, stat_proc_input)
166 ENDIF
167 ENDIF
168
169ELSE ! IF (stat_proc_input /= stat_proc) THEN
170 IF ((stat_proc_input == 0 .AND. stat_proc == 1) .OR. &
171 (stat_proc_input == 1 .AND. stat_proc == 0)) THEN
172 CALL l4f_category_log(this%category, l4f_info, &
173 'computing statistically processed data by integration/differentiation '// &
174 t2c(stat_proc_input)//':'//t2c(stat_proc))
175 CALL volgrid6d_compute_stat_proc_metamorph(this, that, stat_proc_input, &
176 stat_proc, clone)
177 ELSE
178 CALL l4f_category_log(this%category, l4f_error, &
179 'statistical processing '//t2c(stat_proc_input)//':'//t2c(stat_proc)// &
180 ' not implemented or does not make sense')
181 CALL raise_error()
182 ENDIF
183
184ENDIF
185
186END SUBROUTINE volgrid6d_compute_stat_proc
187
188
189!> Specialized method for statistically processing a set of data
190!! already processed with the same statistical processing, on a
191!! different time interval. This method performs statistical
192!! processing by aggregation of shorter intervals.
193!!
194!! The output \a that volgrid6d object contains elements from the original volume
195!! \a this satisfying the conditions
196!! - timerange (vol7d_timerange_class::vol7d_timerange::timerange)
197!! of type \a stat_proc (or \a stat_proc_input if provided)
198!! - any p1 (analysis/observation or forecast)
199!! - p2 > 0 (processing interval non null, non instantaneous data)
200!! and equal to a multiplier of \a step
201!!
202!! Output data will have timerange of type \a stat_proc and p2 = \a
203!! step. The supported statistical processing methods (parameter \a
204!! stat_proc) are:
205!!
206!! - 0 average
207!! - 1 accumulation
208!! - 2 maximum
209!! - 3 minimum
210!! - 4 difference
211!! - 200 vectorial mean
212!!
213!! The start of processing period can be computed automatically from
214!! the input intervals as the first possible interval modulo \a step,
215!! or, for a better control, it can be specified explicitely by the
216!! optional argument \a start. Be warned that, in the final volume,
217!! the first reference time will actually be \a start \a + \a step,
218!! since \a start indicates the beginning of first processing
219!! interval, while reference time (for analysis/oservation) is the end
220!! of the interval.
221!!
222!! The purpose of the optional argument \a stat_proc_input is to allow
223!! processing with a certain statistical processing operator a dataset
224!! already processed with a different operator, by specifying the
225!! latter as stat_proc_input; this is useful, for example, if one
226!! wants to compute the monthly average of daily maximum temperatures;
227!! however this has to be used with care since the resulting data
228!! volume will not carry all the information about the processing
229!! which has been done, in the previous case, for example, the
230!! temperatures will simply look like monthly average temperatures.
231SUBROUTINE volgrid6d_recompute_stat_proc_agg(this, that, stat_proc, &
232 step, start, full_steps, frac_valid, clone, stat_proc_input)
233TYPE(volgrid6d),INTENT(inout) :: this !< volume providing data to be recomputed, it is not modified by the method, apart from performing a \a volgrid6d_alloc_vol on it
234TYPE(volgrid6d),INTENT(out) :: that !< output volume which will contain the recomputed data
235INTEGER,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
236TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed
237TYPE(datetime),INTENT(in),OPTIONAL :: start !< start of statistical processing interval
238LOGICAL,INTENT(in),OPTIONAL :: full_steps !< if \a .TRUE. and \a start is not provided, apply processing only on intervals starting at a forecast time or a reference time modulo \a step
239REAL,INTENT(in),OPTIONAL :: frac_valid !< minimum fraction of valid data required for considering acceptable a recomputed value, default=1.
240LOGICAL, INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a this to \a that
241INTEGER,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
242
243INTEGER :: tri
244INTEGER i, j, n, n1, ndtr, i3, i6
245TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
246INTEGER,POINTER :: dtratio(:)
247REAL :: lfrac_valid
248LOGICAL :: lclone
249REAL,POINTER :: voldatiin(:,:), voldatiout(:,:)
250
251
252NULLIFY(voldatiin, voldatiout)
253IF (PRESENT(stat_proc_input)) THEN
254 tri = stat_proc_input
255ELSE
256 tri = stat_proc
257ENDIF
258IF (PRESENT(frac_valid)) THEN
259 lfrac_valid = frac_valid
260ELSE
261 lfrac_valid = 1.0
262ENDIF
263
264CALL init(that)
265! be safe
266CALL volgrid6d_alloc_vol(this)
267
268! when volume is not decoded it is better to clone anyway to avoid
269! overwriting fields
270lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati)
271! initialise the output volume
272CALL init(that, griddim=this%griddim, time_definition=this%time_definition)
273CALL volgrid6d_alloc(that, dim=this%griddim%dim, ntimerange=1, &
274 nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.)
275that%level = this%level
276that%var = this%var
277
278CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
279 step, this%time_definition, that%time, that%timerange, map_ttr, &
280 dtratio=dtratio, start=start, full_steps=full_steps)
281
282CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati))
283
284do_otimerange: DO j = 1, SIZE(that%timerange)
285 do_otime: DO i = 1, SIZE(that%time)
286
287 DO n1 = 1, SIZE(dtratio)
288 IF (dtratio(n1) <= 0) cycle ! safety check
289
290 DO i6 = 1, SIZE(this%var)
291 DO i3 = 1, SIZE(this%level)
292 CALL volgrid_get_vol_2d(that, i3, i, j, i6, voldatiout)
293 ndtr = 0
294 DO n = 1, map_ttr(i,j)%arraysize
295 IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
296 ndtr = ndtr + 1
297 CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(n)%it, &
298 map_ttr(i,j)%array(n)%itr, i6, voldatiin)
299
300 IF (ndtr == 1) THEN
301 voldatiout = voldatiin
302 IF (lclone) THEN
303 CALL copy(this%gaid(i3, map_ttr(i,j)%array(n)%it,&
304 map_ttr(i,j)%array(n)%itr,i6), that%gaid(i3,i,j,i6))
305 ELSE
306 that%gaid(i3,i,j,i6) = this%gaid(i3, map_ttr(i,j)%array(n)%it, &
307 map_ttr(i,j)%array(n)%itr,i6)
308 ENDIF
309
310 ELSE ! second or more time
311 SELECT CASE(stat_proc)
312 CASE (0, 200, 1, 4) ! average, vectorial mean, accumulation, difference
313 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
314 voldatiout(:,:) = voldatiout(:,:) + voldatiin(:,:)
315 ELSEWHERE
316 voldatiout(:,:) = rmiss
317 END WHERE
318 CASE(2) ! maximum
319 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
320 voldatiout(:,:) = max(voldatiout(:,:), voldatiin(:,:))
321 ELSEWHERE
322 voldatiout(:,:) = rmiss
323 END WHERE
324 CASE(3) ! minimum
325 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
326 voldatiout(:,:) = min(voldatiout(:,:), voldatiin(:,:))
327 ELSEWHERE
328 voldatiout(:,:) = rmiss
329 END WHERE
330 END SELECT
331
332 ENDIF ! first time
333 ENDIF ! dtratio(n1)
334 ENDDO ! ttr
335
336#ifdef DEBUG
337 CALL l4f_log(l4f_debug, &
338 'compute_stat_proc_agg, ndtr/dtratio/frac_valid: '// &
339 t2c(ndtr)//'/'//t2c(dtratio(n1))//'/'//t2c(lfrac_valid))
340#endif
341 IF (ndtr > 0) THEN ! why this condition was not here before?
342 IF (real(ndtr)/real(dtratio(n1)) >= lfrac_valid) THEN ! success
343 IF (stat_proc == 0) THEN ! average
344 WHERE(c_e(voldatiout(:,:)))
345 voldatiout(:,:) = voldatiout(:,:)/ndtr
346 END WHERE
347 ENDIF
348 CALL volgrid_set_vol_2d(that, i3, i, j, i6, voldatiout)
349#ifdef DEBUG
350 CALL l4f_log(l4f_debug, &
351 'compute_stat_proc_agg, coding lev/t/tr/var: '// &
352 t2c(i3)//'/'//t2c(i)//'/'//t2c(j)//'/'//t2c(i6))
353#endif
354 ELSE
355! must nullify the output gaid here, otherwise an incomplete field will be output
356 IF (lclone) THEN
357 CALL delete(that%gaid(i3,i,j,i6))
358 ELSE
359 CALL init(that%gaid(i3,i,j,i6)) ! grid_id lacks a nullify method
360 ENDIF
361#ifdef DEBUG
362 CALL l4f_log(l4f_debug, &
363 'compute_stat_proc_agg, skipping lev/t/tr/var: '// &
364 t2c(i3)//'/'//t2c(i)//'/'//t2c(j)//'/'//t2c(i6))
365#endif
366 ENDIF
367 ENDIF ! ndtr > 0
368
369 ENDDO ! level
370 ENDDO ! var
371 ENDDO ! dtratio
372 CALL delete(map_ttr(i,j))
373 ENDDO do_otime
374ENDDO do_otimerange
375
376DEALLOCATE(dtratio, map_ttr)
377
378END SUBROUTINE volgrid6d_recompute_stat_proc_agg
379
380
381!> Method for statistically processing a set of instantaneous data.
382!! This method performs statistical processing by aggregation of
383!! instantaneous data.
384!!
385!! The output \a that volgrid6d object contains elements from the original volume
386!! \a this satisfying the conditions
387!! - timerange (vol7d_timerange_class::vol7d_timerange::timerange)
388!! of type 254 (instantaeous)
389!! - any p1 (analysis/observation or forecast)
390!! - p2 = 0 (processing interval null, instantaneous data)
391!!
392!! Output data will have timerange of type \a stat_proc, and p2 = \a
393!! step. The supported statistical processing methods (parameter \a
394!! stat_proc) are:
395!!
396!! - 0 average
397!! - 2 maximum
398!! - 3 minimum
399!! - 4 difference
400!!
401!! A maximum distance in time for input valid data can be assigned
402!! with the optional argument \a max_step, in order to filter datasets
403!! with too long "holes".
404SUBROUTINE volgrid6d_compute_stat_proc_agg(this, that, stat_proc, &
405 step, start, full_steps, max_step, clone)
406TYPE(volgrid6d),INTENT(inout) :: this !< volume providing data to be recomputed, it is not modified by the method, apart from performing a \a volgrid6d_alloc_vol on it
407TYPE(volgrid6d),INTENT(out) :: that !< output volume which will contain the recomputed data
408INTEGER,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
409TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed
410TYPE(datetime),INTENT(in),OPTIONAL :: start !< start of statistical processing interval
411LOGICAL,INTENT(in),OPTIONAL :: full_steps !< if \a .TRUE. and \a start is not provided, apply processing only on intervals starting at a forecast time or a reference time modulo \a step
412TYPE(timedelta),INTENT(in),OPTIONAL :: max_step !< maximum allowed distance in time between two contiguougs valid data within an interval, for the interval to be eligible for statistical processing
413LOGICAL , INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a this to \a that
414
415INTEGER :: tri
416INTEGER i, j, n, ninp, i3, i6
417TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
418TYPE(timedelta) :: lmax_step
419LOGICAL :: lclone
420REAL,POINTER :: voldatiin(:,:), voldatiout(:,:)
421
422
423NULLIFY(voldatiin, voldatiout)
424tri = 254
425IF (PRESENT(max_step)) THEN
426 lmax_step = max_step
427ELSE
428 lmax_step = timedelta_max
429ENDIF
430
431CALL init(that)
432! be safe
433CALL volgrid6d_alloc_vol(this)
434
435! when volume is not decoded it is better to clone anyway to avoid
436! overwriting fields
437lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati)
438! initialise the output volume
439CALL init(that, griddim=this%griddim, time_definition=this%time_definition)
440CALL volgrid6d_alloc(that, dim=this%griddim%dim, ntimerange=1, &
441 nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.)
442that%level = this%level
443that%var = this%var
444
445CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
446 step, this%time_definition, that%time, that%timerange, map_ttr, &
447 start=start, full_steps=full_steps)
448
449CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati))
450
451do_otimerange: DO j = 1, SIZE(that%timerange)
452 do_otime: DO i = 1, SIZE(that%time)
453 ninp = map_ttr(i,j)%arraysize
454 IF (ninp <= 0) cycle do_otime
455
456 IF (stat_proc == 4) THEN ! check validity for difference
457 IF (map_ttr(i,j)%array(1)%extra_info /= 1 .OR. &
458 map_ttr(i,j)%array(ninp)%extra_info /= 2) THEN
459 CALL delete(map_ttr(i,j))
460 cycle do_otime
461 ENDIF
462 ELSE
463! check validity condition (missing values in volume are not accounted for)
464 DO n = 2, ninp
465 IF (map_ttr(i,j)%array(n)%time - map_ttr(i,j)%array(n-1)%time > &
466 lmax_step) THEN
467 CALL delete(map_ttr(i,j))
468 cycle do_otime
469 ENDIF
470 ENDDO
471 ENDIF
472
473 DO i6 = 1, SIZE(this%var)
474 DO i3 = 1, SIZE(this%level)
475 CALL volgrid_get_vol_2d(that, i3, i, j, i6, voldatiout)
476
477 IF (stat_proc == 4) THEN ! special treatment for difference
478 IF (lclone) THEN
479 CALL copy(this%gaid(i3, map_ttr(i,j)%array(1)%it,&
480 map_ttr(i,j)%array(1)%itr,i6), that%gaid(i3,i,j,i6))
481 ELSE
482 that%gaid(i3,i,j,i6) = this%gaid(i3, map_ttr(i,j)%array(1)%it, &
483 map_ttr(i,j)%array(1)%itr,i6)
484 ENDIF
485! improve the next workflow?
486 CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(ninp)%it, &
487 map_ttr(i,j)%array(ninp)%itr, i6, voldatiin)
488 voldatiout = voldatiin
489 CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(1)%it, &
490 map_ttr(i,j)%array(1)%itr, i6, voldatiin)
491
492 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
493 voldatiout(:,:) = voldatiout(:,:) - voldatiin(:,:)
494 ELSEWHERE
495 voldatiout(:,:) = rmiss
496 END WHERE
497
498 ELSE ! other stat_proc
499 DO n = 1, ninp
500 CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(n)%it, &
501 map_ttr(i,j)%array(n)%itr, i6, voldatiin)
502
503 IF (n == 1) THEN
504 voldatiout = voldatiin
505 IF (lclone) THEN
506 CALL copy(this%gaid(i3, map_ttr(i,j)%array(n)%it,&
507 map_ttr(i,j)%array(n)%itr,i6), that%gaid(i3,i,j,i6))
508 ELSE
509 that%gaid(i3,i,j,i6) = this%gaid(i3, map_ttr(i,j)%array(n)%it, &
510 map_ttr(i,j)%array(n)%itr,i6)
511 ENDIF
512
513 ELSE ! second or more time
514 SELECT CASE(stat_proc)
515 CASE (0, 1) ! average, accumulation
516 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
517 voldatiout(:,:) = voldatiout(:,:) + voldatiin(:,:)
518 ELSEWHERE
519 voldatiout(:,:) = rmiss
520 END WHERE
521 CASE(2) ! maximum
522 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
523 voldatiout(:,:) = max(voldatiout(:,:), voldatiin(:,:))
524 ELSEWHERE
525 voldatiout(:,:) = rmiss
526 END WHERE
527 CASE(3) ! minimum
528 WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:)))
529 voldatiout(:,:) = min(voldatiout(:,:), voldatiin(:,:))
530 ELSEWHERE
531 voldatiout(:,:) = rmiss
532 END WHERE
533 END SELECT
534
535 ENDIF ! first time
536 ENDDO
537 IF (stat_proc == 0) THEN ! average
538 WHERE(c_e(voldatiout(:,:)))
539 voldatiout(:,:) = voldatiout(:,:)/ninp
540 END WHERE
541 ENDIF
542 ENDIF
543 CALL volgrid_set_vol_2d(that, i3, i, j, i6, voldatiout)
544 ENDDO ! level
545 ENDDO ! var
546 CALL delete(map_ttr(i,j))
547 ENDDO do_otime
548ENDDO do_otimerange
549
550DEALLOCATE(map_ttr)
551
552
553END SUBROUTINE volgrid6d_compute_stat_proc_agg
554
555
556!> Specialized method for statistically processing a set of data
557!! already processed with the same statistical processing, on a
558!! different time interval. This method performs statistical
559!! processing by difference of different intervals. Data with both
560!! analysis/observation or forecast timerange are processed.
561!!
562!! The output \a that volgrid6d object contains elements from the original volume
563!! \a this satisfying the conditions
564!! - timerange (vol7d_timerange_class::vol7d_timerange::timerange)
565!! of type \a stat_proc
566!! - any p1 (analysis/observation or forecast)
567!! - p2 &gt; 0 (processing interval non null, non instantaneous data)
568!! and equal to a multiplier of \a step if \a full_steps is \c .TRUE.
569!!
570!! Output data will have timerange of type \a stat_proc and p2 = \a
571!! step. The supported statistical processing methods (parameter \a
572!! stat_proc) are:
573!!
574!! - 0 average
575!! - 1 accumulation
576!! - 4 difference
577!!
578!! Input volume may have any value of \a this%time_definition, and
579!! that value will be conserved in the output volume.
580SUBROUTINE volgrid6d_recompute_stat_proc_diff(this, that, stat_proc, step, full_steps, start, clone)
581TYPE(volgrid6d),INTENT(inout) :: this !< volume providing data to be recomputed, it is not modified by the method, apart from performing a \a volgrid6d_alloc_vol on it
582TYPE(volgrid6d),INTENT(out) :: that !< output volume which will contain the recomputed data
583INTEGER,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
584TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed
585LOGICAL,INTENT(in),OPTIONAL :: full_steps !< if provided and \a .TRUE., process only data having processing interval (p2) equal to a multiple of \a step
586TYPE(datetime),INTENT(in),OPTIONAL :: start !< if provided, together with \a full_steps, processes data on intervals starting at \a start +- an integer amount of \a step intervals
587LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a this to \a that
588INTEGER :: i3, i4, i6, i, j, k, l, nitr, steps
589INTEGER,ALLOCATABLE :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
590REAL,POINTER :: voldatiin1(:,:), voldatiin2(:,:), voldatiout(:,:)
591!LOGICAL,POINTER :: mask_timerange(:)
592LOGICAL :: lclone
593TYPE(vol7d_var),ALLOCATABLE :: varbufr(:)
594
595
596! be safe
597CALL volgrid6d_alloc_vol(this)
598! when volume is not decoded it is better to clone anyway to avoid
599! overwriting fields
600lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati)
601! initialise the output volume
602CALL init(that, griddim=this%griddim, time_definition=this%time_definition)
603CALL volgrid6d_alloc(that, dim=this%griddim%dim, &
604 nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.)
605that%level = this%level
606that%var = this%var
607
608! compute length of cumulation step in seconds
609CALL getval(step, asec=steps)
610
611! compute the statistical processing relations, output time and
612! timerange are defined here
613CALL recompute_stat_proc_diff_common(this%time, this%timerange, stat_proc, step, &
614 that%time, that%timerange, map_tr, f, keep_tr, &
615 this%time_definition, full_steps, start)
616nitr = SIZE(f)
617
618! complete the definition of the output volume
619CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati))
620! allocate workspace once
621IF (.NOT.ASSOCIATED(that%voldati)) THEN
622 ALLOCATE(voldatiin1(this%griddim%dim%nx, this%griddim%dim%ny), &
623 voldatiin2(this%griddim%dim%nx, this%griddim%dim%ny), &
624 voldatiout(this%griddim%dim%nx, this%griddim%dim%ny))
625ENDIF
626
627! copy the timeranges already satisfying the requested step, if any
628DO i4 = 1, SIZE(this%time)
629 DO i = 1, nitr
630 IF (c_e(keep_tr(i, i4, 2))) THEN
631 l = keep_tr(i, i4, 1)
632 k = keep_tr(i, i4, 2)
633#ifdef DEBUG
634 CALL l4f_category_log(this%category, l4f_debug, &
635 'volgrid6d_recompute_stat_proc_diff, good timerange: '//t2c(f(i))// &
636 '->'//t2c(k))
637#endif
638 DO i6 = 1, SIZE(this%var)
639 DO i3 = 1, SIZE(this%level)
640 IF (c_e(this%gaid(i3,i4,f(i),i6))) THEN
641 IF (lclone) THEN
642 CALL copy(this%gaid(i3,i4,f(i),i6), that%gaid(i3,l,k,i6))
643 ELSE
644 that%gaid(i3,l,k,i6) = this%gaid(i3,i4,f(i),i6)
645 ENDIF
646 IF (ASSOCIATED(that%voldati)) THEN
647 that%voldati(:,:,i3,l,k,i6) = this%voldati(:,:,i3,i4,f(i),i6)
648 ELSE
649 CALL volgrid_get_vol_2d(this, i3, i4, f(i), i6, voldatiout)
650 CALL volgrid_set_vol_2d(that, i3, l, k, i6, voldatiout)
651 ENDIF
652 ENDIF
653 ENDDO
654 ENDDO
655 ENDIF
656 ENDDO
657ENDDO
658
659! varbufr required for setting posdef, optimize with an array
660ALLOCATE(varbufr(SIZE(this%var)))
661DO i6 = 1, SIZE(this%var)
662 varbufr(i6) = convert(this%var(i6))
663ENDDO
664! compute statistical processing
665DO l = 1, SIZE(this%time)
666 DO k = 1, nitr
667 DO j = 1, SIZE(this%time)
668 DO i = 1, nitr
669 IF (c_e(map_tr(i,j,k,l,1))) THEN
670 DO i6 = 1, SIZE(this%var)
671 DO i3 = 1, SIZE(this%level)
672
673 IF (c_e(this%gaid(i3,j,f(i),i6)) .AND. &
674 c_e(this%gaid(i3,l,f(k),i6))) THEN
675! take the gaid from the second time/timerange contributing to the
676! result (l,f(k))
677 IF (lclone) THEN
678 CALL copy(this%gaid(i3,l,f(k),i6), &
679 that%gaid(i3,map_tr(i,j,k,l,1),map_tr(i,j,k,l,2),i6))
680 ELSE
681 that%gaid(i3,map_tr(i,j,k,l,1),map_tr(i,j,k,l,2),i6) = &
682 this%gaid(i3,l,f(k),i6)
683 ENDIF
684
685! get/set 2d sections API is used
686 CALL volgrid_get_vol_2d(this, i3, l, f(k), i6, voldatiin1)
687 CALL volgrid_get_vol_2d(this, i3, j, f(i), i6, voldatiin2)
688 IF (ASSOCIATED(that%voldati)) &
689 CALL volgrid_get_vol_2d(that, i3, &
690 map_tr(i,j,k,l,1), map_tr(i,j,k,l,2), i6, voldatiout)
691
692 IF (stat_proc == 0) THEN ! average
693 WHERE(c_e(voldatiin1(:,:)) .AND. c_e(voldatiin2(:,:)))
694 voldatiout(:,:) = &
695 (voldatiin1(:,:)*this%timerange(f(k))%p2 - &
696 voldatiin2(:,:)*this%timerange(f(i))%p2)/ &
697 steps
698 ELSEWHERE
699 voldatiout(:,:) = rmiss
700 END WHERE
701 ELSE IF (stat_proc == 1 .OR. stat_proc == 4) THEN ! acc, diff
702 WHERE(c_e(voldatiin1(:,:)) .AND. c_e(voldatiin2(:,:)))
703 voldatiout(:,:) = voldatiin1(:,:) - voldatiin2(:,:)
704 ELSEWHERE
705 voldatiout(:,:) = rmiss
706 END WHERE
707 IF (stat_proc == 1) THEN
708 CALL vol7d_var_features_posdef_apply(varbufr(i6), voldatiout)
709 ENDIF
710 ENDIF
711
712 CALL volgrid_set_vol_2d(that, i3, &
713 map_tr(i,j,k,l,1), map_tr(i,j,k,l,2), i6, voldatiout)
714
715 ENDIF
716 ENDDO
717 ENDDO
718 ENDIF
719 ENDDO
720 ENDDO
721 ENDDO
722ENDDO
723
724IF (.NOT.ASSOCIATED(that%voldati)) THEN
725 DEALLOCATE(voldatiin1, voldatiin2, voldatiout)
726ENDIF
727
728END SUBROUTINE volgrid6d_recompute_stat_proc_diff
729
730
731!> Specialized method for statistically processing a set of data
732!! by integration/differentiation.
733!! This method performs statistical processing by integrating
734!! (accumulating) in time values representing time-average rates or
735!! fluxes, (\a stat_proc_input=0 \a stat_proc=1) or by transforming a
736!! time-integrated (accumulated) value in a time-average rate or flux
737!! (\a stat_proc_input=1 \a stat_proc=0). Analysis/observation or
738!! forecast timeranges are processed. The only operation performed is
739!! respectively multiplying or dividing the values by the length of
740!! the time interval in seconds.
741!!
742!! The output \a that volgrid6d object contains elements from the
743!! original volume \a this satisfying the conditions
744!! - timerange (vol7d_timerange_class::vol7d_timerange::timerange)
745!! of type \a stat_proc_input (0 or 1)
746!! - any p1 (analysis/observation or forecast)
747!! - p2 &gt; 0 (processing interval non null, non instantaneous data)
748!!
749!! Output data will have timerange of type \a stat_proc (1 or 0) and
750!! p1 and p2 equal to the corresponding input values. The supported
751!! statistical processing methods (parameter \a stat_proc) are:
752!!
753!! - 0 average
754!! - 1 accumulation
755!!
756!! Input volume may have any value of \a this%time_definition, and
757!! that value will be conserved in the output volume.
758SUBROUTINE volgrid6d_compute_stat_proc_metamorph(this, that, stat_proc_input, stat_proc, clone)
759TYPE(volgrid6d),INTENT(inout) :: this !< volume providing data to be recomputed, it is not modified by the method, apart from performing a \a volgrid6d_alloc_vol on it
760TYPE(volgrid6d),INTENT(out) :: that !< output volume which will contain the recomputed data
761INTEGER,INTENT(in) :: stat_proc_input !< type of statistical processing of data that has to be processed (from grib2 table), only data having timerange of this type will be processed, the actual statistical processing performed and which will appear in the output volume, is however determined by \a stat_proc argument
762INTEGER,INTENT(in) :: stat_proc !< type of statistical processing to be recomputed (from grib2 table), data in output volume \a that will have a timerange of this type
763LOGICAL , INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a this to \a that
764
765INTEGER :: j, i3, i4, i6
766INTEGER,POINTER :: map_tr(:)
767REAL,POINTER :: voldatiin(:,:), voldatiout(:,:)
768REAL,ALLOCATABLE :: int_ratio(:)
769LOGICAL :: lclone
770
771NULLIFY(voldatiin, voldatiout)
772
773! be safe
774CALL volgrid6d_alloc_vol(this)
775! when volume is not decoded it is better to clone anyway to avoid
776! overwriting fields
777lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati)
778
779IF (.NOT.((stat_proc_input == 0 .AND. stat_proc == 1) .OR. &
780 (stat_proc_input == 1 .AND. stat_proc == 0))) THEN
781
782 CALL l4f_category_log(this%category, l4f_warn, &
783 'compute_stat_proc_metamorph, can only be applied to average->accumulated timerange and viceversa')
784! return an empty volume, without signaling error
785 CALL init(that)
786 CALL volgrid6d_alloc_vol(that)
787 RETURN
788ENDIF
789
790! initialise the output volume
791CALL init(that, griddim=this%griddim, time_definition=this%time_definition)
792CALL volgrid6d_alloc(that, dim=this%griddim%dim, ntime=SIZE(this%time), &
793 nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.)
794that%time = this%time
795that%level = this%level
796that%var = this%var
797
798CALL compute_stat_proc_metamorph_common(stat_proc_input, this%timerange, stat_proc, &
799 that%timerange, map_tr)
800
801! complete the definition of the output volume
802CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati))
803
804IF (stat_proc == 0) THEN ! average -> integral
805 int_ratio = 1./real(that%timerange(:)%p2)
806ELSE ! cumulation
807 int_ratio = real(that%timerange(:)%p2)
808ENDIF
809
810DO i6 = 1, SIZE(this%var)
811 DO j = 1, SIZE(map_tr)
812 DO i4 = 1, SIZE(that%time)
813 DO i3 = 1, SIZE(this%level)
814
815 IF (lclone) THEN
816 CALL copy(this%gaid(i3,i4,map_tr(j),i6), that%gaid(i3,i4,j,i6))
817 ELSE
818 that%gaid(i3,i4,map_tr(j),i6) = this%gaid(i3,i4,j,i6)
819 ENDIF
820 CALL volgrid_get_vol_2d(this, i3, i4, map_tr(j), i6, voldatiin)
821 CALL volgrid_get_vol_2d(that, i3, i4, j, i6, voldatiout)
822 WHERE (c_e(voldatiin))
823 voldatiout = voldatiin*int_ratio(j)
824 ELSEWHERE
825 voldatiout = rmiss
826 END WHERE
827 CALL volgrid_set_vol_2d(that, i3, i4, j, i6, voldatiout)
828 ENDDO
829 ENDDO
830 ENDDO
831ENDDO
832
833
834END SUBROUTINE volgrid6d_compute_stat_proc_metamorph
835
836!> Method for building a volume containing the vertical coordinate as
837!! a variable.
838!! This method produces a volgrid6d volume, derived from \a this,
839!! containing a single variable, horizontally constant, on the same
840!! input levels, which describes the vertical coordinate in the form
841!! of a physical variable. The grid, time and timerange metadata are
842!! the same as for the original volume. Only a single vertical level
843!! type, the one matching the \a level argument, is converted to a
844!! variable. The \a level argument can also indicate the layer between
845!! two surfaces of the same type, in that case the variable
846!! representing the vertical coordinate will be set to the value of
847!! the midpoint between the two layers. If something goes wrong,
848!! e.g. no level matches \a level argument or the level canot be
849!! converted to a physical value, an empty volume is returned.
850SUBROUTINE volgrid6d_compute_vert_coord_var(this, level, volgrid_lev)
851TYPE(volgrid6d),INTENT(in) :: this !< volume with the vertical levels
852TYPE(vol7d_level),INTENT(in) :: level !< vertical level to be converted to variable, only the type(s) of level are used not the value(s)
853TYPE(volgrid6d),INTENT(out) :: volgrid_lev !< output volume with the variable describing the vertical coordinate
854
855INTEGER :: nlev, i, ii, iii, iiii
856TYPE(grid_id) :: out_gaid
857LOGICAL,ALLOCATABLE :: levmask(:)
858TYPE(volgrid6d_var) :: lev_var
859
860CALL init(volgrid_lev) ! initialise to null
861IF (.NOT.ASSOCIATED(this%gaid)) THEN
862 CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: input volume not allocated')
863 RETURN
864ENDIF
865! if layer, both surfaces must be of the same type
866IF (c_e(level%level2) .AND. level%level1 /= level%level2) THEN
867 CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: requested (mixed) layer type not valid')
868 RETURN
869ENDIF
870
871! look for valid levels to be converted to vars
872ALLOCATE(levmask(SIZE(this%level)))
873levmask = this%level%level1 == level%level1 .AND. &
874 this%level%level2 == level%level2 .AND. c_e(this%level%l1)
875IF (c_e(level%level2)) levmask = levmask .AND. c_e(this%level%l2)
876nlev = count(levmask)
877IF (nlev == 0) THEN
878 CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: requested level type not available')
879 RETURN
880ENDIF
881
882out_gaid = grid_id_new()
883gaidloop: DO i=1 ,SIZE(this%gaid,1)
884 DO ii=1 ,SIZE(this%gaid,2)
885 DO iii=1 ,SIZE(this%gaid,3)
886 DO iiii=1 ,SIZE(this%gaid,4)
887 IF (c_e(this%gaid(i,ii,iii,iiii))) THEN ! conserve first valid gaid
888 CALL copy(this%gaid(i,ii,iii,iiii), out_gaid)
889 EXIT gaidloop
890 ENDIF
891 ENDDO
892 ENDDO
893 ENDDO
894ENDDO gaidloop
895
896! look for variable corresponding to level
897lev_var = convert(vol7d_var_new(btable=vol7d_level_to_var(level)), &
898 grid_id_template=out_gaid)
899IF (.NOT.c_e(lev_var)) THEN
900 CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: no variable corresponds to requested level type')
901 RETURN
902ENDIF
903
904! prepare output volume
905CALL init(volgrid_lev, griddim=this%griddim, &
906 time_definition=this%time_definition) !, categoryappend=categoryappend)
907CALL volgrid6d_alloc(volgrid_lev, ntime=SIZE(this%time), nlevel=nlev, &
908 ntimerange=SIZE(this%timerange), nvar=1)
909! fill metadata
910volgrid_lev%time = this%time
911volgrid_lev%level = pack(this%level, mask=levmask)
912volgrid_lev%timerange = this%timerange
913volgrid_lev%var(1) = lev_var
914
915CALL volgrid6d_alloc_vol(volgrid_lev, decode=.true.)
916! fill data
917DO i = 1, nlev
918 IF (c_e(level%level2)) THEN
919 volgrid_lev%voldati(:,:,i,:,:,:) = real(volgrid_lev%level(i)%l1 + &
920 volgrid_lev%level(i)%l2)* &
921 vol7d_level_to_var_factor(volgrid_lev%level(i))/2.
922 ELSE
923 volgrid_lev%voldati(:,:,i,:,:,:) = real(volgrid_lev%level(i)%l1)* &
924 vol7d_level_to_var_factor(volgrid_lev%level(i))
925 ENDIF
926ENDDO
927! fill gaid for subsequent export
928IF (c_e(out_gaid)) THEN
929 DO i=1 ,SIZE(volgrid_lev%gaid,1)
930 DO ii=1 ,SIZE(volgrid_lev%gaid,2)
931 DO iii=1 ,SIZE(volgrid_lev%gaid,3)
932 DO iiii=1 ,SIZE(volgrid_lev%gaid,4)
933 CALL copy(out_gaid, volgrid_lev%gaid(i,ii,iii,iiii))
934 ENDDO
935 ENDDO
936 ENDDO
937 ENDDO
938 CALL delete(out_gaid)
939ENDIF
940
941END SUBROUTINE volgrid6d_compute_vert_coord_var
942
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.
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.