volgrid6d_class.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
USE vol7d_level_class
USE vol7d_timerange_class
USE vol7d_class
USE grid_transform_class
USE geo_proj_class
USE grid_class
USE grid_dim_class
USE datetime_class
USE volgrid6d_var_class
USE log4fortran
USE array_utilities
USE grid_id_class
USE gridinfo_class
USE optional_values
!USE file_utilities
#ifdef HAVE_DBALLE
USE vol7d_dballe_class
#endif
IMPLICIT NONE

character (len=255),parameter:: subcategory="volgrid6d_class"

type volgrid6d
  type(griddim_def) :: griddim
  TYPE(datetime),pointer :: time(:)
  TYPE(vol7d_timerange),pointer :: timerange(:)
  TYPE(vol7d_level),pointer :: level(:)
  TYPE(volgrid6d_var),pointer :: var(:)
  TYPE(grid_id),POINTER :: gaid(:,:,:,:)
  REAL,POINTER :: voldati(:,:,:,:,:,:)
  INTEGER :: time_definition
  TYPE(l4f_handle) :: category
end type volgrid6d

INTERFACE init
  MODULE PROCEDURE volgrid6d_init
END INTERFACE

INTERFACE delete
  MODULE PROCEDURE volgrid6d_delete, volgrid6dv_delete
END INTERFACE

INTERFACE import
  MODULE PROCEDURE volgrid6d_read_from_file
  MODULE PROCEDURE import_from_gridinfo, import_from_gridinfovv, &
   volgrid6d_import_from_file
END INTERFACE

INTERFACE export
  MODULE PROCEDURE volgrid6d_write_on_file
  MODULE PROCEDURE export_to_gridinfo, export_to_gridinfov, export_to_gridinfovv,&
   volgrid6d_export_to_file
END INTERFACE

! methods for computing transformations through an initialised
! grid_transform object, probably too low level to be interfaced
INTERFACE compute
  MODULE PROCEDURE volgrid6d_transform_compute, volgrid6d_v7d_transform_compute,&
   v7d_volgrid6d_transform_compute, v7d_v7d_transform_compute
END INTERFACE

INTERFACE transform
  MODULE PROCEDURE volgrid6d_transform, volgrid6dv_transform,&
   volgrid6d_v7d_transform, volgrid6dv_v7d_transform, v7d_volgrid6d_transform, &
   v7d_v7d_transform
END INTERFACE

INTERFACE wind_rot
  MODULE PROCEDURE vg6d_wind_rot
END INTERFACE

INTERFACE wind_unrot
  MODULE PROCEDURE vg6d_wind_unrot
END INTERFACE

INTERFACE display
  MODULE PROCEDURE display_volgrid6d,display_volgrid6dv
END INTERFACE

INTERFACE rounding
  MODULE PROCEDURE vg6d_rounding, vg6dv_rounding
END INTERFACE

private

PUBLIC volgrid6d,init,delete,export,import,compute,transform, &
 wind_rot,wind_unrot,vg6d_c2a,display,volgrid6d_alloc,volgrid6d_alloc_vol, &
 volgrid_get_vol_2d, volgrid_set_vol_2d, volgrid_get_vol_3d, volgrid_set_vol_3d
PUBLIC rounding, vg6d_reduce

CONTAINS


SUBROUTINE volgrid6d_init(this, griddim, time_definition, categoryappend)
TYPE(volgrid6d) :: this
TYPE(griddim_def),OPTIONAL :: griddim
INTEGER,INTENT(IN),OPTIONAL :: time_definition
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

character(len=512) :: a_name

if (present(categoryappend))then
   call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
else
   call l4f_launcher(a_name,a_name_append=trim(subcategory))
endif
this%category=l4f_category_get_handle(a_name)

#ifdef DEBUG
call l4f_category_log(this%category,l4f_debug,"init")
#endif

call init(this%griddim)

if (present(griddim))then
  call copy (griddim,this%griddim)
end if

if(present(time_definition)) then
  this%time_definition = time_definition
else
  this%time_definition = 0 !default to reference time
end if

nullify (this%time,this%timerange,this%level,this%var)
nullify (this%gaid,this%voldati)          

END SUBROUTINE volgrid6d_init


SUBROUTINE volgrid6d_alloc(this, dim, ntime, nlevel, ntimerange, nvar, ini)
TYPE(volgrid6d),INTENT(inout) :: this
TYPE(grid_dim),INTENT(in),OPTIONAL :: dim
INTEGER,INTENT(in),OPTIONAL :: ntime
INTEGER,INTENT(in),OPTIONAL :: nlevel
INTEGER,INTENT(in),OPTIONAL :: ntimerange
INTEGER,INTENT(in),OPTIONAL :: nvar
LOGICAL,INTENT(in),OPTIONAL :: ini

INTEGER :: i, stallo
LOGICAL :: linit

#ifdef DEBUG
call l4f_category_log(this%category,l4f_debug,"alloc")
#endif

IF (PRESENT(ini)) THEN
  linit = ini
ELSE
  linit = .false.
ENDIF


if (present(dim)) call copy (dim,this%griddim%dim)


IF (PRESENT(ntime)) THEN
  IF (ntime >= 0) THEN
    IF (ASSOCIATED(this%time)) DEALLOCATE(this%time)
#ifdef DEBUG
    call l4f_category_log(this%category,l4f_debug,"alloc ntime "//to_char(ntime))
#endif
    ALLOCATE(this%time(ntime),stat=stallo)
    if (stallo /=0)then
      call l4f_category_log(this%category,l4f_fatal,"allocating memory")
      CALL raise_fatal_error()
    end if
    IF (linit) THEN
      DO i = 1, ntime
        this%time(i) = datetime_miss
!        CALL init(this%time(i)) ! senza argomento inizializza a zero non missing
                                 ! baco di datetime?
      ENDDO
    ENDIF
  ENDIF
ENDIF
IF (PRESENT(nlevel)) THEN
  IF (nlevel >= 0) THEN
    IF (ASSOCIATED(this%level)) DEALLOCATE(this%level)
#ifdef DEBUG
    call l4f_category_log(this%category,l4f_debug,"alloc nlevel "//to_char(nlevel))
#endif
    ALLOCATE(this%level(nlevel),stat=stallo)
    if (stallo /=0)then
      call l4f_category_log(this%category,l4f_fatal,"allocating memory")
      CALL raise_fatal_error()
    end if
    IF (linit) THEN
      DO i = 1, nlevel
        CALL init(this%level(i))
      ENDDO
    ENDIF
  ENDIF
ENDIF
IF (PRESENT(ntimerange)) THEN
  IF (ntimerange >= 0) THEN
    IF (ASSOCIATED(this%timerange)) DEALLOCATE(this%timerange)
#ifdef DEBUG
    call l4f_category_log(this%category,l4f_debug,"alloc ntimerange "//to_char(ntimerange))
#endif
    ALLOCATE(this%timerange(ntimerange),stat=stallo)
    if (stallo /=0)then
      call l4f_category_log(this%category,l4f_fatal,"allocating memory")
      CALL raise_fatal_error()
    end if
    IF (linit) THEN
      DO i = 1, ntimerange
        CALL init(this%timerange(i))
      ENDDO
    ENDIF
  ENDIF
ENDIF
IF (PRESENT(nvar)) THEN
  IF (nvar >= 0) THEN
    IF (ASSOCIATED(this%var)) DEALLOCATE(this%var)
#ifdef DEBUG
    call l4f_category_log(this%category,l4f_debug,"alloc nvar "//to_char(nvar))
#endif
    ALLOCATE(this%var(nvar),stat=stallo)
    if (stallo /=0)then
      call l4f_category_log(this%category,l4f_fatal,"allocating memory")
      CALL raise_fatal_error()
    end if
    IF (linit) THEN
      DO i = 1, nvar
        CALL init(this%var(i))
      ENDDO
    ENDIF
  ENDIF
ENDIF

end SUBROUTINE volgrid6d_alloc


SUBROUTINE volgrid6d_alloc_vol(this, ini, inivol, decode)
TYPE(volgrid6d),INTENT(inout) :: this
LOGICAL,INTENT(in),OPTIONAL :: ini
LOGICAL,INTENT(in),OPTIONAL :: inivol
LOGICAL,INTENT(in),OPTIONAL :: decode

INTEGER :: stallo
LOGICAL :: linivol

#ifdef DEBUG
call l4f_category_log(this%category,l4f_debug,"start alloc_vol")
#endif

IF (PRESENT(inivol)) THEN ! opposite condition, cannot use optio_log
  linivol = inivol
ELSE
  linivol = .true.
ENDIF

IF (this%griddim%dim%nx > 0 .AND. this%griddim%dim%ny > 0) THEN
! allocate dimension descriptors to a minimal size for having a
! non-null volume
  IF (.NOT. ASSOCIATED(this%var)) CALL volgrid6d_alloc(this, nvar=1, ini=ini)
  IF (.NOT. ASSOCIATED(this%time)) CALL volgrid6d_alloc(this, ntime=1, ini=ini)
  IF (.NOT. ASSOCIATED(this%level)) CALL volgrid6d_alloc(this, nlevel=1, ini=ini)
  IF (.NOT. ASSOCIATED(this%timerange)) CALL volgrid6d_alloc(this, ntimerange=1, ini=ini)
  
  IF (optio_log(decode)) THEN
    IF (.NOT.ASSOCIATED(this%voldati)) THEN
#ifdef DEBUG
      CALL l4f_category_log(this%category,l4f_debug,"alloc voldati volume")
#endif

      ALLOCATE(this%voldati(this%griddim%dim%nx,this%griddim%dim%ny,&
       SIZE(this%level), SIZE(this%time), &
       SIZE(this%timerange), SIZE(this%var)),stat=stallo)
      IF (stallo /= 0)THEN
        CALL l4f_category_log(this%category,l4f_fatal,"allocating memory")
        CALL raise_fatal_error()
      ENDIF
  
! this is not really needed if we can check other flags for a full
! field missing values
      IF (linivol) this%voldati = rmiss
      this%voldati = rmiss
    ENDIF
  ENDIF

  IF (.NOT.ASSOCIATED(this%gaid)) THEN
#ifdef DEBUG
    CALL l4f_category_log(this%category,l4f_debug,"alloc gaid volume")
#endif
    ALLOCATE(this%gaid(SIZE(this%level), SIZE(this%time), &
     SIZE(this%timerange), SIZE(this%var)),stat=stallo)
    IF (stallo /= 0)THEN
      CALL l4f_category_log(this%category,l4f_fatal,"allocating memory")
      CALL raise_fatal_error()
    ENDIF

    IF (linivol) THEN
!!$    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)
!!$            this%gaid(i,ii,iii,iiii) = grid_id_new() ! optimize?
!!$          ENDDO
!!$        ENDDO
!!$      ENDDO
!!$    ENDDO

      this%gaid = grid_id_new()
    ENDIF
  ENDIF
  
ELSE
  CALL l4f_category_log(this%category,l4f_fatal,'volgrid6d_alloc_vol: &
   &trying to allocate a volume with an invalid or unspecified horizontal grid')
  CALL raise_fatal_error()
ENDIF

END SUBROUTINE volgrid6d_alloc_vol


SUBROUTINE volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
TYPE(volgrid6d),INTENT(in) :: this
INTEGER,INTENT(in) :: ilevel
INTEGER,INTENT(in) :: itime
INTEGER,INTENT(in) :: itimerange
INTEGER,INTENT(in) :: ivar
REAL,POINTER :: voldati(:,:)

IF (ASSOCIATED(this%voldati)) THEN
  voldati => this%voldati(:,:,ilevel,itime,itimerange,ivar)
  RETURN
ELSE
  IF (.NOT.ASSOCIATED(voldati)) THEN
    ALLOCATE(voldati(this%griddim%dim%nx,this%griddim%dim%ny))
  ENDIF
  CALL grid_id_decode_data(this%gaid(ilevel,itime,itimerange,ivar), voldati)
ENDIF

END SUBROUTINE volgrid_get_vol_2d


SUBROUTINE volgrid_get_vol_3d(this, itime, itimerange, ivar, voldati, zlist)
TYPE(volgrid6d),INTENT(in) :: this
INTEGER,INTENT(in) :: itime
INTEGER,INTENT(in) :: itimerange
INTEGER,INTENT(in) :: ivar
REAL,POINTER :: voldati(:,:,:)
LOGICAL,OPTIONAL,INTENT(out),ALLOCATABLE :: zlist(:)

INTEGER :: ilevel
LOGICAL :: done

IF (ASSOCIATED(this%voldati)) THEN
  voldati => this%voldati(:,:,:,itime,itimerange,ivar)
  IF (PRESENT(zlist)) zlist(:) = c_e(this%gaid(:,itime,itimerange,ivar))
  RETURN
ELSE
  IF (.NOT.ASSOCIATED(voldati)) THEN
    ALLOCATE(voldati(this%griddim%dim%nx,this%griddim%dim%ny,SIZE(this%level)))
  ENDIF
  IF (PRESENT(zlist)) THEN
    ALLOCATE(zlist(SIZE(this%level)))
  ENDIF
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP MASTER
  DO ilevel = 1, SIZE(this%level)
!$OMP TASK FIRSTPRIVATE(ilevel),PRIVATE(done)
    CALL grid_id_decode_data(this%gaid(ilevel,itime,itimerange,ivar), &
     voldati(:,:,ilevel),done)
    IF (PRESENT(zlist)) zlist(ilevel) = done
!$OMP END TASK
  ENDDO
!$OMP END MASTER
!$OMP END PARALLEL
ENDIF

END SUBROUTINE volgrid_get_vol_3d


SUBROUTINE volgrid_set_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
TYPE(volgrid6d),INTENT(inout) :: this
INTEGER,INTENT(in) :: ilevel
INTEGER,INTENT(in) :: itime
INTEGER,INTENT(in) :: itimerange
INTEGER,INTENT(in) :: ivar
REAL,INTENT(in) :: voldati(:,:)

IF (ASSOCIATED(this%voldati)) THEN
  RETURN
ELSE
  CALL grid_id_encode_data(this%gaid(ilevel,itime,itimerange,ivar), voldati)
ENDIF

END SUBROUTINE volgrid_set_vol_2d


SUBROUTINE volgrid_set_vol_3d(this, itime, itimerange, ivar, voldati)
TYPE(volgrid6d),INTENT(inout) :: this
INTEGER,INTENT(in) :: itime
INTEGER,INTENT(in) :: itimerange
INTEGER,INTENT(in) :: ivar
REAL,INTENT(in) :: voldati(:,:,:)

INTEGER :: ilevel

IF (ASSOCIATED(this%voldati)) THEN
  RETURN
ELSE
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP MASTER
  DO ilevel = 1, SIZE(this%level)
!$OMP TASK FIRSTPRIVATE(ilevel)
    CALL grid_id_encode_data(this%gaid(ilevel,itime,itimerange,ivar), &
     voldati(:,:,ilevel))
!$OMP END TASK
  ENDDO
!$OMP END MASTER
!$OMP END PARALLEL
ENDIF

END SUBROUTINE volgrid_set_vol_3d


SUBROUTINE volgrid6d_delete(this)
TYPE(volgrid6d),INTENT(inout) :: this

INTEGER :: i, ii, iii, iiii

#ifdef DEBUG
call l4f_category_log(this%category,l4f_debug,"delete")
#endif

if (associated(this%gaid))then

  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)
          CALL delete(this%gaid(i,ii,iii,iiii))
        ENDDO
      ENDDO
    ENDDO
  ENDDO
  DEALLOCATE(this%gaid)

end if

call delete(this%griddim)

!  call delete(this%time)
!  call delete(this%timerange)
!  call delete(this%level)
!  call delete(this%var)

if (associated( this%time )) deallocate(this%time)
if (associated( this%timerange )) deallocate(this%timerange)
if (associated( this%level )) deallocate(this%level)
if (associated( this%var )) deallocate(this%var)

if (associated(this%voldati))deallocate(this%voldati)


                                !chiudo il logger
call l4f_category_delete(this%category)
  
END SUBROUTINE volgrid6d_delete


subroutine volgrid6d_write_on_file (this,unit,description,filename,filename_auto)

TYPE(volgrid6d),INTENT(IN) :: this
integer,optional,intent(inout) :: unit
character(len=*),intent(in),optional :: filename
character(len=*),intent(out),optional :: filename_auto
character(len=*),INTENT(IN),optional :: description

integer :: lunit
character(len=254) :: ldescription,arg,lfilename
integer :: ntime, ntimerange, nlevel, nvar
integer :: tarray(8)
logical :: opened,exist

#ifdef DEBUG
call l4f_category_log(this%category,l4f_debug,"write on file")
#endif

ntime=0
ntimerange=0
nlevel=0
nvar=0

!call idate(im,id,iy)
call date_and_time(values=tarray)
call getarg(0,arg)

if (present(description))then
  ldescription=description
else
  ldescription="Volgrid6d generated by: "//trim(arg)
end if

if (.not. present(unit))then
  lunit=getunit()
else
  if (unit==0)then
    lunit=getunit()
    unit=lunit
  else
    lunit=unit
  end if
end if

lfilename=trim(arg)//".vg6d"
if (index(arg,'/',back=.true.) > 0) lfilename=lfilename(index(arg,'/',back=.true.)+1 : )

if (present(filename))then
  if (filename /= "")then
    lfilename=filename
  end if
end if

if (present(filename_auto))filename_auto=lfilename


inquire(unit=lunit,opened=opened)
if (.not. opened) then 
  inquire(file=lfilename,exist=exist)
  if (exist) CALL raise_error('file exist; cannot open new file')
  if (.not.exist) open (unit=lunit,file=lfilename,form="UNFORMATTED")
  !print *, "opened: ",lfilename
end if

if (associated(this%time)) ntime=size(this%time)
if (associated(this%timerange)) ntimerange=size(this%timerange)
if (associated(this%level)) nlevel=size(this%level)
if (associated(this%var)) nvar=size(this%var)


write(unit=lunit)ldescription
write(unit=lunit)tarray

call write_unit( this%griddim,lunit)
write(unit=lunit) ntime, ntimerange, nlevel, nvar

!! prime 4 dimensioni
if (associated(this%time))      call write_unit(this%time, lunit)
if (associated(this%level))     write(unit=lunit)this%level
if (associated(this%timerange)) write(unit=lunit)this%timerange
if (associated(this%var))       write(unit=lunit)this%var    


!! Volumi di valori dati

if (associated(this%voldati))     write(unit=lunit)this%voldati

if (.not. present(unit)) close(unit=lunit)

end subroutine volgrid6d_write_on_file


subroutine volgrid6d_read_from_file (this,unit,filename,description,tarray,filename_auto)

TYPE(volgrid6d),INTENT(OUT) :: this
integer,intent(inout),optional :: unit
character(len=*),INTENT(in),optional :: filename
character(len=*),intent(out),optional :: filename_auto
character(len=*),INTENT(out),optional :: description
integer,intent(out),optional :: tarray(8)

integer :: ntime, ntimerange, nlevel, nvar

character(len=254) :: ldescription,lfilename,arg
integer :: ltarray(8),lunit
logical :: opened,exist

#ifdef DEBUG
call l4f_category_log(this%category,l4f_debug,"read from file")
#endif

call getarg(0,arg)

if (.not. present(unit))then
  lunit=getunit()
else
  if (unit==0)then
    lunit=getunit()
    unit=lunit
  else
    lunit=unit
  end if
end if

lfilename=trim(arg)//".vg6d"
if (index(arg,'/',back=.true.) > 0) lfilename=lfilename(index(arg,'/',back=.true.)+1 : )

if (present(filename))then
  if (filename /= "")then
    lfilename=filename
  end if
end if

if (present(filename_auto))filename_auto=lfilename


inquire(unit=lunit,opened=opened)
if (.not. opened) then 
  inquire(file=lfilename,exist=exist)
  IF (.NOT. exist) CALL raise_fatal_error('file '//trim(lfilename)//' does not exist, cannot open')
  open (unit=lunit,file=lfilename,form="UNFORMATTED")
end if

read(unit=lunit)ldescription
read(unit=lunit)ltarray

call l4f_log(l4f_info,"Info: reading volgrid6d from file: "//trim(lfilename))
call l4f_log(l4f_info,"Info: description: "//trim(ldescription))
!call l4f_log("Info: written on ",ltarray)

if (present(description))description=ldescription
if (present(tarray))tarray=ltarray


call read_unit( this%griddim,lunit)
read(unit=lunit) ntime, ntimerange, nlevel, nvar


call volgrid6d_alloc (this, &
 ntime=ntime, ntimerange=ntimerange, nlevel=nlevel, nvar=nvar)

call volgrid6d_alloc_vol (this)

if (associated(this%time))      call read_unit(this%time, lunit)
if (associated(this%level))     read(unit=lunit)this%level
if (associated(this%timerange)) read(unit=lunit)this%timerange
if (associated(this%var))       read(unit=lunit)this%var    


!! Volumi di valori 

if (associated(this%voldati))     read(unit=lunit)this%voldati

if (.not. present(unit)) close(unit=lunit)

end subroutine volgrid6d_read_from_file


SUBROUTINE import_from_gridinfo(this, gridinfo, force, dup_mode, clone, &
 isanavar)
TYPE(volgrid6d),INTENT(inout) :: this
TYPE(gridinfo_def),INTENT(in) :: gridinfo
LOGICAL,INTENT(in),OPTIONAL :: force
INTEGER,INTENT(in),OPTIONAL :: dup_mode
LOGICAL , INTENT(in),OPTIONAL :: clone
LOGICAL,INTENT(IN),OPTIONAL :: isanavar

CHARACTER(len=255) :: type
INTEGER :: itime0, itimerange0, itime1, itimerange1, itime, itimerange, &
 ilevel, ivar, ldup_mode
LOGICAL :: dup
TYPE(datetime) :: correctedtime
TYPE(vol7d_timerange) :: correctedtimerange
REAL,ALLOCATABLE :: tmpgrid(:,:)

IF (PRESENT(dup_mode)) THEN
  ldup_mode = dup_mode
ELSE
  ldup_mode = 0
ENDIF

call get_val(this%griddim,proj_type=type)

#ifdef DEBUG
call l4f_category_log(this%category,l4f_debug,"import_from_gridinfo: "//trim(type))
#endif

if (.not. c_e(type))then
  call copy(gridinfo%griddim, this%griddim)
! ho gia` fatto init, altrimenti non potrei fare la get_val(this%griddim)
! per cui meglio non ripetere
!   call init(this,gridinfo%griddim,categoryappend)
  CALL volgrid6d_alloc_vol(this, ini=.true.) ! decode?

else if (.not. (this%griddim == gridinfo%griddim ))then

  CALL l4f_category_log(this%category,l4f_error, &
   "volgrid and gridinfo grid type or size are different, gridinfo rejected")
  CALL raise_error()
  RETURN

end if

! Cerco gli indici del campo da inserire, se non trovo metto nel primo missing
ilevel = index(this%level, gridinfo%level)
IF (ilevel == 0 .AND. optio_log(force)) THEN
  ilevel = index(this%level, vol7d_level_miss)
  IF (ilevel /= 0) this%level(ilevel) = gridinfo%level
ENDIF

IF (ilevel == 0) THEN
  CALL l4f_category_log(this%category,l4f_error, &
   "volgrid6d: level not valid for volume, gridinfo rejected")
  CALL raise_error()
  RETURN
ENDIF

IF (optio_log(isanavar)) THEN ! assign to all times and timeranges
  itime0 = 1
  itime1 = SIZE(this%time)
  itimerange0 = 1
  itimerange1 = SIZE(this%timerange)
ELSE ! usual case
  correctedtime = gridinfo%time
  IF (this%time_definition == 1 .OR. this%time_definition == 2) correctedtime = correctedtime + &
   timedelta_new(sec=gridinfo%timerange%p1)
  itime0 = index(this%time, correctedtime)
  IF (itime0 == 0 .AND. optio_log(force)) THEN
    itime0 = index(this%time, datetime_miss)
    IF (itime0 /= 0) this%time(itime0) = correctedtime
  ENDIF
  IF (itime0 == 0) THEN
    CALL l4f_category_log(this%category,l4f_error, &
     "volgrid6d: time not valid for volume, gridinfo rejected")
    CALL raise_error()
    RETURN
  ENDIF
  itime1 = itime0

  correctedtimerange = gridinfo%timerange
  IF (this%time_definition == 2) correctedtimerange%p1 = 0
  itimerange0 = index(this%timerange, correctedtimerange)
  IF (itimerange0 == 0 .AND. optio_log(force)) THEN
    itimerange0 = index(this%timerange, vol7d_timerange_miss)
    IF (itimerange0 /= 0) this%timerange(itimerange0) = gridinfo%timerange
  ENDIF
  IF (itimerange0 == 0) THEN
    CALL l4f_category_log(this%category,l4f_error, &
     "volgrid6d: timerange not valid for volume, gridinfo rejected")
    CALL raise_error()
    RETURN
  ENDIF
  itimerange1 = itimerange0
ENDIF

ivar = index(this%var, gridinfo%var)
IF (ivar == 0 .AND. optio_log(force)) THEN
  ivar = index(this%var, volgrid6d_var_miss)
  IF (ivar /= 0) this%var(ivar) = gridinfo%var
ENDIF
IF (ivar == 0) THEN
  CALL l4f_category_log(this%category,l4f_error, &
   "volgrid6d: var not valid for volume, gridinfo rejected")
  CALL raise_error()
  RETURN
ENDIF

DO itimerange = itimerange0, itimerange1
  DO itime = itime0, itime1
    IF (ASSOCIATED(this%gaid)) THEN
      dup = .false.
      IF (c_e(this%gaid(ilevel,itime,itimerange,ivar))) THEN
        dup = .true.
        CALL l4f_category_log(this%category,l4f_warn,"gaid exist: grib duplicated")
! avoid memory leaks
        IF (optio_log(clone)) CALL delete(this%gaid(ilevel,itime,itimerange,ivar))
      ENDIF

      IF (optio_log(clone)) THEN
        CALL copy(gridinfo%gaid, this%gaid(ilevel,itime,itimerange,ivar))
#ifdef DEBUG
        CALL l4f_category_log(this%category,l4f_debug,"cloning to a new gaid")
#endif
      ELSE
        this%gaid(ilevel,itime,itimerange,ivar) = gridinfo%gaid
      ENDIF

      IF (ASSOCIATED(this%voldati))THEN
        IF (.NOT.dup .OR. ldup_mode == 0) THEN
          this%voldati(:,:,ilevel,itime,itimerange,ivar) = decode_gridinfo(gridinfo)
        ELSE IF (ldup_mode == 1) THEN
          tmpgrid = decode_gridinfo(gridinfo) ! f2003 automatic allocation
          WHERE(c_e(tmpgrid))
            this%voldati(:,:,ilevel,itime,itimerange,ivar) = tmpgrid(:,:)
          END WHERE
        ELSE IF (ldup_mode == 2) THEN
          WHERE(.NOT.c_e(this%voldati(:,:,ilevel,itime,itimerange,ivar)))
            this%voldati(:,:,ilevel,itime,itimerange,ivar) = decode_gridinfo(gridinfo)
          END WHERE
        ENDIF
      ENDIF

    ELSE
      CALL l4f_category_log(this%category,l4f_error, &
       "gaid not allocated, you probably need to call volgrid6d_alloc_vol first")
      CALL raise_error()
      RETURN
    ENDIF
  ENDDO
ENDDO


END SUBROUTINE import_from_gridinfo


SUBROUTINE export_to_gridinfo(this, gridinfo, itime, itimerange, ilevel, ivar, &
 gaid_template, clone)
TYPE(volgrid6d),INTENT(in) :: this
TYPE(gridinfo_def),INTENT(inout) :: gridinfo
INTEGER :: itime
INTEGER :: itimerange
INTEGER :: ilevel
INTEGER :: ivar
TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template
LOGICAL,INTENT(in),OPTIONAL :: clone

TYPE(grid_id) :: gaid
LOGICAL :: usetemplate
REAL,POINTER :: voldati(:,:)
TYPE(datetime) :: correctedtime

#ifdef DEBUG
CALL l4f_category_log(this%category,l4f_debug,"export_to_gridinfo")
#endif

IF (.NOT.c_e(this%gaid(ilevel,itime,itimerange,ivar))) THEN
#ifdef DEBUG
  CALL l4f_category_log(this%category,l4f_debug,"empty gaid found, skipping export")
#endif
  RETURN
ENDIF

usetemplate = .false.
IF (PRESENT(gaid_template)) THEN
  CALL copy(gaid_template, gaid)
#ifdef DEBUG
  CALL l4f_category_log(this%category,l4f_debug,"template cloned to a new gaid")
#endif
  usetemplate = c_e(gaid)
ENDIF

IF (.NOT.usetemplate) THEN
  IF (optio_log(clone)) THEN
    CALL copy(this%gaid(ilevel,itime,itimerange,ivar), gaid)
#ifdef DEBUG
    CALL l4f_category_log(this%category,l4f_debug,"original gaid cloned to a new one")
#endif
  ELSE
    gaid = this%gaid(ilevel,itime,itimerange,ivar)
  ENDIF
ENDIF

IF (this%time_definition == 1 .OR. this%time_definition == 2) THEN
  correctedtime = this%time(itime) - &
   timedelta_new(sec=this%timerange(itimerange)%p1)
ELSE
  correctedtime = this%time(itime)
ENDIF

CALL init(gridinfo,gaid, this%griddim, correctedtime, this%timerange(itimerange), &
 this%level(ilevel), this%var(ivar))

! reset the gridinfo, bad but necessary at this point for encoding the field
CALL export(gridinfo%griddim, gridinfo%gaid)
! encode the field
IF (ASSOCIATED(this%voldati)) THEN
  CALL encode_gridinfo(gridinfo, this%voldati(:,:,ilevel,itime,itimerange,ivar))
ELSE IF (usetemplate) THEN ! field must be forced into template in this case
  NULLIFY(voldati)
  CALL volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
  CALL encode_gridinfo(gridinfo, voldati)
  DEALLOCATE(voldati)
ENDIF

END SUBROUTINE export_to_gridinfo


SUBROUTINE import_from_gridinfovv(this, gridinfov, dup_mode, clone, decode, &
 time_definition, anavar, categoryappend)
TYPE(volgrid6d),POINTER :: this(:)
TYPE(arrayof_gridinfo),INTENT(in) :: gridinfov
INTEGER,INTENT(in),OPTIONAL :: dup_mode
LOGICAL , INTENT(in),OPTIONAL :: clone
LOGICAL,INTENT(in),OPTIONAL :: decode
INTEGER,INTENT(IN),OPTIONAL :: time_definition
CHARACTER(len=*),INTENT(IN),OPTIONAL :: anavar(:)
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

INTEGER :: i, j, stallo
INTEGER :: ngrid, ntime, ntimerange, nlevel, nvar, ltime_definition
TYPE(l4f_handle) :: category
CHARACTER(len=512) :: a_name
TYPE(datetime),ALLOCATABLE :: correctedtime(:)
LOGICAL,ALLOCATABLE :: isanavar(:)
TYPE(vol7d_var) :: lvar
TYPE(vol7d_timerange),ALLOCATABLE :: correctedtimerange(:)

! category temporanea (altrimenti non possiamo loggare)
if (present(categoryappend))then
  call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
else
  call l4f_launcher(a_name,a_name_append=trim(subcategory))
endif
category=l4f_category_get_handle(a_name)

#ifdef DEBUG
call l4f_category_log(category,l4f_debug,"start import_from_gridinfovv")
#endif

IF (PRESENT(time_definition)) THEN
  ltime_definition = max(min(time_definition, 2), 0)
ELSE
  ltime_definition = 0
ENDIF

ngrid=count_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim,back=.true.)
CALL l4f_category_log(category,l4f_info, t2c(ngrid)// &
 ' different grid definition(s) found in input data')

ALLOCATE(this(ngrid),stat=stallo)
IF (stallo /= 0)THEN
  CALL l4f_category_log(category,l4f_fatal,"allocating memory")
  CALL raise_fatal_error()
ENDIF
DO i = 1, ngrid
  IF (PRESENT(categoryappend))THEN
    CALL init(this(i), time_definition=ltime_definition, categoryappend=trim(categoryappend)//"-vol"//t2c(i))
  ELSE
    CALL init(this(i), time_definition=ltime_definition, categoryappend="vol"//t2c(i))
  ENDIF
ENDDO

this(:)%griddim=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim, &
 ngrid, back=.true.)

! mark elements as ana variables (time-independent)
ALLOCATE(isanavar(gridinfov%arraysize))
isanavar(:) = .false.
IF (PRESENT(anavar)) THEN
  DO i = 1, gridinfov%arraysize
    DO j = 1, SIZE(anavar)
      lvar = convert(gridinfov%array(i)%var)
      IF (lvar%btable == anavar(j)) THEN
        isanavar(i) = .true.
        EXIT
      ENDIF
    ENDDO
  ENDDO
  CALL l4f_category_log(category,l4f_info,t2c(count(isanavar))//'/'// &
   t2c(gridinfov%arraysize)//' constant-data messages found in input data')
ENDIF

IF (ltime_definition == 1 .OR. ltime_definition == 2) THEN ! verification time
  ALLOCATE(correctedtime(gridinfov%arraysize))
  correctedtime(:) = gridinfov%array(1:gridinfov%arraysize)%time
  DO i = 1, gridinfov%arraysize
    correctedtime(i) = correctedtime(i) + &
     timedelta_new(sec=gridinfov%array(i)%timerange%p1)
  ENDDO
ENDIF
IF (ltime_definition == 2) THEN ! set all to analysis
  ALLOCATE(correctedtimerange(gridinfov%arraysize))
  correctedtimerange(:) = gridinfov%array(1:gridinfov%arraysize)%timerange
  correctedtimerange(:)%p1 = 0
ENDIF

DO i = 1, ngrid
  IF (PRESENT(anavar)) THEN
    j = count((this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
     .AND. .NOT.isanavar(:))
    IF (j <= 0) THEN
      CALL l4f_category_log(category, l4f_fatal, 'grid n.'//t2c(i)// &
       ' has only constant data, this is not allowed')
      CALL l4f_category_log(category, l4f_fatal, 'please check anavar argument')
      CALL raise_fatal_error()
    ENDIF
  ENDIF
  IF (ltime_definition == 1 .OR. ltime_definition == 2) THEN ! verification time
    ntime = count_distinct(correctedtime, &
     mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
     .AND. .NOT.isanavar(:), back=.true.)
  ELSE
    ntime = count_distinct(gridinfov%array(1:gridinfov%arraysize)%time, &
     mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
     .AND. .NOT.isanavar(:), back=.true.)
  ENDIF
  IF (ltime_definition == 2) THEN ! set all to analysis
    ntimerange = count_distinct(correctedtimerange, &
     mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
     .AND. .NOT.isanavar(:), back=.true.)
  ELSE
    ntimerange = count_distinct(gridinfov%array(1:gridinfov%arraysize)%timerange, &
     mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
     .AND. .NOT.isanavar(:), back=.true.)
  ENDIF
  nlevel = count_distinct(gridinfov%array(1:gridinfov%arraysize)%level, &
   mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
   back=.true.)
  nvar = count_distinct(gridinfov%array(1:gridinfov%arraysize)%var, &
   mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
   back=.true.)

#ifdef DEBUG
  CALL l4f_category_log(this(i)%category,l4f_debug,"alloc volgrid6d index: "//t2c(i))
#endif

  CALL volgrid6d_alloc(this(i),this(i)%griddim%dim,ntime=ntime, &
   ntimerange=ntimerange,nlevel=nlevel,nvar=nvar)

  IF (ltime_definition == 1 .OR. ltime_definition == 2) THEN ! verification time
    this(i)%time = pack_distinct(correctedtime, ntime, &
     mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
     .AND. .NOT.isanavar(:), back=.true.)
  ELSE
    this(i)%time = pack_distinct(gridinfov%array(1:gridinfov%arraysize)%time, ntime, &
     mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
     .AND. .NOT.isanavar(:), back=.true.)
  ENDIF
  CALL sort(this(i)%time)

  IF (ltime_definition == 2) THEN ! set all to analysis
    this(i)%timerange = pack_distinct(correctedtimerange, ntimerange, &
     mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
     .AND. .NOT.isanavar(:), back=.true.)
  ELSE
    this(i)%timerange = pack_distinct(gridinfov%array(1:gridinfov%arraysize)%timerange, &
     ntimerange, mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
     .AND. .NOT.isanavar(:), back=.true.)
  ENDIF
  CALL sort(this(i)%timerange)

  this(i)%level=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%level, &
   nlevel,mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
   back=.true.)
  CALL sort(this(i)%level)

  this(i)%var=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%var, nvar, &
   mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
   back=.true.)

#ifdef DEBUG
  CALL l4f_category_log(this(i)%category,l4f_debug,"alloc_vol volgrid6d index: "//t2c(i))
#endif
  CALL volgrid6d_alloc_vol(this(i), decode=decode)

ENDDO

IF (ltime_definition == 1 .OR. ltime_definition == 2) DEALLOCATE(correctedtime)
IF (ltime_definition == 2) DEALLOCATE(correctedtimerange)

DO i = 1, gridinfov%arraysize

#ifdef DEBUG
  CALL l4f_category_log(category,l4f_debug,"import from gridinfov index: "//t2c(i))
  CALL l4f_category_log(category,l4f_info, &
   "to volgrid6d index: "//t2c(index(this%griddim, gridinfov%array(i)%griddim)))
#endif

  CALL import(this(index(this%griddim, gridinfov%array(i)%griddim)), &
   gridinfov%array(i), dup_mode=dup_mode, clone=clone, isanavar=isanavar(i))

ENDDO

!chiudo il logger temporaneo
CALL l4f_category_delete(category)

END SUBROUTINE import_from_gridinfovv


SUBROUTINE export_to_gridinfov(this, gridinfov, gaid_template, clone)
TYPE(volgrid6d),INTENT(inout) :: this
TYPE(arrayof_gridinfo),INTENT(inout) :: gridinfov
TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template
LOGICAL,INTENT(in),OPTIONAL :: clone

INTEGER :: i ,itime, itimerange, ilevel, ivar
INTEGER :: ntime, ntimerange, nlevel, nvar
TYPE(gridinfo_def) :: gridinfol

#ifdef DEBUG
CALL l4f_category_log(this%category,l4f_debug,"start export_to_gridinfov")
#endif

! this is necessary in order not to repeatedly and uselessly copy the
! same array of coordinates for each 2d grid during export, the
! side-effect is that the computed projection in this is lost
CALL dealloc(this%griddim%dim)

i=0
ntime=size(this%time)
ntimerange=size(this%timerange)
nlevel=size(this%level)
nvar=size(this%var)

DO itime=1,ntime
  DO itimerange=1,ntimerange
    DO ilevel=1,nlevel
      DO ivar=1,nvar

        CALL init(gridinfol)
        CALL export(this, gridinfol, itime, itimerange, ilevel, ivar, &
         gaid_template=gaid_template, clone=clone)
        IF (c_e(gridinfol%gaid)) THEN
          CALL insert(gridinfov, gridinfol)
        ELSE
          CALL delete(gridinfol)
        ENDIF
        
      ENDDO
    ENDDO
  ENDDO
ENDDO

END SUBROUTINE export_to_gridinfov


SUBROUTINE export_to_gridinfovv(this, gridinfov, gaid_template, clone)
!, &
! categoryappend)
TYPE(volgrid6d),INTENT(inout) :: this(:)
TYPE(arrayof_gridinfo),INTENT(inout) :: gridinfov
TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template
LOGICAL,INTENT(in),OPTIONAL :: clone

INTEGER :: i

DO i = 1, SIZE(this)
#ifdef DEBUG
  CALL l4f_category_log(this(i)%category,l4f_debug, &
   "export_to_gridinfovv grid index: "//t2c(i))
#endif
  CALL export(this(i), gridinfov, gaid_template=gaid_template, clone=clone)
ENDDO

END SUBROUTINE export_to_gridinfovv


SUBROUTINE volgrid6d_import_from_file(this, filename, dup_mode, decode, &
 time_definition, anavar, categoryappend)
TYPE(volgrid6d),POINTER :: this(:)
CHARACTER(len=*),INTENT(in) :: filename
INTEGER,INTENT(in),OPTIONAL :: dup_mode
LOGICAL,INTENT(in),OPTIONAL :: decode
INTEGER,INTENT(IN),OPTIONAL :: time_definition
CHARACTER(len=*),INTENT(IN),OPTIONAL :: anavar(:)
character(len=*),INTENT(in),OPTIONAL :: categoryappend

TYPE(arrayof_gridinfo) :: gridinfo
TYPE(l4f_handle) :: category
CHARACTER(len=512) :: a_name

NULLIFY(this)

IF (PRESENT(categoryappend))THEN
  CALL l4f_launcher(a_name,a_name_append= &
   trim(subcategory)//"."//trim(categoryappend))
ELSE
  CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
ENDIF
category=l4f_category_get_handle(a_name)

CALL import(gridinfo, filename=filename, categoryappend=categoryappend)
  
IF (gridinfo%arraysize > 0) THEN

  CALL import(this, gridinfo, dup_mode=dup_mode, clone=.true., decode=decode, &
   time_definition=time_definition, anavar=anavar, &
   categoryappend=categoryappend)

  CALL l4f_category_log(category,l4f_info,"deleting gridinfo")
  CALL delete(gridinfo)

ELSE
  CALL l4f_category_log(category,l4f_info,"file does not contain gridded data")
ENDIF

! close logger
CALL l4f_category_delete(category)

END SUBROUTINE volgrid6d_import_from_file


SUBROUTINE volgrid6d_export_to_file(this, filename, gaid_template, categoryappend)
TYPE(volgrid6d) :: this(:)
CHARACTER(len=*),INTENT(in) :: filename
TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template
character(len=*),INTENT(in),OPTIONAL :: categoryappend

TYPE(arrayof_gridinfo) :: gridinfo
TYPE(l4f_handle) :: category
CHARACTER(len=512) :: a_name

IF (PRESENT(categoryappend)) THEN
  CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
ELSE
  CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
ENDIF
category=l4f_category_get_handle(a_name)

#ifdef DEBUG
CALL l4f_category_log(category,l4f_debug,"start export to file")
#endif

CALL l4f_category_log(category,l4f_info,"writing volgrid6d to grib file: "//trim(filename))

!IF (ASSOCIATED(this)) THEN
  CALL export(this, gridinfo, gaid_template=gaid_template, clone=.true.)
  IF (gridinfo%arraysize > 0) THEN
    CALL export(gridinfo, filename)
    CALL delete(gridinfo)
  ENDIF
!ELSE
!  CALL l4f_category_log(category,L4F_INFO,"volume volgrid6d is not associated")
!ENDIF

! close logger
CALL l4f_category_delete(category)

END SUBROUTINE volgrid6d_export_to_file


SUBROUTINE volgrid6dv_delete(this)
TYPE(volgrid6d),POINTER :: this(:)

INTEGER :: i

IF (ASSOCIATED(this)) THEN
  DO i = 1, SIZE(this)
#ifdef DEBUG
    CALL l4f_category_log(this(i)%category,l4f_debug, &
     "delete volgrid6d vector index: "//trim(to_char(i)))
#endif
    CALL delete(this(i))
  ENDDO
  DEALLOCATE(this)
ENDIF

END SUBROUTINE volgrid6dv_delete


! Internal method for performing grid to grid computations
SUBROUTINE volgrid6d_transform_compute(this, volgrid6d_in, volgrid6d_out, &
 lev_out, var_coord_vol, clone)
TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per il grigliato
type(volgrid6d), INTENT(in) :: volgrid6d_in ! oggetto da trasformare
type(volgrid6d), INTENT(inout) :: volgrid6d_out ! oggetto trasformato; deve essere completo (init, alloc, alloc_vol)
TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) ! vol7d_level object defining target vertical grid, for vertical interpolations
INTEGER,INTENT(in),OPTIONAL :: var_coord_vol ! index of variable defining vertical coordinate values in input volume
LOGICAL,INTENT(in),OPTIONAL :: clone ! se fornito e \c .TRUE., clona i gaid da volgrid6d_in a volgrid6d_out

INTEGER :: ntime, ntimerange, inlevel, onlevel, nvar, &
 itime, itimerange, ilevel, ivar, levshift, levused, lvar_coord_vol, spos
REAL,POINTER :: voldatiin(:,:,:), voldatiout(:,:,:), coord_3d_in(:,:,:)
TYPE(vol7d_level) :: output_levtype
LOGICAL,ALLOCATABLE :: zlist(:)


#ifdef DEBUG
call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_transform_compute")
#endif

ntime=0
ntimerange=0
inlevel=0
onlevel=0
nvar=0
lvar_coord_vol = optio_i(var_coord_vol)

if (associated(volgrid6d_in%time))then
  ntime=size(volgrid6d_in%time)
  volgrid6d_out%time=volgrid6d_in%time
end if

if (associated(volgrid6d_in%timerange))then
  ntimerange=size(volgrid6d_in%timerange)
  volgrid6d_out%timerange=volgrid6d_in%timerange
end if

IF (ASSOCIATED(volgrid6d_in%level))THEN
  inlevel=SIZE(volgrid6d_in%level)
ENDIF
IF (PRESENT(lev_out)) THEN
  onlevel=SIZE(lev_out)
  volgrid6d_out%level=lev_out
ELSE IF (ASSOCIATED(volgrid6d_in%level))THEN
  onlevel=SIZE(volgrid6d_in%level)
  volgrid6d_out%level=volgrid6d_in%level
ENDIF

if (associated(volgrid6d_in%var))then
  nvar=size(volgrid6d_in%var)
  volgrid6d_out%var=volgrid6d_in%var
end if
! allocate once for speed
IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
  ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, &
   inlevel))
ENDIF
IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
  ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, &
   onlevel))
ENDIF

CALL get_val(this, levshift=levshift, levused=levused)
spos = imiss
IF (c_e(lvar_coord_vol)) THEN
  CALL get_val(this%trans, output_levtype=output_levtype)
  IF (output_levtype%level1 == 103 .OR. output_levtype%level1 == 108) THEN
! here, unlike below, level (1,101) is not accepted because this is
! likely for pressure difference with surface 108, not 103
    spos = firsttrue(volgrid6d_in%level(:) == vol7d_level_new(1))
    IF (spos == 0) THEN
      CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
       'output level '//t2c(output_levtype%level1)// &
       ' requested, but height/press of surface not provided in volume')
    ENDIF
    IF (.NOT.c_e(levshift) .AND. .NOT.c_e(levused)) THEN
      CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
       'internal inconsistence, levshift and levused undefined when they should be')
    ENDIF
  ENDIF
ENDIF

DO ivar=1,nvar
!  IF (c_e(var_coord_vol)) THEN
!    IF (ivar == var_coord_vol) CYCLE ! skip coordinate variable in output
!  ENDIF
  DO itimerange=1,ntimerange
    DO itime=1,ntime
! skip empty columns where possible, improve
      IF (c_e(levshift) .AND. c_e(levused)) THEN
        IF (.NOT.any(c_e( &
         volgrid6d_in%gaid(levshift+1:levshift+levused,itime,itimerange,ivar) &
         ))) cycle
      ENDIF
      DO ilevel = 1, min(inlevel,onlevel)
! if present gaid copy it
        IF (c_e(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)) .AND. .NOT. &
         c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))) THEN

          IF (optio_log(clone)) THEN
            CALL copy(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar),&
             volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
#ifdef DEBUG
            CALL l4f_category_log(volgrid6d_in%category,l4f_debug, &
             "cloning gaid, level "//t2c(ilevel))
#endif
          ELSE
            volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = &
             volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)
          ENDIF
        ENDIF
      ENDDO
! if out level are more, we have to clone anyway
      DO ilevel = min(inlevel,onlevel) + 1, onlevel
        IF (c_e(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar)) .AND. .NOT. &
         c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))) then

          CALL copy(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar),&
           volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
#ifdef DEBUG
          CALL l4f_category_log(volgrid6d_in%category,l4f_debug, &
           "forced cloning gaid, level "//t2c(inlevel)//"->"//t2c(ilevel))
#endif
        ENDIF
      ENDDO

      IF (c_e(lvar_coord_vol)) THEN
        NULLIFY(coord_3d_in)
        CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, lvar_coord_vol, &
         coord_3d_in)
        IF (c_e(spos)) THEN ! compute difference wrt surface coordinate
          IF (spos == 0) THEN ! error condition, set all to missing and goodnight
            coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss
          ELSE
            DO ilevel = levshift+1, levshift+levused
              WHERE(c_e(coord_3d_in(:,:,ilevel)) .AND. c_e(coord_3d_in(:,:,spos)))
                coord_3d_in(:,:,ilevel) = coord_3d_in(:,:,ilevel) - &
                 coord_3d_in(:,:,spos)
              ELSEWHERE
                coord_3d_in(:,:,ilevel) = rmiss
              END WHERE
            ENDDO
          ENDIF
        ENDIF
      ENDIF
      CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, &
       voldatiin, zlist)
      IF (ASSOCIATED(volgrid6d_out%voldati)) & ! improve!!!!
       CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
       voldatiout)
      IF (c_e(lvar_coord_vol)) THEN ! this is a vertint case, zlist is useless
        CALL compute(this, voldatiin, voldatiout, var=convert(volgrid6d_in%var(ivar)), &
         coord_3d_in=coord_3d_in(:,:,levshift+1:levshift+levused)) ! subset coord_3d_in
      ELSE
        CALL compute(this, voldatiin, voldatiout, var=convert(volgrid6d_in%var(ivar)), zlist=zlist)
      ENDIF
      CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
       voldatiout)
    ENDDO
  ENDDO
ENDDO

IF (c_e(lvar_coord_vol)) THEN
  DEALLOCATE(coord_3d_in)
ENDIF
IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
  DEALLOCATE(voldatiin)
ENDIF
IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
  DEALLOCATE(voldatiout)
ENDIF


END SUBROUTINE volgrid6d_transform_compute


SUBROUTINE volgrid6d_transform(this, griddim, volgrid6d_in, volgrid6d_out, &
 lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend)
TYPE(transform_def),INTENT(in) :: this
TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim
TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in
TYPE(volgrid6d),INTENT(out) :: volgrid6d_out
TYPE(vol7d_level),INTENT(in),OPTIONAL,TARGET :: lev_out(:)
TYPE(volgrid6d),INTENT(in),OPTIONAL :: volgrid6d_coord_in
REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
REAL,INTENT(in),OPTIONAL :: maskbounds(:)
LOGICAL,INTENT(in),OPTIONAL :: clone
LOGICAL,INTENT(in),OPTIONAL :: decode
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

TYPE(grid_transform) :: grid_trans
TYPE(vol7d_level),POINTER :: llev_out(:)
TYPE(vol7d_level) :: input_levtype, output_levtype
TYPE(vol7d_var) :: vcoord_var
INTEGER :: i, k, ntime, ntimerange, nlevel, nvar, var_coord_in, var_coord_vol, &
 cf_out, nxc, nyc, nxi, nyi, i3, i4, i5, i6, &
 ulstart, ulend, spos
REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
TYPE(geo_proj) :: proj_in, proj_out
CHARACTER(len=80) :: trans_type
LOGICAL :: ldecode
LOGICAL,ALLOCATABLE :: mask_in(:)

#ifdef DEBUG
call l4f_category_log(volgrid6d_in%category, l4f_debug, "start volgrid6d_transform")
#endif

ntime=0
ntimerange=0
nlevel=0
nvar=0

if (associated(volgrid6d_in%time)) ntime=size(volgrid6d_in%time)
if (associated(volgrid6d_in%timerange)) ntimerange=size(volgrid6d_in%timerange)
if (associated(volgrid6d_in%level)) nlevel=size(volgrid6d_in%level)
if (associated(volgrid6d_in%var)) nvar=size(volgrid6d_in%var)

IF (ntime == 0 .OR. ntimerange == 0 .OR. nlevel == 0 .OR. nvar == 0) THEN
  CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
   "trying to transform an incomplete volgrid6d object, ntime="//t2c(ntime)// &
   ' ntimerange='//t2c(ntimerange)//' nlevel='//t2c(nlevel)//' nvar='//t2c(nvar))
  CALL init(volgrid6d_out) ! initialize to empty
  CALL raise_error()
  RETURN
ENDIF

CALL get_val(this, trans_type=trans_type)

! store desired output component flag and unrotate if necessary
cf_out = imiss
IF (PRESENT(griddim) .AND. (trans_type == 'inter' .OR. trans_type == 'boxinter' &
 .OR. trans_type == 'stencilinter')) THEN ! improve condition!!
  CALL get_val(volgrid6d_in%griddim, proj=proj_in)
  CALL get_val(griddim, component_flag=cf_out, proj=proj_out)
! if different projections wind components must be referred to geographical system
  IF (proj_in /= proj_out) CALL vg6d_wind_unrot(volgrid6d_in)
ELSE IF (PRESENT(griddim)) THEN ! just get component_flag, the rest is rubbish
  CALL get_val(griddim, component_flag=cf_out)
ENDIF


var_coord_in = imiss
var_coord_vol = imiss
IF (trans_type == 'vertint') THEN
  IF (PRESENT(lev_out)) THEN

! if volgrid6d_coord_in provided and allocated, check that it fits
    IF (PRESENT(volgrid6d_coord_in)) THEN
      IF (ASSOCIATED(volgrid6d_coord_in%voldati)) THEN

! strictly 1 time and 1 timerange
        IF (SIZE(volgrid6d_coord_in%voldati,4) /= 1 .OR. &
         SIZE(volgrid6d_coord_in%voldati,5) /= 1) THEN
          CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
           'volume providing constant input vertical coordinate must have &
           &only 1 time and 1 timerange')
          CALL init(volgrid6d_out) ! initialize to empty
          CALL raise_error()
          RETURN
        ENDIF

! search for variable providing vertical coordinate
        CALL get_val(this, output_levtype=output_levtype)
        vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
        IF (.NOT.c_e(vcoord_var)) THEN
          CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
           'requested output level type '//t2c(output_levtype%level1)// &
           ' does not correspond to any known physical variable for &
           &providing vertical coordinate')
          CALL init(volgrid6d_out) ! initialize to empty
          CALL raise_error()
          RETURN
        ENDIF

        DO i = 1, SIZE(volgrid6d_coord_in%var)
          IF (convert(volgrid6d_coord_in%var(i)) == vcoord_var) THEN
            var_coord_in = i
            EXIT
          ENDIF
        ENDDO

        IF (.NOT.c_e(var_coord_in)) THEN
          CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
           'volume providing constant input vertical coordinate contains no &
           &variables matching output level type '//t2c(output_levtype%level1))
          CALL init(volgrid6d_out) ! initialize to empty
          CALL raise_error()
          RETURN
        ENDIF
        CALL l4f_category_log(volgrid6d_in%category, l4f_info, &
         'Coordinate for vertint found in coord volume at position '// &
         t2c(var_coord_in))

! check horizontal grid
! this is too strict (component flag and so on)
!        IF (volgrid6d_coord_in%griddim /= volgrid6d_in%griddim) THEN
        CALL get_val(volgrid6d_coord_in%griddim, nx=nxc, ny=nyc)
        CALL get_val(volgrid6d_in%griddim, nx=nxi, ny=nyi)
        IF (nxc /= nxi .OR. nyc /= nyi) THEN
          CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
           'volume providing constant input vertical coordinate must have &
           &the same grid as the input')
          CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
           'vertical coordinate: '//t2c(nxc)//'x'//t2c(nyc)// &
           ', input volume: '//t2c(nxi)//'x'//t2c(nyi))
          CALL init(volgrid6d_out) ! initialize to empty
          CALL raise_error()
          RETURN
        ENDIF

! check vertical coordinate system
        CALL get_val(this, input_levtype=input_levtype)
        mask_in = & ! implicit allocation
         (volgrid6d_coord_in%level(:)%level1 == input_levtype%level1) .AND. &
         (volgrid6d_coord_in%level(:)%level2 == input_levtype%level2)
        ulstart = firsttrue(mask_in)
        ulend = lasttrue(mask_in)
        IF (ulstart == 0 .OR. ulend == 0) THEN
          CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
           'coordinate file does not contain levels of type '// &
           t2c(input_levtype%level1)//'/'//t2c(input_levtype%level2)// &
           ' specified for input data')
          CALL init(volgrid6d_out) ! initialize to empty
          CALL raise_error()
          RETURN
        ENDIF

        coord_3d_in = volgrid6d_coord_in%voldati(:,:,ulstart:ulend,1,1,var_coord_in) ! implicit allocation
! special case
        IF (output_levtype%level1 == 103 .OR. &
         output_levtype%level1 == 108) THEN ! surface coordinate needed
! allow coding of orography level as in Icon
          spos = firsttrue(volgrid6d_coord_in%level(:) == vol7d_level_new(level1=1) .OR. &
           volgrid6d_coord_in%level(:) == vol7d_level_new(level1=1, level2=101))
          IF (spos == 0) THEN
            CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
             'output level '//t2c(output_levtype%level1)// &
             ' requested, but height/press of surface not provided in coordinate file')
            CALL init(volgrid6d_out) ! initialize to empty
            CALL raise_error()
            RETURN
          ENDIF
          DO k = 1, SIZE(coord_3d_in,3)
            WHERE(c_e(coord_3d_in(:,:,k)) .AND. &
             c_e(volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in)))
              coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - &
               volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in)
            ELSEWHERE
              coord_3d_in(:,:,k) = rmiss
            END WHERE
          ENDDO
        ENDIF

      ENDIF
    ENDIF

    IF (.NOT.c_e(var_coord_in)) THEN ! search for coordinate within volume
! search for variable providing vertical coordinate
      CALL get_val(this, output_levtype=output_levtype)
      vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
      IF (c_e(vcoord_var)) THEN
        DO i = 1, SIZE(volgrid6d_in%var)
          IF (convert(volgrid6d_in%var(i)) == vcoord_var) THEN
            var_coord_vol = i
            EXIT
          ENDIF
        ENDDO

        IF (c_e(var_coord_vol)) THEN
          CALL l4f_category_log(volgrid6d_in%category, l4f_info, &
           'Coordinate for vertint found in input volume at position '// &
           t2c(var_coord_vol))
        ENDIF

      ENDIF
    ENDIF

    CALL init(volgrid6d_out, griddim=volgrid6d_in%griddim, &
     time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend)
    IF (c_e(var_coord_in)) THEN
      CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, &
       coord_3d_in=coord_3d_in, categoryappend=categoryappend)
    ELSE
      CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, &
       categoryappend=categoryappend)
    ENDIF

    CALL get_val(grid_trans, output_level_auto=llev_out) ! get levels if auto-generated
    IF (.NOT.ASSOCIATED(llev_out)) llev_out => lev_out
    nlevel = SIZE(llev_out)
  ELSE
    CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
     'volgrid6d_transform: vertint requested but lev_out not provided')
    CALL init(volgrid6d_out) ! initialize to empty
    CALL raise_error()
    RETURN
  ENDIF

ELSE
  CALL init(volgrid6d_out, griddim=griddim, &
   time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend)
  CALL init(grid_trans, this, in=volgrid6d_in%griddim, out=volgrid6d_out%griddim, &
   maskgrid=maskgrid, maskbounds=maskbounds, categoryappend=categoryappend)
ENDIF


IF (c_e(grid_trans)) THEN ! transformation is valid

  CALL volgrid6d_alloc(volgrid6d_out, ntime=ntime, nlevel=nlevel, &
   ntimerange=ntimerange, nvar=nvar)

  IF (PRESENT(decode)) THEN ! explicitly set decode status
    ldecode = decode
  ELSE ! take it from input
    ldecode = ASSOCIATED(volgrid6d_in%voldati)
  ENDIF
! force decode if gaid is readonly
  decode_loop: DO i6 = 1,nvar
    DO i5 = 1, ntimerange
      DO i4 = 1, ntime
        DO i3 = 1, nlevel
          IF (c_e(volgrid6d_in%gaid(i3,i4,i5,i6))) THEN
            ldecode = ldecode .OR. grid_id_readonly(volgrid6d_in%gaid(i3,i4,i5,i6))
            EXIT decode_loop
          ENDIF
        ENDDO
      ENDDO
    ENDDO
  ENDDO decode_loop

  IF (PRESENT(decode)) THEN
    IF (ldecode.NEQV.decode) THEN
      CALL l4f_category_log(volgrid6d_in%category, l4f_warn, &
     'volgrid6d_transform: decode status forced to .TRUE. because driver does not allow copy')
    ENDIF
  ENDIF

  CALL volgrid6d_alloc_vol(volgrid6d_out, decode=ldecode)

!ensure unproj was called
!call griddim_unproj(volgrid6d_out%griddim)

  IF (trans_type == 'vertint') THEN
#ifdef DEBUG
    CALL l4f_category_log(volgrid6d_in%category, l4f_debug, &
     "volgrid6d_transform: vertint to "//t2c(nlevel)//" levels")
#endif
    CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, lev_out=llev_out, &
     var_coord_vol=var_coord_vol, clone=clone)
  ELSE
    CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, clone=clone)
  ENDIF

  IF (cf_out == 0) THEN ! unrotated components are desired
    CALL wind_unrot(volgrid6d_out) ! unrotate if necessary
  ELSE IF (cf_out == 1) THEN ! rotated components are desired
    CALL wind_rot(volgrid6d_out) ! rotate if necessary
  ENDIF

ELSE
! should log with grid_trans%category, but it is private
  CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
   'volgrid6d_transform: transformation not valid')
  CALL raise_error()
ENDIF

CALL delete (grid_trans)

END SUBROUTINE volgrid6d_transform


SUBROUTINE volgrid6dv_transform(this, griddim, volgrid6d_in, volgrid6d_out, &
 lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend)
TYPE(transform_def),INTENT(in) :: this
TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim
TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in(:)
TYPE(volgrid6d),POINTER :: volgrid6d_out(:)
TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:)
TYPE(volgrid6d),INTENT(in),OPTIONAL :: volgrid6d_coord_in
REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
REAL,INTENT(in),OPTIONAL :: maskbounds(:)
LOGICAL,INTENT(in),OPTIONAL :: clone
LOGICAL,INTENT(in),OPTIONAL :: decode
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

INTEGER :: i, stallo


allocate(volgrid6d_out(size(volgrid6d_in)),stat=stallo)
if (stallo /= 0)then
  call l4f_log(l4f_fatal,"allocating memory")
  call raise_fatal_error()
end if

do i=1,size(volgrid6d_in)
  call transform(this, griddim, volgrid6d_in(i), volgrid6d_out(i), &
   lev_out=lev_out, volgrid6d_coord_in=volgrid6d_coord_in, &
   maskgrid=maskgrid, maskbounds=maskbounds, &
   clone=clone, decode=decode, categoryappend=categoryappend)
end do

END SUBROUTINE volgrid6dv_transform


! Internal method for performing grid to sparse point computations
SUBROUTINE volgrid6d_v7d_transform_compute(this, volgrid6d_in, vol7d_out, &
 networkname, noconvert)
TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per grigliato
type(volgrid6d), INTENT(in) :: volgrid6d_in ! oggetto da trasformare
type(vol7d), INTENT(inout) :: vol7d_out ! oggetto trasformato
CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname ! imposta il network in vol7d_out (default='generic')
LOGICAL,OPTIONAL,INTENT(in) :: noconvert

INTEGER :: nntime, nana, ntime, ntimerange, nlevel, nvar, stallo
INTEGER :: itime, itimerange, ivar, inetwork
REAL,ALLOCATABLE :: voldatir_out(:,:,:)
TYPE(conv_func),POINTER :: c_func(:)
TYPE(datetime),ALLOCATABLE :: validitytime(:,:)
REAL,POINTER :: voldatiin(:,:,:)

#ifdef DEBUG
call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_v7d_transform_compute")
#endif

ntime=0
ntimerange=0
nlevel=0
nvar=0
NULLIFY(c_func)

if (present(networkname))then
  call init(vol7d_out%network(1),name=networkname)
else
  call init(vol7d_out%network(1),name='generic')
end if

if (associated(volgrid6d_in%timerange))then
  ntimerange=size(volgrid6d_in%timerange)
  vol7d_out%timerange=volgrid6d_in%timerange
end if

if (associated(volgrid6d_in%time))then
  ntime=size(volgrid6d_in%time)

  if (vol7d_out%time_definition == volgrid6d_in%time_definition) then

                                ! i time sono definiti uguali: assegno
    vol7d_out%time=volgrid6d_in%time

  else
                                ! converto reference in validity
    allocate (validitytime(ntime,ntimerange),stat=stallo)
    if (stallo /=0)then
      call l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory")
      call raise_fatal_error()
    end if

    do itime=1,ntime
      do itimerange=1,ntimerange
        if (vol7d_out%time_definition > volgrid6d_in%time_definition) then
          validitytime(itime,itimerange) = &
           volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
        else
          validitytime(itime,itimerange) = &
           volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
        end if
      end do
    end do

    nntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.true.)
    vol7d_out%time=pack_distinct(reshape(validitytime,(/ntime*ntimerange/)), nntime,back=.true.)

  end if
end if

IF (ASSOCIATED(volgrid6d_in%level)) THEN
  nlevel = SIZE(volgrid6d_in%level)
  vol7d_out%level=volgrid6d_in%level
ENDIF

IF (ASSOCIATED(volgrid6d_in%var)) THEN
  nvar = SIZE(volgrid6d_in%var)
  IF (.NOT. optio_log(noconvert)) THEN
    CALL vargrib2varbufr(volgrid6d_in%var, vol7d_out%dativar%r, c_func)
  ENDIF
ENDIF

nana = SIZE(vol7d_out%ana)

! allocate once for speed
IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
  ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, &
   nlevel))
ENDIF

ALLOCATE(voldatir_out(nana,1,nlevel),stat=stallo)
IF (stallo /= 0) THEN
  CALL l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory")
  CALL raise_fatal_error()
ENDIF

inetwork=1
do itime=1,ntime
  do itimerange=1,ntimerange
!    do ilevel=1,nlevel
      do ivar=1,nvar
        
                                !non è chiaro se questa sezione è utile o no
                                !ossia il compute sotto sembra prevedere voldatir_out solo in out
!!$        if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
!!$          voldatir_out=reshape(vol7d_out%voldatir(:,itime,ilevel,itimerange,ivar,inetwork),(/nana,1/))
!!$        else
!!$          voldatir_out=reshape(vol7d_out%voldatir(:,index(vol7d_out%time,validitytime(itime,ilevel)),ilevel,itimerange,ivar,inetwork),(/nana,1/))
!!$        end if

        CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, &
         voldatiin)

        CALL compute(this, voldatiin, voldatir_out, vol7d_out%dativar%r(ivar))

        if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
          vol7d_out%voldatir(:,itime,:,itimerange,ivar,inetwork) = &
           voldatir_out(:,1,:)
        else
          vol7d_out%voldatir(:,index(vol7d_out%time,validitytime(itime,itimerange)),:,itimerange,ivar,inetwork)=&
           reshape(voldatir_out,(/nana,nlevel/))
        end if

! 1 indice della dimensione "anagrafica"
! 2 indice della dimensione "tempo"
! 3 indice della dimensione "livello verticale"
! 4 indice della dimensione "intervallo temporale"
! 5 indice della dimensione "variabile"
! 6 indice della dimensione "rete"

      end do
!    end do
  end do
end do

deallocate(voldatir_out)
IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
  DEALLOCATE(voldatiin)
ENDIF
if (allocated(validitytime)) deallocate(validitytime)

! Rescale valid data according to variable conversion table
IF (ASSOCIATED(c_func)) THEN
  DO ivar = 1, nvar
    CALL compute(c_func(ivar), vol7d_out%voldatir(:,:,:,:,ivar,:))
  ENDDO
  DEALLOCATE(c_func)
ENDIF

end SUBROUTINE volgrid6d_v7d_transform_compute


SUBROUTINE volgrid6d_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, &
 maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend)
TYPE(transform_def),INTENT(in) :: this
TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in
TYPE(vol7d),INTENT(out) :: vol7d_out
TYPE(vol7d),INTENT(in),OPTIONAL :: v7d
REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
REAL,INTENT(in),OPTIONAL :: maskbounds(:)
CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname
LOGICAL,OPTIONAL,INTENT(in) :: noconvert
PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

type(grid_transform) :: grid_trans
INTEGER :: ntime, ntimerange, nlevel, nvar, nana, time_definition, nnetwork, stallo
INTEGER :: itime, itimerange, inetwork
TYPE(datetime),ALLOCATABLE :: validitytime(:,:)
INTEGER,ALLOCATABLE :: point_index(:)
TYPE(vol7d) :: v7d_locana

#ifdef DEBUG
call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_v7d_transform")
#endif

call vg6d_wind_unrot(volgrid6d_in)

ntime=0
ntimerange=0
nlevel=0
nvar=0
nnetwork=1

call get_val(this,time_definition=time_definition)
if (.not. c_e(time_definition)) then
  time_definition=1  ! default to validity time
endif

IF (PRESENT(v7d)) THEN
  CALL vol7d_copy(v7d, v7d_locana)
ELSE
  CALL init(v7d_locana)
ENDIF

if (associated(volgrid6d_in%timerange)) ntimerange=size(volgrid6d_in%timerange)

if (associated(volgrid6d_in%time)) then

  ntime=size(volgrid6d_in%time)
  
  if (time_definition /= volgrid6d_in%time_definition) then
    
                                ! converto reference in validity
    allocate (validitytime(ntime,ntimerange),stat=stallo)
    if (stallo /=0)then
      call l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory")
      call raise_fatal_error()
    end if

    do itime=1,ntime
      do itimerange=1,ntimerange
        if (time_definition > volgrid6d_in%time_definition) then
          validitytime(itime,itimerange) = &
           volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
        else
          validitytime(itime,itimerange) = &
           volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
        end if
      end do
    end do

    ntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.true.)
    deallocate (validitytime)
  
  end if
end if


if (associated(volgrid6d_in%level)) nlevel=size(volgrid6d_in%level)
if (associated(volgrid6d_in%var)) nvar=size(volgrid6d_in%var)

CALL init(grid_trans, this, volgrid6d_in%griddim, v7d_locana, &
 maskgrid=maskgrid, maskbounds=maskbounds, find_index=find_index, &
 categoryappend=categoryappend)
CALL init (vol7d_out,time_definition=time_definition)

IF (c_e(grid_trans)) THEN

  nana=SIZE(v7d_locana%ana)
  CALL vol7d_alloc(vol7d_out, nana=nana, ntime=ntime, nlevel=nlevel, &
   ntimerange=ntimerange, ndativarr=nvar, nnetwork=nnetwork)
  vol7d_out%ana = v7d_locana%ana

  CALL get_val(grid_trans, output_point_index=point_index)
  IF (ALLOCATED(point_index)) THEN
! check that size(point_index) == nana?
    CALL vol7d_alloc(vol7d_out, nanavari=1)
    CALL init(vol7d_out%anavar%i(1), 'B01192')
  ENDIF

  CALL vol7d_alloc_vol(vol7d_out)

  IF (ALLOCATED(point_index)) THEN
    DO inetwork = 1, nnetwork
      vol7d_out%volanai(:,1,inetwork) = point_index(:)
    ENDDO
  ENDIF
  CALL compute(grid_trans, volgrid6d_in, vol7d_out, networkname, noconvert)
ELSE
  CALL l4f_log(l4f_error, 'vg6d_v7d_transform: transformation not valid')
  CALL raise_error()
ENDIF

CALL delete(grid_trans)

#ifdef HAVE_DBALLE
! set variables to a conformal state
CALL vol7d_dballe_set_var_du(vol7d_out)
#endif

CALL delete(v7d_locana)

END SUBROUTINE volgrid6d_v7d_transform


SUBROUTINE volgrid6dv_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, &
 maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend)
TYPE(transform_def),INTENT(in) :: this
TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in(:)
TYPE(vol7d),INTENT(out) :: vol7d_out
TYPE(vol7d),INTENT(in),OPTIONAL :: v7d
REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
REAL,INTENT(in),OPTIONAL :: maskbounds(:)
CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname
LOGICAL,OPTIONAL,INTENT(in) :: noconvert
PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

integer :: i
TYPE(vol7d) :: v7dtmp


CALL init(v7dtmp)
CALL init(vol7d_out)

DO i=1,SIZE(volgrid6d_in)
  CALL transform(this, volgrid6d_in(i), v7dtmp, v7d=v7d, &
   maskgrid=maskgrid, maskbounds=maskbounds, &
   networkname=networkname, noconvert=noconvert, find_index=find_index, &
   categoryappend=categoryappend)
  CALL vol7d_append(vol7d_out, v7dtmp)
ENDDO

END SUBROUTINE volgrid6dv_v7d_transform


! Internal method for performing sparse point to grid computations
SUBROUTINE v7d_volgrid6d_transform_compute(this, vol7d_in, volgrid6d_out, networkname, gaid_template)
TYPE(grid_transform),INTENT(in) :: this ! object specifying the specific transformation
type(vol7d), INTENT(in) :: vol7d_in ! object to be transformed
type(volgrid6d), INTENT(inout) :: volgrid6d_out ! transformed object
CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname ! select the network to be processed from the \a vol7d input object, default the first network
TYPE(grid_id),OPTIONAL,INTENT(in) :: gaid_template ! the template (typically grib_api) to be associated with output data, it also helps in improving variable conversion

integer :: nana, ntime, ntimerange, nlevel, nvar
INTEGER :: ilevel, itime, itimerange, ivar, inetwork

REAL,POINTER :: voldatiout(:,:,:)
type(vol7d_network) :: network
TYPE(conv_func), pointer :: c_func(:)
!TODO category sarebbe da prendere da vol7d
#ifdef DEBUG
CALL l4f_category_log(volgrid6d_out%category, l4f_debug, &
 'start v7d_volgrid6d_transform_compute')
#endif

ntime=0
ntimerange=0
nlevel=0
nvar=0

IF (PRESENT(networkname)) THEN
  CALL init(network,name=networkname)
  inetwork = index(vol7d_in%network,network)
  IF (inetwork <= 0) THEN
    CALL l4f_category_log(volgrid6d_out%category,l4f_warn, &
     'network '//trim(networkname)//' not found, first network will be transformed')
    inetwork = 1
  ENDIF
ELSE
  inetwork = 1
ENDIF

! no time_definition conversion implemented, output will be the same as input
if (associated(vol7d_in%time))then
  ntime=size(vol7d_in%time)
  volgrid6d_out%time=vol7d_in%time
end if

if (associated(vol7d_in%timerange))then
  ntimerange=size(vol7d_in%timerange)
  volgrid6d_out%timerange=vol7d_in%timerange
end if

if (associated(vol7d_in%level))then
  nlevel=size(vol7d_in%level)
  volgrid6d_out%level=vol7d_in%level
end if

if (associated(vol7d_in%dativar%r))then
  nvar=size(vol7d_in%dativar%r)
  CALL varbufr2vargrib(vol7d_in%dativar%r, volgrid6d_out%var, c_func, gaid_template)
end if

nana=SIZE(vol7d_in%voldatir, 1)
! allocate once for speed
IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
  ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, &
   nlevel))
ENDIF

DO ivar=1,nvar
  DO itimerange=1,ntimerange
    DO itime=1,ntime

! clone the gaid template where I have data
      IF (PRESENT(gaid_template)) THEN
        DO ilevel = 1, nlevel
          IF (any(c_e(vol7d_in%voldatir(:,itime,ilevel,itimerange,ivar,inetwork)))) THEN
            CALL copy(gaid_template, volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
          ELSE
            volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = grid_id_new()
          ENDIF
        ENDDO
      ENDIF

! get data
      IF (ASSOCIATED(volgrid6d_out%voldati)) & ! improve!!!!
       CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
       voldatiout)
! do the interpolation
      CALL compute(this, &
       vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), voldatiout, &
       vol7d_in%dativar%r(ivar))
! rescale valid data according to variable conversion table
      IF (ASSOCIATED(c_func)) THEN
        CALL compute(c_func(ivar), voldatiout(:,:,:))
      ENDIF
! put data
      CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
       voldatiout)

    ENDDO
  ENDDO
ENDDO

IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
  DEALLOCATE(voldatiout)
ENDIF
IF (ASSOCIATED(c_func)) THEN
  DEALLOCATE(c_func)
ENDIF

END SUBROUTINE v7d_volgrid6d_transform_compute


SUBROUTINE v7d_volgrid6d_transform(this, griddim, vol7d_in, volgrid6d_out, &
 networkname, gaid_template, categoryappend)
TYPE(transform_def),INTENT(in) :: this
TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim
! TODO ripristinare intent(in) dopo le opportune modifiche in grid_class.F90
TYPE(vol7d),INTENT(inout) :: vol7d_in
TYPE(volgrid6d),INTENT(out) :: volgrid6d_out
CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname
TYPE(grid_id),OPTIONAL,INTENT(in) :: gaid_template
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

type(grid_transform) :: grid_trans
integer :: ntime, ntimerange, nlevel, nvar


!TODO la category sarebbe da prendere da vol7d
!call l4f_category_log(vol7d_out%category,L4F_DEBUG,"start volgrid6d_transform")

CALL vol7d_alloc_vol(vol7d_in) ! be safe
ntime=SIZE(vol7d_in%time)
ntimerange=SIZE(vol7d_in%timerange)
nlevel=SIZE(vol7d_in%level)
nvar=0
if (associated(vol7d_in%dativar%r)) nvar=size(vol7d_in%dativar%r)

IF (nvar <= 0) THEN ! use vol7d category once it will be implemented
  CALL l4f_log(l4f_error, &
   "trying to transform a vol7d object incomplete or without real variables")
  CALL init(volgrid6d_out) ! initialize to empty
  CALL raise_error()
  RETURN
ENDIF

CALL init(grid_trans, this, vol7d_in, griddim, categoryappend=categoryappend)
CALL init(volgrid6d_out, griddim, time_definition=vol7d_in%time_definition, & 
 categoryappend=categoryappend)

IF (c_e(grid_trans)) THEN

  CALL volgrid6d_alloc(volgrid6d_out, griddim%dim, ntime=ntime, nlevel=nlevel, &
   ntimerange=ntimerange, nvar=nvar)
! can I avoid decode=.TRUE. here, using gaid_template?
  CALL volgrid6d_alloc_vol(volgrid6d_out, decode=.true.)

  CALL compute(grid_trans, vol7d_in, volgrid6d_out, networkname, gaid_template)

  CALL vg6d_wind_rot(volgrid6d_out)
ELSE
! should log with grid_trans%category, but it is private
  CALL l4f_log(l4f_error, 'v7d_vg6d_transform: transformation not valid')
  CALL raise_error()
ENDIF

CALL delete(grid_trans)

END SUBROUTINE v7d_volgrid6d_transform


! Internal method for performing sparse point to sparse point computations
SUBROUTINE v7d_v7d_transform_compute(this, vol7d_in, vol7d_out, lev_out, &
 var_coord_vol)
TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per grigliato
type(vol7d), INTENT(in) :: vol7d_in ! oggetto da trasformare
type(vol7d), INTENT(inout) :: vol7d_out ! oggetto trasformato
TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) ! vol7d_level object defining target vertical grid, for vertical interpolations
INTEGER,INTENT(in),OPTIONAL :: var_coord_vol ! index of variable defining vertical coordinate values in input volume

INTEGER :: itime, itimerange, ilevel, ivar, inetwork, &
 levshift, levused, lvar_coord_vol, spos
REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
TYPE(vol7d_level) :: output_levtype

lvar_coord_vol = optio_i(var_coord_vol)
vol7d_out%time(:) = vol7d_in%time(:)
vol7d_out%timerange(:) = vol7d_in%timerange(:)
IF (PRESENT(lev_out)) THEN
  vol7d_out%level(:) = lev_out(:)
ELSE
  vol7d_out%level(:) = vol7d_in%level(:)
ENDIF
vol7d_out%network(:) = vol7d_in%network(:)
IF (ASSOCIATED(vol7d_in%dativar%r)) THEN ! work only when real vars are available
  vol7d_out%dativar%r(:) = vol7d_in%dativar%r(:)

  CALL get_val(this, levshift=levshift, levused=levused)
  spos = imiss
  IF (c_e(lvar_coord_vol)) THEN
    CALL get_val(this%trans, output_levtype=output_levtype)
    IF (output_levtype%level1 == 103 .OR. output_levtype%level1 == 108) THEN
      spos = firsttrue(vol7d_in%level(:) == vol7d_level_new(1))
      IF (spos == 0) THEN
        CALL l4f_log(l4f_error, &
         'output level '//t2c(output_levtype%level1)// &
         ' requested, but height/press of surface not provided in volume')
      ENDIF
      IF (.NOT.c_e(levshift) .AND. .NOT.c_e(levused)) THEN
        CALL l4f_log(l4f_error, &
         'internal inconsistence, levshift and levused undefined when they should be')
      ENDIF
      ALLOCATE(coord_3d_in(SIZE(vol7d_in%ana),1,SIZE(vol7d_in%level)))
    ENDIF

  ENDIF

  DO inetwork = 1, SIZE(vol7d_in%network)
    DO ivar = 1, SIZE(vol7d_in%dativar%r)
      DO itimerange = 1, SIZE(vol7d_in%timerange)
        DO itime = 1, SIZE(vol7d_in%time)

! dirty trick to make voldatir look like a 2d-array of shape (nana,1)
          IF (c_e(lvar_coord_vol)) THEN
            IF (c_e(spos)) THEN ! compute difference wrt surface coordinate
              IF (spos == 0) THEN ! error condition, set all to missing and goodnight
                coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss
              ELSE
                DO ilevel = levshift+1, levshift+levused
                  WHERE(c_e(vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork)) .AND. &
                  c_e(vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork))) 
                    coord_3d_in(:,:,ilevel) = vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork) - &
                     vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork)
                  ELSEWHERE
                    coord_3d_in(:,:,ilevel) = rmiss
                  END WHERE
                ENDDO
              ENDIF
              CALL compute(this, &
               vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
               vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
               var=vol7d_in%dativar%r(ivar), &
               coord_3d_in=coord_3d_in)
            ELSE
              CALL compute(this, &
               vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
               vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
               var=vol7d_in%dativar%r(ivar), &
               coord_3d_in=vol7d_in%voldatir(:,itime:itime,:,itimerange, &
               lvar_coord_vol,inetwork))
            ENDIF
          ELSE
            CALL compute(this, &
             vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
             vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
             var=vol7d_in%dativar%r(ivar))
          ENDIF
        ENDDO
      ENDDO
    ENDDO
  ENDDO

ENDIF

END SUBROUTINE v7d_v7d_transform_compute


SUBROUTINE v7d_v7d_transform(this, vol7d_in, vol7d_out, v7d, maskbounds, &
 lev_out, vol7d_coord_in, categoryappend)
TYPE(transform_def),INTENT(in) :: this
TYPE(vol7d),INTENT(inout) :: vol7d_in
TYPE(vol7d),INTENT(out) :: vol7d_out
TYPE(vol7d),INTENT(in),OPTIONAL :: v7d
REAL,INTENT(in),OPTIONAL :: maskbounds(:)
TYPE(vol7d_level),INTENT(in),OPTIONAL,TARGET :: lev_out(:)
TYPE(vol7d),INTENT(in),OPTIONAL :: vol7d_coord_in
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

INTEGER :: nvar, inetwork
TYPE(grid_transform) :: grid_trans
TYPE(vol7d_level),POINTER :: llev_out(:)
TYPE(vol7d_level) :: input_levtype, output_levtype
TYPE(vol7d_var) :: vcoord_var
REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
INTEGER :: var_coord_in, var_coord_vol, i, k, ulstart, ulend, spos
INTEGER,ALLOCATABLE :: point_index(:)
TYPE(vol7d) :: v7d_locana, vol7d_tmpana
CHARACTER(len=80) :: trans_type
LOGICAL,ALLOCATABLE :: mask_in(:), point_mask(:)

CALL vol7d_alloc_vol(vol7d_in) ! be safe
nvar=0
IF (ASSOCIATED(vol7d_in%dativar%r)) nvar=SIZE(vol7d_in%dativar%r)

CALL init(v7d_locana)
IF (PRESENT(v7d)) v7d_locana = v7d
CALL init(vol7d_out, time_definition=vol7d_in%time_definition)

CALL get_val(this, trans_type=trans_type)

var_coord_vol = imiss
IF (trans_type == 'vertint') THEN

  IF (PRESENT(lev_out)) THEN

! if vol7d_coord_in provided and allocated, check that it fits
    var_coord_in = -1
    IF (PRESENT(vol7d_coord_in)) THEN
      IF (ASSOCIATED(vol7d_coord_in%voldatir) .AND. &
       ASSOCIATED(vol7d_coord_in%dativar%r)) THEN

! strictly 1 time, 1 timerange and 1 network
        IF (SIZE(vol7d_coord_in%voldatir,2) /= 1 .OR. &
         SIZE(vol7d_coord_in%voldatir,4) /= 1 .OR. &
         SIZE(vol7d_coord_in%voldatir,6) /= 1) THEN
          CALL l4f_log(l4f_error, &
           'volume providing constant input vertical coordinate must have &
           &only 1 time, 1 timerange and 1 network')
          CALL raise_error()
          RETURN
        ENDIF

! search for variable providing vertical coordinate
        CALL get_val(this, output_levtype=output_levtype)
        vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
        IF (.NOT.c_e(vcoord_var)) THEN
          CALL l4f_log(l4f_error, &
           'requested output level type '//t2c(output_levtype%level1)// &
           ' does not correspond to any known physical variable for &
           &providing vertical coordinate')
          CALL raise_error()
          RETURN
        ENDIF

        var_coord_in = index(vol7d_coord_in%dativar%r, vcoord_var)

        IF (var_coord_in <= 0) THEN
          CALL l4f_log(l4f_error, &
           'volume providing constant input vertical coordinate contains no &
           &real variables matching output level type '//t2c(output_levtype%level1))
          CALL raise_error()
          RETURN
        ENDIF
        CALL l4f_log(l4f_info, &
         'Coordinate for vertint found in coord volume at position '// &
         t2c(var_coord_in))

! check vertical coordinate system
        CALL get_val(this, input_levtype=input_levtype)
        mask_in = & ! implicit allocation
         (vol7d_coord_in%level(:)%level1 == input_levtype%level1) .AND. &
         (vol7d_coord_in%level(:)%level2 == input_levtype%level2)
        ulstart = firsttrue(mask_in)
        ulend = lasttrue(mask_in)
        IF (ulstart == 0 .OR. ulend == 0) THEN
          CALL l4f_log(l4f_error, &
           'coordinate file does not contain levels of type '// &
           t2c(input_levtype%level1)//'/'//t2c(input_levtype%level2)// &
           ' specified for input data')
          CALL raise_error()
          RETURN
        ENDIF

        coord_3d_in = vol7d_coord_in%voldatir(:,1:1,ulstart:ulend,1,var_coord_in,1) ! dirty 1:1, implicit allocation
! special case
        IF (output_levtype%level1 == 103 &
         .OR. output_levtype%level1 == 108) THEN ! surface coordinate needed
          spos = firsttrue(vol7d_coord_in%level(:) == vol7d_level_new(1))
          IF (spos == 0) THEN
            CALL l4f_log(l4f_error, &
             'output level '//t2c(output_levtype%level1)// &
             ' requested, but height/press of surface not provided in coordinate file')
            CALL raise_error()
            RETURN
          ENDIF
          DO k = 1, SIZE(coord_3d_in,3)
            WHERE(c_e(coord_3d_in(:,:,k)) .AND. &
             c_e(vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1)))
              coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - &
               vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1)
            ELSEWHERE
              coord_3d_in(:,:,k) = rmiss
            END WHERE
          ENDDO
        ENDIF

      ENDIF
    ENDIF

    IF (var_coord_in <= 0) THEN ! search for coordinate within volume
! search for variable providing vertical coordinate
      CALL get_val(this, output_levtype=output_levtype)
      vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
      IF (c_e(vcoord_var)) THEN
        DO i = 1, SIZE(vol7d_in%dativar%r)
          IF (vol7d_in%dativar%r(i) == vcoord_var) THEN
            var_coord_vol = i
            EXIT
          ENDIF
        ENDDO

        IF (c_e(var_coord_vol)) THEN
          CALL l4f_log(l4f_info, &
           'Coordinate for vertint found in input volume at position '// &
           t2c(var_coord_vol))
        ENDIF

      ENDIF
    ENDIF

    IF (var_coord_in > 0) THEN
      CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, &
       coord_3d_in=coord_3d_in, categoryappend=categoryappend)
    ELSE
      CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, &
       categoryappend=categoryappend)
    ENDIF

    CALL get_val(grid_trans, output_level_auto=llev_out) ! get levels if auto-generated
    IF (.NOT.associated(llev_out)) llev_out => lev_out

    IF (c_e(grid_trans)) then! .AND. nvar > 0) THEN

      CALL vol7d_alloc(vol7d_out, nana=SIZE(vol7d_in%ana), &
       ntime=SIZE(vol7d_in%time), ntimerange=SIZE(vol7d_in%timerange), &
       nlevel=SIZE(llev_out), nnetwork=SIZE(vol7d_in%network), ndativarr=nvar)
      vol7d_out%ana(:) = vol7d_in%ana(:)

      CALL vol7d_alloc_vol(vol7d_out)

! no need to check c_e(var_coord_vol) here since the presence of
! this%coord_3d_in (external) has precedence over coord_3d_in internal
! in grid_transform_compute
      CALL compute(grid_trans, vol7d_in, vol7d_out, llev_out, &
       var_coord_vol=var_coord_vol)
    ELSE
      CALL l4f_log(l4f_error, 'v7d_v7d_transform: transformation not valid')
      CALL raise_error()
    ENDIF
  ELSE
    CALL l4f_log(l4f_error, &
     'v7d_v7d_transform: vertint requested but lev_out not provided')
    CALL raise_error()
  ENDIF

ELSE

  CALL init(grid_trans, this, vol7d_in, v7d_locana, maskbounds=maskbounds, &
   categoryappend=categoryappend)
! if this init is successful, I am sure that v7d_locana%ana is associated

  IF (c_e(grid_trans)) then! .AND. nvar > 0) THEN

    CALL vol7d_alloc(vol7d_out, nana=SIZE(v7d_locana%ana), &
     ntime=SIZE(vol7d_in%time), ntimerange=SIZE(vol7d_in%timerange), &
     nlevel=SIZE(vol7d_in%level), nnetwork=SIZE(vol7d_in%network), ndativarr=nvar)
    vol7d_out%ana = v7d_locana%ana

    CALL get_val(grid_trans, point_mask=point_mask, output_point_index=point_index)

    IF (ALLOCATED(point_index)) THEN
      CALL vol7d_alloc(vol7d_out, nanavari=1)
      CALL init(vol7d_out%anavar%i(1), 'B01192')
    ENDIF

    CALL vol7d_alloc_vol(vol7d_out)

    IF (ALLOCATED(point_index)) THEN
      DO inetwork = 1, SIZE(vol7d_in%network)
        vol7d_out%volanai(:,1,inetwork) = point_index(:)
      ENDDO
    ENDIF
    CALL compute(grid_trans, vol7d_in, vol7d_out)

    IF (ALLOCATED(point_mask)) THEN ! keep full ana
      IF (SIZE(point_mask) /= SIZE(vol7d_in%ana)) THEN
        CALL l4f_log(l4f_warn, &
         'v7d_v7d_transform: inconsistency in point size: '//t2c(SIZE(point_mask)) &
         //':'//t2c(SIZE(vol7d_in%ana)))
      ELSE
#ifdef DEBUG
        CALL l4f_log(l4f_debug, 'v7d_v7d_transform: merging ana from in to out')
#endif
        CALL vol7d_copy(vol7d_in, vol7d_tmpana, &
         lana=point_mask, lnetwork=(/.true./), &
         ltime=(/.false./), ltimerange=(/.false./), llevel=(/.false./))
        CALL vol7d_append(vol7d_out, vol7d_tmpana)
      ENDIF
    ENDIF

  ELSE
    CALL l4f_log(l4f_error, 'v7d_v7d_transform: transformation not valid')
    CALL raise_error()
  ENDIF

ENDIF

CALL delete (grid_trans)
IF (.NOT. PRESENT(v7d)) CALL delete(v7d_locana)

END SUBROUTINE v7d_v7d_transform


subroutine vg6d_wind_unrot(this)
type(volgrid6d) :: this

integer :: component_flag

call get_val(this%griddim,component_flag=component_flag)

if (component_flag == 1) then
  call l4f_category_log(this%category,l4f_info, &
   "unrotating vector components")
  call vg6d_wind__un_rot(this,.false.)
  call set_val(this%griddim,component_flag=0)
else
  call l4f_category_log(this%category,l4f_info, &
   "no need to unrotate vector components")
end if

end subroutine vg6d_wind_unrot


subroutine vg6d_wind_rot(this)
type(volgrid6d) :: this

integer :: component_flag

call get_val(this%griddim,component_flag=component_flag)

if (component_flag == 0) then
  call l4f_category_log(this%category,l4f_info, &
   "rotating vector components")
  call vg6d_wind__un_rot(this,.true.)
  call set_val(this%griddim,component_flag=1)
else
  call l4f_category_log(this%category,l4f_info, &
   "no need to rotate vector components")
end if

end subroutine vg6d_wind_rot


! Generic UnRotate the wind components.
SUBROUTINE vg6d_wind__un_rot(this,rot)
TYPE(volgrid6d) :: this ! object containing wind to be (un)rotated
LOGICAL :: rot ! if .true. rotate else unrotate

INTEGER :: i, j, k, l, a11, a12, a21, a22, stallo
double precision,pointer :: rot_mat(:,:,:)
real,allocatable :: tmp_arr(:,:)
REAL,POINTER :: voldatiu(:,:), voldativ(:,:)
INTEGER,POINTER :: iu(:), iv(:)

IF (.NOT.ASSOCIATED(this%var)) THEN
  CALL l4f_category_log(this%category, l4f_error, &
   "trying to unrotate an incomplete volgrid6d object")
  CALL raise_fatal_error()
!  RETURN
ENDIF

CALL volgrid6d_var_hor_comp_index(this%var, iu, iv)
IF (.NOT.ASSOCIATED(iu)) THEN
  CALL l4f_category_log(this%category,l4f_error, &
   "unrotation impossible")
  CALL raise_fatal_error()
!  RETURN
ENDIF

! Temporary workspace
ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo)
IF (stallo /= 0) THEN
  CALL l4f_category_log(this%category, l4f_fatal, "allocating memory")
  CALL raise_fatal_error()
ENDIF
! allocate once for speed
IF (.NOT.ASSOCIATED(this%voldati)) THEN
  ALLOCATE(voldatiu(this%griddim%dim%nx, this%griddim%dim%ny), &
   voldativ(this%griddim%dim%nx, this%griddim%dim%ny))
ENDIF

CALL griddim_unproj(this%griddim)
CALL wind_unrot(this%griddim, rot_mat)

a11=1
if (rot)then
  a12=2
  a21=3
else
  a12=3
  a21=2
end if
a22=4

DO l = 1, SIZE(iu)
  DO k = 1, SIZE(this%timerange)
    DO j = 1, SIZE(this%time)
      DO i = 1, SIZE(this%level)
! get data
        CALL volgrid_get_vol_2d(this, i, j, k, iu(l), voldatiu)
        CALL volgrid_get_vol_2d(this, i, j, k, iv(l), voldativ)
! convert units forward
!        CALL compute(conv_fwd(iu(l)), voldatiu)
!        CALL compute(conv_fwd(iv(l)), voldativ)

! multiply wind components by rotation matrix
        WHERE(voldatiu /= rmiss .AND. voldativ /= rmiss)
          tmp_arr(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a11) + &
           voldativ(:,:)*rot_mat(:,:,a12))
          voldativ(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a21) + &
           voldativ(:,:)*rot_mat(:,:,a22))
          voldatiu(:,:) = tmp_arr(:,:)
        END WHERE
! convert units backward
!        CALL uncompute(conv_fwd(iu(l)), voldatiu)
!        CALL uncompute(conv_fwd(iv(l)), voldativ)
! put data
        CALL volgrid_set_vol_2d(this, i, j, k, iu(l), voldatiu)
        CALL volgrid_set_vol_2d(this, i, j, k, iv(l), voldativ)
      ENDDO
    ENDDO
  ENDDO
ENDDO

IF (.NOT.ASSOCIATED(this%voldati)) THEN
  DEALLOCATE(voldatiu, voldativ)
ENDIF
DEALLOCATE(rot_mat, tmp_arr, iu, iv)

END SUBROUTINE vg6d_wind__un_rot


!!$ try to understand the problem:
!!$
!!$ case:
!!$
!!$ 1) we have only one volume: we have to provide the direction of shift
!!$                           compute H and traslate on it
!!$ 2) we have two volumes:
!!$      1) volume U and volume V: compute H and traslate on it
!!$      2) volume U/V and volume H : translate U/V on H
!!$ 3) we have tree volumes: translate U and V on H
!!$
!!$ strange cases:
!!$ 1) do not have U in volume U
!!$ 2) do not have V in volume V
!!$ 3) we have others variables more than U and V in volumes U e V
!!$
!!$
!!$ so the steps are:
!!$ 1) find the volumes 
!!$ 2) define or compute H grid
!!$ 3) trasform the volumes in H 

!!$ N.B.
!!$ case 1) for only one vg6d (U or V) is not managed, but
!!$ the not pubblic subroutines will work but you have to know what you want to do


subroutine vg6d_c2a (this)

TYPE(volgrid6d),INTENT(inout)  :: this(:)

integer :: ngrid,igrid,jgrid,ugrid,vgrid,tgrid
doubleprecision :: xmin, xmax, ymin, ymax
doubleprecision :: xmin_t, xmax_t, ymin_t, ymax_t
doubleprecision :: step_lon_t,step_lat_t
character(len=80) :: type_t,type
TYPE(griddim_def):: griddim_t

ngrid=size(this)

do igrid=1,ngrid

  call init(griddim_t)

  call get_val(this(igrid)%griddim,xmin=xmin_t, xmax=xmax_t, ymin=ymin_t, ymax=ymax_t,proj_type=type_t)
  step_lon_t=(xmax_t-xmin_t)/dble(this(igrid)%griddim%dim%nx-1)
  step_lat_t=(ymax_t-ymin_t)/dble(this(igrid)%griddim%dim%ny-1)

  do jgrid=1,ngrid

    ugrid=imiss
    vgrid=imiss
    tgrid=imiss

#ifdef DEBUG
    call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: search U/V/T points:"//to_char(igrid)//to_char(jgrid))
    call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test delta: "//to_char(step_lon_t)//to_char(step_lat_t))
#endif
    
    if (this(igrid)%griddim == this(jgrid)%griddim ) cycle

    if (this(igrid)%griddim%dim%nx == this(jgrid)%griddim%dim%nx .and. &
     this(igrid)%griddim%dim%ny == this(jgrid)%griddim%dim%ny ) then

      call get_val(this(jgrid)%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,proj_type=type)
      
      if (type_t /= type )cycle

#ifdef DEBUG
      call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test U "//&
       to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))

      call l4f_category_log(this(igrid)%category,l4f_debug,"diff coordinate lon"//&
       to_char(abs(xmin - (xmin_t+step_lon_t/2.d0)))//&
       to_char(abs(xmax - (xmax_t+step_lon_t/2.d0))))
      call l4f_category_log(this(igrid)%category,l4f_debug,"diff coordinate lat"//&
       to_char(abs(ymin - (ymin_t+step_lat_t/2.d0)))//&
       to_char(abs(ymax - (ymax_t+step_lat_t/2.d0))))
#endif

      if ( abs(xmin - (xmin_t+step_lon_t/2.d0)) < 1.d-3 .and. abs(xmax - (xmax_t+step_lon_t/2.d0)) < 1.d-3 ) then
        if ( abs(ymin - ymin_t) < 1.d-3 .and. abs(ymax - ymax_t) < 1.d-3 ) then

#ifdef DEBUG
          call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found U")
#endif
          ugrid=jgrid
          tgrid=igrid

        end if
      end if

#ifdef DEBUG
      call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test V "//&
       to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))
#endif

      if ( abs(ymin - (ymin_t+step_lat_t/2.d0)) < 1.d-3 .and. abs(ymax - (ymax_t+step_lat_t/2.d0)) < 1.d-3 ) then
        if ( abs(xmin - xmin_t) < 1.d-3 .and. abs(xmax - xmax_t) < 1.d-3  ) then
          
#ifdef DEBUG
          call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found V")
#endif
          vgrid=jgrid
          tgrid=igrid

        end if
      end if


      ! test if we have U and V

#ifdef DEBUG
      call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test U and V"//&
       to_char(xmin_t)//to_char(xmax_t)//to_char(ymin_t)//to_char(ymax_t)//&
       to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))

      call l4f_category_log(this(igrid)%category,l4f_debug,"UV diff coordinate lon"//&
       to_char(abs(xmin_t - xmin)-step_lon_t/2.d0)//&
       to_char(abs(xmax_t - xmax)-step_lon_t/2.d0))
      call l4f_category_log(this(igrid)%category,l4f_debug,"UV diff coordinate lat"//&
       to_char(abs(ymin_t - ymin) -step_lat_t/2.d0)//&
       to_char(abs(ymax_t - ymax)-step_lat_t/2.d0))
#endif
      
      if ( abs(ymin - (ymin_t+step_lat_t/2.d0)) < 2.d-3 .and. abs(ymax - (ymax_t+step_lat_t/2.d0)) < 2.d-3 ) then
        if ( abs(xmin_t - (xmin+step_lon_t/2.d0)) < 2.d-3 .and. abs(xmax_t - (xmax+step_lon_t/2.d0)) < 2.d-3 ) then
          
#ifdef DEBUG
          call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found U and V case up and right")
#endif

          vgrid=jgrid
          ugrid=igrid

          call init(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin_t, ymax=ymax_t)

        end if
      end if
    end if

                                ! abbiamo almeno due volumi: riportiamo U e/o V su T
    if (c_e(ugrid)) then
#ifdef DEBUG
      call l4f_category_log(this(igrid)%category,l4f_debug,"C grid U points: interpolate U on T "//to_char(tgrid)//to_char(ugrid))
#endif
      if (c_e(tgrid))then
        call vg6d_c2a_grid(this(ugrid),this(tgrid)%griddim,cgrid=1)
      else
        call vg6d_c2a_grid(this(ugrid),griddim_t,cgrid=1)
      end if
      call vg6d_c2a_mat(this(ugrid),cgrid=1)
    end if
    
    if (c_e(vgrid)) then
#ifdef DEBUG
      call l4f_category_log(this(igrid)%category,l4f_debug,"C grid V points: interpolate V on T "//to_char(tgrid)//to_char(vgrid))
#endif
      if (c_e(tgrid))then
        call vg6d_c2a_grid(this(vgrid),this(tgrid)%griddim,cgrid=2)
      else
        call vg6d_c2a_grid(this(vgrid),griddim_t,cgrid=2)
      end if
      call vg6d_c2a_mat(this(vgrid),cgrid=2)
    end if

  end do

  call delete(griddim_t)

end do
  

end subroutine vg6d_c2a


! Convert C grid to A grid
subroutine vg6d_c2a_grid(this,griddim_t,cgrid)

type(volgrid6d),intent(inout) :: this ! object containing fields to be translated (U or V points)
type(griddim_def),intent(in),optional :: griddim_t ! object containing grid of T points
integer,intent(in) :: cgrid ! in C grid (Arakawa) we have 0=T,1=U,2=V points

doubleprecision :: xmin, xmax, ymin, ymax
doubleprecision :: step_lon,step_lat


if (present(griddim_t)) then

 call get_val(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
 call set_val(this%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
! improve grid definition precision
 CALL griddim_setsteps(this%griddim)

else

  select case (cgrid)

  case(0)

#ifdef DEBUG
    call l4f_category_log(this%category,l4f_debug,"C grid: T points, nothing to do")
#endif
    return

  case (1)

#ifdef DEBUG
    call l4f_category_log(this%category,l4f_debug,"C grid: U points, we need interpolation")
#endif

    call get_val(this%griddim, xmin=xmin, xmax=xmax)
    step_lon=(xmax-xmin)/dble(this%griddim%dim%nx-1)
    xmin=xmin-step_lon/2.d0
    xmax=xmax-step_lon/2.d0
    call set_val(this%griddim, xmin=xmin, xmax=xmax)
! improve grid definition precision
    CALL griddim_setsteps(this%griddim)
    
  case (2)

#ifdef DEBUG
    call l4f_category_log(this%category,l4f_debug,"C grid: V points, we need interpolation")
#endif

    call get_val(this%griddim, ymin=ymin, ymax=ymax)
    step_lat=(ymax-ymin)/dble(this%griddim%dim%ny-1)
    ymin=ymin-step_lat/2.d0
    ymax=ymax-step_lat/2.d0
    call set_val(this%griddim, ymin=ymin, ymax=ymax)
! improve grid definition precision
    CALL griddim_setsteps(this%griddim)
    
  case default

    call l4f_category_log(this%category,l4f_fatal,"C grid type not known")
    call raise_fatal_error ()

  end select

end if


call griddim_unproj(this%griddim)


end subroutine vg6d_c2a_grid

! Convert C grid to A grid
subroutine vg6d_c2a_mat(this,cgrid)

type(volgrid6d),intent(inout) :: this ! object containing fields to be translated
integer,intent(in) :: cgrid ! in C grid (Arakawa) we have 0=T,1=U,2=V points

INTEGER :: i, j, k, iv, stallo
REAL,ALLOCATABLE :: tmp_arr(:,:)
REAL,POINTER :: voldatiuv(:,:)


IF (cgrid == 0) RETURN ! nothing to do
IF (cgrid /= 1 .AND. cgrid /= 2) THEN ! wrong cgrid
  CALL l4f_category_log(this%category,l4f_fatal,"C grid type "// &
   trim(to_char(cgrid))//" not known")
  CALL raise_error()
  RETURN
ENDIF

! Temporary workspace
ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo)
if (stallo /=0)then
  call l4f_log(l4f_fatal,"allocating memory")
  call raise_fatal_error()
end if

! allocate once for speed
IF (.NOT.ASSOCIATED(this%voldati)) THEN
  ALLOCATE(voldatiuv(this%griddim%dim%nx, this%griddim%dim%ny), stat=stallo)
  IF (stallo /= 0) THEN
    CALL l4f_log(l4f_fatal,"allocating memory")
    CALL raise_fatal_error()
  ENDIF
ENDIF

IF (cgrid == 1) THEN ! U points to H points
  DO iv = 1, SIZE(this%var)
    DO k = 1, SIZE(this%timerange)
      DO j = 1, SIZE(this%time)
        DO i = 1, SIZE(this%level)
          tmp_arr(:,:) = rmiss
          CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv)

! West boundary extrapolation
          WHERE(voldatiuv(1,:) /= rmiss .AND. voldatiuv(2,:) /= rmiss)
            tmp_arr(1,:) = voldatiuv(1,:) - (voldatiuv(2,:) - voldatiuv(1,:)) / 2.
          END WHERE

! Rest of the field
          WHERE(voldatiuv(1:this%griddim%dim%nx-1,:) /= rmiss .AND. &
           voldatiuv(2:this%griddim%dim%nx,:) /= rmiss)
            tmp_arr(2:this%griddim%dim%nx,:) = &
             (voldatiuv(1:this%griddim%dim%nx-1,:) + &
             voldatiuv(2:this%griddim%dim%nx,:)) / 2.
          END WHERE

          voldatiuv(:,:) = tmp_arr
          CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv)
        ENDDO
      ENDDO
    ENDDO
  ENDDO

ELSE IF (cgrid == 2) THEN ! V points to H points
  DO iv = 1, SIZE(this%var)
    DO k = 1, SIZE(this%timerange)
      DO j = 1, SIZE(this%time)
        DO i = 1, SIZE(this%level)
          tmp_arr(:,:) = rmiss
          CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv)

! South boundary extrapolation
          WHERE(voldatiuv(:,1) /= rmiss .AND. voldatiuv(:,2) /= rmiss)
            tmp_arr(:,1) = voldatiuv(:,1) - (voldatiuv(:,2) - voldatiuv(:,1)) / 2.
          END WHERE

! Rest of the field
          WHERE(voldatiuv(:,1:this%griddim%dim%ny-1) /= rmiss .AND. &
           voldatiuv(:,2:this%griddim%dim%ny) /= rmiss)
            tmp_arr(:,2:this%griddim%dim%ny) = &
             (voldatiuv(:,1:this%griddim%dim%ny-1) + &
             voldatiuv(:,2:this%griddim%dim%ny)) / 2.
          END WHERE

          voldatiuv(:,:) = tmp_arr
          CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv)
        ENDDO
      ENDDO
    ENDDO
  ENDDO
ENDIF ! tertium non datur

IF (.NOT.ASSOCIATED(this%voldati)) THEN
  DEALLOCATE(voldatiuv)
ENDIF
DEALLOCATE (tmp_arr)

end subroutine vg6d_c2a_mat


subroutine display_volgrid6d (this)

TYPE(volgrid6d),intent(in) :: this
integer :: i

#ifdef DEBUG
call l4f_category_log(this%category,l4f_debug,"ora mostro gridinfo " )
#endif

print*,"----------------------- volgrid6d display ---------------------"
call display(this%griddim)

IF (ASSOCIATED(this%time))then
  print*,"---- time vector ----"
  print*,"elements=",size(this%time)
  do i=1, size(this%time)
    call display(this%time(i))
  end do
end IF

IF (ASSOCIATED(this%timerange))then
  print*,"---- timerange vector ----"
  print*,"elements=",size(this%timerange)
  do i=1, size(this%timerange)
    call display(this%timerange(i))
  end do
end IF

IF (ASSOCIATED(this%level))then
  print*,"---- level vector ----"
  print*,"elements=",size(this%level)
  do i=1, size(this%level)
    call display(this%level(i))
  end do
end IF

IF (ASSOCIATED(this%var))then
  print*,"---- var vector ----"
  print*,"elements=",size(this%var)
  do i=1, size(this%var)
    call display(this%var(i))
  end do
end IF

IF (ASSOCIATED(this%gaid))then
  print*,"---- gaid vector (present mask only) ----"
  print*,"elements=",shape(this%gaid)
  print*,c_e(reshape(this%gaid,(/SIZE(this%gaid)/)))
end IF

print*,"--------------------------------------------------------------"


end subroutine display_volgrid6d


subroutine display_volgrid6dv (this)

TYPE(volgrid6d),intent(in) :: this(:)
integer :: i

print*,"----------------------- volgrid6d  vector ---------------------"

print*,"elements=",size(this)

do i=1, size(this)

#ifdef DEBUG
  call l4f_category_log(this(i)%category,l4f_debug,"ora mostro il vettore volgrid6d" )
#endif

  call display(this(i))

end do
print*,"--------------------------------------------------------------"

end subroutine display_volgrid6dv


subroutine vg6dv_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge)
type(volgrid6d),intent(in) :: vg6din(:)
type(volgrid6d),intent(out),pointer :: vg6dout(:)
type(vol7d_level),intent(in),optional :: level(:)
type(vol7d_timerange),intent(in),optional :: timerange(:)
logical,intent(in),optional :: merge
logical,intent(in),optional :: nostatproc

integer :: i

allocate(vg6dout(size(vg6din)))

do i = 1, size(vg6din)
  call vg6d_rounding(vg6din(i),vg6dout(i),level,timerange,nostatproc,merge)
end do

end subroutine vg6dv_rounding

subroutine vg6d_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge)
type(volgrid6d),intent(in) :: vg6din
type(volgrid6d),intent(out) :: vg6dout
type(vol7d_level),intent(in),optional :: level(:)
type(vol7d_timerange),intent(in),optional :: timerange(:)
logical,intent(in),optional :: merge
logical,intent(in),optional :: nostatproc

integer :: ilevel,itimerange
type(vol7d_level) :: roundlevel(size(vg6din%level)) 
type(vol7d_timerange) :: roundtimerange(size(vg6din%timerange))

roundlevel=vg6din%level

if (present(level))then
  do ilevel = 1, size(vg6din%level)
    if ((any(vg6din%level(ilevel) .almosteq. level))) then
      roundlevel(ilevel)=level(1)
    end if
  end do
end if

roundtimerange=vg6din%timerange

if (present(timerange))then
  do itimerange = 1, size(vg6din%timerange)
    if ((any(vg6din%timerange(itimerange) .almosteq. timerange))) then
      roundtimerange(itimerange)=timerange(1)
    end if
  end do
end if

!set istantaneous values everywere
!preserve p1 for forecast time
if (optio_log(nostatproc)) then
  roundtimerange(:)%timerange=254
  roundtimerange(:)%p2=0
end if


call vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge)

end subroutine vg6d_rounding

subroutine vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge)
type(volgrid6d),intent(in) :: vg6din
type(volgrid6d),intent(out) :: vg6dout
type(vol7d_level),intent(in) :: roundlevel(:)
type(vol7d_timerange),intent(in) :: roundtimerange(:)
logical,intent(in),optional :: merge

integer :: nlevel,ntime,ntimerange,nvar,ilevel,itimerange,ivar,indl,indt,itime,nx,ny
real,allocatable :: vol2d(:,:)

nx=vg6din%griddim%dim%nx
ny=vg6din%griddim%dim%ny
nlevel=count_distinct(roundlevel,back=.true.)
ntime=size(vg6din%time)
ntimerange=count_distinct(roundtimerange,back=.true.)
nvar=size(vg6din%var)

call init(vg6dout, vg6din%griddim, vg6din%time_definition, categoryappend="generated by vg6d_reduce")
call volgrid6d_alloc(vg6dout, vg6din%griddim%dim, ntime, nlevel, ntimerange, nvar)

if ( ASSOCIATED(vg6din%voldati) .or. optio_log(merge)) then          
  call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.true.)
  allocate(vol2d(nx,ny))
else
  call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.false.)
end if

vg6dout%time=vg6din%time
vg6dout%var=vg6din%var
vg6dout%timerange=pack_distinct(roundtimerange,ntimerange,back=.true.)
vg6dout%level=pack_distinct(roundlevel,nlevel,back=.true.)
! sort modified dimensions
CALL sort(vg6dout%timerange)
CALL sort(vg6dout%level)

do ilevel=1,size(vg6din%level)
  indl=index(vg6dout%level,roundlevel(ilevel))
  do itimerange=1,size(vg6din%timerange)
    indt=index(vg6dout%timerange,roundtimerange(itimerange))
    do ivar=1, nvar
      do itime=1,ntime

        if ( ASSOCIATED(vg6din%voldati)) then
          vol2d=vg6din%voldati(:,:,ilevel,itime,itimerange,ivar)
        end if

        if (optio_log(merge)) then          

          if ( .not. ASSOCIATED(vg6din%voldati)) then
            CALL grid_id_decode_data(vg6din%gaid(ilevel,itime,itimerange,ivar), vol2d)
          end if

          !! merge present data point by point
          where (.not. c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar))) 
          
            vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d
          
          end where
        else if ( ASSOCIATED(vg6din%voldati)) then
          if (.not. any(c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar))))then
            vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d
          end if
        end if
      
        if (c_e(vg6din%gaid(ilevel,itime,itimerange,ivar)).and. .not. c_e(vg6dout%gaid(indl,itime,indt,ivar)))then
          call copy (vg6din%gaid(ilevel,itime,itimerange,ivar), vg6dout%gaid(indl,itime,indt,ivar))
        end if
      end do
    end do
  end do
end do

if ( ASSOCIATED(vg6din%voldati) .or. optio_log(merge)) then          
  deallocate(vol2d)
end if

end subroutine vg6d_reduce


end module volgrid6d_class