volgrid6d_compute_stat_proc volgrid6d_compute_stat_proc_agg volgrid6d_compute_stat_proc_metamorph volgrid6d_compute_vert_coord_var volgrid6d_recompute_stat_proc_agg volgrid6d_recompute_stat_proc_diff

volgrid6d_compute_vert_coord_var() subroutine volgrid6d_compute_vert_coord_var ( type(volgrid6d), intent(in) this, type(vol7d_level), intent(in) level, type(volgrid6d), intent(out) volgrid_lev ) Method for building a volume containing the vertical coordinate as a variable. This method produces a volgrid6d volume, derived from this, containing a single variable, horizontally constant, on the same input levels, which describes the vertical coordinate in the form of a physical variable. The grid, time and timerange metadata are the same as for the original volume. Only a single vertical level type, the one matching the level argument, is converted to a variable. The level argument can also indicate the layer between two surfaces of the same type, in that case the variable representing the vertical coordinate will be set to the value of the midpoint between the two layers. If something goes wrong, e.g. no level matches level argument or the level canot be converted to a physical value, an empty volume is returned. Parametri [in] this volume with the vertical levels [in] level vertical level to be converted to variable, only the type(s) of level are used not the value(s) [out] volgrid_lev output volume with the variable describing the vertical coordinate Definizione alla linea 1041 del file volgrid6d_class_compute.F90. ! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it> ! authors: ! Davide Cesari <dcesari@arpa.emr.it> ! Paolo Patruno <ppatruno@arpa.emr.it> ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! You should have received a copy of the GNU General Public License ! along with this program. If not, see <http://www.gnu.org/licenses/>. #include "config.h" !> Extension of volgrid6d_class with methods for performing simple !! statistical operations on entire volumes of data. This module !! includes methods for performing operations such as (re)averaging !! and (re)cumulation in time of entire volgrid6d_class::volgrid6d !! objects. !! !! \ingroup volgrid6d MODULE volgrid6d_class_compute USE datetime_class USE volgrid6d_var_class USE volgrid6d_class USE grid_id_class USE stat_proc_engine USE simple_stat IMPLICIT NONE CONTAINS !> General-purpose method for computing a statistical processing on !! data in a volgrid6d object already processed with the same statistical !! processing, on a different time interval specified by \a step and !! \a start. This method tries to apply all the suitable specialized !! statistical processing methods according to the input and output !! statistical processing requested. The argument \a stat_proc_input !! determines which data will be processed, while the \a stat_proc !! argument determines the type of statistical process to be applied !! and which will be owned by output data. !! !! The possible combinations are: !! !! - \a stat_proc_input = 254 !! - \a stat_proc = 0 average instantaneous observations !! - \a stat_proc = 2 compute maximum of instantaneous observations !! - \a stat_proc = 3 compute minimum of instantaneous observations !! - \a stat_proc = 4 compute difference of instantaneous observations !! !! processing is computed on longer time intervals by aggregation, !! see the description of volgrid6d_compute_stat_proc_agg() !! !! - \a stat_proc_input = * !! - \a stat_proc = 254 consider statistically processed values as !! instantaneous without any extra processing !! !! see the description of volgrid6d_decompute_stat_proc() !! !! - \a stat_proc_input = 0, 1, 2, 3, 4, 200 !! - \a stat_proc = \a stat_proc_input recompute input data on !! different intervals !! !! the same statistical processing is applied to obtain data !! processed on a different interval, either longer, by !! aggregation, or shorter, by differences, see the description of !! volgrid6d_recompute_stat_proc_agg() and !! volgrid6d_recompute_stat_proc_diff() respectively; it is also !! possible to provide \a stat_proc_input \a /= \a stat_proc, but !! it has to be used with care. !! !! - \a stat_proc_input = 0 !! - \a stat_proc = 1 !! !! a time-averaged rate or flux is transformed into a !! time-integrated value (sometimes called accumulated) on the same !! interval by multiplying the values by the length of the time !! interval in seconds, keeping constant all the rest, including !! the variable; the unit of the variable implicitly changes !! accordingly, this is supported officially in grib2 standard, in !! the other cases it is a forcing of the standards. !! !! - \a stat_proc_input = 1 !! - \a stat_proc = 0 !! !! a time-integrated value (sometimes called accumulated) is !! transformed into a time-averaged rate or flux on the same !! interval by dividing the values by the length of the time !! interval in seconds, see also the previous description of the !! opposite computation. !! !! If a particular statistical processing cannot be performed on the !! input data, the program continues with a warning and, if requested, !! the input data is passed over to the volume specified by the \a !! other argument, in order to allow continuation of processing. All !! the other parameters are passed over to the specifical statistical !! processing methods and are documented there. SUBROUTINE volgrid6d_compute_stat_proc(this, that, stat_proc_input, stat_proc, & step, start, full_steps, frac_valid, max_step, weighted, clone) TYPE(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 TYPE(volgrid6d),INTENT(out) :: that !< output volume which will contain the recomputed data INTEGER,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 INTEGER,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 TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed TYPE(datetime),INTENT(in),OPTIONAL :: start !< start of statistical processing interval LOGICAL,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 REAL,INTENT(in),OPTIONAL :: frac_valid !< minimum fraction of valid data required for considering acceptable a recomputed value, default=1. TYPE(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 LOGICAL,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 LOGICAL , INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a this to \a that INTEGER :: dtmax, dtstep IF (stat_proc_input == 254) THEN CALL l4f_category_log(this%category, l4f_info, & 'computing statistical processing by aggregation '//& trim(to_char(stat_proc_input))//':'//trim(to_char(stat_proc))) CALL volgrid6d_compute_stat_proc_agg(this, that, stat_proc, & step, start, full_steps, max_step, clone) ELSE IF (stat_proc == 254) THEN CALL l4f_category_log(this%category, l4f_error, & 'statistical processing to instantaneous data not implemented for gridded fields') CALL raise_error() ELSE IF (stat_proc_input == stat_proc .OR. & (stat_proc == 0 .OR. stat_proc == 2 .OR. stat_proc == 3)) THEN ! avg, min and max can be computed from any input, with care IF (count(this%timerange(:)%timerange == stat_proc_input) == 0) THEN CALL l4f_category_log(this%category, l4f_warn, & 'no timeranges of the desired statistical processing type '//t2c(stat_proc)//' available') ! return an empty volume, without signaling error CALL init(that) CALL volgrid6d_alloc_vol(that) ELSE ! euristically determine whether aggregation or difference is more suitable dtmax = maxval(this%timerange(:)%p2, & mask=(this%timerange(:)%timerange == stat_proc)) CALL getval(step, asec=dtstep) #ifdef DEBUG CALL l4f_category_log(this%category, l4f_debug, & 'stat_proc='//t2c(stat_proc)//' dtmax='//t2c(dtmax)//' dtstep='//t2c(dtstep)) #endif IF (dtstep <= dtmax) THEN CALL l4f_category_log(this%category, l4f_info, & 'recomputing statistically processed data by difference '// & t2c(stat_proc_input)//':'//t2c(stat_proc)) CALL volgrid6d_recompute_stat_proc_diff(this, that, stat_proc, step, & full_steps, start, clone) ELSE CALL l4f_category_log(this%category, l4f_info, & 'recomputing statistically processed data by aggregation '// & t2c(stat_proc_input)//':'//t2c(stat_proc)) CALL volgrid6d_recompute_stat_proc_agg(this, that, stat_proc, step, start, & full_steps, frac_valid, clone, stat_proc_input) ENDIF ENDIF ELSE ! IF (stat_proc_input /= stat_proc) THEN IF ((stat_proc_input == 0 .AND. stat_proc == 1) .OR. & (stat_proc_input == 1 .AND. stat_proc == 0)) THEN CALL l4f_category_log(this%category, l4f_info, & 'computing statistically processed data by integration/differentiation '// & t2c(stat_proc_input)//':'//t2c(stat_proc)) CALL volgrid6d_compute_stat_proc_metamorph(this, that, stat_proc_input, & stat_proc, clone) ELSE CALL l4f_category_log(this%category, l4f_error, & 'statistical processing '//t2c(stat_proc_input)//':'//t2c(stat_proc)// & ' not implemented or does not make sense') CALL raise_error() ENDIF ENDIF END SUBROUTINE volgrid6d_compute_stat_proc !> Specialized method for statistically processing a set of data !! already processed with the same statistical processing, on a !! different time interval. This method performs statistical !! processing by aggregation of shorter intervals. !! !! The output \a that volgrid6d object contains elements from the original volume !! \a this satisfying the conditions !! - timerange (vol7d_timerange_class::vol7d_timerange::timerange) !! of type \a stat_proc (or \a stat_proc_input if provided) !! - any p1 (analysis/observation or forecast) !! - p2 > 0 (processing interval non null, non instantaneous data) !! and equal to a multiplier of \a step !! !! Output data will have timerange of type \a stat_proc and p2 = \a !! step. The supported statistical processing methods (parameter \a !! stat_proc) are: !! !! - 0 average !! - 1 accumulation !! - 2 maximum !! - 3 minimum !! - 4 difference !! - 200 vectorial mean !! !! The start of processing period can be computed automatically from !! the input intervals as the first possible interval modulo \a step, !! or, for a better control, it can be specified explicitely by the !! optional argument \a start. Be warned that, in the final volume, !! the first reference time will actually be \a start \a + \a step, !! since \a start indicates the beginning of first processing !! interval, while reference time (for analysis/oservation) is the end !! of the interval. !! !! The purpose of the optional argument \a stat_proc_input is to allow !! processing with a certain statistical processing operator a dataset !! already processed with a different operator, by specifying the !! latter as stat_proc_input; this is useful, for example, if one !! wants to compute the monthly average of daily maximum temperatures; !! however this has to be used with care since the resulting data !! volume will not carry all the information about the processing !! which has been done, in the previous case, for example, the !! temperatures will simply look like monthly average temperatures. SUBROUTINE volgrid6d_recompute_stat_proc_agg(this, that, stat_proc, & step, start, full_steps, frac_valid, clone, stat_proc_input) TYPE(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 TYPE(volgrid6d),INTENT(out) :: that !< output volume which will contain the recomputed data 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 TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed TYPE(datetime),INTENT(in),OPTIONAL :: start !< start of statistical processing interval LOGICAL,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 REAL,INTENT(in),OPTIONAL :: frac_valid !< minimum fraction of valid data required for considering acceptable a recomputed value, default=1. LOGICAL, INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a this to \a that 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 INTEGER :: tri INTEGER i, j, n, n1, ndtr, i3, i6 TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:) INTEGER,POINTER :: dtratio(:) REAL :: lfrac_valid LOGICAL :: lclone REAL,POINTER :: voldatiin(:,:), voldatiout(:,:) NULLIFY(voldatiin, voldatiout) IF (PRESENT(stat_proc_input)) THEN tri = stat_proc_input ELSE tri = stat_proc ENDIF IF (PRESENT(frac_valid)) THEN lfrac_valid = frac_valid ELSE lfrac_valid = 1.0 ENDIF CALL init(that) ! be safe CALL volgrid6d_alloc_vol(this) ! when volume is not decoded it is better to clone anyway to avoid ! overwriting fields lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati) ! initialise the output volume CALL init(that, griddim=this%griddim, time_definition=this%time_definition) CALL volgrid6d_alloc(that, dim=this%griddim%dim, ntimerange=1, & nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.) that%level = this%level that%var = this%var CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, & step, this%time_definition, that%time, that%timerange, map_ttr, & dtratio=dtratio, start=start, full_steps=full_steps) CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati)) do_otimerange: DO j = 1, SIZE(that%timerange) do_otime: DO i = 1, SIZE(that%time) DO n1 = 1, SIZE(dtratio) IF (dtratio(n1) <= 0) cycle ! safety check DO i6 = 1, SIZE(this%var) DO i3 = 1, SIZE(this%level) CALL volgrid_get_vol_2d(that, i3, i, j, i6, voldatiout) ndtr = 0 DO n = 1, map_ttr(i,j)%arraysize IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN ndtr = ndtr + 1 CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(n)%it, & map_ttr(i,j)%array(n)%itr, i6, voldatiin) IF (ndtr == 1) THEN voldatiout = voldatiin IF (lclone) THEN CALL copy(this%gaid(i3, map_ttr(i,j)%array(n)%it,& map_ttr(i,j)%array(n)%itr,i6), that%gaid(i3,i,j,i6)) ELSE that%gaid(i3,i,j,i6) = this%gaid(i3, map_ttr(i,j)%array(n)%it, & map_ttr(i,j)%array(n)%itr,i6) ENDIF ELSE ! second or more time SELECT CASE(stat_proc) CASE (0, 200, 1, 4) ! average, vectorial mean, accumulation, difference WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:))) voldatiout(:,:) = voldatiout(:,:) + voldatiin(:,:) ELSEWHERE voldatiout(:,:) = rmiss END WHERE CASE(2) ! maximum WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:))) voldatiout(:,:) = max(voldatiout(:,:), voldatiin(:,:)) ELSEWHERE voldatiout(:,:) = rmiss END WHERE CASE(3) ! minimum WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:))) voldatiout(:,:) = min(voldatiout(:,:), voldatiin(:,:)) ELSEWHERE voldatiout(:,:) = rmiss END WHERE END SELECT ENDIF ! first time ENDIF ! dtratio(n1) ENDDO ! ttr #ifdef DEBUG CALL l4f_log(l4f_debug, & 'compute_stat_proc_agg, ndtr/dtratio/frac_valid: '// & t2c(ndtr)//'/'//t2c(dtratio(n1))//'/'//t2c(lfrac_valid)) #endif IF (ndtr > 0) THEN ! why this condition was not here before? IF (real(ndtr)/real(dtratio(n1)) >= lfrac_valid) THEN ! success IF (stat_proc == 0) THEN ! average WHERE(c_e(voldatiout(:,:))) voldatiout(:,:) = voldatiout(:,:)/ndtr END WHERE ENDIF CALL volgrid_set_vol_2d(that, i3, i, j, i6, voldatiout) #ifdef DEBUG CALL l4f_log(l4f_debug, & 'compute_stat_proc_agg, coding lev/t/tr/var: '// & t2c(i3)//'/'//t2c(i)//'/'//t2c(j)//'/'//t2c(i6)) #endif ELSE ! must nullify the output gaid here, otherwise an incomplete field will be output IF (lclone) THEN CALL delete(that%gaid(i3,i,j,i6)) ELSE CALL init(that%gaid(i3,i,j,i6)) ! grid_id lacks a nullify method ENDIF #ifdef DEBUG CALL l4f_log(l4f_debug, & 'compute_stat_proc_agg, skipping lev/t/tr/var: '// & t2c(i3)//'/'//t2c(i)//'/'//t2c(j)//'/'//t2c(i6)) #endif ENDIF ENDIF ! ndtr > 0 ENDDO ! level ENDDO ! var ENDDO ! dtratio CALL delete(map_ttr(i,j)) ENDDO do_otime ENDDO do_otimerange DEALLOCATE(dtratio, map_ttr) END SUBROUTINE volgrid6d_recompute_stat_proc_agg !> Method for statistically processing a set of instantaneous data. !! This method performs statistical processing by aggregation of !! instantaneous data. !! !! The output \a that volgrid6d object contains elements from the original volume !! \a this satisfying the conditions !! - timerange (vol7d_timerange_class::vol7d_timerange::timerange) !! of type 254 (instantaeous) !! - any p1 (analysis/observation or forecast) !! - p2 = 0 (processing interval null, instantaneous data) !! !! Output data will have timerange of type \a stat_proc, and p2 = \a !! step. The supported statistical processing methods (parameter \a !! stat_proc) are: !! !! - 0 average !! - 2 maximum !! - 3 minimum !! - 4 difference !! !! A maximum distance in time for input valid data can be assigned !! with the optional argument \a max_step, in order to filter datasets !! with too long "holes". SUBROUTINE volgrid6d_compute_stat_proc_agg(this, that, stat_proc, & step, start, full_steps, max_step, clone) TYPE(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 TYPE(volgrid6d),INTENT(out) :: that !< output volume which will contain the recomputed data 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 TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed TYPE(datetime),INTENT(in),OPTIONAL :: start !< start of statistical processing interval LOGICAL,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 TYPE(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 LOGICAL , INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a this to \a that INTEGER :: tri INTEGER i, j, n, ninp, i3, i6 TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:) TYPE(timedelta) :: lmax_step LOGICAL :: lclone REAL,POINTER :: voldatiin(:,:), voldatiout(:,:) NULLIFY(voldatiin, voldatiout) tri = 254 IF (PRESENT(max_step)) THEN lmax_step = max_step ELSE lmax_step = timedelta_max ENDIF CALL init(that) ! be safe CALL volgrid6d_alloc_vol(this) ! when volume is not decoded it is better to clone anyway to avoid ! overwriting fields lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati) ! initialise the output volume CALL init(that, griddim=this%griddim, time_definition=this%time_definition) CALL volgrid6d_alloc(that, dim=this%griddim%dim, ntimerange=1, & nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.) that%level = this%level that%var = this%var CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, & step, this%time_definition, that%time, that%timerange, map_ttr, & start=start, full_steps=full_steps) CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati)) do_otimerange: DO j = 1, SIZE(that%timerange) do_otime: DO i = 1, SIZE(that%time) ninp = map_ttr(i,j)%arraysize IF (ninp <= 0) cycle do_otime IF (stat_proc == 4) THEN ! check validity for difference IF (map_ttr(i,j)%array(1)%extra_info /= 1 .OR. & map_ttr(i,j)%array(ninp)%extra_info /= 2) THEN CALL delete(map_ttr(i,j)) cycle do_otime ENDIF ELSE ! check validity condition (missing values in volume are not accounted for) DO n = 2, ninp IF (map_ttr(i,j)%array(n)%time - map_ttr(i,j)%array(n-1)%time > & lmax_step) THEN CALL delete(map_ttr(i,j)) cycle do_otime ENDIF ENDDO ENDIF DO i6 = 1, SIZE(this%var) DO i3 = 1, SIZE(this%level) CALL volgrid_get_vol_2d(that, i3, i, j, i6, voldatiout) IF (stat_proc == 4) THEN ! special treatment for difference IF (lclone) THEN CALL copy(this%gaid(i3, map_ttr(i,j)%array(1)%it,& map_ttr(i,j)%array(1)%itr,i6), that%gaid(i3,i,j,i6)) ELSE that%gaid(i3,i,j,i6) = this%gaid(i3, map_ttr(i,j)%array(1)%it, & map_ttr(i,j)%array(1)%itr,i6) ENDIF ! improve the next workflow? CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(ninp)%it, & map_ttr(i,j)%array(ninp)%itr, i6, voldatiin) voldatiout = voldatiin CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(1)%it, & map_ttr(i,j)%array(1)%itr, i6, voldatiin) WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:))) voldatiout(:,:) = voldatiout(:,:) - voldatiin(:,:) ELSEWHERE voldatiout(:,:) = rmiss END WHERE ELSE ! other stat_proc DO n = 1, ninp CALL volgrid_get_vol_2d(this, i3, map_ttr(i,j)%array(n)%it, & map_ttr(i,j)%array(n)%itr, i6, voldatiin) IF (n == 1) THEN voldatiout = voldatiin IF (lclone) THEN CALL copy(this%gaid(i3, map_ttr(i,j)%array(n)%it,& map_ttr(i,j)%array(n)%itr,i6), that%gaid(i3,i,j,i6)) ELSE that%gaid(i3,i,j,i6) = this%gaid(i3, map_ttr(i,j)%array(n)%it, & map_ttr(i,j)%array(n)%itr,i6) ENDIF ELSE ! second or more time SELECT CASE(stat_proc) CASE (0, 1) ! average, accumulation WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:))) voldatiout(:,:) = voldatiout(:,:) + voldatiin(:,:) ELSEWHERE voldatiout(:,:) = rmiss END WHERE CASE(2) ! maximum WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:))) voldatiout(:,:) = max(voldatiout(:,:), voldatiin(:,:)) ELSEWHERE voldatiout(:,:) = rmiss END WHERE CASE(3) ! minimum WHERE(c_e(voldatiin(:,:)) .AND. c_e(voldatiout(:,:))) voldatiout(:,:) = min(voldatiout(:,:), voldatiin(:,:)) ELSEWHERE voldatiout(:,:) = rmiss END WHERE END SELECT ENDIF ! first time ENDDO IF (stat_proc == 0) THEN ! average WHERE(c_e(voldatiout(:,:))) voldatiout(:,:) = voldatiout(:,:)/ninp END WHERE ENDIF ENDIF CALL volgrid_set_vol_2d(that, i3, i, j, i6, voldatiout) ENDDO ! level ENDDO ! var CALL delete(map_ttr(i,j)) ENDDO do_otime ENDDO do_otimerange DEALLOCATE(map_ttr) END SUBROUTINE volgrid6d_compute_stat_proc_agg !> Specialized method for statistically processing a set of data !! already processed with the same statistical processing, on a !! different time interval. This method performs statistical !! processing by difference of different intervals. Data with both !! analysis/observation or forecast timerange are processed. !! !! The output \a that volgrid6d object contains elements from the original volume !! \a this satisfying the conditions !! - timerange (vol7d_timerange_class::vol7d_timerange::timerange) !! of type \a stat_proc !! - any p1 (analysis/observation or forecast) !! - p2 &gt; 0 (processing interval non null, non instantaneous data) !! and equal to a multiplier of \a step if \a full_steps is \c .TRUE. !! !! Output data will have timerange of type \a stat_proc and p2 = \a !! step. The supported statistical processing methods (parameter \a !! stat_proc) are: !! !! - 0 average !! - 1 accumulation !! - 4 difference !! !! Input volume may have any value of \a this%time_definition, and !! that value will be conserved in the output volume. SUBROUTINE volgrid6d_recompute_stat_proc_diff(this, that, stat_proc, step, full_steps, start, clone) TYPE(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 TYPE(volgrid6d),INTENT(out) :: that !< output volume which will contain the recomputed data 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 TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed LOGICAL,INTENT(in),OPTIONAL :: full_steps !< if provided and \a .TRUE., process only data having processing interval (p2) equal to a multiple of \a step TYPE(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 LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a this to \a that INTEGER :: i3, i4, i6, i, j, k, l, nitr, steps INTEGER,ALLOCATABLE :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:) REAL,POINTER :: voldatiin1(:,:), voldatiin2(:,:), voldatiout(:,:) !LOGICAL,POINTER :: mask_timerange(:) LOGICAL :: lclone TYPE(vol7d_var),ALLOCATABLE :: varbufr(:) ! be safe CALL volgrid6d_alloc_vol(this) ! when volume is not decoded it is better to clone anyway to avoid ! overwriting fields lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati) ! initialise the output volume CALL init(that, griddim=this%griddim, time_definition=this%time_definition) CALL volgrid6d_alloc(that, dim=this%griddim%dim, & nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.) that%level = this%level that%var = this%var ! compute length of cumulation step in seconds CALL getval(step, asec=steps) ! compute the statistical processing relations, output time and ! timerange are defined here CALL recompute_stat_proc_diff_common(this%time, this%timerange, stat_proc, step, & that%time, that%timerange, map_tr, f, keep_tr, & this%time_definition, full_steps, start) nitr = SIZE(f) ! complete the definition of the output volume CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati)) ! allocate workspace once IF (.NOT.ASSOCIATED(that%voldati)) THEN ALLOCATE(voldatiin1(this%griddim%dim%nx, this%griddim%dim%ny), & voldatiin2(this%griddim%dim%nx, this%griddim%dim%ny), & voldatiout(this%griddim%dim%nx, this%griddim%dim%ny)) ENDIF ! copy the timeranges already satisfying the requested step, if any DO i4 = 1, SIZE(this%time) DO i = 1, nitr IF (c_e(keep_tr(i, i4, 2))) THEN l = keep_tr(i, i4, 1) k = keep_tr(i, i4, 2) #ifdef DEBUG CALL l4f_category_log(this%category, l4f_debug, & 'volgrid6d_recompute_stat_proc_diff, good timerange: '//t2c(f(i))// & '->'//t2c(k)) #endif DO i6 = 1, SIZE(this%var) DO i3 = 1, SIZE(this%level) IF (c_e(this%gaid(i3,i4,f(i),i6))) THEN IF (lclone) THEN CALL copy(this%gaid(i3,i4,f(i),i6), that%gaid(i3,l,k,i6)) ELSE that%gaid(i3,l,k,i6) = this%gaid(i3,i4,f(i),i6) ENDIF IF (ASSOCIATED(that%voldati)) THEN that%voldati(:,:,i3,l,k,i6) = this%voldati(:,:,i3,i4,f(i),i6) ELSE CALL volgrid_get_vol_2d(this, i3, i4, f(i), i6, voldatiout) CALL volgrid_set_vol_2d(that, i3, l, k, i6, voldatiout) ENDIF ENDIF ENDDO ENDDO ENDIF ENDDO ENDDO ! varbufr required for setting posdef, optimize with an array IF (stat_proc == 1) THEN ALLOCATE(varbufr(SIZE(this%var))) DO i6 = 1, SIZE(this%var) varbufr(i6) = convert(this%var(i6)) CALL vol7d_var_features_set(varbufr(i6)) ENDDO ENDIF ! compute statistical processing DO l = 1, SIZE(this%time) DO k = 1, nitr DO j = 1, SIZE(this%time) DO i = 1, nitr IF (c_e(map_tr(i,j,k,l,1))) THEN DO i6 = 1, SIZE(this%var) DO i3 = 1, SIZE(this%level) IF (c_e(this%gaid(i3,j,f(i),i6)) .AND. & c_e(this%gaid(i3,l,f(k),i6))) THEN ! take the gaid from the second time/timerange contributing to the ! result (l,f(k)) IF (lclone) THEN CALL copy(this%gaid(i3,l,f(k),i6), & that%gaid(i3,map_tr(i,j,k,l,1),map_tr(i,j,k,l,2),i6)) ELSE that%gaid(i3,map_tr(i,j,k,l,1),map_tr(i,j,k,l,2),i6) = & this%gaid(i3,l,f(k),i6) ENDIF ! get/set 2d sections API is used CALL volgrid_get_vol_2d(this, i3, l, f(k), i6, voldatiin1) CALL volgrid_get_vol_2d(this, i3, j, f(i), i6, voldatiin2) IF (ASSOCIATED(that%voldati)) & CALL volgrid_get_vol_2d(that, i3, & map_tr(i,j,k,l,1), map_tr(i,j,k,l,2), i6, voldatiout) IF (stat_proc == 0) THEN ! average WHERE(c_e(voldatiin1(:,:)) .AND. c_e(voldatiin2(:,:))) voldatiout(:,:) = & (voldatiin1(:,:)*this%timerange(f(k))%p2 - & voldatiin2(:,:)*this%timerange(f(i))%p2)/ & steps ELSEWHERE voldatiout(:,:) = rmiss END WHERE ELSE IF (stat_proc == 1 .OR. stat_proc == 4) THEN ! acc, diff WHERE(c_e(voldatiin1(:,:)) .AND. c_e(voldatiin2(:,:))) voldatiout(:,:) = voldatiin1(:,:) - voldatiin2(:,:) ELSEWHERE voldatiout(:,:) = rmiss END WHERE IF (stat_proc == 1) THEN CALL vol7d_var_features_posdef_apply(varbufr(i6), voldatiout) ENDIF ENDIF CALL volgrid_set_vol_2d(that, i3, & map_tr(i,j,k,l,1), map_tr(i,j,k,l,2), i6, voldatiout) ENDIF ENDDO ENDDO ENDIF ENDDO ENDDO ENDDO ENDDO IF (.NOT.ASSOCIATED(that%voldati)) THEN DEALLOCATE(voldatiin1, voldatiin2, voldatiout) ENDIF END SUBROUTINE volgrid6d_recompute_stat_proc_diff !> Specialized method for statistically processing a set of data !! by integration/differentiation. !! This method performs statistical processing by integrating !! (accumulating) in time values representing time-average rates or !! fluxes, (\a stat_proc_input=0 \a stat_proc=1) or by transforming a !! time-integrated (accumulated) value in a time-average rate or flux !! (\a stat_proc_input=1 \a stat_proc=0). Analysis/observation or !! forecast timeranges are processed. The only operation performed is !! respectively multiplying or dividing the values by the length of !! the time interval in seconds. !! !! The output \a that volgrid6d object contains elements from the !! original volume \a this satisfying the conditions !! - timerange (vol7d_timerange_class::vol7d_timerange::timerange) !! of type \a stat_proc_input (0 or 1) !! - any p1 (analysis/observation or forecast) !! - p2 &gt; 0 (processing interval non null, non instantaneous data) !! !! Output data will have timerange of type \a stat_proc (1 or 0) and !! p1 and p2 equal to the corresponding input values. The supported !! statistical processing methods (parameter \a stat_proc) are: !! !! - 0 average !! - 1 accumulation !! !! Input volume may have any value of \a this%time_definition, and !! that value will be conserved in the output volume. SUBROUTINE volgrid6d_compute_stat_proc_metamorph(this, that, stat_proc_input, stat_proc, clone) TYPE(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 TYPE(volgrid6d),INTENT(out) :: that !< output volume which will contain the recomputed data INTEGER,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 INTEGER,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 LOGICAL , INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a this to \a that INTEGER :: j, i3, i4, i6 INTEGER,POINTER :: map_tr(:) REAL,POINTER :: voldatiin(:,:), voldatiout(:,:) REAL,ALLOCATABLE :: int_ratio(:) LOGICAL :: lclone NULLIFY(voldatiin, voldatiout) ! be safe CALL volgrid6d_alloc_vol(this) ! when volume is not decoded it is better to clone anyway to avoid ! overwriting fields lclone = optio_log(clone) .OR. .NOT.ASSOCIATED(this%voldati) IF (.NOT.((stat_proc_input == 0 .AND. stat_proc == 1) .OR. & (stat_proc_input == 1 .AND. stat_proc == 0))) THEN CALL l4f_category_log(this%category, l4f_warn, & 'compute_stat_proc_metamorph, can only be applied to average->accumulated timerange and viceversa') ! return an empty volume, without signaling error CALL init(that) CALL volgrid6d_alloc_vol(that) RETURN ENDIF ! initialise the output volume CALL init(that, griddim=this%griddim, time_definition=this%time_definition) CALL volgrid6d_alloc(that, dim=this%griddim%dim, ntime=SIZE(this%time), & nlevel=SIZE(this%level), nvar=SIZE(this%var), ini=.false.) that%time = this%time that%level = this%level that%var = this%var CALL compute_stat_proc_metamorph_common(stat_proc_input, this%timerange, stat_proc, & that%timerange, map_tr) ! complete the definition of the output volume CALL volgrid6d_alloc_vol(that, decode=ASSOCIATED(this%voldati)) IF (stat_proc == 0) THEN ! average -> integral int_ratio = 1./real(that%timerange(:)%p2) ELSE ! cumulation int_ratio = real(that%timerange(:)%p2) ENDIF DO i6 = 1, SIZE(this%var) DO j = 1, SIZE(map_tr) DO i4 = 1, SIZE(that%time) DO i3 = 1, SIZE(this%level) IF (lclone) THEN CALL copy(this%gaid(i3,i4,map_tr(j),i6), that%gaid(i3,i4,j,i6)) ELSE that%gaid(i3,i4,map_tr(j),i6) = this%gaid(i3,i4,j,i6) ENDIF CALL volgrid_get_vol_2d(this, i3, i4, map_tr(j), i6, voldatiin) CALL volgrid_get_vol_2d(that, i3, i4, j, i6, voldatiout) WHERE (c_e(voldatiin)) voldatiout = voldatiin*int_ratio(j) ELSEWHERE voldatiout = rmiss END WHERE CALL volgrid_set_vol_2d(that, i3, i4, j, i6, voldatiout) ENDDO ENDDO ENDDO ENDDO END SUBROUTINE volgrid6d_compute_stat_proc_metamorph !> Method for building a volume containing the vertical coordinate as !! a variable. !! This method produces a volgrid6d volume, derived from \a this, !! containing a single variable, horizontally constant, on the same !! input levels, which describes the vertical coordinate in the form !! of a physical variable. The grid, time and timerange metadata are !! the same as for the original volume. Only a single vertical level !! type, the one matching the \a level argument, is converted to a !! variable. The \a level argument can also indicate the layer between !! two surfaces of the same type, in that case the variable !! representing the vertical coordinate will be set to the value of !! the midpoint between the two layers. If something goes wrong, !! e.g. no level matches \a level argument or the level canot be !! converted to a physical value, an empty volume is returned. SUBROUTINE volgrid6d_compute_vert_coord_var(this, level, volgrid_lev) TYPE(volgrid6d),INTENT(in) :: this !< volume with the vertical levels TYPE(vol7d_level),INTENT(in) :: level !< vertical level to be converted to variable, only the type(s) of level are used not the value(s) TYPE(volgrid6d),INTENT(out) :: volgrid_lev !< output volume with the variable describing the vertical coordinate INTEGER :: nlev, i, ii, iii, iiii TYPE(grid_id) :: out_gaid LOGICAL,ALLOCATABLE :: levmask(:) TYPE(volgrid6d_var) :: lev_var CALL init(volgrid_lev) ! initialise to null IF (.NOT.ASSOCIATED(this%gaid)) THEN CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: input volume not allocated') RETURN ENDIF ! if layer, both surfaces must be of the same type IF (c_e(level%level2) .AND. level%level1 /= level%level2) THEN CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: requested (mixed) layer type not valid') RETURN ENDIF ! look for valid levels to be converted to vars ALLOCATE(levmask(SIZE(this%level))) levmask = this%level%level1 == level%level1 .AND. & this%level%level2 == level%level2 .AND. c_e(this%level%l1) IF (c_e(level%level2)) levmask = levmask .AND. c_e(this%level%l2) nlev = count(levmask) IF (nlev == 0) THEN CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: requested level type not available') RETURN ENDIF out_gaid = grid_id_new() gaidloop: DO i=1 ,SIZE(this%gaid,1) DO ii=1 ,SIZE(this%gaid,2) DO iii=1 ,SIZE(this%gaid,3) DO iiii=1 ,SIZE(this%gaid,4) IF (c_e(this%gaid(i,ii,iii,iiii))) THEN ! conserve first valid gaid CALL copy(this%gaid(i,ii,iii,iiii), out_gaid) EXIT gaidloop ENDIF ENDDO ENDDO ENDDO ENDDO gaidloop ! look for variable corresponding to level lev_var = convert(vol7d_var_new(btable=vol7d_level_to_var(level)), & grid_id_template=out_gaid) IF (.NOT.c_e(lev_var)) THEN CALL l4f_log(l4f_error, 'volgrid6d_compute_vert_coord_var: no variable corresponds to requested level type') RETURN ENDIF ! prepare output volume CALL init(volgrid_lev, griddim=this%griddim, & time_definition=this%time_definition) !, categoryappend=categoryappend) CALL volgrid6d_alloc(volgrid_lev, ntime=SIZE(this%time), nlevel=nlev, & ntimerange=SIZE(this%timerange), nvar=1) ! fill metadata volgrid_lev%time = this%time volgrid_lev%level = pack(this%level, mask=levmask) volgrid_lev%timerange = this%timerange volgrid_lev%var(1) = lev_var CALL volgrid6d_alloc_vol(volgrid_lev, decode=.true.) ! fill data DO i = 1, nlev IF (c_e(level%level2)) THEN volgrid_lev%voldati(:,:,i,:,:,:) = real(volgrid_lev%level(i)%l1 + & volgrid_lev%level(i)%l2)* & vol7d_level_to_var_factor(volgrid_lev%level(i))/2. ELSE volgrid_lev%voldati(:,:,i,:,:,:) = real(volgrid_lev%level(i)%l1)* & vol7d_level_to_var_factor(volgrid_lev%level(i)) ENDIF ENDDO ! fill gaid for subsequent export IF (c_e(out_gaid)) THEN DO i=1 ,SIZE(volgrid_lev%gaid,1) DO ii=1 ,SIZE(volgrid_lev%gaid,2) DO iii=1 ,SIZE(volgrid_lev%gaid,3) DO iiii=1 ,SIZE(volgrid_lev%gaid,4) CALL copy(out_gaid, volgrid_lev%gaid(i,ii,iii,iiii)) ENDDO ENDDO ENDDO ENDDO CALL delete(out_gaid) ENDIF END SUBROUTINE volgrid6d_compute_vert_coord_var END MODULE volgrid6d_class_compute