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"

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