arrayof_ttr_mapper_append arrayof_ttr_mapper_append_unique arrayof_ttr_mapper_delete arrayof_ttr_mapper_insert arrayof_ttr_mapper_insert_array arrayof_ttr_mapper_insert_sorted arrayof_ttr_mapper_insert_unique arrayof_ttr_mapper_packarray arrayof_ttr_mapper_remove count_distinct_sorted_ttr_mapper count_distinct_ttr_mapper index_sorted_ttr_mapper index_ttr_mapper inssor_ttr_mapper map_distinct_ttr_mapper map_inv_distinct_ttr_mapper pack_distinct_sorted_ttr_mapper pack_distinct_ttr_mapper sort_ttr_mapper

inssor_ttr_mapper() subroutine inssor_ttr_mapper ( type(ttr_mapper), dimension (:), intent(inout) xdont ) Sorts into increasing order (Insertion sort) Sorts XDONT into increasing order (Insertion sort) This subroutine uses insertion sort. It does not use any work array and is faster when XDONT is of very small size (< 20), or already almost sorted, so it is used in a final pass when the partial quicksorting has left a sequence of small subsets and that sorting is only necessary within each subset to complete the process. Michel Olagnon - Apr. 2000 Definizione alla linea 1626 del file stat_proc_engine.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" !> This module contains functions that are only for internal use of !! the library. It should not be used by user procedures because it is !! subject to change !! \ingroup vol7d MODULE stat_proc_engine USE datetime_class USE vol7d_class IMPLICIT NONE ! per ogni elemento i,j dell'output, definire n elementi di input ad ! esso contribuenti (arrayof_ttr_mapper) con le seguenti informazioni TYPE ttr_mapper INTEGER :: it=imiss ! k INTEGER :: itr=imiss ! l INTEGER :: extra_info=imiss ! dtratio for intervals, extreme for point in time TYPE(datetime) :: time=datetime_miss ! time for point in time END TYPE ttr_mapper INTERFACE OPERATOR (==) MODULE PROCEDURE ttr_mapper_eq END INTERFACE INTERFACE OPERATOR (/=) MODULE PROCEDURE ttr_mapper_ne END INTERFACE INTERFACE OPERATOR (>) MODULE PROCEDURE ttr_mapper_gt END INTERFACE INTERFACE OPERATOR (<) MODULE PROCEDURE ttr_mapper_lt END INTERFACE INTERFACE OPERATOR (>=) MODULE PROCEDURE ttr_mapper_ge END INTERFACE INTERFACE OPERATOR (<=) MODULE PROCEDURE ttr_mapper_le END INTERFACE #undef VOL7D_POLY_TYPE #undef VOL7D_POLY_TYPES #undef ENABLE_SORT #define VOL7D_POLY_TYPE TYPE(ttr_mapper) #define VOL7D_POLY_TYPES _ttr_mapper #define ENABLE_SORT #include "array_utilities_pre.F90" #define ARRAYOF_ORIGTYPE TYPE(ttr_mapper) #define ARRAYOF_TYPE arrayof_ttr_mapper #define ARRAYOF_ORIGEQ 1 #define ARRAYOF_ORIGGT 1 #include "arrayof_pre.F90" ! from arrayof TYPE,PRIVATE :: ttr_decoder TYPE(datetime) :: pstart=datetime_miss TYPE(datetime) :: pend=datetime_miss TYPE(datetime) :: reftime=datetime_miss LOGICAL :: isana=.false. END TYPE ttr_decoder CONTAINS ! simplified comparisons only on time ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res) TYPE(ttr_mapper),INTENT(IN) :: this, that LOGICAL :: res res = this%time == that%time END FUNCTION ttr_mapper_eq ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res) TYPE(ttr_mapper),INTENT(IN) :: this, that LOGICAL :: res res = this%time /= that%time END FUNCTION ttr_mapper_ne ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res) TYPE(ttr_mapper),INTENT(IN) :: this, that LOGICAL :: res res = this%time > that%time END FUNCTION ttr_mapper_gt ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res) TYPE(ttr_mapper),INTENT(IN) :: this, that LOGICAL :: res res = this%time < that%time END FUNCTION ttr_mapper_lt ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res) TYPE(ttr_mapper),INTENT(IN) :: this, that LOGICAL :: res res = this%time >= that%time END FUNCTION ttr_mapper_ge ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res) TYPE(ttr_mapper),INTENT(IN) :: this, that LOGICAL :: res res = this%time <= that%time END FUNCTION ttr_mapper_le #include "arrayof_post.F90" #include "array_utilities_inc.F90" ! common operations for statistical processing by differences SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, & otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, & start) TYPE(datetime),INTENT(in) :: itime(:) TYPE(vol7d_timerange),INTENT(in) :: itimerange(:) INTEGER,INTENT(in) :: stat_proc TYPE(timedelta),INTENT(in) :: step TYPE(datetime),POINTER :: otime(:) TYPE(vol7d_timerange),POINTER :: otimerange(:) INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:) INTEGER :: nitr LOGICAL,ALLOCATABLE :: mask_timerange(:) INTEGER,INTENT(in) :: time_definition LOGICAL,INTENT(in),OPTIONAL :: full_steps TYPE(datetime),INTENT(in),OPTIONAL :: start INTEGER :: i, j, k, l, dirtyrep INTEGER :: steps LOGICAL :: lfull_steps, useful TYPE(datetime) :: lstart, tmptime ! pstart1, pstart2, pend1, pend2, reftime1, reftime2 TYPE(vol7d_timerange) :: tmptimerange TYPE(arrayof_datetime) :: a_otime TYPE(arrayof_vol7d_timerange) :: a_otimerange TYPE(timedelta) :: start_delta TYPE(ttr_decoder),ALLOCATABLE :: ttr_d(:,:) ! compute length of cumulation step in seconds CALL getval(step, asec=steps) lstart = datetime_miss IF (PRESENT(start)) lstart = start lfull_steps = optio_log(full_steps) ! create a mask of suitable timeranges ALLOCATE(mask_timerange(SIZE(itimerange))) mask_timerange(:) = itimerange(:)%timerange == stat_proc & .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss & .AND. itimerange(:)%p1 >= 0 & .AND. itimerange(:)%p2 > 0 IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges mask_timerange(:) = mask_timerange(:) .AND. & (itimerange(:)%p1 == 0 .OR. & ! all analyses mod(itimerange(:)%p1, steps) == 0 .OR. & ! end time is mod step mod(itimerange(:)%p1 - itimerange(:)%p2, steps) == 0) ! start time is mod step ENDIF ! mask_timerange includes all candidate timeranges nitr = count(mask_timerange) ALLOCATE(f(nitr)) j = 1 DO i = 1, nitr DO WHILE(.NOT.mask_timerange(j)) j = j + 1 ENDDO f(i) = j j = j + 1 ENDDO ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2)) map_tr(:,:,:,:,:) = imiss ! precompute ttr ALLOCATE(ttr_d(nitr,SIZE(itime))) DO l = 1, SIZE(itime) DO k = 1, nitr CALL ttr_decoder_decode(itime(l), itimerange(f(k)), & time_definition, ttr_d(k,l)) IF (ttr_d(k,l)%reftime == ttr_d(k,l)%pend) ttr_d(k,l)%isana = .true. ENDDO ENDDO ! now we have to evaluate time/timerage pairs which do not need processing ALLOCATE(keep_tr(nitr, SIZE(itime), 2)) CALL compute_keep_tr() ! scan through all possible combinations of time and timerange DO dirtyrep = 1, 2 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors CALL packarray(a_otime) CALL packarray(a_otimerange) CALL sort(a_otime%array) CALL sort(a_otimerange%array) ENDIF DO l = 1, SIZE(itime) DO k = 1, nitr ! CALL time_timerange_get_period(itime(l), itimerange(f(k)), & ! time_definition, pstart2, pend2, reftime2) DO j = 1, SIZE(itime) DO i = 1, nitr useful = .false. ! CALL time_timerange_get_period(itime(j), itimerange(f(i)), & ! time_definition, pstart1, pend1, reftime1) ! tmptimerange = vol7d_timerange_new(timerange=stat_proc) ! IF (ttr_d(k,l)%reftime == ttr_d(k,l)%pend .AND. ttr_d(i,j)%reftime == ttr_d(i,j)%pend) THEN ! analysis IF (ttr_d(k,l)%isana .AND. ttr_d(i,j)%isana) THEN ! analysis IF (ttr_d(k,l)%pstart == ttr_d(i,j)%pstart .AND. ttr_d(k,l)%pend > ttr_d(i,j)%pend) THEN ! =-| CALL time_timerange_set_period(tmptime, tmptimerange, & time_definition, stat_proc, ttr_d(i,j)%pend, ttr_d(k,l)%pend, ttr_d(k,l)%reftime) IF (lfull_steps) THEN IF (mod(ttr_d(k,l)%reftime, step) == timedelta_0) THEN useful = .true. ENDIF ELSE useful = .true. ENDIF ! keep only data after lstart IF (c_e(lstart)) THEN IF (lstart > ttr_d(i,j)%pend) useful = .false. ENDIF ELSE IF (ttr_d(k,l)%pstart < ttr_d(i,j)%pstart .AND. ttr_d(k,l)%pend == ttr_d(i,j)%pend) THEN ! -=| CALL time_timerange_set_period(tmptime, tmptimerange, & time_definition, stat_proc, ttr_d(k,l)%pstart, ttr_d(i,j)%pstart, ttr_d(i,j)%pstart) IF (lfull_steps) THEN IF (mod(ttr_d(i,j)%pstart, step) == timedelta_0) THEN useful = .true. ENDIF ELSE useful = .true. ENDIF ! keep only data after lstart IF (c_e(lstart)) THEN IF (lstart > ttr_d(k,l)%pstart) useful = .false. ENDIF ENDIF ELSE IF (ttr_d(k,l)%reftime == ttr_d(i,j)%reftime) THEN ! forecast, same reftime IF (lfull_steps) THEN IF (c_e(lstart)) THEN ! lstart shifts the interval for computing modulo step, this does not ! remove data before lstart but just shifts the phase start_delta = lstart-ttr_d(k,l)%reftime ELSE start_delta = timedelta_0 ENDIF ENDIF IF (ttr_d(k,l)%pstart == ttr_d(i,j)%pstart .AND. ttr_d(k,l)%pend > ttr_d(i,j)%pend) THEN ! |=- CALL time_timerange_set_period(tmptime, tmptimerange, & time_definition, stat_proc, ttr_d(i,j)%pend, ttr_d(k,l)%pend, ttr_d(k,l)%reftime) IF (lfull_steps) THEN IF (mod(ttr_d(k,l)%pend-ttr_d(k,l)%reftime-start_delta, step) == timedelta_0) THEN useful = .true. ENDIF ELSE useful = .true. ENDIF ! keep only data after lstart IF (c_e(lstart)) THEN IF (lstart > ttr_d(i,j)%pend) useful = .false. ENDIF ELSE IF (ttr_d(k,l)%pstart < ttr_d(i,j)%pstart .AND. ttr_d(k,l)%pend == ttr_d(i,j)%pend) THEN ! |-= CALL time_timerange_set_period(tmptime, tmptimerange, & time_definition, stat_proc, ttr_d(k,l)%pstart, ttr_d(i,j)%pstart, ttr_d(k,l)%reftime) IF (lfull_steps) THEN IF (mod(ttr_d(i,j)%pstart-ttr_d(k,l)%reftime-start_delta, step) == timedelta_0) THEN useful = .true. ENDIF ELSE useful = .true. ENDIF ! keep only data after lstart IF (c_e(lstart)) THEN IF (lstart > ttr_d(k,l)%pstart) useful = .false. ENDIF ENDIF ENDIF useful = useful .AND. tmptime /= datetime_miss .AND. & tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps IF (useful) THEN ! add a_otime, a_otimerange map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime) map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange) ENDIF ENDDO ENDDO ENDDO ENDDO ENDDO ! we have to repeat the computation with sorted arrays CALL compute_keep_tr() otime => a_otime%array otimerange => a_otimerange%array ! delete local objects keeping the contents CALL delete(a_otime, nodealloc=.true.) CALL delete(a_otimerange, nodealloc=.true.) #ifdef DEBUG CALL l4f_log(l4f_debug, & 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr,1)))//', '// & t2c((SIZE(map_tr,2)))//', '// & t2c((SIZE(map_tr,3)))//', '// & t2c((SIZE(map_tr,4)))) CALL l4f_log(l4f_debug, & 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr))/2)//', '// & t2c(count(c_e(map_tr))/2)) CALL l4f_log(l4f_debug, & 'recompute_stat_proc_diff, nitr: '//t2c(nitr)) CALL l4f_log(l4f_debug, & 'recompute_stat_proc_diff, good timeranges: '//t2c(count(c_e(keep_tr))/2)) CALL l4f_log(l4f_debug, & 'recompute_stat_proc_diff, output times: '//t2c(SIZE(otime))) CALL l4f_log(l4f_debug, & 'recompute_stat_proc_diff, output timeranges: '//t2c(SIZE(otimerange))) #endif CONTAINS SUBROUTINE compute_keep_tr() INTEGER :: start_deltas keep_tr(:,:,:) = imiss DO l = 1, SIZE(itime) itrloop: DO k = 1, nitr IF (itimerange(f(k))%p2 == steps) THEN ! CALL time_timerange_get_period(itime(l), itimerange(f(k)), & ! time_definition, pstart2, pend2, reftime2) useful = .false. ! keep only data after lstart IF (c_e(lstart)) THEN IF (lstart > ttr_d(k,l)%pstart) cycle itrloop ENDIF IF (ttr_d(k,l)%reftime == ttr_d(k,l)%pend) THEN ! analysis IF (c_e(lstart)) THEN ! in analysis mode start wins over full_steps IF (mod(ttr_d(k,l)%reftime-lstart, step) == timedelta_0) THEN useful = .true. ENDIF ELSE IF (lfull_steps) THEN IF (mod(ttr_d(k,l)%reftime, step) == timedelta_0) THEN useful = .true. ENDIF ELSE useful = .true. ENDIF ELSE ! forecast IF (lfull_steps) THEN ! same as above for start_delta, but in seconds and not in timerange/timedelta units IF (c_e(lstart)) THEN start_deltas = timedelta_getamsec(lstart-ttr_d(k,l)%reftime)/1000_int_ll ELSE start_deltas = 0 ENDIF IF (mod(itimerange(f(k))%p1 - start_deltas, steps) == 0) THEN useful = .true. ENDIF ELSE useful = .true. ENDIF ENDIF IF (useful) THEN ! CALL time_timerange_set_period(tmptime, tmptimerange, & ! time_definition, stat_proc, pstart2, pend2, reftime2) keep_tr(k,l,1) = append_unique(a_otime, itime(l)) keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k))) ENDIF ENDIF ENDDO itrloop ENDDO END SUBROUTINE compute_keep_tr END SUBROUTINE recompute_stat_proc_diff_common ! common operations for statistical processing by metamorphosis SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, & otimerange, map_tr) INTEGER,INTENT(in) :: istat_proc TYPE(vol7d_timerange),INTENT(in) :: itimerange(:) INTEGER,INTENT(in) :: ostat_proc TYPE(vol7d_timerange),POINTER :: otimerange(:) INTEGER,POINTER :: map_tr(:) INTEGER :: i LOGICAL :: tr_mask(SIZE(itimerange)) IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume ALLOCATE(otimerange(0), map_tr(0)) RETURN ENDIF ! useful timeranges tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss & .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0 ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask))) otimerange = pack(itimerange, mask=tr_mask) otimerange(:)%timerange = ostat_proc map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask) END SUBROUTINE compute_stat_proc_metamorph_common ! common operations for statistical processing by aggregation SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, & step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps) TYPE(datetime),INTENT(in) :: itime(:) TYPE(vol7d_timerange),INTENT(in) :: itimerange(:) INTEGER,INTENT(in) :: stat_proc INTEGER,INTENT(in) :: tri TYPE(timedelta),INTENT(in) :: step INTEGER,INTENT(in) :: time_definition TYPE(datetime),POINTER :: otime(:) TYPE(vol7d_timerange),POINTER :: otimerange(:) TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:) INTEGER,POINTER,OPTIONAL :: dtratio(:) TYPE(datetime),INTENT(in),OPTIONAL :: start LOGICAL,INTENT(in),OPTIONAL :: full_steps INTEGER :: i, j, k, l, na, nf, n INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps INTEGER(kind=int_ll) :: stepms, mstepms LOGICAL :: lforecast TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime TYPE(arrayof_datetime) :: a_otime TYPE(arrayof_vol7d_timerange) :: a_otimerange TYPE(arrayof_integer) :: a_dtratio LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!! TYPE(ttr_mapper) :: lmapper CHARACTER(len=8) :: env_var LOGICAL :: climat_behavior ! enable bad behavior for climat database (only for instantaneous data) env_var = '' CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var) climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio) ! compute length of cumulation step in seconds CALL getval(timedelta_depop(step), asec=steps) ! create a mask of suitable timeranges ALLOCATE(mask_timerange(SIZE(itimerange))) mask_timerange(:) = itimerange(:)%timerange == tri & .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0 IF (PRESENT(dtratio)) THEN WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0) mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0 ELSEWHERE mask_timerange(:) = .false. END WHERE ELSE ! instantaneous mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0 ENDIF #ifdef DEBUG CALL l4f_log(l4f_debug, & '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// & t2c(count(mask_timerange))) #endif ! euristically determine whether we are dealing with an ! analysis/observation or a forecast dataset na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0) nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0) lforecast = nf >= na #ifdef DEBUG CALL l4f_log(l4f_debug, & 'recompute_stat_proc_agg, na: '//t2c(na)//', nf: '//t2c(nf)) #endif ! keep only timeranges of one type (really necessary?) IF (lforecast) THEN CALL l4f_log(l4f_info, & 'recompute_stat_proc_agg, processing in forecast mode') ELSE mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0 CALL l4f_log(l4f_info, & 'recompute_stat_proc_agg, processing in analysis mode') ENDIF #ifdef DEBUG CALL l4f_log(l4f_debug, & '(re)compute_stat_proc_agg, number of useful timeranges: '// & t2c(count(mask_timerange))) #endif IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume ALLOCATE(otime(0), otimerange(0), map_ttr(0,0)) IF (PRESENT(dtratio)) ALLOCATE(dtratio(0)) RETURN ENDIF ! determine start and end of processing period, should work with p2==0 lstart = datetime_miss IF (PRESENT(start)) lstart = start lend = itime(SIZE(itime)) ! compute some quantities maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange) maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange) minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange) IF (time_definition == 0) THEN ! reference time lend = lend + timedelta_new(sec=maxp1) ENDIF ! extend interval at the end in order to include all the data in case ! frac_valid<1; must use < and not <= in the DO WHILE loops some lines ! below in order to exclude the last full interval which would be empty lend = lend + step IF (lstart == datetime_miss) THEN ! autodetect lstart = itime(1) ! if autodetected, adjust to obtain real absolute start of data IF (time_definition == 0) THEN ! reference time lstart = lstart + timedelta_new(sec=minp1mp2) ELSE ! verification time ! go back to start of longest processing interval lstart = lstart - timedelta_new(sec=maxp2) ENDIF ! apply full_steps in analysis mode and when start is not specified ! (start by itself in analysis mode implies full_steps with respect to ! start instead of absolute full steps), todo in forecast mode IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN lstart = lstart - (mod(lstart, step)) ! round to step, (should be MODULO, not MOD) ENDIF ENDIF #ifdef DEBUG CALL l4f_log(l4f_debug, & 'recompute_stat_proc_agg, processing period: '//t2c(lstart)//' - '//t2c(lend)) #endif ! create output time and timerange lists IF (lforecast) THEN ! forecast mode IF (time_definition == 0) THEN ! reference time CALL insert(a_otime, itime) ! should I limit to elements itime >= lstart? ! apply start shift to timerange, not to reftime ! why did we use itime(SIZE(itime)) (last ref time)? ! CALL getval(lstart-itime(SIZE(itime)), asec=dstart) CALL getval(lstart-itime(1), asec=dstart) ! allow starting before first reftime but restrict dtstart to -steps ! dstart = MAX(0, dstart) IF (dstart < 0) dstart = mod(dstart, steps) DO p1 = steps + dstart, maxp1, steps CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps)) ENDDO ELSE ! verification time ! verification time in forecast mode is the ugliest case, because we ! don't know a priori how many different (thus incompatible) reference ! times we have, so some assumption of regularity has to be made. For ! this reason msteps, the minimum step between two times, is ! computed. We choose to compute it as a difference between itime ! elements, it could be also computed as difference of itimerange%p1 ! elements. But what if it is not modulo steps? mstepms = steps*1000_int_ll DO i = 2, SIZE(itime) CALL getval(itime(i)-itime(i-1), amsec=stepms) IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN msteps = stepms/1000_int_ll IF (mod(steps, msteps) == 0) mstepms = stepms ENDIF ENDDO msteps = mstepms/1000_int_ll tmptime = lstart + step DO WHILE(tmptime < lend) ! < since lend has been extended CALL insert_unique(a_otime, tmptime) tmptime = tmptime + step ENDDO ! in next loop, we used initially steps instead of msteps (see comment ! above), this gave much less combinations of time/timeranges with ! possible empty output; we start from msteps instead of steps in ! order to include partial initial intervals in case frac_valid<1; ! however, as a gemeral rule, for aggregation of forecasts it's better ! to use reference time DO p1 = msteps, maxp1, msteps CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps)) ENDDO ENDIF ELSE ! analysis mode tmptime = lstart + step DO WHILE(tmptime < lend) ! < since lend has been extended CALL insert_unique(a_otime, tmptime) tmptime = tmptime + step ENDDO CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps)) ENDIF CALL packarray(a_otime) CALL packarray(a_otimerange) otime => a_otime%array otimerange => a_otimerange%array CALL sort(otime) CALL sort(otimerange) ! delete local objects keeping the contents CALL delete(a_otime, nodealloc=.true.) CALL delete(a_otimerange, nodealloc=.true.) #ifdef DEBUG CALL l4f_log(l4f_debug, & 'recompute_stat_proc_agg, output time and timerange: '//& t2c(SIZE(otime))//', '//t2c(size(otimerange))) #endif IF (PRESENT(dtratio)) THEN ! count the possible i/o interval ratios DO k = 1, SIZE(itimerange) IF (itimerange(k)%p2 /= 0) & CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer ENDDO CALL packarray(a_dtratio) dtratio => a_dtratio%array CALL sort(dtratio) ! delete local object keeping the contents CALL delete(a_dtratio, nodealloc=.true.) #ifdef DEBUG CALL l4f_log(l4f_debug, & 'recompute_stat_proc_agg, found '//t2c(size(dtratio))// & ' possible aggregation ratios, from '// & t2c(dtratio(1))//' to '//t2c(dtratio(SIZE(dtratio)))) #endif ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange))) do_itimerange1: DO l = 1, SIZE(itimerange) IF (.NOT.mask_timerange(l)) cycle do_itimerange1 do_itime1: DO k = 1, SIZE(itime) CALL time_timerange_get_period(itime(k), itimerange(l), & time_definition, pstart1, pend1, reftime1) do_otimerange1: DO j = 1, SIZE(otimerange) do_otime1: DO i = 1, SIZE(otime) CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, & time_definition, pstart2, pend2, reftime2) IF (lforecast) THEN IF (reftime1 /= reftime2) cycle do_otime1 ENDIF IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. & mod(pstart1-pstart2, pend1-pstart1) == timedelta_0) THEN ! useful lmapper%it = k lmapper%itr = l lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer n = append(map_ttr(i,j), lmapper) cycle do_itime1 ! can contribute only to a single interval ENDIF ENDDO do_otime1 ENDDO do_otimerange1 ENDDO do_itime1 ENDDO do_itimerange1 ELSE ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange))) do_itimerange2: DO l = 1, SIZE(itimerange) IF (.NOT.mask_timerange(l)) cycle do_itimerange2 do_itime2: DO k = 1, SIZE(itime) CALL time_timerange_get_period(itime(k), itimerange(l), & time_definition, pstart1, pend1, reftime1) do_otimerange2: DO j = 1, SIZE(otimerange) do_otime2: DO i = 1, SIZE(otime) CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, & time_definition, pstart2, pend2, reftime2) IF (lforecast) THEN IF (reftime1 /= reftime2) cycle do_otime2 ENDIF IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful lmapper%it = k lmapper%itr = l IF (pstart1 == pstart2) THEN lmapper%extra_info = 1 ! start of interval ELSE IF (pend1 == pend2) THEN lmapper%extra_info = 2 ! end of interval ELSE lmapper%extra_info = imiss ENDIF lmapper%time = pstart1 ! = pend1, order by time? n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.) ! no CYCLE, a single input can contribute to multiple output intervals ENDIF ENDDO do_otime2 ENDDO do_otimerange2 ENDDO do_itime2 ENDDO do_itimerange2 ENDIF END SUBROUTINE recompute_stat_proc_agg_common SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, & max_step, weights) TYPE(datetime),INTENT(in) :: vertime(:) TYPE(datetime),INTENT(in) :: pstart TYPE(datetime),INTENT(in) :: pend LOGICAL,INTENT(in) :: time_mask(:) TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:) INTEGER :: i, nt TYPE(datetime),ALLOCATABLE :: lvertime(:) TYPE(datetime) :: half, nexthalf INTEGER(kind=int_ll) :: dt, tdt nt = count(time_mask) ALLOCATE(lvertime(nt)) lvertime = pack(vertime, mask=time_mask) IF (PRESENT(max_step)) THEN ! new way ! max_step = timedelta_0 ! DO i = 1, nt - 1 ! IF (lvertime(i+1) - lvertime(i) > max_step) & ! max_step = lvertime(i+1) - lvertime(i) ! ENDDO ! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart ! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt) ! old way IF (nt == 1) THEN max_step = pend - pstart ELSE half = lvertime(1) + (lvertime(2) - lvertime(1))/2 max_step = half - pstart DO i = 2, nt - 1 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2 IF (nexthalf - half > max_step) max_step = nexthalf - half half = nexthalf ENDDO IF (pend - half > max_step) max_step = pend - half ENDIF ENDIF IF (PRESENT(weights)) THEN IF (nt == 1) THEN weights(1) = 1.0 ELSE CALL getval(pend - pstart, amsec=tdt) half = lvertime(1) + (lvertime(2) - lvertime(1))/2 CALL getval(half - pstart, amsec=dt) weights(1) = dble(dt)/dble(tdt) DO i = 2, nt - 1 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2 CALL getval(nexthalf - half, amsec=dt) weights(i) = dble(dt)/dble(tdt) half = nexthalf ENDDO CALL getval(pend - half, amsec=dt) weights(nt) = dble(dt)/dble(tdt) ENDIF ENDIF END SUBROUTINE compute_stat_proc_agg_sw ! get start of period, end of period and reference time from time, ! timerange, according to time_definition. SUBROUTINE time_timerange_get_period(time, timerange, time_definition, & pstart, pend, reftime) TYPE(datetime),INTENT(in) :: time TYPE(vol7d_timerange),INTENT(in) :: timerange INTEGER,INTENT(in) :: time_definition TYPE(datetime),INTENT(out) :: reftime TYPE(datetime),INTENT(out) :: pstart TYPE(datetime),INTENT(out) :: pend TYPE(timedelta) :: p1, p2 p1 = timedelta_new(sec=timerange%p1) ! end of period p2 = timedelta_new(sec=timerange%p2) ! length of period IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. & ! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. & timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow? pstart = datetime_miss pend = datetime_miss reftime = datetime_miss RETURN ENDIF IF (time_definition == 0) THEN ! time == reference time reftime = time pend = time + p1 pstart = pend - p2 ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time pend = time pstart = time - p2 reftime = time - p1 ELSE pstart = datetime_miss pend = datetime_miss reftime = datetime_miss ENDIF END SUBROUTINE time_timerange_get_period ! get start of period, end of period and reference time from time, ! timerange, according to time_definition. SUBROUTINE ttr_decoder_decode(time, timerange, time_definition, decoded) TYPE(datetime),INTENT(in) :: time TYPE(vol7d_timerange),INTENT(in) :: timerange INTEGER,INTENT(in) :: time_definition TYPE(ttr_decoder),INTENT(out) :: decoded TYPE(timedelta) :: p1, p2 p1 = timedelta_new(sec=timerange%p1) ! end of period p2 = timedelta_new(sec=timerange%p2) ! length of period IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. & ! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. & timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow? decoded%pstart = datetime_miss decoded%pend = datetime_miss decoded%reftime = datetime_miss RETURN ENDIF IF (time_definition == 0) THEN ! time == reference time decoded%reftime = time decoded%pend = time + p1 decoded%pstart = decoded%pend - p2 ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time decoded%pend = time decoded%pstart = time - p2 decoded%reftime = time - p1 ELSE decoded%pstart = datetime_miss decoded%pend = datetime_miss decoded%reftime = datetime_miss ENDIF END SUBROUTINE ttr_decoder_decode ! get start of period, end of period and reference time from time, ! timerange, according to time_definition. step is taken separately ! from timerange and can be "popular" SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, & pstart, pend, reftime) TYPE(datetime),INTENT(in) :: time TYPE(vol7d_timerange),INTENT(in) :: timerange TYPE(timedelta),INTENT(in) :: step INTEGER,INTENT(in) :: time_definition TYPE(datetime),INTENT(out) :: reftime TYPE(datetime),INTENT(out) :: pstart TYPE(datetime),INTENT(out) :: pend TYPE(timedelta) :: p1 p1 = timedelta_new(sec=timerange%p1) ! end of period IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. & ! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. & timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow? pstart = datetime_miss pend = datetime_miss reftime = datetime_miss RETURN ENDIF IF (time_definition == 0) THEN ! time == reference time reftime = time pend = time + p1 pstart = pend - step ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time pend = time pstart = time - step reftime = time - p1 ELSE pstart = datetime_miss pend = datetime_miss reftime = datetime_miss ENDIF END SUBROUTINE time_timerange_get_period_pop ! get start of period, end of period and reference time from time, ! timerange, according to time_definition. step is taken separately ! from timerange and can be "popular" SUBROUTINE ttr_decoder_decode_pop(time, timerange, step, time_definition, decoded) TYPE(datetime),INTENT(in) :: time TYPE(vol7d_timerange),INTENT(in) :: timerange TYPE(timedelta),INTENT(in) :: step INTEGER,INTENT(in) :: time_definition TYPE(ttr_decoder),INTENT(out) :: decoded TYPE(timedelta) :: p1 p1 = timedelta_new(sec=timerange%p1) ! end of period IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. & ! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. & timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow? decoded%pstart = datetime_miss decoded%pend = datetime_miss decoded%reftime = datetime_miss RETURN ENDIF IF (time_definition == 0) THEN ! time == reference time decoded%reftime = time decoded%pend = time + p1 decoded%pstart = decoded%pend - step ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time decoded%pend = time decoded%pstart = time - step decoded%reftime = time - p1 ELSE decoded%pstart = datetime_miss decoded%pend = datetime_miss decoded%reftime = datetime_miss ENDIF END SUBROUTINE ttr_decoder_decode_pop ! set time, timerange%p1, timerange%p2 according to pstart, pend, ! reftime and time_definition. SUBROUTINE time_timerange_set_period(time, timerange, time_definition, & stat_proc, pstart, pend, reftime) TYPE(datetime),INTENT(out) :: time TYPE(vol7d_timerange),INTENT(out) :: timerange INTEGER,INTENT(in) :: time_definition INTEGER,INTENT(in) :: stat_proc TYPE(datetime),INTENT(in) :: pstart TYPE(datetime),INTENT(in) :: pend TYPE(datetime),INTENT(in) :: reftime TYPE(timedelta) :: p1, p2 INTEGER(kind=int_ll) :: dmsec IF (time_definition == 0) THEN ! time == reference time time = reftime p1 = pend - reftime p2 = pend - pstart ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time time = pend p1 = pend - reftime p2 = pend - pstart ELSE time = datetime_miss ENDIF timerange%timerange = stat_proc IF (time /= datetime_miss) THEN CALL getval(p1, amsec=dmsec) ! end of period timerange%p1 = int(dmsec/1000_int_ll) CALL getval(p2, amsec=dmsec) ! length of period timerange%p2 = int(dmsec/1000_int_ll) ELSE timerange%p1 = imiss timerange%p2 = imiss ENDIF END SUBROUTINE time_timerange_set_period END MODULE stat_proc_engine