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"

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

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
TYPE(volgrid6d),INTENT(out) :: that
INTEGER,INTENT(in) :: stat_proc_input
INTEGER,INTENT(in) :: stat_proc
TYPE(timedelta),INTENT(in) :: step
TYPE(datetime),INTENT(in),OPTIONAL :: start
LOGICAL,INTENT(in),OPTIONAL :: full_steps
REAL,INTENT(in),OPTIONAL :: frac_valid
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
LOGICAL , INTENT(in),OPTIONAL :: clone

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


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
TYPE(volgrid6d),INTENT(out) :: that
INTEGER,INTENT(in) :: stat_proc
TYPE(timedelta),INTENT(in) :: step
TYPE(datetime),INTENT(in),OPTIONAL :: start
LOGICAL,INTENT(in),OPTIONAL :: full_steps
REAL,INTENT(in),OPTIONAL :: frac_valid
LOGICAL, INTENT(in),OPTIONAL :: clone
INTEGER,INTENT(in),OPTIONAL :: stat_proc_input

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


SUBROUTINE volgrid6d_compute_stat_proc_agg(this, that, stat_proc, &
 step, start, full_steps, max_step, clone)
TYPE(volgrid6d),INTENT(inout) :: this
TYPE(volgrid6d),INTENT(out) :: that
INTEGER,INTENT(in) :: stat_proc
TYPE(timedelta),INTENT(in) :: step
TYPE(datetime),INTENT(in),OPTIONAL :: start
LOGICAL,INTENT(in),OPTIONAL :: full_steps
TYPE(timedelta),INTENT(in),OPTIONAL :: max_step
LOGICAL , INTENT(in),OPTIONAL :: clone

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


SUBROUTINE volgrid6d_recompute_stat_proc_diff(this, that, stat_proc, step, full_steps, start, clone)
TYPE(volgrid6d),INTENT(inout) :: this
TYPE(volgrid6d),INTENT(out) :: that
INTEGER,INTENT(in) :: stat_proc
TYPE(timedelta),INTENT(in) :: step
LOGICAL,INTENT(in),OPTIONAL :: full_steps
TYPE(datetime),INTENT(in),OPTIONAL :: start
LOGICAL,INTENT(in),OPTIONAL :: clone
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


SUBROUTINE volgrid6d_compute_stat_proc_metamorph(this, that, stat_proc_input, stat_proc, clone)
TYPE(volgrid6d),INTENT(inout) :: this
TYPE(volgrid6d),INTENT(out) :: that
INTEGER,INTENT(in) :: stat_proc_input
INTEGER,INTENT(in) :: stat_proc
LOGICAL , INTENT(in),OPTIONAL :: clone

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

SUBROUTINE volgrid6d_compute_vert_coord_var(this, level, volgrid_lev)
TYPE(volgrid6d),INTENT(in) :: this
TYPE(vol7d_level),INTENT(in) :: level
TYPE(volgrid6d),INTENT(out) :: volgrid_lev

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