grid_id_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_id_class
#ifdef HAVE_LIBGRIBAPI
USE grib_api
#endif
#ifdef HAVE_LIBGDAL
USE gdal
#endif
USE missing_values
USE optional_values
USE char_utilities
USE file_utilities
USE log4fortran
USE err_handling
IMPLICIT NONE

INTEGER,PARAMETER :: grid_id_no_driver = 0
INTEGER,PARAMETER :: grid_id_grib_api = 1
INTEGER,PARAMETER :: grid_id_gdal = 2

#if defined HAVE_LIBGRIBAPI
INTEGER,PARAMETER :: grid_id_default = grid_id_grib_api 
#elif defined HAVE_LIBGDAL
INTEGER,PARAMETER :: grid_id_default = grid_id_gdal 
#else
INTEGER,PARAMETER :: grid_id_default = grid_id_no_driver 
#endif

CHARACTER(len=12),PARAMETER :: driverlist(0:2) = &
 (/'no_driver   ','grib_api    ','gdal        '/)

#ifdef HAVE_LIBGDAL

TYPE gdal_file_id_options
  DOUBLE PRECISION :: xmin=dmiss 
  DOUBLE PRECISION :: ymin=dmiss 
  DOUBLE PRECISION :: xmax=dmiss 
  DOUBLE PRECISION :: ymax=dmiss 
END TYPE gdal_file_id_options
#endif

TYPE grid_file_id
PRIVATE
#ifdef HAVE_LIBGRIBAPI
INTEGER :: gaid=imiss
#endif
#ifdef HAVE_LIBGDAL
TYPE(gdaldataseth) :: gdalid
INTEGER :: nlastband=0
TYPE(gdal_file_id_options) :: gdal_options
TYPE(grid_file_id),POINTER :: file_id_copy=>null()
#endif
INTEGER :: driver=grid_id_default
END TYPE grid_file_id


TYPE grid_id
PRIVATE
INTEGER :: nodriverid=imiss
#ifdef HAVE_LIBGRIBAPI
INTEGER :: gaid=imiss
#endif
#ifdef HAVE_LIBGDAL
TYPE(gdalrasterbandh) :: gdalid
TYPE(grid_file_id),POINTER :: file_id=>null()
#endif
INTEGER :: driver=grid_id_default
END TYPE grid_id

INTERFACE init
  MODULE PROCEDURE grid_file_id_init, grid_id_init
END INTERFACE

INTERFACE delete
  MODULE PROCEDURE grid_file_id_delete, grid_id_delete
END INTERFACE

INTERFACE copy
  MODULE PROCEDURE grid_id_copy
END INTERFACE

INTERFACE export
  MODULE PROCEDURE grid_id_export
END INTERFACE

!INTERFACE ASSIGNMENT(=)
!  MODULE PROCEDURE &
!#ifdef HAVE_LIBGRIBAPI
!   grid_id_assign_int, &
!#endif
!#ifdef HAVE_LIBGDAL
!   grid_id_assign_dal, &
!this%gdalid = that%gdalid
!#endif
!  grid_id_assign_grid_id
!END INTERFACE


INTERFACE c_e
  MODULE PROCEDURE grid_id_c_e, grid_id_c_e_v, grid_file_id_c_e, grid_file_id_c_e_v
END INTERFACE

INTERFACE display
  MODULE PROCEDURE grid_id_display
END INTERFACE

PRIVATE grid_file_id_delete, grid_id_delete, grid_id_copy, &
 grid_id_c_e, grid_file_id_c_e, grid_id_c_e_v, grid_file_id_c_e_v, grid_id_display

CONTAINS


SUBROUTINE grid_file_id_init(this, filename, mode, driver, from_grid_id)
TYPE(grid_file_id),INTENT(out) :: this ! object to initialise
CHARACTER(len=*),INTENT(in) :: filename ! name of file containing gridded data, in the format [driver:]pathname
CHARACTER(len=*),INTENT(in) :: mode ! access mode for file, 'r' or 'w'
INTEGER,INTENT(in),OPTIONAL :: driver ! select the driver that will be associated to the grid_file_id created, use the constants \a grid_id_notype, \a grid_id_grib_api, \a grid_id_gdal
TYPE(grid_id),INTENT(in),OPTIONAL :: from_grid_id ! select the driver as the one associated to the provided grid_id object

this = grid_file_id_new(filename, mode, driver, from_grid_id)

END SUBROUTINE grid_file_id_init


FUNCTION grid_file_id_new(filename, mode, driver, from_grid_id) RESULT(this)
CHARACTER(len=*),INTENT(in) :: filename
CHARACTER(len=*),INTENT(in) :: mode
INTEGER,INTENT(in),OPTIONAL :: driver
TYPE(grid_id),INTENT(in),OPTIONAL :: from_grid_id
TYPE(grid_file_id) :: this

INTEGER :: n, ier, nf
#ifdef HAVE_LIBGDAL
INTEGER :: imode
#endif
TYPE(csv_record) :: driveropts
CHARACTER(len=12) :: drivername

#ifdef HAVE_LIBGDAL
CALL gdalnullify(this%gdalid)
#endif

IF (filename == '' .OR. .NOT.c_e(filename)) RETURN

n = index(filename,':')
IF (n > 1) THEN ! override with driver from filename
  CALL init(driveropts, filename(:n-1), nfield=nf)
  CALL csv_record_getfield(driveropts, drivername)
#ifdef HAVE_LIBGRIBAPI
  IF (drivername == 'grib_api') THEN
    this%driver = grid_id_grib_api
  ENDIF
#endif
#ifdef HAVE_LIBGDAL
  IF (drivername == 'gdal') THEN
    IF (nf > 4) THEN
      this%driver = grid_id_gdal
      CALL csv_record_getfield(driveropts, this%gdal_options%xmin)
      CALL csv_record_getfield(driveropts, this%gdal_options%ymin)
      CALL csv_record_getfield(driveropts, this%gdal_options%xmax)
      CALL csv_record_getfield(driveropts, this%gdal_options%ymax)
      ! set extreme values if missing, magnitude defined empirically
      ! because of integer overflow in fortrangis
      IF (.NOT.c_e(this%gdal_options%xmin)) this%gdal_options%xmin = -1.0d6
      IF (.NOT.c_e(this%gdal_options%ymin)) this%gdal_options%ymin = -1.0d6
      IF (.NOT.c_e(this%gdal_options%xmax)) this%gdal_options%xmax = 1.0d6
      IF (.NOT.c_e(this%gdal_options%ymax)) this%gdal_options%ymax = 1.0d6
    ELSE
      CALL l4f_log(l4f_error, 'gdal driver requires 4 extra arguments (bounding box)')
      CALL raise_error()
    ENDIF
  ENDIF
#endif
  CALL delete(driveropts)
ENDIF
IF (PRESENT(driver)) THEN ! override with driver
  this%driver = driver
ENDIF
IF (PRESENT(from_grid_id)) THEN ! override with from_grid_id
  this%driver = from_grid_id%driver
ENDIF

#ifdef HAVE_LIBGRIBAPI
IF (this%driver == grid_id_grib_api) THEN
  CALL grib_open_file(this%gaid, filename(n+1:), trim(mode), ier)
  IF (ier /= grib_success) this%gaid = imiss
ENDIF
#endif
#ifdef HAVE_LIBGDAL
IF (this%driver == grid_id_gdal) THEN
  IF (mode(1:1) == 'w') THEN
    imode = ga_update
  ELSE
    imode = ga_readonly
  ENDIF
  CALL gdalallregister()
  this%gdalid = gdalopen(trim(filename(n+1:))//c_null_char, imode)
! dirty trick, with gdal I have to keep a copy of the file_id, memory leak
  ALLOCATE(this%file_id_copy)
  this%file_id_copy = this
ENDIF
#endif

END FUNCTION grid_file_id_new


FUNCTION grid_file_id_count(this) RESULT(count)
TYPE(grid_file_id),INTENT(in) :: this
INTEGER :: count

INTEGER :: ier

count = 0
#ifdef HAVE_LIBGRIBAPI
IF (this%driver == grid_id_grib_api) THEN
  IF (c_e(this%gaid)) THEN
    CALL grib_count_in_file(this%gaid, count, ier)
    IF (ier /= grib_success) count = 0
  ENDIF
ENDIF
#endif
#ifdef HAVE_LIBGDAL
IF (this%driver == grid_id_gdal) THEN
  IF (gdalassociated(this%gdalid)) THEN
    count = gdalgetrastercount(this%gdalid)
  ENDIF
ENDIF
#endif

END FUNCTION grid_file_id_count


SUBROUTINE grid_file_id_delete(this)
TYPE(grid_file_id),INTENT(inout) :: this

#ifdef HAVE_LIBGRIBAPI
IF (this%driver == grid_id_grib_api) THEN
  IF (c_e(this%gaid)) CALL grib_close_file(this%gaid)
ENDIF
this%gaid = imiss
#endif
#ifdef HAVE_LIBGDAL
IF (this%driver == grid_id_gdal) THEN
! dirty trick, with gdal I have to keep the file open
!  IF (gdalassociated(this%gdalid)) CALL gdalclose(this%gdalid)
  this%nlastband = 0
ENDIF
CALL gdalnullify(this%gdalid)
#endif

this%driver = imiss

END SUBROUTINE grid_file_id_delete


! Function to check whether a \a grid_file_id object has been correctly associated
! to the requested file. It returns \a .FALSE. if the file has not
! been opened correctly or if the object has been initialized as empty.
FUNCTION grid_file_id_c_e(this)
TYPE(grid_file_id),INTENT(in) :: this ! object to be checked
LOGICAL :: grid_file_id_c_e

grid_file_id_c_e = .false.

#ifdef HAVE_LIBGRIBAPI
IF (this%driver == grid_id_grib_api) THEN
  grid_file_id_c_e = c_e(this%gaid)
ENDIF
#endif
#ifdef HAVE_LIBGDAL
IF (this%driver == grid_id_gdal) THEN
  grid_file_id_c_e = gdalassociated(this%gdalid)
ENDIF
#endif

END FUNCTION grid_file_id_c_e


! Function to check whether a \a grid_file_id object has been correctly associated
! to the requested file. It returns \a .FALSE. if the file has not
! been opened correctly or if the object has been initialized as empty.
FUNCTION grid_file_id_c_e_v(this)
TYPE(grid_file_id),INTENT(in) :: this(:) ! object to be checked
LOGICAL :: grid_file_id_c_e_v(SIZE(this))

INTEGER :: i

DO i = 1, SIZE(this)
  grid_file_id_c_e_v(i) = c_e(this(i))
ENDDO

END FUNCTION grid_file_id_c_e_v


SUBROUTINE grid_id_init(this, from_grid_file_id, grib_api_template, grib_api_id)
TYPE(grid_id),INTENT(out) :: this ! object to be initialized
TYPE(grid_file_id),INTENT(inout),OPTIONAL :: from_grid_file_id ! file object from which grid object has to be created
CHARACTER(len=*),INTENT(in),OPTIONAL :: grib_api_template ! grib_api template file from which grid_object has to be created
INTEGER,INTENT(in),OPTIONAL :: grib_api_id ! grib_api id obtained directly from a \a grib_get subroutine call

this = grid_id_new(from_grid_file_id, grib_api_template, grib_api_id)

END SUBROUTINE grid_id_init


FUNCTION grid_id_new(from_grid_file_id, grib_api_template, grib_api_id, &
 no_driver_id) RESULT(this)
TYPE(grid_file_id),INTENT(inout),OPTIONAL,TARGET :: from_grid_file_id
CHARACTER(len=*),INTENT(in),OPTIONAL :: grib_api_template
INTEGER,INTENT(in),OPTIONAL :: grib_api_id
INTEGER,INTENT(in),OPTIONAL :: no_driver_id
TYPE(grid_id) :: this

INTEGER :: ier

#ifdef HAVE_LIBGDAL
CALL gdalnullify(this%gdalid)
#endif

IF (PRESENT(from_grid_file_id)) THEN
  this%driver = from_grid_file_id%driver ! take driver from file_id

#ifdef HAVE_LIBGRIBAPI
  IF (this%driver == grid_id_grib_api) THEN
    IF (c_e(from_grid_file_id%gaid)) THEN
      CALL grib_new_from_file(from_grid_file_id%gaid, this%gaid, ier)
      IF (ier /= grib_success) this%gaid = imiss
    ENDIF
  ENDIF
#endif
#ifdef HAVE_LIBGDAL
  IF (this%driver == grid_id_gdal) THEN
    IF (gdalassociated(from_grid_file_id%gdalid) .AND. &
     ASSOCIATED(from_grid_file_id%file_id_copy)) THEN
      IF (from_grid_file_id%nlastband < &
       gdalgetrastercount(from_grid_file_id%gdalid)) THEN ! anything to read?
        from_grid_file_id%nlastband = from_grid_file_id%nlastband + 1
        this%gdalid = &
         gdalgetrasterband(from_grid_file_id%gdalid, from_grid_file_id%nlastband)
        this%file_id => from_grid_file_id%file_id_copy ! for gdal remember copy of file_id

      ENDIF
    ENDIF
  ENDIF
#endif

#ifdef HAVE_LIBGRIBAPI
ELSE IF (PRESENT(grib_api_template)) THEN
  this%driver = grid_id_grib_api
  CALL grib_new_from_samples(this%gaid, grib_api_template, ier)
  IF (ier /= grib_success) this%gaid = imiss
ELSE IF (PRESENT(grib_api_id)) THEN
  this%driver = grid_id_grib_api
  this%gaid = grib_api_id
#endif
ELSE IF (PRESENT(no_driver_id)) THEN
  this%driver = grid_id_no_driver
  this%nodriverid = no_driver_id
ENDIF

END FUNCTION grid_id_new


SUBROUTINE grid_id_delete(this)
TYPE(grid_id),INTENT(inout) :: this

this%nodriverid = imiss
#ifdef HAVE_LIBGRIBAPI
IF (this%driver == grid_id_grib_api) THEN
  IF (c_e(this%gaid)) CALL grib_release(this%gaid)
ENDIF
this%gaid = imiss
#endif
#ifdef HAVE_LIBGDAL
CALL gdalnullify(this%gdalid)
NULLIFY(this%file_id)
#endif

this%driver = imiss

END SUBROUTINE grid_id_delete


FUNCTION grid_id_readonly(this) RESULT(readonly)
TYPE(grid_id),INTENT(in) :: this
LOGICAL :: readonly

readonly = this%driver /= grid_id_grib_api

END FUNCTION grid_id_readonly


SUBROUTINE grid_id_copy(this, that)
TYPE(grid_id),INTENT(in) :: this
TYPE(grid_id),INTENT(out) :: that

that = this ! start with a shallow copy

#ifdef HAVE_LIBGRIBAPI
IF (this%driver == grid_id_grib_api) THEN
  IF (c_e(this%gaid)) THEN
    that%gaid = -1
    CALL grib_clone(this%gaid, that%gaid)
  ENDIF
ENDIF
#endif
#ifdef HAVE_LIBGDAL
IF (this%driver == grid_id_gdal) THEN
  IF (c_e(this)) THEN
!   that = grid_id_new(no_driver_id=1)
!    that%gdalid = this%gdalid ! better idea?
!    that%file_id => this%file_id
  ENDIF
ENDIF
#endif

END SUBROUTINE grid_id_copy


SUBROUTINE grid_id_export(this, file_id)
TYPE(grid_id),INTENT(inout) :: this
TYPE(grid_file_id),INTENT(in) :: file_id

INTEGER :: ier

IF (c_e(this) .AND. c_e(file_id)) THEN
#ifdef HAVE_LIBGRIBAPI
  IF (this%driver == grid_id_grib_api .AND. file_id%driver == grid_id_grib_api) &
   CALL grib_write(this%gaid, file_id%gaid, ier) ! log ier?
#endif
ENDIF
#ifdef HAVE_LIBGDAL
IF (this%driver == grid_id_gdal .AND. file_id%driver == grid_id_gdal) THEN
  ! not implemented, log?
ENDIF
#endif

END SUBROUTINE grid_id_export


! Function to check whether a \a _file_id object has been correctly associated
! to a grid. It returns \a .FALSE. if the grid has not been correctly
! obtained from the file or if the object has been initialized as
! empty.
FUNCTION grid_id_c_e(this)
TYPE(grid_id),INTENT(in) :: this ! object to be checked
LOGICAL :: grid_id_c_e

grid_id_c_e = .false.

#ifdef HAVE_LIBGRIBAPI
IF (this%driver == grid_id_grib_api) THEN
  grid_id_c_e = c_e(this%gaid)
ENDIF
#endif
#ifdef HAVE_LIBGDAL
IF (this%driver == grid_id_gdal) THEN
  grid_id_c_e = gdalassociated(this%gdalid)
ENDIF
#endif
IF (this%driver == grid_id_no_driver) THEN
  grid_id_c_e = c_e(this%nodriverid)
ENDIF

END FUNCTION grid_id_c_e


! Function to check whether a \a _file_id object has been correctly associated
! to a grid. It returns \a .FALSE. if the grid has not been correctly
! obtained from the file or if the object has been initialized as
! empty.
FUNCTION grid_id_c_e_v(this)
TYPE(grid_id),INTENT(in) :: this(:) ! object to be checked
LOGICAL :: grid_id_c_e_v(size(this))

INTEGER :: i

DO i = 1, SIZE(this)
  grid_id_c_e_v(i) = c_e(this(i))
ENDDO

END FUNCTION grid_id_c_e_v


FUNCTION grid_file_id_get_driver(this) RESULT(driver)
TYPE(grid_file_id),INTENT(in) :: this
CHARACTER(len=LEN(driverlist)) :: driver

IF (this%driver > 0 .AND. this%driver <= SIZE(driverlist)) THEN
  driver = driverlist(this%driver)
ELSE
  driver = driverlist(0)
ENDIF

END FUNCTION grid_file_id_get_driver
  

FUNCTION grid_id_get_driver(this) RESULT(driver)
TYPE(grid_id),INTENT(in) :: this
CHARACTER(len=LEN(driverlist)) :: driver

IF (this%driver > 0 .AND. this%driver <= SIZE(driverlist)) THEN
  driver = driverlist(this%driver)
ELSE
  driver = driverlist(0)
ENDIF

END FUNCTION grid_id_get_driver


SUBROUTINE grid_id_display(this, namespace)
TYPE(grid_id),INTENT(in) :: this
CHARACTER(len=*),OPTIONAL :: namespace

INTEGER :: kiter, iret
CHARACTER(len=255) :: key, value, lnamespace


#ifdef HAVE_LIBGRIBAPI
IF (this%driver == grid_id_grib_api) THEN

  lnamespace = optio_c(namespace,255)
  IF (.NOT.c_e(lnamespace))THEN
    lnamespace = "ls"
  ENDIF

  print*,"GRIB_API namespace:",trim(lnamespace)

  CALL grib_keys_iterator_new(this%gaid, kiter, namespace=trim(lnamespace))

  iter: DO
    CALL grib_keys_iterator_next(kiter, iret) 
  
    IF (iret /= 1) THEN
      EXIT iter
    END IF
  
    CALL grib_keys_iterator_get_name(kiter, key)

    IF (key == 'computeStatistics') cycle

    CALL grib_get(this%gaid, trim(key), value, iret)
    IF (iret == 0)THEN
      print*, trim(key)//' = '//trim(VALUE)
    ELSE
      print*, trim(key)//' = '//"KEY NOT FOUND, namespace :"//trim(lnamespace)//" ( bug ? )"
    ENDIF
  ENDDO iter

  CALL grib_keys_iterator_delete(kiter)

ENDIF

#endif

END SUBROUTINE grid_id_display


#ifdef HAVE_LIBGRIBAPI

FUNCTION grid_file_id_get_gaid(this) RESULT(gaid)
TYPE(grid_file_id),INTENT(in) :: this
INTEGER :: gaid
gaid = this%gaid
END FUNCTION grid_file_id_get_gaid

FUNCTION grid_id_get_gaid(this) RESULT(gaid)
TYPE(grid_id),INTENT(in) :: this
INTEGER :: gaid
gaid = this%gaid
END FUNCTION grid_id_get_gaid
#endif


#ifdef HAVE_LIBGDAL

FUNCTION grid_file_id_get_gdalid(this) RESULT(gdalid)
TYPE(grid_file_id),INTENT(in) :: this
TYPE(gdaldataseth) :: gdalid
gdalid = this%gdalid
END FUNCTION grid_file_id_get_gdalid

FUNCTION grid_id_get_gdalid(this) RESULT(gdalid)
TYPE(grid_id),INTENT(in) :: this
TYPE(gdalrasterbandh) :: gdalid
gdalid = this%gdalid
END FUNCTION grid_id_get_gdalid

FUNCTION grid_id_get_gdal_options(this) RESULT(gdal_options)
TYPE(grid_id),INTENT(in) :: this
TYPE(gdal_file_id_options) :: gdal_options

TYPE(gdal_file_id_options) :: gdal_options_local

IF (ASSOCIATED(this%file_id)) THEN
  gdal_options = this%file_id%gdal_options
ELSE
  gdal_options = gdal_options_local ! empty object
ENDIF

END FUNCTION grid_id_get_gdal_options
#endif


SUBROUTINE grid_id_decode_data(this, field, done)
TYPE(grid_id),INTENT(in) :: this
REAL,INTENT(out) :: field(:,:)
LOGICAL,OPTIONAL,INTENT(out) :: done

LOGICAL :: ldone

ldone = .false.
#ifdef HAVE_LIBGRIBAPI
IF (c_e(this%gaid)) THEN
  CALL grid_id_decode_data_gribapi(this%gaid, field)
  ldone = .true.
ENDIF
#endif
#ifdef HAVE_LIBGDAL
! subarea?
IF (gdalassociated(this%gdalid)) THEN
  CALL grid_id_decode_data_gdal(this%gdalid, field, this%file_id%gdal_options)
  ldone = .true.
ENDIF
#endif
IF (.NOT.ldone) THEN
  field(:,:) = rmiss
ENDIF
IF (PRESENT(done)) done = ldone

END SUBROUTINE grid_id_decode_data


SUBROUTINE grid_id_encode_data(this, field)
TYPE(grid_id),INTENT(inout) :: this
REAL,intent(in) :: field(:,:)

#ifdef HAVE_LIBGRIBAPI
IF (this%driver == grid_id_grib_api) THEN

!  call display(this,"parameter")
!  print *,minval(field,c_e(field)),maxval(field,c_e(field))
!  print *,"-----------------------"

  IF (c_e(this%gaid)) CALL grid_id_encode_data_gribapi(this%gaid, field)
ENDIF
#endif
#ifdef HAVE_LIBGDAL
IF (this%driver == grid_id_gdal) THEN
!gdalid = grid_id_get_gdalid(this%gaid)
  CALL l4f_log(l4f_warn,"export to gdal not implemented" )
! subarea?
ENDIF
#endif

END SUBROUTINE grid_id_encode_data


#ifdef HAVE_LIBGRIBAPI
SUBROUTINE grid_id_decode_data_gribapi(gaid, field)
INTEGER,INTENT(in) :: gaid ! grib_api id
REAL,INTENT(out),TARGET,CONTIGUOUS :: field(:,:) ! array of decoded values

INTEGER :: editionnumber
INTEGER :: alternativerowscanning, &
 iscansnegatively, jscanspositively, jpointsareconsecutive
#ifdef DEBUG
INTEGER :: numberofvalues
#endif
INTEGER :: numberofpoints
REAL :: vector(size(field))
REAL,POINTER :: fieldp(:)
INTEGER :: x1, x2, xs, y1, y2, ys, ord(2), ierr


call grib_get(gaid,'GRIBEditionNumber',editionnumber)

if (editionnumber == 2) then

  CALL grib_get(gaid,'alternativeRowScanning',alternativerowscanning,ierr)
  IF (ierr == grib_success .AND. alternativerowscanning /= 0) THEN
    CALL l4f_log(l4f_error, "grib_api alternativeRowScanning not supported: " &
     //t2c(alternativerowscanning))
    CALL raise_error()
    field(:,:) = rmiss
    RETURN
  ENDIF

else if (editionnumber /= 1) then

  CALL l4f_log(l4f_error, &
   "grib_api GribEditionNumber not supported: "//t2c(editionnumber))
  call raise_error()
  field(:,:) = rmiss
  RETURN

end if

CALL grib_get(gaid,'iScansNegatively',iscansnegatively,ierr)
IF (ierr /= grib_success) iscansnegatively=0
CALL grib_get(gaid,'jScansPositively',jscanspositively,ierr)
IF (ierr /= grib_success) jscanspositively=1
CALL grib_get(gaid,'jPointsAreConsecutive',jpointsareconsecutive,ierr)
IF (ierr /= grib_success) jpointsareconsecutive=0

call grib_get(gaid,'numberOfPoints',numberofpoints)

IF (numberofpoints /= SIZE(field)) THEN
  CALL l4f_log(l4f_error, 'grid_id_decode_data_gribapi numberOfPoints and grid size different')
  CALL l4f_log(l4f_error, 'grid_id_decode_data_gribapi numberOfPoints: ' &
   //t2c(numberofpoints)//', nx,ny: '&
   //t2c(SIZE(field,1))//','//t2c(SIZE(field,2)))
  CALL raise_error()
  field(:,:) = rmiss
  RETURN
ENDIF

! get data values
#ifdef DEBUG
call grib_get(gaid,'numberOfValues',numberofvalues)
call l4f_log(l4f_info,'grib_api number of values: '//to_char(numberofvalues))
call l4f_log(l4f_info,'grib_api number of points: '//to_char(numberofpoints))
#endif

CALL grib_set(gaid,'missingValue',rmiss)

IF (iscansnegatively == 0 .AND. jscanspositively == 0 .AND. jpointsareconsecutive == 0) THEN

  fieldp(1:SIZE(field)) => field
  CALL grib_get(gaid,'values',fieldp)

ELSE

CALL grib_get(gaid,'values',vector)
! suspect bug in grib_api, when all field is missing it is set to zero
! IF (numberOfValues == 0) vector = rmiss

#ifdef DEBUG
CALL l4f_log(l4f_debug, 'grib_api, decoded field in interval: '// &
 t2c(minval(vector,mask=c_e(vector)))//' '//t2c(maxval(vector,mask=c_e(vector))))
CALL l4f_log(l4f_debug, 'grib_api, decoded field with number of missing: '// &
 t2c(count(.NOT.c_e(vector))))

IF (numberofvalues /= count(c_e(vector))) THEN
  CALL l4f_log(l4f_warn, 'grid_id_decode_data_gribapi numberOfValues and valid data count different')
  CALL l4f_log(l4f_warn, 'grid_id_decode_data_gribapi numberOfValues: ' &
   //t2c(numberofvalues)//', valid data: '//t2c(count(c_e(vector))))
!  CALL raise_warning()
ENDIF
#endif

! Transfer data field changing scanning mode to 64
IF (iscansnegatively == 0) THEN
  x1 = 1
  x2 = SIZE(field,1)
  xs = 1
ELSE
  x1 = SIZE(field,1)
  x2 = 1
  xs = -1
ENDIF
IF (jscanspositively == 0) THEN
  y1 = SIZE(field,2)
  y2 = 1
  ys = -1
ELSE
  y1 = 1
  y2 = SIZE(field,2)
  ys = 1
ENDIF

IF (jpointsareconsecutive == 0) THEN
  ord = (/1,2/)
ELSE
  ord = (/2,1/)
ENDIF

field(x1:x2:xs,y1:y2:ys) = reshape(vector, &
 (/SIZE(field,1),SIZE(field,2)/), order=ord)

ENDIF

END SUBROUTINE grid_id_decode_data_gribapi


SUBROUTINE grid_id_encode_data_gribapi(gaid, field)
INTEGER,INTENT(in) :: gaid ! grib_api id
REAL,intent(in) :: field(:,:) ! data array to be encoded

INTEGER :: editionnumber
INTEGER :: alternativerowscanning, iscansnegatively, &
 jscanspositively, jpointsareconsecutive
INTEGER :: x1, x2, xs, y1, y2, ys, ierr

call grib_get(gaid,'GRIBEditionNumber',editionnumber)

if (editionnumber == 2) then

  CALL grib_get(gaid,'alternativeRowScanning',alternativerowscanning,ierr)
  IF (ierr == grib_success .AND. alternativerowscanning /= 0)THEN
    CALL l4f_log(l4f_error, "grib_api alternativeRowScanning not supported: " &
     //trim(to_char(alternativerowscanning)))
    CALL raise_error()
    RETURN
  ENDIF

else if( editionnumber /= 1) then

  call l4f_log(l4f_error, &
   "grib_api GribEditionNumber not supported: "//t2c(editionnumber))
  call raise_error()
  RETURN

end if

CALL grib_get(gaid,'iScansNegatively',iscansnegatively,ierr)
IF (ierr /= grib_success) iscansnegatively=0
CALL grib_get(gaid,'jScansPositively',jscanspositively,ierr)
IF (ierr /= grib_success) jscanspositively=1
CALL grib_get(gaid,'jPointsAreConsecutive',jpointsareconsecutive,ierr)
IF (ierr /= grib_success) jpointsareconsecutive=0

! these grib_sets are alredy done in gridinfo_export, but it seems
! that it is necessary to repeat them here, they can fail with
! unstructured grids, thus ierr
#ifdef DEBUG
CALL l4f_log(l4f_debug, 'grib_api, Ni,Nj:'//t2c(SIZE(field,1))//','//t2c(SIZE(field,2)))
#endif
CALL grib_set(gaid,'Ni',SIZE(field,1), ierr)
CALL grib_set(gaid,'Nj',SIZE(field,2), ierr)

! Transfer data field changing scanning mode from 64
IF (iscansnegatively  == 0) THEN
  x1 = 1
  x2 = SIZE(field,1)
  xs = 1
ELSE
  x1 = SIZE(field,1)
  x2 = 1
  xs = -1
ENDIF
IF (jscanspositively == 0) THEN
  y1 = SIZE(field,2)
  y2 = 1
  ys = -1
ELSE
  y1 = 1
  y2 = SIZE(field,2)
  ys = 1
ENDIF


IF (any(field == rmiss)) THEN

  CALL grib_set(gaid,'missingValue',rmiss)
  IF (editionnumber == 1) THEN
! enable bitmap in grib1
! grib_api 1.9.9 went into an infinite loop with second order packing here
!    CALL grib_set(gaid,'packingType','grid_simple')
! now it's fixed, leaving second order if present
    CALL grib_set(gaid,"bitmapPresent",1)
  ELSE
! enable bitmap in grib2
    CALL grib_set(gaid,"bitMapIndicator",0)
  ENDIF

ELSE

  IF (editionnumber == 1) THEN
! disable bitmap in grib1
    CALL grib_set(gaid,"bitmapPresent",0)
  ELSE
! disable bitmap in grib2
    CALL grib_set(gaid,"bitMapIndicator",255)
  ENDIF

ENDIF

#ifdef DEBUG
CALL l4f_log(l4f_debug, 'grib_api, coding field in interval: '// &
 t2c(minval(field,mask=c_e(field)))//' '//t2c(maxval(field,mask=c_e(field))))
CALL l4f_log(l4f_debug, 'grib_api, coding field with number of missing: '// &
 t2c(count(.NOT.c_e(field))))
CALL l4f_log(l4f_debug, 'grib_api, sizex:'//t2c(x1)//','//t2c(x2)//','//t2c(xs))
CALL l4f_log(l4f_debug, 'grib_api, sizey:'//t2c(y1)//','//t2c(y2)//','//t2c(ys))
#endif
IF (jpointsareconsecutive == 0) THEN
  CALL grib_set(gaid,'values', reshape(field(x1:x2:xs,y1:y2:ys), &
   (/SIZE(field)/)))
ELSE
  CALL grib_set(gaid,'values', reshape(transpose(field(x1:x2:xs,y1:y2:ys)), &
   (/SIZE(field)/)))
ENDIF

END SUBROUTINE grid_id_encode_data_gribapi
#endif


#ifdef HAVE_LIBGDAL
SUBROUTINE grid_id_decode_data_gdal(gdalid, field, gdal_options)
#ifdef F2003_FULL_FEATURES
USE ieee_arithmetic
#endif
TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal id
REAL,INTENT(out) :: field(:,:) ! array of decoded values
TYPE(gdal_file_id_options),INTENT(in) :: gdal_options

TYPE(gdaldataseth) :: hds
REAL(kind=c_double) :: geotrans(6), dummy1, dummy2, dummy3, dummy4
REAL :: gdalmiss
REAL,ALLOCATABLE :: buffer(:,:)
INTEGER :: ix1, iy1, ix2, iy2, ixs, iys, ord(2), ier
INTEGER(kind=c_int) :: nrx, nry
LOGICAL :: must_trans


hds = gdalgetbanddataset(gdalid) ! go back to dataset
ier = gdalgetgeotransform(hds, geotrans)

IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN
! transformation is diagonal, no transposing
  IF (geotrans(2) > 0.0_c_double) THEN
    ix1 = 1 
    ix2 = SIZE(field,1)
    ixs = 1
  ELSE
    ix1 = SIZE(field,1)
    ix2 = 1 
    ixs = -1
  ENDIF
  IF (geotrans(6) > 0.0_c_double) THEN
    iy1 = 1 
    iy2 = SIZE(field,2)
    iys = 1
  ELSE
    iy1 = SIZE(field,2)
    iy2 = 1 
    iys = -1
  ENDIF
  nrx = SIZE(field,1)
  nry = SIZE(field,2)
  ord = (/1,2/)
  must_trans = .false.
!  ALLOCATE(buffer(nrx,nry))

ELSE IF (geotrans(2) == 0.0_c_double .AND. geotrans(6) == 0.0_c_double) THEN
! transformation is anti-diagonal, transposing required this should not happen
  IF (geotrans(3) > 0.0_c_double) THEN
    ix1 = 1 
    ix2 = SIZE(field,1)
    ixs = 1
  ELSE
    ix1 = SIZE(field,1)
    ix2 = 1 
    ixs = -1
  ENDIF
  IF (geotrans(5) > 0.0_c_double) THEN
    iy1 = 1 
    iy2 = SIZE(field,2)
    iys = 1
  ELSE
    iy1 = SIZE(field,2)
    iy2 = 1 
    iys = -1
  ENDIF
  nrx = SIZE(field,2)
  nry = SIZE(field,1)
  ord = (/2,1/)
  must_trans = .true.
!  ALLOCATE(buffer(nry,nrx))

ELSE ! transformation is a rotation, not supported
  CALL l4f_log(l4f_error, 'gdal geotransform is a generic rotation, not supported')
  CALL raise_error()
  field(:,:) = rmiss
  RETURN
ENDIF

! read data from file
CALL gdalrastersimpleread_f(gdalid, gdal_options%xmin, gdal_options%ymin, &
 gdal_options%xmax, gdal_options%ymax, buffer, dummy1, dummy2, dummy3, dummy4)

IF (.NOT.ALLOCATED(buffer)) THEN ! error in read
  CALL l4f_log(l4f_error, 'gdal error in reading with gdal driver')
  CALL raise_error()
  field(:,:) = rmiss
  RETURN
ENDIF

IF (SIZE(buffer) /= (SIZE(field)))THEN
  CALL l4f_log(l4f_error, 'gdal raster band and gridinfo size different')
  CALL l4f_log(l4f_error, 'gdal rasterband: ' &
   //t2c(SIZE(buffer,1))//'X'//t2c(SIZE(buffer,2))//', nx,ny:' &
   //t2c(SIZE(field,ord(1)))//'X'//t2c(SIZE(field,ord(2))))
  CALL raise_error()
  field(:,:) = rmiss
  RETURN
ENDIF

#ifdef F2003_FULL_FEATURES
! aparently gdal datasets may contain NaN
WHERE(ieee_is_nan(buffer))
  buffer = rmiss
END WHERE
#else
WHERE(buffer /= buffer)
  buffer = rmiss
END WHERE
#endif

! set missing value if necessary
gdalmiss = REAL(gdalgetrasternodatavalue(gdalid, ier))
IF (ier /= 0) THEN ! success -> there are missing values
#ifdef DEBUG
  CALL l4f_log(l4f_info, 'gdal missing data value: '//trim(to_char(gdalmiss)))
#endif
  WHERE(buffer(:,:) == gdalmiss)
    buffer(:,:) = rmiss
  END WHERE
ELSE
#ifdef DEBUG
  CALL l4f_log(l4f_info, 'gdal no missing data found in band')
#endif
ENDIF

! reshape the field
IF (must_trans) THEN
  field(ix1:ix2:ixs,iy1:iy2:iys) = transpose(buffer)
ELSE
  field(ix1:ix2:ixs,iy1:iy2:iys) = buffer(:,:)
ENDIF


END SUBROUTINE grid_id_decode_data_gdal
#endif

END MODULE grid_id_class