grid_transform_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 grid_transform_class
USE vol7d_class
USE err_handling
USE geo_proj_class
USE grid_class
USE grid_dim_class
USE optional_values
USE array_utilities
USE georef_coord_class
USE simple_stat
USE log4fortran
IMPLICIT NONE

CHARACTER(len=255),PARAMETER:: subcategory="grid_transform_class"

! information for interpolation aver a rectangular area
TYPE area_info
  double precision :: boxdx ! longitudinal/x extension of the box for box interpolation, default the target x grid step
  double precision :: boxdy ! latitudinal/y extension of the box for box interpolation, default the target y grid step
  double precision :: radius ! radius in gridpoints for stencil interpolation
END TYPE area_info

! information for statistical processing of interpoland data
TYPE stat_info
  DOUBLE PRECISION :: percentile ! percentile [0,100] of the distribution of points in the box to use as interpolated value, if missing, the average is used, if 0.or 100. the MIN() and MAX() are used as a shortcut
END TYPE stat_info

! information for point interval
TYPE interval_info
  REAL :: gt=rmiss ! lower limit of interval, missing for -inf
  REAL :: lt=rmiss ! upper limit of interval, missing for +inf
  LOGICAL :: ge=.true. ! if true >= otherwise >
  LOGICAL :: le=.true. ! if true <= otherwise <
END TYPE interval_info

! rectangle index information
type rect_ind
  INTEGER :: ix ! index of initial point of new grid on x
  INTEGER :: iy ! index of initial point of new grid on y
  INTEGER :: fx ! index of final point of new grid on x
  INTEGER :: fy ! index of final point of new grid on y
end type rect_ind

! rectangle coord information
type rect_coo
  DOUBLEPRECISION ilon ! coordinate of initial point of new grid on x
  DOUBLEPRECISION ilat ! coordinate of initial point of new grid on y
  DOUBLEPRECISION flon ! coordinate of final point of new grid on x
  DOUBLEPRECISION flat ! coordinate of final point of new grid on y
end type rect_coo

! box information
type box_info
  INTEGER :: npx ! number of points along x direction
  INTEGER :: npy ! number of points along y direction
end type box_info

! Vertical interpolation information.
! The input vertical coordinate can be indicated either as the value
! of the vertical level (so that it will be the same on every point
! at a given vertical level), or as the value of a specified variable
! at each point in space (so that it will depend on the horizontal
! position too).
TYPE vertint
!  CHARACTER(len=80) :: sub_type ! subtype of transformation, can be \c 'linear'
  TYPE(vol7d_level) :: input_levtype ! type of vertical level of input data (only type of first and second surface are used, level values are ignored)
  TYPE(vol7d_var) :: input_coordvar ! variable that defines the vertical coordinate in the input volume, if missing, the value of the vertical level is used
  TYPE(vol7d_level) :: output_levtype ! type of vertical level of output data (only type of first and second surface are used, level values are ignored)
END TYPE vertint

TYPE transform_def
  PRIVATE
  CHARACTER(len=80) :: trans_type ! type of transformation, can be \c 'zoom', \c 'boxregrid', \c 'inter', \c 'vertint' ...
  CHARACTER(len=80) :: sub_type ! subtype of transformation, can be \c 'linear'
  LOGICAL :: extrap ! enable elaboration outside data bounding box
  TYPE(rect_ind) :: rect_ind ! rectangle information by index
  TYPE(rect_coo) :: rect_coo ! rectangle information by coordinates
  TYPE(area_info) :: area_info ! 
  TYPE(arrayof_georef_coord_array) :: poly ! polygon information
  TYPE(stat_info) :: stat_info ! 
  TYPE(interval_info) :: interval_info ! 
  TYPE(box_info) :: box_info ! boxregrid specification
  TYPE(vertint) :: vertint ! vertical interpolation specification
  INTEGER :: time_definition ! time definition for interpolating to sparse points
  TYPE(l4f_handle) :: category ! category for log4fortran
END TYPE transform_def


TYPE grid_transform
  PRIVATE
  TYPE(transform_def),PUBLIC :: trans ! type of transformation required
  INTEGER :: innx = imiss
  INTEGER :: inny = imiss
  INTEGER :: innz = imiss
  INTEGER :: outnx = imiss
  INTEGER :: outny = imiss
  INTEGER :: outnz = imiss
  INTEGER :: levshift = imiss
  INTEGER :: levused = imiss
  INTEGER :: iniox,inioy,infox,infoy,outinx,outiny,outfnx,outfny
  INTEGER,POINTER :: inter_index_x(:,:) => null()
  INTEGER,POINTER :: inter_index_y(:,:) => null()
  INTEGER,POINTER :: inter_index_z(:) => null()
  INTEGER,POINTER :: point_index(:,:) => null()
  DOUBLE PRECISION,POINTER :: inter_x(:,:) => null()
  DOUBLE PRECISION,POINTER :: inter_y(:,:) => null()
  DOUBLE PRECISION,POINTER :: inter_xp(:,:) => null()
  DOUBLE PRECISION,POINTER :: inter_yp(:,:) => null()
  DOUBLE PRECISION,POINTER :: inter_zp(:) => null()
  DOUBLE PRECISION,POINTER :: vcoord_in(:) => null()
  DOUBLE PRECISION,POINTER :: vcoord_out(:) => null()
  LOGICAL,POINTER :: point_mask(:,:) => null()
  LOGICAL,POINTER :: stencil(:,:) => null()
  REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
  REAL,ALLOCATABLE :: val_mask(:,:)
  REAL :: val1 = rmiss
  REAL :: val2 = rmiss
  LOGICAL :: recur = .false.
  LOGICAL :: dolog = .false. ! must compute log() of vert coord before vertint

!  type(volgrid6d) :: input_vertcoordvol ! volume which provides the input vertical coordinate if separated from the data volume itself (for vertint) cannot be here because of cross-use, should be an argument of compute
!  type(vol7d_level), pointer :: output_vertlevlist(:) ! list of vertical levels of output data (for vertint) can be here or an argument of compute, how to do?
  TYPE(vol7d_level),POINTER :: output_level_auto(:) => null() ! array of auto-generated levels, stored for successive query
  TYPE(l4f_handle) :: category ! category for log4fortran
  LOGICAL :: valid = .false. ! the transformation has been successfully initialised
  PROCEDURE(basic_find_index),NOPASS,POINTER :: find_index => basic_find_index ! allow a local implementation of find_index
END TYPE grid_transform


INTERFACE init
  MODULE PROCEDURE transform_init, grid_transform_levtype_levtype_init, &
   grid_transform_init, &
   grid_transform_grid_vol7d_init, grid_transform_vol7d_grid_init, &
   grid_transform_vol7d_vol7d_init
END INTERFACE

INTERFACE delete
  MODULE PROCEDURE transform_delete, grid_transform_delete
END INTERFACE

INTERFACE get_val
  MODULE PROCEDURE transform_get_val, grid_transform_get_val
END INTERFACE

INTERFACE compute
  MODULE PROCEDURE grid_transform_compute, grid_transform_v7d_grid_compute
END INTERFACE

INTERFACE c_e
  MODULE PROCEDURE grid_transform_c_e
END INTERFACE

PRIVATE
PUBLIC init, delete, get_val, compute, c_e
PUBLIC transform_def, grid_transform
PUBLIC interval_info, interval_info_new, interval_info_valid, basic_find_index

CONTAINS


FUNCTION interval_info_new(interv_gt, interv_ge, interv_lt, interv_le) RESULT(this)
REAL,INTENT(in),OPTIONAL :: interv_gt
REAL,INTENT(in),OPTIONAL :: interv_ge
REAL,INTENT(in),OPTIONAL :: interv_lt
REAL,INTENT(in),OPTIONAL :: interv_le

TYPE(interval_info) :: this

IF (PRESENT(interv_gt)) THEN
  IF (c_e(interv_gt)) THEN
    this%gt = interv_gt
    this%ge = .false.
  ENDIF
ENDIF
IF (PRESENT(interv_ge)) THEN
  IF (c_e(interv_ge)) THEN
    this%gt = interv_ge
    this%ge = .true.
  ENDIF
ENDIF
IF (PRESENT(interv_lt)) THEN
  IF (c_e(interv_lt)) THEN
    this%lt = interv_lt
    this%le = .false.
  ENDIF
ENDIF
IF (PRESENT(interv_le)) THEN
  IF (c_e(interv_le)) THEN
    this%lt = interv_le
    this%le = .true.
  ENDIF
ENDIF

END FUNCTION interval_info_new

! Private method for testing whether \a val is included in \a this
! interval taking into account all cases, zero, one or two extremes,
! strict or non strict inclusion, empty interval, etc, while no check
! is made for val being missing. Returns 1.0 if val is in interval and
! 0.0 if not.
ELEMENTAL FUNCTION interval_info_valid(this, val)
TYPE(interval_info),INTENT(in) :: this
REAL,INTENT(in) :: val

REAL :: interval_info_valid

interval_info_valid = 1.0

IF (c_e(this%gt)) THEN
  IF (val < this%gt) interval_info_valid = 0.0
  IF (.NOT.this%ge) THEN
    IF (val == this%gt) interval_info_valid = 0.0
  ENDIF
ENDIF
IF (c_e(this%lt)) THEN
  IF (val > this%lt) interval_info_valid = 0.0
  IF (.NOT.this%le) THEN
    IF (val == this%lt) interval_info_valid = 0.0
  ENDIF
ENDIF

END FUNCTION interval_info_valid

SUBROUTINE transform_init(this, trans_type, sub_type, &
 ix, iy, fx, fy, ilon, ilat, flon, flat, &
 npx, npy, boxdx, boxdy, radius, poly, percentile, &
 interv_gt, interv_ge, interv_lt, interv_le, &
 extrap, time_definition, &
 input_levtype, input_coordvar, output_levtype, categoryappend)
TYPE(transform_def),INTENT(out) :: this
CHARACTER(len=*) :: trans_type
CHARACTER(len=*) :: sub_type
INTEGER,INTENT(in),OPTIONAL :: ix
INTEGER,INTENT(in),OPTIONAL :: iy
INTEGER,INTENT(in),OPTIONAL :: fx
INTEGER,INTENT(in),OPTIONAL :: fy
DOUBLEPRECISION,INTENT(in),OPTIONAL :: ilon
DOUBLEPRECISION,INTENT(in),OPTIONAL :: ilat
DOUBLEPRECISION,INTENT(in),OPTIONAL :: flon
DOUBLEPRECISION,INTENT(in),OPTIONAL :: flat
INTEGER,INTENT(IN),OPTIONAL :: npx
INTEGER,INTENT(IN),OPTIONAL :: npy
DOUBLEPRECISION,INTENT(in),OPTIONAL :: boxdx
DOUBLEPRECISION,INTENT(in),OPTIONAL :: boxdy
DOUBLEPRECISION,INTENT(in),OPTIONAL :: radius
TYPE(arrayof_georef_coord_array),OPTIONAL :: poly
DOUBLEPRECISION,INTENT(in),OPTIONAL :: percentile
REAL,INTENT(in),OPTIONAL :: interv_gt
REAL,INTENT(in),OPTIONAL :: interv_ge
REAL,INTENT(in),OPTIONAL :: interv_lt
REAL,INTENT(in),OPTIONAL :: interv_le
LOGICAL,INTENT(IN),OPTIONAL :: extrap
INTEGER,INTENT(IN),OPTIONAL :: time_definition
TYPE(vol7d_level),INTENT(IN),OPTIONAL :: input_levtype
TYPE(vol7d_var),INTENT(IN),OPTIONAL :: input_coordvar
TYPE(vol7d_level),INTENT(IN),OPTIONAL :: output_levtype
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)

this%trans_type = trans_type
this%sub_type = sub_type

CALL optio(extrap,this%extrap)

call optio(ix,this%rect_ind%ix)
call optio(iy,this%rect_ind%iy)
call optio(fx,this%rect_ind%fx)
call optio(fy,this%rect_ind%fy)

call optio(ilon,this%rect_coo%ilon)
call optio(ilat,this%rect_coo%ilat)
call optio(flon,this%rect_coo%flon)
call optio(flat,this%rect_coo%flat)

CALL optio(boxdx,this%area_info%boxdx)
CALL optio(boxdy,this%area_info%boxdy)
CALL optio(radius,this%area_info%radius)
IF (PRESENT(poly)) this%poly = poly
CALL optio(percentile,this%stat_info%percentile)

this%interval_info = interval_info_new(interv_gt, interv_ge, interv_lt, interv_le)

CALL optio(npx,this%box_info%npx)
CALL optio(npy,this%box_info%npy)

IF (PRESENT(input_levtype)) THEN
  this%vertint%input_levtype = input_levtype
ELSE
  this%vertint%input_levtype = vol7d_level_miss
ENDIF
IF (PRESENT(input_coordvar)) THEN
  this%vertint%input_coordvar = input_coordvar
ELSE
  this%vertint%input_coordvar = vol7d_var_miss
ENDIF
IF (PRESENT(output_levtype)) THEN
  this%vertint%output_levtype = output_levtype
ELSE
  this%vertint%output_levtype = vol7d_level_miss
ENDIF

call optio(time_definition,this%time_definition)
if (c_e(this%time_definition) .and. &
 (this%time_definition < 0 .OR. this%time_definition > 2))THEN
  call l4f_category_log(this%category,l4f_error,"Error time_definition invalid: "//to_char(this%time_definition))
  call raise_fatal_error()
end if


IF (this%trans_type == 'zoom') THEN

  IF (this%sub_type == 'coord' .OR. this%sub_type == 'projcoord')THEN

    if (c_e(this%rect_coo%ilon) .and. c_e(this%rect_coo%ilat) .and. &
        c_e(this%rect_coo%flon) .and. c_e(this%rect_coo%flat)) then ! coordinates given
    
!check
      if ( this%rect_coo%ilon > this%rect_coo%flon .or. &
       this%rect_coo%ilat > this%rect_coo%flat ) then

        call l4f_category_log(this%category,l4f_error, &
         "invalid zoom coordinates: ")
        call l4f_category_log(this%category,l4f_error, &
         trim(to_char(this%rect_coo%ilon))//'/'// &
         trim(to_char(this%rect_coo%flon)))
        call l4f_category_log(this%category,l4f_error, &
         trim(to_char(this%rect_coo%ilat))//'/'// &
         trim(to_char(this%rect_coo%flat)))
        call raise_fatal_error()
      end if

    else

      call l4f_category_log(this%category,l4f_error,"zoom: coord parameters missing")
      call raise_fatal_error()
        
    end if

  else if (this%sub_type == 'coordbb')then

    if (c_e(this%rect_coo%ilon) .and. c_e(this%rect_coo%ilat) .and. &
        c_e(this%rect_coo%flon) .and. c_e(this%rect_coo%flat)) then ! coordinates given
    else

      call l4f_category_log(this%category,l4f_error,"zoom: coordbb parameters missing")
      call raise_fatal_error()
        
    end if

  else if (this%sub_type == 'index')then

    IF (c_e(this%rect_ind%ix) .AND. c_e(this%rect_ind%iy) .AND. &
     c_e(this%rect_ind%fx) .AND. c_e(this%rect_ind%fy)) THEN

! check
      IF (this%rect_ind%ix > this%rect_ind%fx .OR. &
       this%rect_ind%iy > this%rect_ind%fy) THEN

        CALL l4f_category_log(this%category,l4f_error,'invalid zoom indices: ')
        CALL l4f_category_log(this%category,l4f_error, &
         trim(to_char(this%rect_ind%ix))//'/'// &
         trim(to_char(this%rect_ind%fx)))
        CALL l4f_category_log(this%category,l4f_error, &
         trim(to_char(this%rect_ind%iy))//'/'// &
         trim(to_char(this%rect_ind%fy)))

        CALL raise_fatal_error()
      ENDIF

    ELSE

      CALL l4f_category_log(this%category,l4f_error,&
       'zoom: index parameters ix, iy, fx, fy not provided')
      CALL raise_fatal_error()

    ENDIF

  ELSE
    CALL sub_type_error()
    RETURN
  END IF

ELSE IF (this%trans_type == 'inter' .OR. this%trans_type == 'intersearch') THEN

  IF (this%sub_type == 'near' .OR. this%sub_type == 'bilin' .OR. &
   this%sub_type == 'linear' .OR. this%sub_type == 'shapiro_near') THEN
! nothing to do here
  ELSE
    CALL sub_type_error()
    RETURN
  ENDIF

ELSE IF (this%trans_type == 'boxregrid' .OR. this%trans_type == 'boxinter' .OR. &
 this%trans_type == 'polyinter' .OR. this%trans_type == 'maskinter' .OR. &
 this%trans_type == 'stencilinter') THEN

  IF (this%trans_type == 'boxregrid') THEN
    IF (c_e(this%box_info%npx) .AND. c_e(this%box_info%npy)) THEN
      IF (this%box_info%npx <= 0 .OR. this%box_info%npy <= 0 ) THEN
        CALL l4f_category_log(this%category,l4f_error,'boxregrid: invalid parameters  '//&
         trim(to_char(this%box_info%npx))//' '//trim(to_char(this%box_info%npy)))
        CALL raise_fatal_error()
      ENDIF
    ELSE
      CALL l4f_category_log(this%category,l4f_error, &
       'boxregrid: parameters npx, npy missing')
      CALL raise_fatal_error()
    ENDIF
  ENDIF

  IF (this%trans_type == 'polyinter') THEN
    IF (this%poly%arraysize <= 0) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       "polyinter: poly parameter missing or empty")
      CALL raise_fatal_error()
    ENDIF
  ENDIF

  IF (this%trans_type == 'stencilinter') THEN
    IF (.NOT.c_e(this%area_info%radius)) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       "stencilinter: radius parameter missing")
      CALL raise_fatal_error()
    ENDIF
  ENDIF

  IF (this%sub_type == 'average' .OR. this%sub_type == 'stddev' &
   .OR. this%sub_type == 'stddevnm1') THEN
    this%stat_info%percentile = rmiss
  ELSE IF (this%sub_type == 'max') THEN
    this%stat_info%percentile = 101.
  ELSE IF (this%sub_type == 'min') THEN
    this%stat_info%percentile = -1.
  ELSE IF (this%sub_type == 'percentile') THEN
    IF (.NOT.c_e(this%stat_info%percentile)) THEN
      CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
       ':percentile: percentile value not provided')
      CALL raise_fatal_error()
    ELSE IF (this%stat_info%percentile >= 100.) THEN
      this%sub_type = 'max'
    ELSE IF (this%stat_info%percentile <= 0.) THEN
      this%sub_type = 'min'
    ENDIF
  ELSE IF (this%sub_type == 'frequency') THEN
    IF (.NOT.c_e(this%interval_info%gt) .AND. .NOT.c_e(this%interval_info%gt)) THEN
      CALL l4f_category_log(this%category,l4f_error,trim(this%trans_type)// &
       ':frequency: lower and/or upper limit not provided')
      CALL raise_fatal_error()
    ENDIF
  ELSE
    CALL sub_type_error()
    RETURN
  ENDIF

ELSE IF (this%trans_type == 'maskgen')THEN

  IF (this%sub_type == 'poly') THEN

    IF (this%poly%arraysize <= 0) THEN
      CALL l4f_category_log(this%category,l4f_error,"maskgen:poly poly parameter missing or empty")
      CALL raise_fatal_error()
    ENDIF

  ELSE IF (this%sub_type == 'grid') THEN
! nothing to do for now

  ELSE
    CALL sub_type_error()
    RETURN
  ENDIF

ELSE IF (this%trans_type == 'vertint') THEN

  IF (this%vertint%input_levtype == vol7d_level_miss) THEN
    CALL l4f_category_log(this%category,l4f_error, &
     'vertint parameter input_levtype not provided')
    CALL raise_fatal_error()
  ENDIF

  IF (this%vertint%output_levtype == vol7d_level_miss) THEN
    CALL l4f_category_log(this%category,l4f_error, &
     'vertint parameter output_levtype not provided')
    CALL raise_fatal_error()
  ENDIF

  IF (this%sub_type == 'linear' .OR. this%sub_type == 'linearsparse') THEN
! nothing to do here
  ELSE
    CALL sub_type_error()
    RETURN
  ENDIF

ELSE IF (this%trans_type == 'metamorphosis') THEN

  IF (this%sub_type == 'all') THEN
! nothing to do here
  ELSE IF (this%sub_type == 'coordbb')THEN

    IF (c_e(this%rect_coo%ilon) .AND. c_e(this%rect_coo%ilat) .AND. &
     c_e(this%rect_coo%flon) .AND. c_e(this%rect_coo%flat)) THEN ! coordinates given
    ELSE

      CALL l4f_category_log(this%category,l4f_error,"metamorphosis: coordbb parameters missing")
      CALL raise_fatal_error()
        
    ENDIF

  ELSE IF (this%sub_type == 'poly')THEN

    IF (this%poly%arraysize <= 0) THEN
      CALL l4f_category_log(this%category,l4f_error,"metamorphosis:poly: poly parameter missing or empty")
      CALL raise_fatal_error()
    ENDIF

  ELSE IF (this%sub_type == 'mask' .OR. this%sub_type == 'maskvalid' .OR. &
   this%sub_type == 'maskinvalid' .OR. this%sub_type == 'setinvalidto' .OR. &
   this%sub_type == 'settoinvalid' .OR. this%sub_type == 'lemaskinvalid' .OR. &
   this%sub_type == 'ltmaskinvalid' .OR. this%sub_type == 'gemaskinvalid' .OR. &
   this%sub_type == 'gtmaskinvalid' .OR. this%sub_type == 'index') THEN
! nothing to do here
  ELSE
    CALL sub_type_error()
    RETURN
  ENDIF

ELSE
  CALL trans_type_error()
  RETURN
ENDIF

CONTAINS

SUBROUTINE sub_type_error()

CALL l4f_category_log(this%category, l4f_error, trim(this%trans_type) &
 //': sub_type '//trim(this%sub_type)//' is not defined')
CALL raise_fatal_error()

END SUBROUTINE sub_type_error

SUBROUTINE trans_type_error()

CALL l4f_category_log(this%category, l4f_error, 'trans_type '//this%trans_type &
 //' is not defined')
CALL raise_fatal_error()

END SUBROUTINE trans_type_error


END SUBROUTINE transform_init


SUBROUTINE transform_delete(this)
TYPE(transform_def),INTENT(inout) :: this

this%trans_type=cmiss
this%sub_type=cmiss

this%rect_ind%ix=imiss
this%rect_ind%iy=imiss
this%rect_ind%fx=imiss
this%rect_ind%fy=imiss

this%rect_coo%ilon=dmiss
this%rect_coo%ilat=dmiss
this%rect_coo%flon=dmiss
this%rect_coo%flat=dmiss

this%box_info%npx=imiss
this%box_info%npy=imiss

this%extrap=.false.

!chiudo il logger
CALL l4f_category_delete(this%category)

END SUBROUTINE transform_delete


SUBROUTINE transform_get_val(this, time_definition, trans_type, sub_type, &
 input_levtype, output_levtype)
type(transform_def),intent(in) :: this
INTEGER,INTENT(out),OPTIONAL :: time_definition
CHARACTER(len=*),INTENT(out),OPTIONAL :: trans_type
CHARACTER(len=*),INTENT(out),OPTIONAL :: sub_type
TYPE(vol7d_level),INTENT(out),OPTIONAL :: input_levtype

TYPE(vol7d_level),INTENT(out),OPTIONAL :: output_levtype


IF (PRESENT(time_definition)) time_definition=this%time_definition
IF (PRESENT(trans_type)) trans_type = this%trans_type
IF (PRESENT(sub_type)) sub_type = this%sub_type
IF (PRESENT(input_levtype)) input_levtype = this%vertint%input_levtype
IF (PRESENT(output_levtype)) output_levtype = this%vertint%output_levtype


END SUBROUTINE transform_get_val


SUBROUTINE grid_transform_levtype_levtype_init(this, trans, lev_in, lev_out, &
 coord_3d_in, categoryappend)
TYPE(grid_transform),INTENT(out) :: this
TYPE(transform_def),INTENT(in) :: trans
TYPE(vol7d_level),INTENT(in) :: lev_in(:)
TYPE(vol7d_level),INTENT(in) :: lev_out(:)
REAL,INTENT(inout),OPTIONAL,ALLOCATABLE :: coord_3d_in(:,:,:)
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

DOUBLE PRECISION :: coord_in(SIZE(lev_in))
DOUBLE PRECISION,ALLOCATABLE :: coord_out(:)
LOGICAL :: mask_in(SIZE(lev_in))
LOGICAL,ALLOCATABLE :: mask_out(:)
LOGICAL :: dolog
INTEGER :: i, j, icache, inused, istart, iend, ostart, oend


CALL grid_transform_init_common(this, trans, categoryappend)
#ifdef DEBUG
CALL l4f_category_log(this%category, l4f_debug, "grid_transform vertint")
#endif

IF (this%trans%trans_type == 'vertint') THEN

  IF (c_e(trans%vertint%input_levtype%level2) .AND. &
   trans%vertint%input_levtype%level1 /= trans%vertint%input_levtype%level2) THEN
    CALL l4f_category_log(this%category, l4f_error, &
     'vertint: input upper and lower surface must be of the same type, '// &
     t2c(trans%vertint%input_levtype%level1)//'/='// &
     t2c(trans%vertint%input_levtype%level2))
    CALL raise_error()
    RETURN
  ENDIF
  IF (c_e(trans%vertint%output_levtype%level2) .AND. &
   trans%vertint%output_levtype%level1 /= trans%vertint%output_levtype%level2) THEN
    CALL l4f_category_log(this%category, l4f_error, &
     'vertint: output upper and lower surface must be of the same type'// &
     t2c(trans%vertint%output_levtype%level1)//'/='// &
     t2c(trans%vertint%output_levtype%level2))
    CALL raise_error()
    RETURN
  ENDIF

  mask_in(:) = (lev_in(:)%level1 == trans%vertint%input_levtype%level1) .AND. &
   (lev_in(:)%level2 == trans%vertint%input_levtype%level2)
  CALL make_vert_coord(lev_in, mask_in, coord_in, dolog)
  this%innz = SIZE(lev_in)
  istart = firsttrue(mask_in)
  iend = lasttrue(mask_in)
  inused = iend - istart + 1
  IF (inused /= count(mask_in)) THEN
    CALL l4f_category_log(this%category, l4f_error, &
     'grid_transform_levtype_levtype_init: input levels badly sorted '//&
     t2c(inused)//'/'//t2c(count(mask_in)))
    CALL raise_error()
    RETURN
  ENDIF
  this%levshift = istart-1
  this%levused = inused

  IF (trans%vertint%input_levtype%level1 /= trans%vertint%output_levtype%level1) THEN
#ifdef DEBUG
    CALL l4f_category_log(this%category, l4f_debug, &
     'vertint: different input and output level types '// &
     t2c(trans%vertint%input_levtype%level1)//' '// &
     t2c(trans%vertint%output_levtype%level1))
#endif

    ALLOCATE(mask_out(SIZE(lev_out)), this%vcoord_out(SIZE(lev_out)))
    mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
     (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
    CALL make_vert_coord(lev_out, mask_out, this%vcoord_out, dolog)
    this%outnz = SIZE(mask_out)
    DEALLOCATE(mask_out)

    IF (.NOT.PRESENT(coord_3d_in)) THEN
      CALL l4f_category_log(this%category, l4f_warn, &
       'vertint: different input and output level types &
       &and no coord_3d_in, expecting vert. coord. in volume')
      this%dolog = dolog ! a little bit dirty, I must compute log later
    ELSE
      IF (SIZE(coord_3d_in,3) /= inused) THEN
        CALL l4f_category_log(this%category, l4f_error, &
         'vertint: vertical size of coord_3d_in (vertical coordinate) &
         &different from number of input levels suitable for interpolation')
        CALL l4f_category_log(this%category, l4f_error, &
         'coord_3d_in: '//t2c(SIZE(coord_3d_in,3))// &
         ', input levels for interpolation: '//t2c(inused))
        CALL raise_error()
        RETURN
      ENDIF

      CALL move_alloc(coord_3d_in, this%coord_3d_in) ! steal allocation
      IF (dolog) THEN
        WHERE(c_e(this%coord_3d_in) .AND. this%coord_3d_in > 0.0)
          this%coord_3d_in = log(this%coord_3d_in)
        ELSE WHERE
          this%coord_3d_in = rmiss
        END WHERE
      ENDIF
    ENDIF

    this%valid = .true. ! warning, no check of subtype

  ELSE
! here we assume that valid levels are contiguous and ordered

#ifdef DEBUG
    CALL l4f_category_log(this%category, l4f_debug, &
     'vertint: equal input and output level types '// &
     t2c(trans%vertint%input_levtype%level1))
#endif

    IF (SIZE(lev_out) > 0) THEN ! output level list provided
      ALLOCATE(mask_out(SIZE(lev_out)), coord_out(SIZE(lev_out)))
      mask_out(:) = (lev_out(:)%level1 == trans%vertint%output_levtype%level1) .AND. &
       (lev_out(:)%level2 == trans%vertint%output_levtype%level2)
      CALL make_vert_coord(lev_out, mask_out, coord_out, dolog)

    ELSE ! output level list not provided, try to autogenerate
      IF (c_e(trans%vertint%input_levtype%level2) .AND. &
       .NOT.c_e(trans%vertint%output_levtype%level2)) THEN ! full -> half
        IF (trans%vertint%output_levtype%level1 == 105 .OR. &
         trans%vertint%output_levtype%level1 == 150) THEN
          ALLOCATE(this%output_level_auto(inused-1))
          CALL l4f_category_log(this%category,l4f_info, &
           'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
           //'/'//t2c(iend-istart)//' output levels (f->h)')
          DO i = istart, iend - 1
            CALL init(this%output_level_auto(i-istart+1), &
             trans%vertint%input_levtype%level1, lev_in(i)%l2)
          ENDDO
        ELSE
          CALL l4f_category_log(this%category, l4f_error, &
           'grid_transform_levtype_levtype_init: automatic generation of output levels &
           &available only for hybrid levels')
          CALL raise_error()
          RETURN
        ENDIF
      ELSE IF (.NOT.c_e(trans%vertint%input_levtype%level2) .AND. &
       c_e(trans%vertint%output_levtype%level2)) THEN ! half -> full
        ALLOCATE(this%output_level_auto(inused-1))
        IF (trans%vertint%output_levtype%level1 == 105 .OR. &
         trans%vertint%output_levtype%level1 == 150) THEN
          CALL l4f_category_log(this%category,l4f_info, &
           'grid_transform_levtype_levtype_init: autogenerating '//t2c(inused-1) &
           //'/'//t2c(iend-istart)//' output levels (h->f)')
          DO i = istart, iend - 1
            CALL init(this%output_level_auto(i-istart+1), trans%vertint%input_levtype%level1, &
             lev_in(i)%l1, trans%vertint%input_levtype%level1, &
             lev_in(i)%l1+1)
          ENDDO
        ELSE
          CALL l4f_category_log(this%category, l4f_error, &
           'grid_transform_levtype_levtype_init: automatic generation of output levels &
           &available only for hybrid levels')
          CALL raise_error()
          RETURN
        ENDIF
      ELSE
        CALL l4f_category_log(this%category, l4f_error, &
         'grid_transform_levtype_levtype_init: strange situation'// &
         to_char(c_e(trans%vertint%input_levtype%level2))//' '// &
         to_char(c_e(trans%vertint%output_levtype%level2)))
        CALL raise_error()
        RETURN
      ENDIF
      ALLOCATE(coord_out(inused-1), mask_out(inused-1))
      mask_out(:) = .true.
      CALL make_vert_coord(this%output_level_auto, mask_out, coord_out, dolog)
    ENDIF

    this%outnz = SIZE(mask_out)
    ostart = firsttrue(mask_out)
    oend = lasttrue(mask_out)

! set valid = .FALSE. here?
    IF (istart == 0) THEN
      CALL l4f_category_log(this%category, l4f_warn, &
       'grid_transform_levtype_levtype_init: &
       &input contains no vertical levels of type ('// &
       trim(to_char(trans%vertint%input_levtype%level1))//','// &
       trim(to_char(trans%vertint%input_levtype%level2))// &
       ') suitable for interpolation')
      RETURN
!      iend = -1 ! for loops
    ELSE IF (istart == iend) THEN
      CALL l4f_category_log(this%category, l4f_warn, &
       'grid_transform_levtype_levtype_init: &
       &input contains only 1 vertical level of type ('// &
       trim(to_char(trans%vertint%input_levtype%level1))//','// &
       trim(to_char(trans%vertint%input_levtype%level2))// &
       ') suitable for interpolation')
    ENDIF
    IF (ostart == 0) THEN
      CALL l4f_category_log(this%category, l4f_warn, &
       'grid_transform_levtype_levtype_init: &
       &output contains no vertical levels of type ('// &
       trim(to_char(trans%vertint%output_levtype%level1))//','// &
       trim(to_char(trans%vertint%output_levtype%level2))// &
       ') suitable for interpolation')
      RETURN
!      oend = -1 ! for loops
    ENDIF

! end of code common for all vertint subtypes
    IF (this%trans%sub_type == 'linear') THEN

      ALLOCATE(this%inter_index_z(this%outnz), this%inter_zp(this%outnz))
      this%inter_index_z(:) = imiss
      this%inter_zp(:) = dmiss
      IF (this%trans%extrap .AND. istart > 0) THEN
        WHERE(mask_out)
! extrapolate down by default
          this%inter_index_z(:) = istart
          this%inter_zp(:) = 1.0d0
        ENDWHERE
      ENDIF

      icache = istart + 1
      outlev: DO j = ostart, oend
        inlev: DO i = icache, iend
          IF (coord_in(i) >= coord_out(j)) THEN
            IF (coord_out(j) >= coord_in(i-1)) THEN
              this%inter_index_z(j) = i - 1
              this%inter_zp(j) = (coord_out(j)-coord_in(i-1)) / &
               (coord_in(i)-coord_in(i-1)) ! weight for (i)
              icache = i ! speedup next j iteration
            ENDIF
            cycle outlev ! found or extrapolated down
          ENDIF
        ENDDO inlev
! if I'm here I must extrapolate up
        IF (this%trans%extrap .AND. iend > 1) THEN
          this%inter_index_z(j) = iend - 1
          this%inter_zp(j) = 0.0d0
        ENDIF
      ENDDO outlev

      DEALLOCATE(coord_out, mask_out)
      this%valid = .true.

    ELSE IF (this%trans%sub_type == 'linearsparse') THEN
! just store vertical coordinates, dirty work is done later
      ALLOCATE(this%vcoord_in(this%levused),  this%vcoord_out(SIZE(coord_out)))
      this%vcoord_in(:) = coord_in(this%levshift+1:this%levshift+this%levused)
      this%vcoord_out(:) = coord_out(:)
      DEALLOCATE(coord_out, mask_out)
      this%valid = .true.

    ENDIF

  ENDIF ! levels are different

!ELSE IF (this%trans%trans_type == 'verttrans') THEN

ENDIF

END SUBROUTINE grid_transform_levtype_levtype_init


! internal subroutine for computing vertical coordinate values, for
! pressure-based coordinates the logarithm is computed
SUBROUTINE make_vert_coord(lev, mask, coord, dolog)
TYPE(vol7d_level),INTENT(in) :: lev(:)
LOGICAL,INTENT(inout) :: mask(:)
DOUBLE PRECISION,INTENT(out) :: coord(:)
LOGICAL,INTENT(out) :: dolog

INTEGER :: k
DOUBLE PRECISION :: fact

dolog = .false.
k = firsttrue(mask)
IF (k <= 0) RETURN
coord(:) = dmiss

IF (any(lev(k)%level1 == height_level)) THEN ! improve with a conversion table somewhere
  fact = 1.0d-3
ELSE IF (any(lev(k)%level1 == thermo_level)) THEN
  fact = 1.0d-1
ELSE IF (any(lev(k)%level1 == sigma_level)) THEN
  fact = 1.0d-4
ELSE
  fact = 1.0d0
ENDIF

IF (c_e(lev(k)%level2) .AND. lev(k)%level1 == lev(k)%level2) THEN ! layer between 2 levels
  IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN ! pressure, compute log
    WHERE(mask(:) .AND. lev(:)%l1 > 0 .AND. lev(:)%l2 > 0)
      coord(:) = (log(dble(lev(:)%l1)*fact) + log(dble(lev(:)%l2)*fact))*0.5d0
    END WHERE
    dolog = .true.
  ELSE
    WHERE(mask(:))
      coord(:) = (lev(:)%l1 + lev(:)%l2)*fact*0.5d0
    END WHERE
  ENDIF
ELSE ! half level
  IF (lev(k)%level1 == 100 .OR. lev(k)%level1 == 108) THEN ! pressure, compute log
    WHERE(mask(:) .AND. lev(:)%l1 > 0)
      coord(:) = log(dble(lev(:)%l1)*fact)
    END WHERE
    dolog = .true.
  ELSE
    WHERE(mask(:))
      coord(:) = lev(:)%l1*fact
    END WHERE
  ENDIF
ENDIF

! refine mask
mask(:) = mask(:) .AND. c_e(coord(:))

END SUBROUTINE make_vert_coord


SUBROUTINE grid_transform_init(this, trans, in, out, maskgrid, maskbounds, &
 categoryappend)
TYPE(grid_transform),INTENT(out) :: this
TYPE(transform_def),INTENT(in) :: trans
TYPE(griddim_def),INTENT(inout) :: in
TYPE(griddim_def),INTENT(inout) :: out
REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
REAL,INTENT(in),OPTIONAL :: maskbounds(:)
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

INTEGER :: nx, ny, i, j, ix, iy, n, nm, nr, cf_i, cf_o, nprev, &
 xnmin, xnmax, ynmin, ynmax
DOUBLE PRECISION :: xmin, xmax, ymin, ymax, steplon, steplat, &
 xmin_new, ymin_new, ellips_smaj_axis, ellips_flatt, r2
TYPE(geo_proj) :: proj_in, proj_out
TYPE(georef_coord) :: point
LOGICAL,ALLOCATABLE :: point_mask(:,:)
TYPE(griddim_def) :: lin, lout


CALL grid_transform_init_common(this, trans, categoryappend)
#ifdef DEBUG
CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-vg6d")
#endif

! output ellipsoid has to be the same as for the input (improve)
CALL get_val(in, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
CALL set_val(out, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)

IF (this%trans%trans_type == 'zoom') THEN

  IF (this%trans%sub_type == 'coord') THEN

    CALL griddim_zoom_coord(in, &
     this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
     this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
     this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
     this%trans%rect_ind%fx, this%trans%rect_ind%fy)

  ELSE IF (this%trans%sub_type == 'projcoord') THEN

    CALL griddim_zoom_projcoord(in, &
     this%trans%rect_coo%ilon, this%trans%rect_coo%ilat,&
     this%trans%rect_coo%flon, this%trans%rect_coo%flat,&
     this%trans%rect_ind%ix, this%trans%rect_ind%iy, &
     this%trans%rect_ind%fx, this%trans%rect_ind%fy)

  ELSE IF (this%trans%sub_type == 'coordbb') THEN

! compute coordinates of input grid in geo system
    CALL copy(in, lin)
    CALL unproj(lin)
    CALL get_val(lin, nx=nx, ny=ny)

    ALLOCATE(point_mask(nx,ny))
    point_mask(:,:) = .false.

! mark points falling into requested bounding-box
    DO j = 1, ny
      DO i = 1, nx
!        IF (geo_coord_inside_rectang())
        IF (lin%dim%lon(i,j) > this%trans%rect_coo%ilon .AND. &
         lin%dim%lon(i,j) < this%trans%rect_coo%flon .AND. &
         lin%dim%lat(i,j) > this%trans%rect_coo%ilat .AND. &
         lin%dim%lat(i,j) < this%trans%rect_coo%flat) THEN ! improve!
          point_mask(i,j) = .true.
        ENDIF
      ENDDO
    ENDDO

! determine cut indices keeping all points which fall inside b-b
    DO i = 1, nx
      IF (any(point_mask(i,:))) EXIT
    ENDDO
    this%trans%rect_ind%ix = i
    DO i = nx, this%trans%rect_ind%ix, -1
      IF (any(point_mask(i,:))) EXIT
    ENDDO
    this%trans%rect_ind%fx = i

    DO j = 1, ny
      IF (any(point_mask(:,j))) EXIT
    ENDDO
    this%trans%rect_ind%iy = j
    DO j = ny, this%trans%rect_ind%iy, -1
      IF (any(point_mask(:,j))) EXIT
    ENDDO
    this%trans%rect_ind%fy = j

    DEALLOCATE(point_mask)

    IF (this%trans%rect_ind%ix > this%trans%rect_ind%fx .OR. &
     this%trans%rect_ind%iy > this%trans%rect_ind%fy) THEN

      CALL l4f_category_log(this%category,l4f_error, &
       "zoom coordbb: no points inside bounding box "//&
       trim(to_char(this%trans%rect_coo%ilon))//","// &
       trim(to_char(this%trans%rect_coo%flon))//","// &
       trim(to_char(this%trans%rect_coo%ilat))//","// &
       trim(to_char(this%trans%rect_coo%flat)))
      CALL raise_error()
      RETURN

    ENDIF
    CALL delete(lin)
  ENDIF
! to do in all zoom cases

  CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
   ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)

! old indices
  this%iniox = min(max(this%trans%rect_ind%ix,1),nx) ! iox
  this%inioy = min(max(this%trans%rect_ind%iy,1),ny) ! ioy
  this%infox = max(min(this%trans%rect_ind%fx,nx),1) ! fox
  this%infoy = max(min(this%trans%rect_ind%fy,ny),1) ! foy
! new indices
  this%outinx = min(max(2-this%trans%rect_ind%ix,1),nx) ! inx
  this%outiny = min(max(2-this%trans%rect_ind%iy,1),ny) ! iny
  this%outfnx = min(this%trans%rect_ind%fx,nx)-this%trans%rect_ind%ix+1 ! fnx
  this%outfny = min(this%trans%rect_ind%fy,ny)-this%trans%rect_ind%iy+1 ! fny

  xmin=xmin+steplon*(this%trans%rect_ind%ix-1)
  ymin=ymin+steplat*(this%trans%rect_ind%iy-1)
  xmax=xmax+steplon*(this%trans%rect_ind%fx-nx)
  ymax=ymax+steplat*(this%trans%rect_ind%fy-ny)

  CALL copy(in, out)
! deallocate coordinates if allocated because they will change
  CALL dealloc(out%dim)

  out%dim%nx = this%trans%rect_ind%fx - this%trans%rect_ind%ix + 1 ! newx
  out%dim%ny = this%trans%rect_ind%fy - this%trans%rect_ind%iy + 1 ! newy
  this%outnx = out%dim%nx
  this%outny = out%dim%ny
  this%innx = nx
  this%inny = ny

  CALL set_val(out, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)

  this%valid = .true. ! warning, no check of subtype
  
ELSE IF (this%trans%trans_type == 'boxregrid') THEN

  CALL get_val(in, nx=nx, ny=ny, xmin=xmin, xmax=xmax, &
   ymin=ymin, ymax=ymax, dx=steplon, dy=steplat)

  this%innx = nx
  this%inny = ny

! new grid
  xmin_new = xmin + (this%trans%box_info%npx - 1)*0.5d0*steplon
  ymin_new = ymin + (this%trans%box_info%npy - 1)*0.5d0*steplat

  CALL copy(in, out)
! deallocate coordinates if allocated because they will change
  CALL dealloc(out%dim)

  out%dim%nx = nx/this%trans%box_info%npx
  out%dim%ny = ny/this%trans%box_info%npy
  this%outnx = out%dim%nx
  this%outny = out%dim%ny
  steplon = steplon*this%trans%box_info%npx
  steplat = steplat*this%trans%box_info%npy

  CALL set_val(out, xmin=xmin_new, ymin=ymin_new, &
   xmax=xmin_new + dble(out%dim%nx-1)*steplon, dx=steplon, &
   ymax=ymin_new + dble(out%dim%ny-1)*steplat, dy=steplat)

  this%valid = .true. ! warning, no check of subtype

ELSE IF (this%trans%trans_type == 'inter' .OR. this%trans%trans_type == 'intersearch') THEN

  CALL outgrid_setup() ! common setup for grid-generating methods

  IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin'&
    .OR. this%trans%sub_type == 'shapiro_near') THEN

    CALL get_val(in, nx=this%innx, ny=this%inny, &
     xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
    CALL get_val(out, nx=this%outnx, ny=this%outny)

    ALLOCATE(this%inter_index_x(this%outnx,this%outny), &
     this%inter_index_y(this%outnx,this%outny))
    CALL copy(out, lout)
    CALL unproj(lout)

    IF (this%trans%sub_type == 'bilin') THEN
      CALL this%find_index(in, .false., &
       this%innx, this%inny, xmin, xmax, ymin, ymax, &
       lout%dim%lon, lout%dim%lat, this%trans%extrap, &
       this%inter_index_x, this%inter_index_y)

      ALLOCATE(this%inter_x(this%innx,this%inny), &
       this%inter_y(this%innx,this%inny))
      ALLOCATE(this%inter_xp(this%outnx,this%outny), &
       this%inter_yp(this%outnx,this%outny))

! compute coordinates of input grid
      CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
! compute coordinates of output grid in input system
      CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)

    ELSE ! near, shapiro_near
      CALL this%find_index(in, .true., &
       this%innx, this%inny, xmin, xmax, ymin, ymax, &
       lout%dim%lon, lout%dim%lat, this%trans%extrap, &
       this%inter_index_x, this%inter_index_y)

      IF (this%trans%trans_type == 'intersearch') THEN ! replicate code above
        ALLOCATE(this%inter_x(this%innx,this%inny), &
         this%inter_y(this%innx,this%inny))
        ALLOCATE(this%inter_xp(this%outnx,this%outny), &
         this%inter_yp(this%outnx,this%outny))

! compute coordinates of input grid
        CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
! compute coordinates of output grid in input system
        CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
      ENDIF
    ENDIF

    CALL delete(lout)
    this%valid = .true.
  ENDIF

ELSE IF (this%trans%trans_type == 'boxinter') THEN

  CALL outgrid_setup() ! common setup for grid-generating methods
  CALL get_val(in, nx=this%innx, ny=this%inny)
  CALL get_val(out, nx=this%outnx, ny=this%outny, &
   xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
! TODO now box size is ignored
! if box size not provided, use the actual grid step
  IF (.NOT.c_e(this%trans%area_info%boxdx)) &
   CALL get_val(out, dx=this%trans%area_info%boxdx)
  IF (.NOT.c_e(this%trans%area_info%boxdy)) &
   CALL get_val(out, dx=this%trans%area_info%boxdy)
! half size is actually needed
  this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
  this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
! unlike before, here index arrays must have the shape of input grid
  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
   this%inter_index_y(this%innx,this%inny))

! compute coordinates of input grid in geo system
  CALL copy(in, lin)
  CALL unproj(lin)
! use find_index in the opposite way, here extrap does not make sense
  CALL this%find_index(out, .true., &
   this%outnx, this%outny, xmin, xmax, ymin, ymax, &
   lin%dim%lon, lin%dim%lat, .false., &
   this%inter_index_x, this%inter_index_y)

  CALL delete(lin)
  this%valid = .true. ! warning, no check of subtype

ELSE IF (this%trans%trans_type == 'stencilinter') THEN

  CALL outgrid_setup() ! common setup for grid-generating methods
! from inter:near
  CALL get_val(in, nx=this%innx, ny=this%inny, &
   xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
  CALL get_val(out, nx=this%outnx, ny=this%outny)

  ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
   this%inter_index_y(this%outnx,this%outny))
  CALL copy(out, lout)
  CALL unproj(lout)
  CALL this%find_index(in, .true., &
   this%innx, this%inny, xmin, xmax, ymin, ymax, &
   lout%dim%lon, lout%dim%lat, this%trans%extrap, &
   this%inter_index_x, this%inter_index_y)

! define the stencil mask
  nr = int(this%trans%area_info%radius) ! integer radius
  n = nr*2+1 ! n. of points
  nm = nr + 1 ! becomes index of center
  r2 = this%trans%area_info%radius**2
  ALLOCATE(this%stencil(n,n))
  this%stencil(:,:) = .true.
  DO iy = 1, n
    DO ix = 1, n
      IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
    ENDDO
  ENDDO

! set to missing the elements of inter_index too close to the domain
! borders
  xnmin = nr + 1
  xnmax = this%innx - nr
  ynmin = nr + 1
  ynmax = this%inny - nr
  DO iy = 1, this%outny
    DO ix = 1, this%outnx
      IF (this%inter_index_x(ix,iy) < xnmin .OR. &
       this%inter_index_x(ix,iy) > xnmax .OR. &
       this%inter_index_y(ix,iy) < ynmin .OR. &
       this%inter_index_y(ix,iy) > ynmax) THEN
        this%inter_index_x(ix,iy) = imiss
        this%inter_index_y(ix,iy) = imiss
      ENDIF
    ENDDO
  ENDDO

#ifdef DEBUG
  CALL l4f_category_log(this%category, l4f_debug, &
   'stencilinter: stencil size '//t2c(n*n))
  CALL l4f_category_log(this%category, l4f_debug, &
   'stencilinter: stencil points '//t2c(count(this%stencil)))
#endif

  CALL delete(lout)
  this%valid = .true. ! warning, no check of subtype

ELSE IF (this%trans%trans_type == 'maskgen') THEN

  IF (this%trans%sub_type == 'poly') THEN

    CALL copy(in, out)
    CALL get_val(in, nx=this%innx, ny=this%inny)
    this%outnx = this%innx
    this%outny = this%inny

! unlike before, here index arrays must have the shape of input grid
    ALLOCATE(this%inter_index_x(this%innx,this%inny), &
     this%inter_index_y(this%innx,this%inny))
    this%inter_index_x(:,:) = imiss
    this%inter_index_y(:,:) = 1

! compute coordinates of input grid in geo system
    CALL unproj(out) ! should be unproj(lin)

    nprev = 1
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev)
    DO j = 1, this%inny
      inside_x: DO i = 1, this%innx
        point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))

        DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
          IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
            this%inter_index_x(i,j) = n
            nprev = n
            cycle inside_x
          ENDIF
        ENDDO
        DO n = nprev-1, 1, -1 ! test the other polygons
          IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
            this%inter_index_x(i,j) = n
            nprev = n
            cycle inside_x
          ENDIF
        ENDDO

!     CALL delete(point) ! speedup
      ENDDO inside_x
    ENDDO
!$OMP END PARALLEL

  ELSE IF (this%trans%sub_type == 'grid') THEN
! here out(put grid) is abused for indicating the box-generating grid
! but the real output grid is the input grid
    CALL copy(out, lout) ! save out for local use
    CALL delete(out) ! needed before copy
    CALL copy(in, out)
    CALL get_val(in, nx=this%innx, ny=this%inny)
    this%outnx = this%innx
    this%outny = this%inny
    CALL get_val(lout, nx=nx, ny=ny, &
     xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)

! unlike before, here index arrays must have the shape of input grid
    ALLOCATE(this%inter_index_x(this%innx,this%inny), &
     this%inter_index_y(this%innx,this%inny))

! compute coordinates of input/output grid in geo system
    CALL unproj(out)

! use find_index in the opposite way, here extrap does not make sense
    CALL this%find_index(lout, .true., &
     nx, ny, xmin, xmax, ymin, ymax, &
     out%dim%lon, out%dim%lat, .false., &
     this%inter_index_x, this%inter_index_y)
! transform indices to 1-d for mask generation
    WHERE(c_e(this%inter_index_x(:,:)))
      this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
       (this%inter_index_y(:,:)-1)*nx
    END WHERE

    CALL delete(lout)
  ENDIF

  this%valid = .true.

ELSE IF (this%trans%trans_type == 'polyinter') THEN

! this is the only difference wrt maskgen:poly
   this%recur = .true. ! grid-to-grid polyinter is done in two steps!

  CALL copy(in, out)
  CALL get_val(in, nx=this%innx, ny=this%inny)
  this%outnx = this%innx
  this%outny = this%inny

! unlike before, here index arrays must have the shape of input grid
  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
   this%inter_index_y(this%innx,this%inny))
  this%inter_index_x(:,:) = imiss
  this%inter_index_y(:,:) = 1

! compute coordinates of input grid in geo system
  CALL unproj(out) ! should be unproj(lin)

  nprev = 1
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev)
  DO j = 1, this%inny
    inside_x_2: DO i = 1, this%innx
      point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))

      DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
        IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
          this%inter_index_x(i,j) = n
          nprev = n
          cycle inside_x_2
        ENDIF
      ENDDO
      DO n = nprev-1, 1, -1 ! test the other polygons
        IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
          this%inter_index_x(i,j) = n
          nprev = n
          cycle inside_x_2
        ENDIF
      ENDDO

!     CALL delete(point) ! speedup
    ENDDO inside_x_2
  ENDDO
!$OMP END PARALLEL

  this%valid = .true. ! warning, no check of subtype

ELSE IF (this%trans%trans_type == 'metamorphosis') THEN

  CALL copy(in, out)
  CALL get_val(in, nx=this%innx, ny=this%inny)
  this%outnx = this%innx
  this%outny = this%inny

  IF (this%trans%sub_type == 'maskvalid' .OR. this%trans%sub_type == 'maskinvalid') THEN

    IF (.NOT.PRESENT(maskgrid)) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init maskgrid argument missing for metamorphosis:'// &
       trim(this%trans%sub_type)//' transformation')
      CALL raise_error()
      RETURN
    ENDIF

    IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init mask non conformal with input field')
      CALL l4f_category_log(this%category,l4f_error, &
       'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
       ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
      CALL raise_error()
      RETURN
    ENDIF

    ALLOCATE(this%point_mask(this%innx,this%inny))

    IF (this%trans%sub_type == 'maskvalid') THEN
! behavior depends on the presence/usability of maskbounds,
! simplified wrt its use in metamorphosis:mask
      IF (.NOT.PRESENT(maskbounds)) THEN
        this%point_mask(:,:) = c_e(maskgrid(:,:))
      ELSE IF (SIZE(maskbounds) < 2) THEN
        this%point_mask(:,:) = c_e(maskgrid(:,:))
      ELSE
        this%point_mask(:,:) = c_e(maskgrid(:,:)) .AND. &
         maskgrid(:,:) > maskbounds(1) .AND. &
         maskgrid(:,:) <= maskbounds(SIZE(maskbounds))
      ENDIF
    ELSE ! reverse condition
      this%point_mask(:,:) = .NOT.c_e(maskgrid(:,:))
    ENDIF

    this%valid = .true.

  ELSE IF (this%trans%sub_type == 'lemaskinvalid' .OR. &
   this%trans%sub_type == 'ltmaskinvalid' .OR. &
   this%trans%sub_type == 'gemaskinvalid' .OR. &
   this%trans%sub_type == 'gtmaskinvalid') THEN
! here i can only store field for computing mask runtime

    this%val_mask = maskgrid
    this%valid = .true.

  ELSE IF (this%trans%sub_type == 'setinvalidto') THEN

    IF (.NOT.PRESENT(maskbounds)) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init maskbounds missing for metamorphosis:'// &
       trim(this%trans%sub_type)//' transformation')
      RETURN
    ELSE IF (SIZE(maskbounds) < 1) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init maskbounds empty for metamorphosis:'// &
       trim(this%trans%sub_type)//' transformation')
      RETURN
    ELSE
      this%val1 = maskbounds(1)
#ifdef DEBUG
      CALL l4f_category_log(this%category, l4f_debug, &
       "grid_transform_init setting invalid data to "//t2c(this%val1))
#endif
    ENDIF

    this%valid = .true.

  ELSE IF (this%trans%sub_type == 'settoinvalid') THEN

    IF (.NOT.PRESENT(maskbounds)) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init maskbounds missing for metamorphosis:'// &
       trim(this%trans%sub_type)//' transformation')
      CALL raise_error()
      RETURN
    ELSE IF (SIZE(maskbounds) < 2) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
       trim(this%trans%sub_type)//' transformation')
      CALL raise_error()
      RETURN
    ELSE
      this%val1 = maskbounds(1)
      this%val2 = maskbounds(SIZE(maskbounds))
#ifdef DEBUG
      CALL l4f_category_log(this%category, l4f_debug, &
       "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//','// &
       t2c(this%val2)//']')
#endif
    ENDIF

    this%valid = .true.

  ENDIF

ENDIF

CONTAINS

! local subroutine to be called by all methods interpolating to a new
! grid, no parameters passed, used as a macro to avoid repeating code
SUBROUTINE outgrid_setup()

! set increments in new grid in order for all the baraque to work
CALL griddim_setsteps(out)
! check component flag
CALL get_val(in, proj=proj_in, component_flag=cf_i)
CALL get_val(out, proj=proj_out, component_flag=cf_o)
IF (proj_in == proj_out) THEN
! same projection: set output component flag equal to input regardless
! of its value
  CALL set_val(out, component_flag=cf_i)
ELSE
! different projection, interpolation possible only with vector data
! referred to geograpical axes
  IF (cf_i == 1) THEN
    CALL l4f_category_log(this%category,l4f_warn, &
     'trying to interpolate a grid with component flag 1 to a grid on a different projection')
    CALL l4f_category_log(this%category,l4f_warn, &
     'vector fields will probably be wrong')
  ELSE
    CALL set_val(out, component_flag=cf_i)
  ENDIF
ENDIF
! rotate the input grid by n*360 degrees to bring it closer to the output grid
CALL griddim_set_central_lon(in, griddim_central_lon(out))

END SUBROUTINE outgrid_setup

END SUBROUTINE grid_transform_init


SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
 maskgrid, maskbounds, find_index, categoryappend)
TYPE(grid_transform),INTENT(out) :: this
TYPE(transform_def),INTENT(in) :: trans
TYPE(griddim_def),INTENT(in) :: in
TYPE(vol7d),INTENT(inout) :: v7d_out
REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
REAL,INTENT(in),OPTIONAL :: maskbounds(:)
PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
 time_definition
DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
DOUBLE PRECISION,ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
REAL,ALLOCATABLE :: lmaskbounds(:)
TYPE(georef_coord) :: point
TYPE(griddim_def) :: lin
!$ INTEGER :: outnx

IF (PRESENT(find_index)) THEN ! move in init_common?
  IF (ASSOCIATED(find_index)) THEN
    this%find_index => find_index
  ENDIF
ENDIF
CALL grid_transform_init_common(this, trans, categoryappend)
#ifdef DEBUG
CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-v7d")
#endif

! used after in some transformations
CALL get_val(trans, time_definition=time_definition)
IF (.NOT. c_e(time_definition)) THEN
  time_definition=1  ! default to validity time
ENDIF

IF (this%trans%trans_type == 'inter') THEN

  IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin' &
   .OR. this%trans%sub_type == 'shapiro_near') THEN

    CALL copy(in, lin)
    CALL get_val(lin, nx=this%innx, ny=this%inny)
    this%outnx = SIZE(v7d_out%ana)
    this%outny = 1

    ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
     this%inter_index_y(this%outnx,this%outny))
    ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))

    CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
! equalize in/out coordinates
    lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
    CALL griddim_set_central_lon(lin, lonref)
    CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)

    IF (this%trans%sub_type == 'bilin') THEN
      CALL this%find_index(lin, .false., &
       this%innx, this%inny, xmin, xmax, ymin, ymax, &
       lon, lat, this%trans%extrap, &
       this%inter_index_x, this%inter_index_y)

      ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
      ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))

      CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
      CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)

    ELSE ! near shapiro_near
      CALL this%find_index(lin, .true., &
       this%innx, this%inny, xmin, xmax, ymin, ymax, &
       lon, lat, this%trans%extrap, &
       this%inter_index_x, this%inter_index_y)

      IF (this%trans%trans_type == 'intersearch') THEN ! replicate code above
        ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
        ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))

        CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
        CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
      ENDIF
    ENDIF

    DEALLOCATE(lon,lat)
    CALL delete(lin)

    this%valid = .true.

  ENDIF

ELSE IF (this%trans%trans_type == 'polyinter') THEN

  CALL get_val(in, nx=this%innx, ny=this%inny)
! unlike before, here index arrays must have the shape of input grid
  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
   this%inter_index_y(this%innx,this%inny))
  this%inter_index_x(:,:) = imiss
  this%inter_index_y(:,:) = 1

! compute coordinates of input grid in geo system
  CALL copy(in, lin)
  CALL unproj(lin)

  this%outnx = this%trans%poly%arraysize
  this%outny = 1
  CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
  CALL init(v7d_out, time_definition=time_definition)
  CALL vol7d_alloc(v7d_out, nana=this%outnx)

! equalize in/out coordinates
  ALLOCATE(lon(this%outnx,1))
  CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
  CALL griddim_set_central_lon(lin, lonref)
  DEALLOCATE(lon)

! setup output point list, equal to average of polygon points
! warning, in case of concave areas points may coincide!
  CALL poly_to_coordinates(this%trans%poly, v7d_out)

  nprev = 1
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(iy, ix, n, point) FIRSTPRIVATE(nprev)
  DO iy = 1, this%inny
    inside_x: DO ix = 1, this%innx
      point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))

      DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
        IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
          this%inter_index_x(ix,iy) = n
          nprev = n
          cycle inside_x
        ENDIF
      ENDDO
      DO n = nprev-1, 1, -1 ! test the other polygons
        IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
          this%inter_index_x(ix,iy) = n
          nprev = n
          cycle inside_x
        ENDIF
      ENDDO

!     CALL delete(point) ! speedup
    ENDDO inside_x
  ENDDO
!$OMP END PARALLEL

#ifdef DEBUG
  DO n = 1, this%trans%poly%arraysize
    CALL l4f_category_log(this%category, l4f_debug, &
     'Polygon: '//t2c(n)//' grid points: '// &
     t2c(count(this%inter_index_x(:,:) == n)))
  ENDDO
#endif

  CALL delete(lin)
  this%valid = .true. ! warning, no check of subtype

ELSE IF (this%trans%trans_type == 'stencilinter') THEN
  
! from inter:near
  CALL copy(in, lin)
  CALL get_val(lin, nx=this%innx, ny=this%inny)
  this%outnx = SIZE(v7d_out%ana)
  this%outny = 1

  ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
   this%inter_index_y(this%outnx,this%outny))
  ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))

  CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
! equalize in/out coordinates
  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
  CALL griddim_set_central_lon(lin, lonref)

  CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)

  CALL this%find_index(lin, .true., &
   this%innx, this%inny, xmin, xmax, ymin, ymax, &
   lon, lat, this%trans%extrap, &
   this%inter_index_x, this%inter_index_y)

! define the stencil mask
  nr = int(this%trans%area_info%radius) ! integer radius
  n = nr*2+1 ! n. of points
  nm = nr + 1 ! becomes index of center
  r2 = this%trans%area_info%radius**2
  ALLOCATE(this%stencil(n,n))
  this%stencil(:,:) = .true.
  DO iy = 1, n
    DO ix = 1, n
      IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
    ENDDO
  ENDDO

! set to missing the elements of inter_index too close to the domain
! borders
  xnmin = nr + 1
  xnmax = this%innx - nr
  ynmin = nr + 1
  ynmax = this%inny - nr
  DO iy = 1, this%outny
    DO ix = 1, this%outnx
      IF (this%inter_index_x(ix,iy) < xnmin .OR. &
       this%inter_index_x(ix,iy) > xnmax .OR. &
       this%inter_index_y(ix,iy) < ynmin .OR. &
       this%inter_index_y(ix,iy) > ynmax) THEN
        this%inter_index_x(ix,iy) = imiss
        this%inter_index_y(ix,iy) = imiss
      ENDIF
    ENDDO
  ENDDO

#ifdef DEBUG
  CALL l4f_category_log(this%category, l4f_debug, &
   'stencilinter: stencil size '//t2c(n*n))
  CALL l4f_category_log(this%category, l4f_debug, &
   'stencilinter: stencil points '//t2c(count(this%stencil)))
#endif

  DEALLOCATE(lon,lat)
  CALL delete(lin)

  this%valid = .true. ! warning, no check of subtype

ELSE IF (this%trans%trans_type == 'maskinter') THEN

  IF (.NOT.PRESENT(maskgrid)) THEN
    CALL l4f_category_log(this%category,l4f_error, &
     'grid_transform_init maskgrid argument missing for maskinter transformation')
    CALL raise_fatal_error()
  ENDIF

  CALL get_val(in, nx=this%innx, ny=this%inny)
  IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
    CALL l4f_category_log(this%category,l4f_error, &
     'grid_transform_init mask non conformal with input field')
    CALL l4f_category_log(this%category,l4f_error, &
     'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
     ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
    CALL raise_fatal_error()
  ENDIF
! unlike before, here index arrays must have the shape of input grid
  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
   this%inter_index_y(this%innx,this%inny))
  this%inter_index_x(:,:) = imiss
  this%inter_index_y(:,:) = 1

! generate the subarea boundaries according to maskgrid and maskbounds
  CALL gen_mask_class()

! compute coordinates of input grid in geo system
  CALL copy(in, lin)
  CALL unproj(lin)

!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(iy, ix, n)
  DO iy = 1, this%inny
    DO ix = 1, this%innx
      IF (c_e(maskgrid(ix,iy))) THEN
        IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
          DO n = nmaskarea, 1, -1
            IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
              this%inter_index_x(ix,iy) = n
              EXIT
            ENDIF
          ENDDO
        ENDIF
      ENDIF
    ENDDO
  ENDDO
!$OMP END PARALLEL

  this%outnx = nmaskarea
  this%outny = 1
  CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
  CALL init(v7d_out, time_definition=time_definition)
  CALL vol7d_alloc(v7d_out, nana=nmaskarea)

! setup output point list, equal to average of polygon points
! warning, in case of concave areas points may coincide!
  DO n = 1, nmaskarea
    CALL init(v7d_out%ana(n), &
     lon=stat_average(pack(lin%dim%lon(:,:), &
     mask=(this%inter_index_x(:,:) == n))), &
     lat=stat_average(pack(lin%dim%lat(:,:), &
     mask=(this%inter_index_x(:,:) == n))))
  ENDDO

  CALL delete(lin)
  this%valid = .true. ! warning, no check of subtype

ELSE IF (this%trans%trans_type == 'metamorphosis') THEN

! common to all metamorphosis subtypes
! compute coordinates of input grid in geo system
  CALL copy(in, lin)
  CALL unproj(lin)

  CALL get_val(in, nx=this%innx, ny=this%inny)
! allocate index array
  ALLOCATE(this%point_index(this%innx,this%inny))
  this%point_index(:,:) = imiss
! setup output coordinates
  IF (this%trans%sub_type /= 'index') THEN ! improve!
    CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
    CALL init(v7d_out, time_definition=time_definition)
  ENDIF

  IF (this%trans%sub_type == 'all' ) THEN

    this%outnx = this%innx*this%inny
    this%outny = 1
    CALL vol7d_alloc(v7d_out, nana=this%outnx)

    n = 0
    DO iy=1,this%inny
      DO ix=1,this%innx
        CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
         lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
        n = n + 1
        this%point_index(ix,iy) = n
      ENDDO
    ENDDO

    this%valid = .true.

  ELSE IF (this%trans%sub_type == 'coordbb' ) THEN

! count and mark points falling into requested bounding-box
    this%outnx = 0
    this%outny = 1
    DO iy = 1, this%inny
      DO ix = 1, this%innx
!        IF (geo_coord_inside_rectang()
        IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
         lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
         lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
         lin%dim%lat(ix,iy) < this%trans%rect_coo%flat) THEN ! improve!
          this%outnx = this%outnx + 1
          this%point_index(ix,iy) = this%outnx
        ENDIF
      ENDDO
    ENDDO

    IF (this%outnx <= 0) THEN
      CALL l4f_category_log(this%category,l4f_warn, &
       "metamorphosis:coordbb: no points inside bounding box "//&
       trim(to_char(this%trans%rect_coo%ilon))//","// &
       trim(to_char(this%trans%rect_coo%flon))//","// &
       trim(to_char(this%trans%rect_coo%ilat))//","// &
       trim(to_char(this%trans%rect_coo%flat)))
    ENDIF

    CALL vol7d_alloc(v7d_out, nana=this%outnx)

! collect coordinates of points falling into requested bounding-box
    n = 0
    DO iy = 1, this%inny
      DO ix = 1, this%innx
        IF (c_e(this%point_index(ix,iy))) THEN
          n = n + 1
          CALL init(v7d_out%ana(n), &
           lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
        ENDIF
      ENDDO
    ENDDO

    this%valid = .true.

  ELSE IF (this%trans%sub_type == 'poly' ) THEN

! count and mark points falling into requested polygon
#ifdef _OPENMP
    outnx = 0
#else
    this%outnx = 0
#endif
    this%outny = 1

! this OMP block has to be checked
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(iy, ix, point, n), REDUCTION(+:outnx)
    DO iy = 1, this%inny
      DO ix = 1, this%innx
        point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
        DO n = 1, this%trans%poly%arraysize
          IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
#ifdef _OPENMP
                outnx = outnx + 1
#else
                this%outnx = this%outnx + 1
#endif
            this%point_index(ix,iy) = n
            EXIT
          ENDIF
        ENDDO
!     CALL delete(point) ! speedup
      ENDDO
    ENDDO
!$OMP END PARALLEL
!$    this%outnx = outnx
    IF (this%outnx <= 0) THEN
      CALL l4f_category_log(this%category,l4f_warn, &
       "metamorphosis:poly: no points inside polygons")
    ENDIF

    CALL vol7d_alloc(v7d_out, nana=this%outnx)
! collect coordinates of points falling into requested polygon
    n = 0
    DO iy = 1, this%inny
      DO ix = 1, this%innx
        IF (c_e(this%point_index(ix,iy))) THEN
          n = n + 1
          CALL init(v7d_out%ana(n), &
           lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
        ENDIF
      ENDDO
    ENDDO

    this%valid = .true.

  ELSE IF (this%trans%sub_type == 'mask' ) THEN

    IF (.NOT.PRESENT(maskgrid)) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
      CALL raise_error()
      RETURN
    ENDIF

    IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init mask non conformal with input field')
      CALL l4f_category_log(this%category,l4f_error, &
       'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
       ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
      CALL raise_error()
      RETURN
    ENDIF

! generate the subarea boundaries according to maskgrid and maskbounds
    CALL gen_mask_class()

#ifdef _OPENMP
    outnx = 0
#else
    this%outnx = 0
#endif
    this%outny = 1

! this OMP block has to be checked
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(iy, ix), REDUCTION(+:outnx)
    DO iy = 1, this%inny
      DO ix = 1, this%innx
        IF (c_e(maskgrid(ix,iy))) THEN
          IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1)) THEN
            DO n = nmaskarea, 1, -1
              IF (maskgrid(ix,iy) > lmaskbounds(n)) THEN
#ifdef _OPENMP
                outnx = outnx + 1
#else
                this%outnx = this%outnx + 1
#endif
                this%point_index(ix,iy) = n
                EXIT
              ENDIF
            ENDDO
          ENDIF
        ENDIF
      ENDDO
    ENDDO
!$OMP END PARALLEL
!$    this%outnx = outnx

    IF (this%outnx <= 0) THEN
      CALL l4f_category_log(this%category,l4f_warn, &
       "grid_transform_init no points inside mask for metamorphosis:mask transformation")
    ENDIF
#ifdef DEBUG
    DO n = 1, nmaskarea
      CALL l4f_category_log(this%category,l4f_info, &
       "points in subarea "//t2c(n)//": "// &
       t2c(count(this%point_index(:,:) == n)))
    ENDDO
#endif
    CALL vol7d_alloc(v7d_out, nana=this%outnx)
! collect coordinates of points falling into requested polygon
    n = 0
    DO iy = 1, this%inny
      DO ix = 1, this%innx
        IF (c_e(this%point_index(ix,iy))) THEN
          n = n + 1
          CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
        ENDIF
      ENDDO
    ENDDO

    this%valid = .true.

  ELSE IF (this%trans%sub_type == 'index') THEN
    this%outnx = SIZE(v7d_out%ana)
    this%outny = 1

! collect points having requested index
! check that v7d_out has just 1 variable B01192 (index)
    IF (any(c_e(lin%dim%lon))) THEN ! here we set also coordinates, find a more efficient way to check
      DO n = 1, SIZE(v7d_out%ana)
        iy = (v7d_out%volanai(n,1,1) - 1)/this%innx + 1
        ix = v7d_out%volanai(n,1,1) - (iy-1)*this%innx
        this%point_index(ix,iy) = v7d_out%volanai(n,1,1)
        CALL init(v7d_out%ana(n), &
         lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy)) ! (re)set coordinates from truth
      ENDDO
    ELSE ! here we trust coordinates from input
      DO n = 1, SIZE(v7d_out%ana)
        iy = (v7d_out%volanai(n,1,1) - 1)/this%innx + 1
        ix = v7d_out%volanai(n,1,1) - (iy-1)*this%innx
        this%point_index(ix,iy) = v7d_out%volanai(n,1,1)
      ENDDO
    ENDIF
    this%valid = .true.
  ENDIF
  CALL delete(lin)
ENDIF

CONTAINS

SUBROUTINE gen_mask_class()
REAL :: startmaskclass, mmin, mmax
INTEGER :: i

IF (PRESENT(maskbounds)) THEN
  nmaskarea = SIZE(maskbounds) - 1
  IF (nmaskarea > 0) THEN
    lmaskbounds = maskbounds ! automatic f2003 allocation
    RETURN
  ELSE IF (nmaskarea == 0) THEN
    CALL l4f_category_log(this%category,l4f_warn, &
     'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
     //', it will be ignored')
    CALL l4f_category_log(this%category,l4f_warn, &
     'at least 2 values are required for maskbounds')
  ENDIF
ENDIF

mmin = minval(maskgrid, mask=c_e(maskgrid))
mmax = maxval(maskgrid, mask=c_e(maskgrid))

nmaskarea = int(mmax - mmin + 1.5)
startmaskclass = mmin - 0.5
! assign limits for each class
ALLOCATE(lmaskbounds(nmaskarea+1))
lmaskbounds(:) = (/(startmaskclass+REAL(i),i=0,nmaskarea)/)
#ifdef DEBUG
CALL l4f_category_log(this%category,l4f_debug, &
 'maskinter, '//t2c(nmaskarea)//' subareas defined, with boundaries:')
DO i = 1, SIZE(lmaskbounds)
  CALL l4f_category_log(this%category,l4f_debug, &
   'maskinter '//t2c(i)//' '//t2c(lmaskbounds(i)))
ENDDO
#endif

END SUBROUTINE gen_mask_class

END SUBROUTINE grid_transform_grid_vol7d_init


SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
TYPE(grid_transform),INTENT(out) :: this
TYPE(transform_def),INTENT(in) :: trans
TYPE(vol7d),INTENT(in) :: v7d_in
TYPE(griddim_def),INTENT(in) :: out
character(len=*),INTENT(in),OPTIONAL :: categoryappend

INTEGER :: nx, ny
DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:)
TYPE(griddim_def) :: lout


CALL grid_transform_init_common(this, trans, categoryappend)
#ifdef DEBUG
CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-vg6d")
#endif

IF (this%trans%trans_type == 'inter') THEN

  IF ( this%trans%sub_type == 'linear' ) THEN
    
    this%innx=SIZE(v7d_in%ana)
    this%inny=1
    ALLOCATE(lon(this%innx,1),lat(this%innx,1))
    ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
    CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))

    CALL copy (out, lout)
! equalize in/out coordinates
    lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
    CALL griddim_set_central_lon(lout, lonref)

    CALL get_val(lout, nx=nx, ny=ny)
    this%outnx=nx
    this%outny=ny
    ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))

    CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
    CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
    CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)

    DEALLOCATE(lon,lat)
    CALL delete(lout)

    this%valid = .true.

  ENDIF

ELSE IF (this%trans%trans_type == 'boxinter') THEN

  this%innx=SIZE(v7d_in%ana)
  this%inny=1
! index arrays must have the shape of input grid
  ALLOCATE(lon(this%innx,1),lat(this%innx,1))
  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
   this%inter_index_y(this%innx,this%inny))
! get coordinates of input grid in geo system
  CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))

  CALL copy (out, lout)
! equalize in/out coordinates
  lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
  CALL griddim_set_central_lon(lout, lonref)

  CALL get_val(lout, nx=this%outnx, ny=this%outny, &
   xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
! TODO now box size is ignored
! if box size not provided, use the actual grid step
  IF (.NOT.c_e(this%trans%area_info%boxdx)) &
   CALL get_val(out, dx=this%trans%area_info%boxdx)
  IF (.NOT.c_e(this%trans%area_info%boxdy)) &
   CALL get_val(out, dx=this%trans%area_info%boxdy)
! half size is actually needed
  this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
  this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0

! use find_index in the opposite way, here extrap does not make sense
  CALL this%find_index(lout, .true., &
   this%outnx, this%outny, xmin, xmax, ymin, ymax, &
   lon, lat, .false., &
   this%inter_index_x, this%inter_index_y)

  DEALLOCATE(lon,lat)
  CALL delete(lout)

  this%valid = .true. ! warning, no check of subtype

ENDIF

END SUBROUTINE grid_transform_vol7d_grid_init


SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
 maskbounds, categoryappend)
TYPE(grid_transform),INTENT(out) :: this
TYPE(transform_def),INTENT(in) :: trans
TYPE(vol7d),INTENT(in) :: v7d_in
TYPE(vol7d),INTENT(inout) :: v7d_out
REAL,INTENT(in),OPTIONAL :: maskbounds(:)
CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend

INTEGER :: i, n
DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
! temporary, improve!!!!
DOUBLE PRECISION :: lon1, lat1
TYPE(georef_coord) :: point


CALL grid_transform_init_common(this, trans, categoryappend)
#ifdef DEBUG
CALL l4f_category_log(this%category, l4f_debug, "grid_transform v7d-v7d")
#endif

IF (this%trans%trans_type == 'inter') THEN

  IF ( this%trans%sub_type == 'linear' ) THEN
    
    CALL vol7d_alloc_vol(v7d_out) ! be safe
    this%outnx=SIZE(v7d_out%ana)
    this%outny=1
  
    this%innx=SIZE(v7d_in%ana)
    this%inny=1
  
    ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
    ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))

    CALL getval(v7d_in%ana(:)%coord,lon=this%inter_xp(:,1),lat=this%inter_yp(:,1))
    CALL getval(v7d_out%ana(:)%coord,lon=this%inter_x(:,1),lat=this%inter_y(:,1))

    this%valid = .true.

  ENDIF

ELSE IF (this%trans%trans_type == 'polyinter') THEN

  this%innx=SIZE(v7d_in%ana)
  this%inny=1
! unlike before, here index arrays must have the shape of input grid
  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
   this%inter_index_y(this%innx,this%inny))
  this%inter_index_x(:,:) = imiss
  this%inter_index_y(:,:) = 1

  DO i = 1, SIZE(v7d_in%ana)
! temporary, improve!!!!
    CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
    point = georef_coord_new(x=lon1, y=lat1)

    DO n = 1, this%trans%poly%arraysize
      IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
        this%inter_index_x(i,1) = n
        EXIT
      ENDIF
    ENDDO
  ENDDO

  this%outnx=this%trans%poly%arraysize
  this%outny=1
  CALL vol7d_alloc(v7d_out, nana=this%outnx)

! setup output point list, equal to average of polygon points
! warning, in case of concave areas points may coincide!
  CALL poly_to_coordinates(this%trans%poly, v7d_out)

  this%valid = .true. ! warning, no check of subtype

ELSE IF (this%trans%trans_type == 'metamorphosis') THEN

! common to all metamorphosis subtypes
  this%innx = SIZE(v7d_in%ana)
  this%inny = 1
! allocate index array
  ALLOCATE(this%point_index(this%innx,this%inny))
  this%point_index(:,:) = imiss

  IF (this%trans%sub_type == 'all' ) THEN

    CALL metamorphosis_all_setup()

  ELSE IF (this%trans%sub_type == 'coordbb' ) THEN

    ALLOCATE(lon(this%innx),lat(this%innx))

! count and mark points falling into requested bounding-box
    this%outnx = 0
    this%outny = 1
    CALL getval(v7d_in%ana(:)%coord,lon=lon,lat=lat)
    DO i = 1, this%innx
!      IF (geo_coord_inside_rectang()
      IF (lon(i) > this%trans%rect_coo%ilon .AND. &
       lon(i) < this%trans%rect_coo%flon .AND. &
       lat(i) > this%trans%rect_coo%ilat .AND. &
       lat(i) < this%trans%rect_coo%flat) THEN ! improve!
        this%outnx = this%outnx + 1
        this%point_index(i,1) = this%outnx
      ENDIF
    ENDDO

    IF (this%outnx <= 0) THEN
      CALL l4f_category_log(this%category,l4f_warn, &
       "metamorphosis:coordbb: no points inside bounding box "//&
       trim(to_char(this%trans%rect_coo%ilon))//","// &
       trim(to_char(this%trans%rect_coo%flon))//","// &
       trim(to_char(this%trans%rect_coo%ilat))//","// &
       trim(to_char(this%trans%rect_coo%flat)))
    ENDIF

    CALL vol7d_alloc(v7d_out, nana=this%outnx)

! collect coordinates of points falling into requested bounding-box
    n = 0
    DO i = 1, this%innx
      IF (c_e(this%point_index(i,1))) THEN
        n = n + 1
        CALL init(v7d_out%ana(n),lon=lon(i),lat=lat(i))
      ENDIF
    ENDDO
    DEALLOCATE(lon, lat)

    this%valid = .true.

  ELSE IF (this%trans%sub_type == 'poly' ) THEN

! count and mark points falling into requested polygon
    this%outnx = 0
    this%outny = 1
    DO i = 1, this%innx
! temporary, improve!!!!
      CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
      point = georef_coord_new(x=lon1, y=lat1)
      DO n = 1, this%trans%poly%arraysize
        IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
          this%outnx = this%outnx + 1
          this%point_index(i,1) = n
          EXIT
        ENDIF
      ENDDO
!     CALL delete(point) ! speedup
    ENDDO

    IF (this%outnx <= 0) THEN
      CALL l4f_category_log(this%category,l4f_warn, &
       "metamorphosis:poly: no points inside polygons")
    ENDIF

    CALL vol7d_alloc(v7d_out, nana=this%outnx)

! collect coordinates of points falling into requested polygon
    n = 0
    DO i = 1, this%innx
      IF (c_e(this%point_index(i,1))) THEN
        n = n + 1
! temporary, improve!!!!
        CALL getval(v7d_in%ana(i)%coord,lon=lon1,lat=lat1)
        CALL init(v7d_out%ana(n),lon=lon1,lat=lat1)
      ENDIF
    ENDDO

    this%valid = .true.

  ELSE IF (this%trans%sub_type == 'setinvalidto' ) THEN

    IF (.NOT.PRESENT(maskbounds)) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init maskbounds missing for metamorphosis:'// &
       trim(this%trans%sub_type)//' transformation')
      RETURN
    ELSE IF (SIZE(maskbounds) < 1) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init maskbounds empty for metamorphosis:'// &
       trim(this%trans%sub_type)//' transformation')
      RETURN
    ELSE
      this%val1 = maskbounds(1)
#ifdef DEBUG
      CALL l4f_category_log(this%category, l4f_debug, &
       "grid_transform_init setting invalid data to "//t2c(this%val1))
#endif
    ENDIF

    CALL metamorphosis_all_setup()

  ELSE IF (this%trans%sub_type == 'settoinvalid' ) THEN

    IF (.NOT.PRESENT(maskbounds)) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init maskbounds missing for metamorphosis:'// &
       trim(this%trans%sub_type)//' transformation')
      CALL raise_error()
      RETURN
    ELSE IF (SIZE(maskbounds) < 2) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
       trim(this%trans%sub_type)//' transformation')
      CALL raise_error()
      RETURN
    ELSE
      this%val1 = maskbounds(1)
      this%val2 = maskbounds(SIZE(maskbounds))
#ifdef DEBUG
      CALL l4f_category_log(this%category, l4f_debug, &
       "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//','// &
       t2c(this%val2)//']')
#endif
    ENDIF

    CALL metamorphosis_all_setup()

  ENDIF
ENDIF

CONTAINS

! common code to metamorphosis transformations conserving the number
! of points
SUBROUTINE metamorphosis_all_setup()

this%outnx = SIZE(v7d_in%ana)
this%outny = 1
this%point_index(:,1) = (/(i,i=1,this%innx)/)
CALL vol7d_alloc(v7d_out, nana=SIZE(v7d_in%ana))
v7d_out%ana = v7d_in%ana

this%valid = .true.

END SUBROUTINE metamorphosis_all_setup

END SUBROUTINE grid_transform_vol7d_vol7d_init


! Private subroutine for performing operations common to all constructors
SUBROUTINE grid_transform_init_common(this, trans, categoryappend)
TYPE(grid_transform),INTENT(inout) :: this
TYPE(transform_def),INTENT(in) :: trans
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,"start init_grid_transform")
#endif

this%trans=trans

END SUBROUTINE grid_transform_init_common

! internal subroutine to correctly initialise the output coordinates
! with polygon centroid coordinates
SUBROUTINE poly_to_coordinates(poly, v7d_out)
TYPE(arrayof_georef_coord_array),intent(in) :: poly
TYPE(vol7d),INTENT(inout) :: v7d_out

INTEGER :: n, sz
DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)

DO n = 1, poly%arraysize
  CALL getval(poly%array(n), x=lon, y=lat)
  sz = min(SIZE(lon), SIZE(lat))
  IF (lon(1) == lon(sz) .AND. lat(1) == lat(sz)) THEN ! closed polygon
    sz = sz - 1
  ENDIF
  CALL init(v7d_out%ana(n), lon=stat_average(lon(1:sz)), lat=stat_average(lat(1:sz)))
ENDDO

END SUBROUTINE poly_to_coordinates

SUBROUTINE grid_transform_delete(this)
TYPE(grid_transform),INTENT(inout) :: this

CALL delete(this%trans)

this%innx=imiss
this%inny=imiss
this%outnx=imiss
this%outny=imiss
this%iniox=imiss
this%inioy=imiss
this%infox=imiss
this%infoy=imiss
this%outinx=imiss
this%outiny=imiss
this%outfnx=imiss
this%outfny=imiss

if (associated(this%inter_index_x)) deallocate (this%inter_index_x)
if (associated(this%inter_index_y)) deallocate (this%inter_index_y)
if (associated(this%inter_index_z)) deallocate (this%inter_index_z)
if (associated(this%point_index)) deallocate (this%point_index)

if (associated(this%inter_x)) deallocate (this%inter_x)
if (associated(this%inter_y)) deallocate (this%inter_y)

if (associated(this%inter_xp)) deallocate (this%inter_xp)
if (associated(this%inter_yp)) deallocate (this%inter_yp)
if (associated(this%inter_zp)) deallocate (this%inter_zp)
if (associated(this%vcoord_in)) deallocate (this%vcoord_in)
if (associated(this%vcoord_out)) deallocate (this%vcoord_out)
if (associated(this%point_mask)) deallocate (this%point_mask)
if (associated(this%stencil)) deallocate (this%stencil)
if (associated(this%output_level_auto)) deallocate (this%output_level_auto)
IF (ALLOCATED(this%coord_3d_in)) DEALLOCATE(this%coord_3d_in)
this%valid = .false.

! close the logger
call l4f_category_delete(this%category)

END SUBROUTINE grid_transform_delete


SUBROUTINE grid_transform_get_val(this, output_level_auto, point_mask, &
 point_index, output_point_index, levshift, levused)
TYPE(grid_transform),INTENT(in) :: this
TYPE(vol7d_level),POINTER,OPTIONAL :: output_level_auto(:)
LOGICAL,INTENT(out),ALLOCATABLE,OPTIONAL :: point_mask(:)
INTEGER,INTENT(out),ALLOCATABLE,OPTIONAL :: point_index(:)
INTEGER,INTENT(out),ALLOCATABLE,OPTIONAL :: output_point_index(:)
INTEGER,INTENT(out),OPTIONAL :: levshift
INTEGER,INTENT(out),OPTIONAL :: levused

INTEGER :: i

IF (PRESENT(output_level_auto)) output_level_auto => this%output_level_auto
IF (PRESENT(point_mask)) THEN
  IF (ASSOCIATED(this%point_index)) THEN
    point_mask = c_e(reshape(this%point_index, (/SIZE(this%point_index)/)))
  ENDIF
ENDIF
IF (PRESENT(point_index)) THEN
  IF (ASSOCIATED(this%point_index)) THEN
    point_index = reshape(this%point_index, (/SIZE(this%point_index)/))
  ENDIF
ENDIF
IF (PRESENT(output_point_index)) THEN
  IF (ASSOCIATED(this%point_index)) THEN
! metamorphosis, index is computed from input origin of output point
    output_point_index = pack(this%point_index(:,:), c_e(this%point_index))
  ELSE IF (this%trans%trans_type == 'polyinter' .OR. &
   this%trans%trans_type == 'maskinter') THEN
! other cases, index is order of output point
    output_point_index = (/(i,i=1,this%outnx)/)
  ENDIF
ENDIF
IF (PRESENT(levshift)) levshift = this%levshift
IF (PRESENT(levused)) levused = this%levused

END SUBROUTINE grid_transform_get_val


FUNCTION grid_transform_c_e(this)
TYPE(grid_transform),INTENT(in) :: this
LOGICAL :: grid_transform_c_e

grid_transform_c_e = this%valid

END FUNCTION grid_transform_c_e


RECURSIVE SUBROUTINE grid_transform_compute(this, field_in, field_out, var, &
 coord_3d_in, zlist)
TYPE(grid_transform),INTENT(in),TARGET :: this
REAL,INTENT(in) :: field_in(:,:,:)
REAL,INTENT(out) :: field_out(:,:,:)
TYPE(vol7d_var),INTENT(in),OPTIONAL :: var
REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:)
LOGICAL,INTENT(in),OPTIONAL :: zlist(:)

INTEGER :: i, j, k, l, m, s, ii, jj, ie, je, n, navg, kk, kkcache, kkup, kkdown, &
 kfound, kfoundin, inused, i1, i2, j1, j2, np, ns, ix, iy
INTEGER,ALLOCATABLE :: nval(:,:)
REAL :: z1,z2,z3,z4,z(4)
DOUBLE PRECISION  :: x1,x3,y1,y3,xp,yp, disttmp, dist
INTEGER :: innx, inny, innz, outnx, outny, outnz, vartype, nearcount
REAL,ALLOCATABLE :: coord_in(:)
LOGICAL,ALLOCATABLE :: mask_in(:)
REAL,ALLOCATABLE :: val_in(:), field_tmp(:,:,:)
REAL,POINTER :: coord_3d_in_act(:,:,:)
TYPE(grid_transform) :: likethis
LOGICAL :: alloc_coord_3d_in_act, nm1, optsearch, farenough
CHARACTER(len=4) :: env_var
LOGICAL,ALLOCATABLE :: lzlist(:)

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

field_out(:,:,:) = rmiss

IF (.NOT.this%valid) THEN
  CALL l4f_category_log(this%category,l4f_error, &
   "refusing to perform a non valid transformation")
  RETURN
ENDIF

IF (this%recur) THEN ! if recursive transformation, recur here and exit
#ifdef DEBUG
  CALL l4f_category_log(this%category,l4f_debug, &
   "recursive transformation in grid_transform_compute")
#endif

  IF (this%trans%trans_type == 'polyinter') THEN
    likethis = this
    likethis%recur = .false.
    likethis%outnx = this%trans%poly%arraysize
    likethis%outny = 1
    ALLOCATE(field_tmp(this%trans%poly%arraysize,1,SIZE(field_out,3)))
    CALL grid_transform_compute(likethis, field_in, field_tmp, var, coord_3d_in, zlist)

    DO k = 1, SIZE(field_out,3)
      DO j = 1, this%inny
        DO i = 1, this%innx
          IF (c_e(this%inter_index_x(i,j))) THEN
            field_out(i,j,k) = field_tmp(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
          ENDIF
        ENDDO
      ENDDO
    ENDDO
    DEALLOCATE(field_tmp)
  ENDIF

  RETURN
ENDIF

vartype = var_ord
IF (PRESENT(var)) THEN
  vartype = vol7d_vartype(var)
ENDIF

innx = SIZE(field_in,1); inny = SIZE(field_in,2); innz = SIZE(field_in,3)
outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)
ALLOCATE(lzlist(innz))
IF (PRESENT(zlist)) THEN
  lzlist(:) = zlist(:)
ELSE
  lzlist(:) = .true.
ENDIF

! check size of field_in, field_out
IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation

  IF (innz /= this%innz) THEN
    CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
     t2c(this%innz)//" /= "//t2c(innz))
    CALL raise_error()
    RETURN
  ENDIF

  IF (outnz /= this%outnz) THEN
    CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
     t2c(this%outnz)//" /= "//t2c(outnz))
    CALL raise_error()
    RETURN
  ENDIF

  IF (innx /= outnx .OR. inny /= outny) THEN
    CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//&
     t2c(innx)//","//t2c(inny)//" /= "//&
     t2c(outnx)//","//t2c(outny))
    CALL raise_error()
    RETURN
  ENDIF

ELSE ! horizontal interpolation

  IF (innx /= this%innx .OR. inny /= this%inny) THEN
    CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
     t2c(this%innx)//","//t2c(this%inny)//" /= "//&
     t2c(innx)//","//t2c(inny))
    CALL raise_error()
    RETURN
  ENDIF

  IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
    CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
     t2c(this%outnx)//","//t2c(this%outny)//" /= "//&
     t2c(outnx)//","//t2c(outny))
    CALL raise_error()
    RETURN
  ENDIF

  IF (innz /= outnz) THEN
    CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//&
     t2c(innz)//" /= "//t2c(outnz))
    CALL raise_error()
    RETURN
  ENDIF

ENDIF

#ifdef DEBUG
call l4f_category_log(this%category,l4f_debug, &
 "start grid_transform_compute "//trim(this%trans%trans_type)//':'// &
 trim(this%trans%sub_type))
#endif

IF (this%trans%trans_type == 'zoom') THEN

  field_out(this%outinx:this%outfnx, &
   this%outiny:this%outfny,:) = &
   field_in(this%iniox:this%infox, &
   this%inioy:this%infoy,:)

ELSE IF (this%trans%trans_type == 'boxregrid') THEN

  IF (this%trans%sub_type == 'average') THEN
    IF (vartype == var_dir360) THEN
      DO k = 1, innz
        IF (.NOT.lzlist(k)) cycle
        jj = 0
        DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
          je = j+this%trans%box_info%npy-1
          jj = jj+1
          ii = 0
          DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
            ie = i+this%trans%box_info%npx-1
            ii = ii+1
            navg = count(field_in(i:ie,j:je,k) /= rmiss)
            IF (navg > 0) THEN
              field_out(ii,jj,k) = find_prevailing_direction(field_in(i:ie,j:je,k), &
               0.0, 360.0, 5.0)
            ENDIF
          ENDDO
        ENDDO
      ENDDO

    ELSE
      DO k = 1, innz
        IF (.NOT.lzlist(k)) cycle
        jj = 0
        DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
          je = j+this%trans%box_info%npy-1
          jj = jj+1
          ii = 0
          DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
            ie = i+this%trans%box_info%npx-1
            ii = ii+1
            navg = count(field_in(i:ie,j:je,k) /= rmiss)
            IF (navg > 0) THEN
              field_out(ii,jj,k) = sum(field_in(i:ie,j:je,k), &
               mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
            ENDIF
          ENDDO
        ENDDO
      ENDDO

    ENDIF

  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
   this%trans%sub_type == 'stddevnm1') THEN

    IF (this%trans%sub_type == 'stddev') THEN
      nm1 = .false.
    ELSE
      nm1 = .true.
    ENDIF

    navg = this%trans%box_info%npx*this%trans%box_info%npy
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      jj = 0
      DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
        je = j+this%trans%box_info%npy-1
        jj = jj+1
        ii = 0
        DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
          ie = i+this%trans%box_info%npx-1
          ii = ii+1
          field_out(ii,jj,k) = stat_stddev( &
           reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
        ENDDO
      ENDDO
    ENDDO

  ELSE IF (this%trans%sub_type == 'max') THEN
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      jj = 0
      DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
        je = j+this%trans%box_info%npy-1
        jj = jj+1
        ii = 0
        DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
          ie = i+this%trans%box_info%npx-1
          ii = ii+1
          navg = count(field_in(i:ie,j:je,k) /= rmiss)
          IF (navg > 0) THEN
            field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
             mask=(field_in(i:ie,j:je,k) /= rmiss))
          ENDIF
        ENDDO
      ENDDO
    ENDDO

  ELSE IF (this%trans%sub_type == 'min') THEN
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      jj = 0
      DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
        je = j+this%trans%box_info%npy-1
        jj = jj+1
        ii = 0
        DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
          ie = i+this%trans%box_info%npx-1
          ii = ii+1
          navg = count(field_in(i:ie,j:je,k) /= rmiss)
          IF (navg > 0) THEN
            field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
             mask=(field_in(i:ie,j:je,k) /= rmiss))
          ENDIF
        ENDDO
      ENDDO
    ENDDO

  ELSE IF (this%trans%sub_type == 'percentile') THEN

    navg = this%trans%box_info%npx*this%trans%box_info%npy
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      jj = 0
      DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
        je = j+this%trans%box_info%npy-1
        jj = jj+1
        ii = 0
        DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
          ie = i+this%trans%box_info%npx-1
          ii = ii+1
          field_out(ii:ii,jj,k) = stat_percentile( &
           reshape(field_in(i:ie,j:je,k),(/navg/)), &
           (/REAL(this%trans%stat_info%percentile)/))
        ENDDO
      ENDDO
    ENDDO

  ELSE IF (this%trans%sub_type == 'frequency') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      jj = 0
      DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
        je = j+this%trans%box_info%npy-1
        jj = jj+1
        ii = 0
        DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
          ie = i+this%trans%box_info%npx-1
          ii = ii+1
          navg = count(field_in(i:ie,j:je,k) /= rmiss)
          IF (navg > 0) THEN
            field_out(ii,jj,k) = sum(interval_info_valid( &
             this%trans%interval_info, field_in(i:ie,j:je,k)), &
             mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
          ENDIF
        ENDDO
      ENDDO
    ENDDO

  ENDIF

ELSE IF (this%trans%trans_type == 'inter') THEN

  IF (this%trans%sub_type == 'near') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%outny
        DO i = 1, this%outnx

          IF (c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
           field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)

        ENDDO
      ENDDO
    ENDDO

  ELSE IF (this%trans%sub_type == 'bilin') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%outny
        DO i = 1, this%outnx

          IF (c_e(this%inter_index_x(i,j))) THEN

            z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
            z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
            z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
            z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)

            IF (c_e(z1) .AND. c_e(z2) .AND. c_e(z3) .AND. c_e(z4)) THEN

              x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
              y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
              x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
              y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)

              xp=this%inter_xp(i,j)
              yp=this%inter_yp(i,j)

              field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)

            ENDIF
          ENDIF

        ENDDO
      ENDDO
    ENDDO
  ELSE IF (this%trans%sub_type == 'shapiro_near') THEN
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%outny
        DO i = 1, this%outnx

          IF (c_e(this%inter_index_x(i,j))) THEN

            IF(this%inter_index_x(i,j)-1>0)THEN          
              z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j)  ,k)
            ELSE
              z(1) = rmiss
            END IF
            IF(this%inter_index_x(i,j)+1<=this%outnx)THEN
              z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)  ,k)
            ELSE
              z(3) = rmiss
            END IF
            IF(this%inter_index_y(i,j)+1<=this%outny)THEN          
              z(2)=field_in(this%inter_index_x(i,j)  ,this%inter_index_y(i,j)+1,k)
            ELSE
              z(2) = rmiss
            END IF
            IF(this%inter_index_y(i,j)-1>0)THEN
              z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
            ELSE
              z(4) = rmiss
            END IF
            field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))

          ENDIF

        ENDDO
      ENDDO
    ENDDO

  ENDIF
ELSE IF (this%trans%trans_type == 'intersearch') THEN

  likethis = this
  likethis%trans%trans_type = 'inter' ! fake type and make a recursive call to compute base field
  CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in, zlist)
  CALL getenv('LIBSIM_DISABLEOPTSEARCH', env_var)
  optsearch = len_trim(env_var) == 0

  DO k = 1, innz
    IF (.NOT.lzlist(k)) cycle
    IF ((.NOT.all(c_e(field_out(:,:,k)))) .AND. (any(c_e(field_in(:,:,k))))) THEN ! must fill some values
      DO j = 1, this%outny
        DO i = 1, this%outnx
          IF (.NOT.c_e(field_out(i,j,k))) THEN
            dist = huge(dist)
            nearcount = 0
            IF (optsearch) THEN ! optimized, error-prone algorithm
              ix = this%inter_index_x(i,j)
              iy = this%inter_index_y(i,j)
              DO s = 0, max(this%innx, this%inny)
                farenough = .true.
                DO m = iy-s, iy+s, max(2*s, 1) ! y loop on upper and lower frames
                  IF (m < 1 .OR. m > this%inny) cycle
                  DO l = max(1, ix-s), min(this%innx, ix+s) ! x loop on upper and lower frames
                    disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
                    IF (c_e(field_in(l,m,k))) THEN
                      IF (disttmp < dist) THEN
                        dist = disttmp
                        field_out(i,j,k) = field_in(l,m,k)
                        nearcount = 1
                      ELSE IF (disttmp == dist) THEN
                        field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
                        nearcount = nearcount + 1
                      ENDIF
                    ENDIF
                    IF (disttmp < dist) farenough = .false.
                  ENDDO
                ENDDO
                DO m = max(1, iy-s+1), min(this%inny, iy+s-1) ! y loop on left and right frames (avoid corners)
                  DO l = ix-s, ix+s, 2*s ! x loop on left and right frames (exchange loops?)
                    IF (l < 1 .OR. l > this%innx) cycle
                    disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
                    IF (c_e(field_in(l,m,k))) THEN
                      IF (disttmp < dist) THEN
                        dist = disttmp
                        field_out(i,j,k) = field_in(l,m,k)
                        nearcount = 1
                      ELSE IF (disttmp == dist) THEN
                        field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
                        nearcount = nearcount + 1
                      ENDIF
                    ENDIF
                    IF (disttmp < dist) farenough = .false.
                  ENDDO
                ENDDO
                IF (s > 0 .AND. farenough) EXIT ! nearest point found, do not trust the same point, in case of bilin it could be not the nearest
              ENDDO
            ELSE ! linear, simple, slow algorithm
              DO m = 1, this%inny
                DO l = 1, this%innx
                  IF (c_e(field_in(l,m,k))) THEN
                    disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
                    IF (disttmp < dist) THEN
                      dist = disttmp
                      field_out(i,j,k) = field_in(l,m,k)
                      nearcount = 1
                    ELSE IF (disttmp == dist) THEN
                      field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
                      nearcount = nearcount + 1
                    ENDIF
                  ENDIF
                ENDDO
              ENDDO
            ENDIF
! average points with same minimum distance
            IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
          ENDIF
        ENDDO
      ENDDO
    ENDIF
  ENDDO

ELSE IF (this%trans%trans_type == 'boxinter' &
 .OR. this%trans%trans_type == 'polyinter' &
 .OR. this%trans%trans_type == 'maskinter') THEN

  IF (this%trans%sub_type == 'average') THEN

    IF (vartype == var_dir360) THEN
      DO k = 1, innz
        IF (.NOT.lzlist(k)) cycle
        DO j = 1, this%outny
          DO i = 1, this%outnx
            field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
             0.0, 360.0, 5.0, &
             mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
          ENDDO
        ENDDO
      ENDDO

    ELSE
      ALLOCATE(nval(this%outnx, this%outny))
      field_out(:,:,:) = 0.0
      DO k = 1, innz
        IF (.NOT.lzlist(k)) cycle
        nval(:,:) = 0
        DO j = 1, this%inny
          DO i = 1, this%innx
            IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
              field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
               field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
               field_in(i,j,k)
              nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
               nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
            ENDIF
          ENDDO
        ENDDO
        WHERE (nval(:,:) /= 0)
          field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
        ELSEWHERE
          field_out(:,:,k) = rmiss
        END WHERE
      ENDDO
      DEALLOCATE(nval)
    ENDIF

  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
   this%trans%sub_type == 'stddevnm1') THEN

    IF (this%trans%sub_type == 'stddev') THEN
      nm1 = .false.
    ELSE
      nm1 = .true.
    ENDIF
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%outny
        DO i = 1, this%outnx
! da paura
          field_out(i:i,j,k) = stat_stddev( &
           reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
           mask=reshape((this%inter_index_x == i .AND. &
           this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)), nm1=nm1)
        ENDDO
      ENDDO
    ENDDO

  ELSE IF (this%trans%sub_type == 'max') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%inny
        DO i = 1, this%innx
          IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
            IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
              field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
               max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
               field_in(i,j,k))
            ELSE
              field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
               field_in(i,j,k)
            ENDIF
          ENDIF
        ENDDO
      ENDDO
    ENDDO


  ELSE IF (this%trans%sub_type == 'min') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%inny
        DO i = 1, this%innx
          IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
            IF (c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k))) THEN
              field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
               min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
               field_in(i,j,k))
            ELSE
              field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
               field_in(i,j,k)
            ENDIF
          ENDIF
        ENDDO
      ENDDO
    ENDDO

  ELSE IF (this%trans%sub_type == 'percentile') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%outny
        DO i = 1, this%outnx
! da paura
          field_out(i:i,j,k) = stat_percentile( &
           reshape(field_in(:,:,k), (/SIZE(field_in(:,:,k))/)), &
           (/REAL(this%trans%stat_info%percentile)/), &
           mask=reshape((this%inter_index_x == i .AND. &
           this%inter_index_y == j), (/SIZE(field_in(:,:,k))/)))
        ENDDO
      ENDDO
    ENDDO

  ELSE IF (this%trans%sub_type == 'frequency') THEN

    ALLOCATE(nval(this%outnx, this%outny))
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      field_out(:,:,k) = 0.0
      nval(:,:) = 0
      DO j = 1, this%inny
        DO i = 1, this%innx
          IF (c_e(this%inter_index_x(i,j)) .AND. c_e(field_in(i,j,k))) THEN
            field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
             field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
             interval_info_valid(this%trans%interval_info, field_in(i,j,k))
            nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
             nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
          ENDIF
        ENDDO
      ENDDO
      WHERE (nval(:,:) /= 0)
        field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
      ELSEWHERE
        field_out(:,:,k) = rmiss
      END WHERE
    ENDDO
    DEALLOCATE(nval)

  ENDIF

ELSE IF (this%trans%trans_type == 'stencilinter') THEN
  np = SIZE(this%stencil,1)/2
  ns = SIZE(this%stencil)

  IF (this%trans%sub_type == 'average') THEN

    IF (vartype == var_dir360) THEN
      DO k = 1, innz
        IF (.NOT.lzlist(k)) cycle
        DO j = 1, this%outny
          DO i = 1, this%outnx
            IF (c_e(this%inter_index_x(i,j))) THEN
              i1 = this%inter_index_x(i,j) - np
              i2 = this%inter_index_x(i,j) + np
              j1 = this%inter_index_y(i,j) - np
              j2 = this%inter_index_y(i,j) + np
              field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
               0.0, 360.0, 5.0, &
               mask=this%stencil(:,:)) ! simpler and equivalent form
!               mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
            ENDIF
          ENDDO
        ENDDO
      ENDDO

    ELSE
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
      DO k = 1, innz
        IF (.NOT.lzlist(k)) cycle
        DO j = 1, this%outny
          DO i = 1, this%outnx
            IF (c_e(this%inter_index_x(i,j))) THEN
              i1 = this%inter_index_x(i,j) - np
              i2 = this%inter_index_x(i,j) + np
              j1 = this%inter_index_y(i,j) - np
              j2 = this%inter_index_y(i,j) + np
              n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
              IF (n > 0) THEN
                field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
                 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
              ENDIF
            ENDIF
          ENDDO
        ENDDO
      ENDDO
!$OMP END PARALLEL
    ENDIF

  ELSE IF (this%trans%sub_type == 'stddev' .OR. &
   this%trans%sub_type == 'stddevnm1') THEN

    IF (this%trans%sub_type == 'stddev') THEN
      nm1 = .false.
    ELSE
      nm1 = .true.
    ENDIF

!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%outny
        DO i = 1, this%outnx
          IF (c_e(this%inter_index_x(i,j))) THEN
            i1 = this%inter_index_x(i,j) - np
            i2 = this%inter_index_x(i,j) + np
            j1 = this%inter_index_y(i,j) - np
            j2 = this%inter_index_y(i,j) + np
! da paura
            field_out(i:i,j,k) = stat_stddev( &
             reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
             mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
             this%stencil(:,:), (/ns/)), nm1=nm1)
          ENDIF
        ENDDO
      ENDDO
    ENDDO
!$OMP END PARALLEL

  ELSE IF (this%trans%sub_type == 'max') THEN

!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%outny
        DO i = 1, this%outnx
          IF (c_e(this%inter_index_x(i,j))) THEN
            i1 = this%inter_index_x(i,j) - np
            i2 = this%inter_index_x(i,j) + np
            j1 = this%inter_index_y(i,j) - np
            j2 = this%inter_index_y(i,j) + np
            n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
            IF (n > 0) THEN
              field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
               mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
            ENDIF
          ENDIF
        ENDDO
      ENDDO
    ENDDO
!$OMP END PARALLEL

  ELSE IF (this%trans%sub_type == 'min') THEN

!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%outny
        DO i = 1, this%outnx
          IF (c_e(this%inter_index_x(i,j))) THEN
            i1 = this%inter_index_x(i,j) - np
            i2 = this%inter_index_x(i,j) + np
            j1 = this%inter_index_y(i,j) - np
            j2 = this%inter_index_y(i,j) + np
            n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
            IF (n > 0) THEN
              field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
               mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
            ENDIF
          ENDIF
        ENDDO
      ENDDO
    ENDDO
!$OMP END PARALLEL

  ELSE IF (this%trans%sub_type == 'percentile') THEN

!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2)
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%outny
        DO i = 1, this%outnx
          IF (c_e(this%inter_index_x(i,j))) THEN
            i1 = this%inter_index_x(i,j) - np
            i2 = this%inter_index_x(i,j) + np
            j1 = this%inter_index_y(i,j) - np
            j2 = this%inter_index_y(i,j) + np
! da paura
            field_out(i:i,j,k) = stat_percentile( &
             reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
             (/REAL(this%trans%stat_info%percentile)/), &
             mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
             this%stencil(:,:), (/ns/)))
          ENDIF
        ENDDO
      ENDDO
    ENDDO
!$OMP END PARALLEL

  ELSE IF (this%trans%sub_type == 'frequency') THEN

!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(i, j, k, i1, i2, j1, j2, n)
    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      DO j = 1, this%outny
        DO i = 1, this%outnx
          IF (c_e(this%inter_index_x(i,j))) THEN
            i1 = this%inter_index_x(i,j) - np
            i2 = this%inter_index_x(i,j) + np
            j1 = this%inter_index_y(i,j) - np
            j2 = this%inter_index_y(i,j) + np
            n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
            IF (n > 0) THEN
              field_out(i,j,k) = sum(interval_info_valid( &
               this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
               mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
            ENDIF
          ENDIF
        ENDDO
      ENDDO
    ENDDO
!$OMP END PARALLEL
    
  ENDIF

ELSE IF (this%trans%trans_type == 'maskgen') THEN

  DO k = 1, innz
    WHERE(c_e(this%inter_index_x(:,:)))
      field_out(:,:,k) = REAL(this%inter_index_x(:,:))
    END WHERE
  ENDDO

ELSE IF (this%trans%trans_type == 'metamorphosis') THEN

  IF (this%trans%sub_type == 'all') THEN

    field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))

  ELSE IF (this%trans%sub_type == 'coordbb' .OR. this%trans%sub_type == 'poly' &
   .OR. this%trans%sub_type == 'mask' .OR. this%trans%sub_type == 'index') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
! this is to sparse-points only, so field_out(:,1,k) is acceptable
      field_out(:,1,k) = pack(field_in(:,:,k), c_e(this%point_index(:,:)))
    ENDDO

  ELSE IF (this%trans%sub_type == 'maskvalid' .OR. &
   this%trans%sub_type == 'maskinvalid') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      WHERE (this%point_mask(:,:)) ! this is almost the same as the previous syntax
        field_out(:,:,k) = field_in(:,:,k) ! check which is more efficient
      END WHERE
    ENDDO

  ELSE IF (this%trans%sub_type == 'lemaskinvalid') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
        field_out(:,:,k) = field_in(:,:,k)
      ELSEWHERE
        field_out(:,:,k) = rmiss
      END WHERE
    ENDDO

  ELSE IF (this%trans%sub_type == 'ltmaskinvalid') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
        field_out(:,:,k) = field_in(:,:,k)
      ELSEWHERE
        field_out(:,:,k) = rmiss
      END WHERE
    ENDDO

  ELSE IF (this%trans%sub_type == 'gemaskinvalid') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
        field_out(:,:,k) = field_in(:,:,k)
      ELSEWHERE
        field_out(:,:,k) = rmiss
      END WHERE
    ENDDO

  ELSE IF (this%trans%sub_type == 'gtmaskinvalid') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle
      WHERE (c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
        field_out(:,:,k) = field_in(:,:,k)
      ELSEWHERE
        field_out(:,:,k) = rmiss
      END WHERE
    ENDDO

  ELSE IF (this%trans%sub_type == 'setinvalidto') THEN

    DO k = 1, innz
      IF (.NOT.lzlist(k)) cycle ! should i skip levels here?
      WHERE (c_e(field_in(:,:,k)))
        field_out(:,:,k) = field_in(:,:,k)
      ELSE WHERE
        field_out(:,:,k) = this%val1
      END WHERE
    ENDDO

  ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
    IF (c_e(this%val1) .AND. c_e(this%val2)) THEN
      WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
       .AND. field_in(:,:,:) <= this%val2)
        field_out(:,:,:) = rmiss
      ELSE WHERE
        field_out(:,:,:) = field_in(:,:,:)
      END WHERE
    ELSE IF (c_e(this%val1)) THEN
      WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
        field_out(:,:,:) = rmiss
      ELSE WHERE
        field_out(:,:,:) = field_in(:,:,:)
      END WHERE
    ELSE IF (c_e(this%val2)) THEN
      WHERE (c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
        field_out(:,:,:) = rmiss
      ELSE WHERE
        field_out(:,:,:) = field_in(:,:,:)
      END WHERE
    ENDIF
  ENDIF

ELSE IF (this%trans%trans_type == 'vertint') THEN

  IF (this%trans%sub_type == 'linear') THEN

    alloc_coord_3d_in_act = .false.
    IF (ASSOCIATED(this%inter_index_z)) THEN

      DO k = 1, outnz
        IF (c_e(this%inter_index_z(k))) THEN
          z1 = REAL(this%inter_zp(k)) ! weight for k+1
          z2 = REAL(1.0D0 - this%inter_zp(k)) ! weight for k
          SELECT CASE(vartype)

          CASE(var_dir360)
            DO j = 1, inny
              DO i = 1, innx
                IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
                 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
                  field_out(i,j,k) = &
                   interp_var_360(field_in(i,j,this%inter_index_z(k)), &
                   field_in(i,j,this%inter_index_z(k)+1), z1, z2)
                ENDIF
              ENDDO
            ENDDO

          CASE(var_press)
            DO j = 1, inny
              DO i = 1, innx
                IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
                 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
                 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
                 field_in(i,j,this%inter_index_z(k)+1) > 0.) THEN
                  field_out(i,j,k) = exp( &
                   log(field_in(i,j,this%inter_index_z(k)))*z2 + &
                   log(field_in(i,j,this%inter_index_z(k)+1))*z1)
                ENDIF
              ENDDO
            ENDDO

          CASE default
            DO j = 1, inny
              DO i = 1, innx
                IF (c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
                 c_e(field_in(i,j,this%inter_index_z(k)+1))) THEN
                  field_out(i,j,k) = &
                   field_in(i,j,this%inter_index_z(k))*z2 + &
                   field_in(i,j,this%inter_index_z(k)+1)*z1
                ENDIF
              ENDDO
            ENDDO

          END SELECT

        ENDIF
      ENDDO

    ELSE ! use coord_3d_in

      IF (ALLOCATED(this%coord_3d_in)) THEN
        coord_3d_in_act => this%coord_3d_in
#ifdef DEBUG
        CALL l4f_category_log(this%category,l4f_debug, &
         "Using external vertical coordinate for vertint")
#endif
      ELSE
        IF (PRESENT(coord_3d_in)) THEN
#ifdef DEBUG
          CALL l4f_category_log(this%category,l4f_debug, &
           "Using internal vertical coordinate for vertint")
#endif
          IF (this%dolog) THEN
            ALLOCATE(coord_3d_in_act(SIZE(coord_3d_in,1), &
             SIZE(coord_3d_in,2), SIZE(coord_3d_in,3)))
            alloc_coord_3d_in_act = .true.
            WHERE (c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
              coord_3d_in_act = log(coord_3d_in)
            ELSEWHERE
              coord_3d_in_act = rmiss
            END WHERE
          ELSE
            coord_3d_in_act => coord_3d_in
          ENDIF
        ELSE
          CALL l4f_category_log(this%category,l4f_error, &
           "Missing internal and external vertical coordinate for vertint")
          CALL raise_error()
          RETURN
        ENDIF
      ENDIF

      inused = SIZE(coord_3d_in_act,3)
      IF (inused < 2) RETURN ! to avoid algorithm failure
      kkcache = 1

      DO k = 1, outnz
        IF (.NOT.c_e(this%vcoord_out(k))) cycle
        DO j = 1, inny
          DO i = 1, innx
            kfound = imiss
            DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
              kkup = kkcache + kk
              kkdown = kkcache - kk + 1

              IF (kkdown >= 1) THEN ! search down
                IF (this%vcoord_out(k) >= &
                 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
                 this%vcoord_out(k) < &
                 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1))) THEN
                  kkcache = kkdown
                  kfoundin = kkcache
                  kfound = kkcache + this%levshift
                  EXIT ! kk
                ENDIF
              ENDIF
              IF (kkup < inused) THEN ! search up
                IF (this%vcoord_out(k) >= &
                 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
                 this%vcoord_out(k) < &
                 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1))) THEN
                  kkcache = kkup
                  kfoundin = kkcache
                  kfound = kkcache + this%levshift
                  EXIT ! kk
                ENDIF
              ENDIF

            ENDDO

            IF (c_e(kfound)) THEN
              IF (c_e(field_in(i,j,kfound)) .AND. c_e(field_in(i,j,kfound+1))) THEN
                z1 = REAL((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
                 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
                z2 = 1.0 - z1
                SELECT CASE(vartype)

                CASE(var_dir360)
                  field_out(i,j,k) = &
                   interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
                CASE(var_press)
                  field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
                   log(field_in(i,j,kfound+1))*z1)

                CASE default
                  field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
                END SELECT

              ENDIF
            ELSE
            ENDIF
          ENDDO ! i
        ENDDO ! j
      ENDDO ! k
      IF (alloc_coord_3d_in_act) DEALLOCATE(coord_3d_in_act)
    ENDIF

  ELSE IF (this%trans%sub_type == 'linearsparse') THEN


    IF (.NOT.ASSOCIATED(this%vcoord_in) .AND. .NOT.PRESENT(coord_3d_in)) THEN
      CALL l4f_category_log(this%category,l4f_error, &
       "linearsparse interpolation, no input vert coord available")
      RETURN
    ENDIF

    ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
    DO j = 1, inny
      DO i = 1, innx

        IF (ASSOCIATED(this%vcoord_in)) THEN
          mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
           .AND. c_e(this%vcoord_in(:))
        ELSE
          mask_in = c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
           .AND. c_e(coord_3d_in(i,j,:))
        ENDIF

        IF (vartype == var_press) THEN
          mask_in(:) = mask_in(:) .AND. &
           (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
        ENDIF
        inused = count(mask_in)
        IF (inused > 1) THEN
          IF (ASSOCIATED(this%vcoord_in)) THEN
            coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
          ELSE
            coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
          ENDIF
          IF (vartype == var_press) THEN
            val_in(1:inused) = log(pack( &
             field_in(i,j,this%levshift+1:this%levshift+this%levused), &
             mask=mask_in))
          ELSE
            val_in(1:inused) = pack( &
             field_in(i,j,this%levshift+1:this%levshift+this%levused), &
             mask=mask_in)
          ENDIF
          kkcache = 1
          DO k = 1, outnz

            kfound = imiss
            DO kk = 1, max(inused-kkcache-1, kkcache) ! +-1
              kkup = kkcache + kk
              kkdown = kkcache - kk + 1

              IF (kkdown >= 1) THEN ! search down
                IF (this%vcoord_out(k) >= &
                 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
                 this%vcoord_out(k) < &
                 max(coord_in(kkdown), coord_in(kkdown+1))) THEN
                  kkcache = kkdown
                  kfoundin = kkcache
                  kfound = kkcache
                  EXIT ! kk
                ENDIF
              ENDIF
              IF (kkup < inused) THEN ! search up
                IF (this%vcoord_out(k) >= &
                 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
                 this%vcoord_out(k) < &
                 max(coord_in(kkup), coord_in(kkup+1))) THEN
                  kkcache = kkup
                  kfoundin = kkcache
                  kfound = kkcache
                  EXIT ! kk
                ENDIF
              ENDIF

            ENDDO

            IF (c_e(kfound)) THEN
              z1 = REAL((this%vcoord_out(k) - coord_in(kfound-1))/ &
               (coord_in(kfound) - coord_in(kfound-1)))
              z2 = 1.0 - z1
              IF (vartype == var_dir360) THEN
                field_out(i,j,k) = &
                 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
              ELSE IF (vartype == var_press) THEN
                field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
              ELSE
                field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
              ENDIF
            ENDIF

          ENDDO

        ENDIF

      ENDDO
    ENDDO
    DEALLOCATE(coord_in,val_in)


  ENDIF

ELSE IF (this%trans%trans_type == '' .OR. this%trans%trans_type == 'none') THEN

  field_out(:,:,:) = field_in(:,:,:)

ENDIF


CONTAINS


! internal function for interpolating directions from 0 to 360 degree
! hope it is inlined by the compiler
FUNCTION interp_var_360(v1, v2, w1, w2)
REAL,INTENT(in) :: v1, v2, w1, w2
REAL :: interp_var_360

REAL :: lv1, lv2

IF (abs(v1 - v2) > 180.) THEN
  IF (v1 > v2) THEN
    lv1 = v1 - 360.
    lv2 = v2
  ELSE
    lv1 = v1
    lv2 = v2 - 360.
  ENDIF
  interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
ELSE
  interp_var_360 = v1*w2 + v2*w1
ENDIF

END FUNCTION interp_var_360


RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask) RESULT(prev)
REAL,INTENT(in) :: v1(:,:)
REAL,INTENT(in) :: l, h, res
LOGICAL,INTENT(in),OPTIONAL :: mask(:,:)
REAL :: prev

REAL :: m
INTEGER :: nh, nl
!REAL,PARAMETER :: res = 1.0

m = (l + h)/2.
IF ((h - l) <= res) THEN
  prev = m
  RETURN
ENDIF

IF (PRESENT(mask)) THEN
  nl = count(v1 >= l .AND. v1 < m .AND. mask)
  nh = count(v1 >= m .AND. v1 < h .AND. mask)
ELSE
  nl = count(v1 >= l .AND. v1 < m)
  nh = count(v1 >= m .AND. v1 < h)
ENDIF
IF (nh > nl) THEN
  prev = find_prevailing_direction(v1, m, h, res)
ELSE IF (nl > nh) THEN
  prev = find_prevailing_direction(v1, l, m, res)
ELSE IF (nl == 0 .AND. nh == 0) THEN
  prev = rmiss
ELSE
  prev = m
ENDIF

END FUNCTION find_prevailing_direction


END SUBROUTINE grid_transform_compute


SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
 coord_3d_in)
TYPE(grid_transform),INTENT(in) :: this
REAL, INTENT(in) :: field_in(:,:)
REAL, INTENT(out):: field_out(:,:,:)
TYPE(vol7d_var),INTENT(in),OPTIONAL :: var
REAL,INTENT(in),OPTIONAL,TARGET :: coord_3d_in(:,:,:)

real,allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
INTEGER :: inn_p, ier, k
INTEGER :: innx, inny, innz, outnx, outny, outnz

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

field_out(:,:,:) = rmiss

IF (.NOT.this%valid) THEN
  call l4f_category_log(this%category,l4f_error, &
   "refusing to perform a non valid transformation")
  call raise_error()
  RETURN
ENDIF

innx = SIZE(field_in,1); inny = 1; innz = SIZE(field_in,2)
outnx = SIZE(field_out,1); outny = SIZE(field_out,2); outnz = SIZE(field_out,3)

! check size of field_in, field_out
IF (this%trans%trans_type == 'vertint') THEN ! vertical interpolation

  IF (innz /= this%innz) THEN
    CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
     t2c(this%innz)//" /= "//t2c(innz))
    CALL raise_error()
    RETURN
  ENDIF

  IF (outnz /= this%outnz) THEN
    CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
     t2c(this%outnz)//" /= "//t2c(outnz))
    CALL raise_error()
    RETURN
  ENDIF

  IF (innx /= outnx .OR. inny /= outny) THEN
    CALL l4f_category_log(this%category,l4f_error,"vertical interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent hor. sizes: "//&
     t2c(innx)//","//t2c(inny)//" /= "//&
     t2c(outnx)//","//t2c(outny))
    CALL raise_error()
    RETURN
  ENDIF

ELSE ! horizontal interpolation

  IF (innx /= this%innx .OR. inny /= this%inny) THEN
    CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent input shape: "//&
     t2c(this%innx)//","//t2c(this%inny)//" /= "//&
     t2c(innx)//","//t2c(inny))
    CALL raise_error()
    RETURN
  ENDIF

  IF (outnx /= this%outnx .OR. outny /= this%outny) THEN
    CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent output shape: "//&
     t2c(this%outnx)//","//t2c(this%outny)//" /= "//&
     t2c(outnx)//","//t2c(outny))
    CALL raise_error()
    RETURN
  ENDIF

  IF (innz /= outnz) THEN
    CALL l4f_category_log(this%category,l4f_error,"horizontal interpolation")
    CALL l4f_category_log(this%category,l4f_error,"inconsistent vert. sizes: "//&
     t2c(innz)//" /= "//t2c(outnz))
    CALL raise_error()
    RETURN
  ENDIF

ENDIF

#ifdef DEBUG
  call l4f_category_log(this%category,l4f_debug, &
   "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//':'// &
   trim(this%trans%sub_type))
#endif

IF (this%trans%trans_type == 'inter') THEN

  IF (this%trans%sub_type == 'linear') THEN

#ifdef HAVE_LIBNGMATH
! optimization, allocate only once with a safe size
    ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
    DO k = 1, innz
      inn_p = count(c_e(field_in(:,k)))

      CALL l4f_category_log(this%category,l4f_info, &
       "Number of sparse data points: "//t2c(inn_p)//','//t2c(SIZE(field_in(:,k))))

      IF (inn_p > 2) THEN

        field_in_p(1:inn_p) = pack(field_in(:,k), c_e(field_in(:,k)))
        x_in_p(1:inn_p) = pack(this%inter_xp(:,1), c_e(field_in(:,k)))
        y_in_p(1:inn_p) = pack(this%inter_yp(:,1), c_e(field_in(:,k)))

        IF (.NOT.this%trans%extrap) THEN
          CALL nnseti('ext', 0) ! 0 = no extrapolation
          CALL nnsetr('nul', rmiss)
        ENDIF

        CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, & ! (1:inn_p) omitted
         this%outnx, this%outny, REAL(this%inter_x(:,1)), & ! no f90 interface
         REAL(this%inter_y(1,:)), field_out(1,1,k), ier)

        IF (ier /= 0) THEN
          CALL l4f_category_log(this%category,l4f_error, &
           "Error code from NCAR natgrids: "//t2c(ier))
          CALL raise_error()
          EXIT
        ENDIF ! exit loop to deallocate
      ELSE

        CALL l4f_category_log(this%category,l4f_info, &
         "insufficient data in gridded region to triangulate")

      ENDIF
    ENDDO
    DEALLOCATE(field_in_p, x_in_p, y_in_p)

#else
    CALL l4f_category_log(this%category,l4f_error, &
     "libsim compiled without NATGRIDD (ngmath ncarg library)")
    CALL raise_error()
    RETURN
#endif

  ENDIF

ELSE IF (this%trans%trans_type == 'boxinter' .OR. &
 this%trans%trans_type == 'polyinter' .OR. &
 this%trans%trans_type == 'vertint' .OR. &
 this%trans%trans_type == 'metamorphosis') THEN ! use the grid-to-grid method

  CALL compute(this, &
   reshape(field_in, (/SIZE(field_in,1), 1, SIZE(field_in,2)/)), field_out, var, &
   coord_3d_in)

ENDIF

END SUBROUTINE grid_transform_v7d_grid_compute


! Bilinear interpolation
! Effettua interpolazione bilineare dati i valori nei punti 1,2,3,4 e
! le coordinate dei punti 1 e 3 oltre a quelle del punto p dove viene
! valutato il campo.
!_____________________________________________________________
!                               disposizione punti
!       4       3
!
!         p
!
!       1       2
! _____________________________________________________________
ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp) RESULT(zp)
REAL,INTENT(in) :: z1,z2,z3,z4 ! Z values on the four points
DOUBLE PRECISION,INTENT(in):: x1,y1 ! coordinate of the lower left point
DOUBLE PRECISION,INTENT(in):: x3,y3 ! coordinate of the upper right point
DOUBLE PRECISION,INTENT(in):: xp,yp ! coordinate of point where interpolate
REAL :: zp

REAL :: p1,p2
REAL :: z5,z6


p2 = REAL((yp-y1)/(y3-y1))
p1 = REAL((xp-x1)/(x3-x1))

z5 = (z4-z1)*p2+z1
z6 = (z3-z2)*p2+z2

zp = (z6-z5)*(p1)+z5
      
END FUNCTION hbilin


! Shapiro filter of order 2
FUNCTION shapiro (z,zp) RESULT(zs)
!! z_smoothed(i,j) = z(i,j) * (1-S) + S * sum(z_vicini)/N 
!!                 = z(i,j) - S/N (N*z(i,j) - sum(z_vicini))
REAL,INTENT(in) :: z(:) ! Z values on the four surrounding points
REAL,INTENT(in) :: zp          ! Z values on the central point
REAL   :: zs                   ! Z smoothed value on the central point
INTEGER:: N

IF(c_e(zp))THEN
  n = count(c_e(z))
  zs = zp - 0.5* ( n*zp - sum(z, c_e(z)) )/n
ELSE
  zs = rmiss
END IF

END FUNCTION shapiro


! Locate index of requested point
SUBROUTINE basic_find_index(this, near, nx, ny, xmin, xmax, ymin, ymax, &
 lon, lat, extrap, index_x, index_y)
TYPE(griddim_def),INTENT(in) :: this ! griddim object (from grid)
LOGICAL,INTENT(in) :: near ! near or bilin interpolation (determine wich point is requested)
INTEGER,INTENT(in) :: nx,ny ! dimension (to grid)
DOUBLE PRECISION,INTENT(in) :: xmin, xmax, ymin, ymax ! extreme coordinate (to grid)
DOUBLE PRECISION,INTENT(in) :: lon(:,:),lat(:,:) ! target coordinate
LOGICAL,INTENT(in) :: extrap ! extrapolate
INTEGER,INTENT(out) :: index_x(:,:),index_y(:,:) ! index of point requested

INTEGER :: lnx, lny
DOUBLE PRECISION :: x(SIZE(lon,1),SIZE(lon,2)),y(SIZE(lon,1),SIZE(lon,2))

IF (near) THEN
  CALL proj(this,lon,lat,x,y)
  index_x = nint((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
  index_y = nint((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
  lnx = nx
  lny = ny
ELSE
  CALL proj(this,lon,lat,x,y)
  index_x = floor((x-xmin)/((xmax-xmin)/dble(nx-1)))+1
  index_y = floor((y-ymin)/((ymax-ymin)/dble(ny-1)))+1
  lnx = nx-1
  lny = ny-1
ENDIF

IF (extrap) THEN ! trim indices outside grid for extrapolation
  index_x = max(index_x, 1)
  index_y = max(index_y, 1)
  index_x = min(index_x, lnx)
  index_y = min(index_y, lny)
ELSE ! nullify indices outside grid
  WHERE(index_x < 1 .OR. index_x > lnx .OR. index_y < 1 .OR. index_y > lny)
    index_x = imiss
    index_y = imiss
  END WHERE
ENDIF

END SUBROUTINE basic_find_index

END MODULE grid_transform_class