ecearth3/sources/tm5mp/base/file_grib_api.F90

!#######################################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status >0) then; TRACEBACK; action; return; end if
!
! trap errors from GRIB_API:
#define IF_GRIB_NOTOK_RETURN(action) if (status/=GRIB_SUCCESS) then; call GRIB_Get_Error_String(status,gol); call goErr; TRACEBACK; action; return; end if
!
!#######################################################################

MODULE FILE_GRIB_API

  USE GO,       only : gol, goPr, goErr
  USE GRIB_API, only : HGRIB => kindOfInt  ! handle kind
  USE GRIB_API, only : GRIB_SUCCESS, GRIB_Get_Error_String

  IMPLICIT NONE
  PRIVATE

  PUBLIC :: TGribFile
  PUBLIC :: Init, Done, Opened  
  PUBLIC :: Get
  PUBLIC :: Check  
  PUBLIC :: ReadRecord

  ! --- const ------------------------------
  
  character(len=*), parameter  ::  mname = 'file_grib_api'  
  
  ! --- types -------------------------------

  TYPE TGribFile
    integer(kind=HGRIB)    ::  fu
    character(len=256)     ::  fname
    ! grib id and return value (end_of_file check)
    integer(kind=HGRIB)    ::  igrib
    integer                ::  iret
    logical                :: opened
  END TYPE TGribFile

  ! --- interfaces ---------------------------

  interface Init
    module procedure gribfile_Init
  end interface

  interface Done
    module procedure gribfile_Done
  end interface
  
  interface Opened
    module procedure gribfile_Opened
  end interface
  
  interface ReadRecord
    module procedure gribfile_ReadRecord
  end interface
  
  interface Get
    module procedure gribfile_Get
  end interface

  interface Check
     module procedure gribfile_Check
  end interface
  

CONTAINS

  ! =================================================================

  ! Open gribfile.
  ! 
  ! USAGE
  !   call Init( gribfile, 'input.gb', 'r'|'w'  )
  !
  ! DESCRIPTION
  !   Interface around routine 'pbOpen'.
  !   In 'gribfile', space is allocated to store grib sections.
  !

  subroutine gribfile_Init( F, file, mode, status )

    use Grib_API, only : Grib_Open_File

    ! --- in/out -----------------------

    type(TGribFile), intent(out)              ::  F
    character(len=*), intent(in)              ::  file
    character(len=*), intent(in)              ::  mode
    integer, intent(out)                      ::  status
    
    ! --- const ------------------------
    
    character(len=*), parameter  ::  rname = mname//'/gribfile_Init'

    ! --- begin ------------------------

    F%fname = file              ! filename
    
    select case ( mode )
    case ( 'r' )

       call grib_open_file ( F%fu, file, 'r' )
       F%opened=.true.

    case ( 'w' )
      
       write (gol,'("writing not yet supported for grib-api")') ; call goErr
       TRACEBACK; status=1; return
       
    case default
      
       write (gol,'("unknown mode `",a,"`")') mode; call goErr
       TRACEBACK; status=1; return
        
    end select
    
    ! no message loaded yet:
    F%igrib = -1
    
  end subroutine gribfile_Init

  ! ===

  ! Close gribfile, clear buffers
  ! 
  ! USAGE
  !   call Done( gribfile )
  !
  ! DESCRIPTION
  !   Interface around routine 'pbClose'.
  !

  subroutine gribfile_Done( F, status )

    use Grib_API, only : GRIB_Release
    use Grib_API, only : Grib_Close_File

    ! --- in/out -----------------------

    type(TGribFile), intent(inout)  ::  F
    integer, intent(inout)          ::  status

    ! --- const ------------------------
    
    character(len=*), parameter  ::  rname = mname//'/gribfile_Done'

    ! --- begin ------------------------

    ! old message loaded ?
    if ( F%igrib > 0 ) then
      ! remove from memory:
      call GRIB_Release( F%igrib, status=status )
      IF_GRIB_NOTOK_RETURN(status=1)
      ! no message loaded yet:
      F%igrib = -1
    end if
       
    ! write(*,'(1x,a,a)') 'Close grib file ',trim(F%fname)
    call grib_close_file ( F%fu, status )
    if ( status /= 0 ) then
       write (gol,'("from closing grib file:")'); call goErr
       write (gol,'("  status : ",i6)') status; call goErr
       write (gol,'("  file   : ",a)') trim(F%fname); call goErr
       TRACEBACK; status=1; return
    else
       F%opened=.false.
    end if

  end subroutine gribfile_Done


  ! ===
  
  
  ! Grib file opened ?
  
  logical function gribfile_Opened( gribfile ) 
  
    ! --- in/out ------------------------------
    
    type(TGribFile), intent(in)     ::  gribfile
    
    ! --- begin -------------------------
    
    !inquire( unit=gribfile%fu, opened=gribfile_Opened )
    gribfile_Opened=gribfile%opened

  end function gribfile_Opened
  

  ! ==================================================================
  
  ! USAGE
  !   call ReadRecord( gribfile [,debug=0|1] [,status=status] )
  !
  ! DESCRIPTION
  !   Reads next grib record (a/k/a grib message) in buffers.
  !
  ! RETURN STATUS
  !    0 : no error
  !    1 : eof
  !    2 : some error
  !   Execution might stop if other errors are encountered.
  !
  subroutine gribfile_ReadRecord( F, status )

    use Grib_API, only : GRIB_END_OF_FILE
    use Grib_API, only : GRIB_SUCCESS
    use Grib_API, only : Grib_New_From_File
    use Grib_API, only : Grib_Get
    use GRIB_API, only : GRIB_Release

    ! --- in/out -------------------------

    type(TGribFile), intent(inout) :: F
    integer,         intent(out)   :: status
    
    ! --- const --------------------------
    
    character(len=*), parameter   ::  rname = mname//'/gribfile_ReadRecord'

    integer                ::  iret
    logical                ::  verbose
    
    ! --- begin --------------------------

    ! old message loaded ?
    if ( F%igrib > 0 ) then
      ! remove from memory:
      call GRIB_Release( F%igrib, status=status )
      IF_GRIB_NOTOK_RETURN(status=1)
      ! no message loaded yet:
      F%igrib = -1
    end if
    
    !     a new grib message is loaded from file
    !     igrib is the grib id to be used in subsequent calls
    !     ret is the return value to check on end_of_file

    call  grib_new_from_file ( F%fu, F%igrib, iret ) 

    if ( iret == GRIB_END_OF_FILE ) then
       write (gol,'("end-of-file is hit before a GRIB product is read")') ; call goErr
       write (gol,'("  grib file : ",a)') F%fname                         ; call goErr
       TRACEBACK; status = 1; return
    end if

    if ( iret /= GRIB_SUCCESS ) then
       write (gol,'("problem reading GRIB data block")') ; call goErr
       write (gol,'("  grib file : ",a)') F%fname        ; call goErr
       call grib_get_error_string(iret,gol)
       call goErr
       TRACEBACK;  status = 2; return
    end if
    
    ! no error ...
    status = 0
    
  end subroutine gribfile_ReadRecord


  ! ============================================================


  ! Extract data from grib sections.
  ! Optional arguments are passed to read parameter id, date, etc.
  ! If the same option is preceded by 'check_', an error
  ! message is written if the grib file does not match
  !
  !  pid   :  parameter identifier (integer number)
  !
  !  level :  index of hybride level
  !
  !  reftime   :  5 element integer array : yy, mm, dd, hh, min
  !  timerange :  4 element ingeger array : unit, val1, val2, indicator
  !

  subroutine gribfile_Get( F, status, &
                           model_id, pid, pidShortName, &
                           levtype, level, nlev, hyb_a, hyb_b, &
                           reftime, timerange, &
                           gridtype, &
                           ll, &
                           lon_first, lon_last, lon_inc, lon_n, &
                           lat_first, lat_last, lat_inc, lat_n, &
                           T, sh, &
                           N, gg  )

    use Grib_API,       only : grib_get, grib_get_size
    use file_grib_base, only : gridtype_ll, gridtype_gg, gridtype_sh, gridtype_red_gg
    use file_grib_base, only : levtype_sfc, levtype_hyb, levtype_land, levtype_hyb2

    ! --- in/out -------------------------

    type(TGribFile), intent(in)      ::  F
    integer,         intent(out)     ::  status

    integer,       intent(out), optional ::  model_id         
    integer,       intent(out), optional ::  pid         
    character(72), intent(out), optional ::  pidShortName         
    integer,       intent(out), optional ::  levtype, level, nlev
    real,          intent(out), optional ::  hyb_a(:), hyb_b(:)      
    integer,       intent(out), optional ::  reftime(5)  
    integer,       intent(out), optional ::  timerange(4)
    integer,       intent(out), optional ::  gridtype
    real,          intent(out), optional ::  lon_first, lon_last, lon_inc
    integer,       intent(out), optional ::  lon_n
    real,          intent(out), optional ::  lat_first, lat_last, lat_inc
    integer,       intent(out), optional ::  lat_n
    real,          intent(out), optional ::  ll(:,:)
    integer,       intent(out), optional ::  T
    complex,       intent(out), optional ::  sh(:)
    integer,       intent(out), optional ::  N
    real,          intent(out), optional ::  gg(:)
   
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/gribfile_Get'
    
    ! --- local --------------------------
    
    integer                ::  grib_reftime(5)
    integer                ::  grib_timerange(4)
    integer                ::  PVPresent, nb_pv, step, iret, ni, nj
    integer                ::  dataDate, dataTime
    integer                ::  numberOfValues 
    real, allocatable      ::  pv(:)
    real, allocatable      ::  values(:)
    character(72)          ::  gridtypestring
    integer :: TypeOfLevel, edition

    ! --- begin --------------------------
    
    ! this grib1/grib2 business is getting messy!
    call grib_get(F%igrib, 'editionNumber', edition)
    
    ! -- model id

    if ( present( model_id ) ) then
       call grib_get(F%igrib, 'generatingProcessIdentifier', model_id )
    end if

    ! --  parameter

    if ( present( pid ) ) then
       call grib_get(F%igrib, 'paramId', pid )
    end if

    ! --  parameter short name

    if ( present( pidShortName ) ) then
       call grib_get(F%igrib, 'shortName', pidShortName )
    end if

    ! --  level

    if ( present(levtype) ) then
       call grib_get(F%igrib, 'levtype', levtype )
       ! make TM5 code works with edition2:
       if ((edition==2).and.(levtype==levtype_hyb2)) then
          levtype=levtype_hyb
       endif
    endif
    
    if ( present(level  ) ) &
         call grib_get(F%igrib, 'level', level )

    ! -- Vertical Coordinate Parameters 

    if ( present(hyb_a) .or. present(hyb_b) ) then

       if ( .not. (present(hyb_a) .and. present(hyb_b)) ) then
          write (gol,'("extract both hybrid params or none")'); call goErr
          TRACEBACK ; status = 1 ; return
       end if
       
       if (edition==1) then
          call grib_get(F%igrib,'indicatorOfTypeOfLevel', TypeOfLevel)
       
          select case (TypeOfLevel)
          case(levtype_hyb)
             ! get size, check, read, and split into 'a' and 'b'
             call grib_get_size(F%igrib,'pv',nb_pv)   
             if ( nb_pv /= size(hyb_a) + size(hyb_b) ) then
                write (gol,'("numbers of hybride parameters do not match:")'); call goErr
                write (gol,'("  expected : ",i6)') size(hyb_a)+size(hyb_b); call goErr
                write (gol,'("  found    : ",i6)') nb_pv; call goErr
                write (gol,'("  file     : ",a)') trim(F%fname); call goErr
                TRACEBACK ; status = 1 ; return
             end if
             allocate(pv(nb_pv))
             call grib_get(F%igrib,'pv',pv) 
             hyb_a = pv(1:size(hyb_a))
             hyb_b = pv(size(hyb_a)+1:nb_pv)
             deallocate(pv)
             !       case(levtype_sfc)
             !          !surface
             !          nlev=1
          case default
             write (gol,'("there are no hybride parameters in the grib message")'); call goErr
             write (gol,'("  file     : ",a)') trim(F%fname); call goErr
             TRACEBACK ; status = 1 ; return
          end select
          
       elseif (edition==2) then
          call grib_get(F%igrib,'PVPresent',PVPresent)
          if (PVPresent == 1) then
             ! get size, check, read, and split into 'a' and 'b'
             call grib_get_size(F%igrib,'pv',nb_pv)   
             if ( nb_pv /= size(hyb_a) + size(hyb_b) ) then
                write (gol,'("numbers of hybride parameters do not match:")'); call goErr
                write (gol,'("  expected : ",i6)') size(hyb_a)+size(hyb_b); call goErr
                write (gol,'("  found    : ",i6)') nb_pv; call goErr
                write (gol,'("  file     : ",a)') trim(F%fname); call goErr
                TRACEBACK ; status = 1 ; return
             end if
             allocate(pv(nb_pv))
             call grib_get(F%igrib,'pv',pv) 
             hyb_a = pv(1:size(hyb_a))
             hyb_b = pv(size(hyb_a)+1:nb_pv)
             deallocate(pv)
          else
             write (gol,'("there are no hybride parameters in the grib message")'); call goErr
             write (gol,'("  file     : ",a)') trim(F%fname); call goErr
             TRACEBACK ; status = 1 ; return
          end if
       end if
  
 !PLS: PVPresent flag is not available anymore. (for ei in grib1 I think)
 !PLS      
 !PLS       ! set PVPresent as an integer 
 !PLS       call grib_get(F%igrib,'PVPresent',PVPresent)
 !PLS       if (PVPresent == 1) then 
 !PLS         ! number of parameters:
 !PLS         call grib_get_size(F%igrib,'pv',nb_pv) 
 !PLS         ! parameters are stored a(0),..,a(L),b(0),..,b(L),
 !PLS         ! thus length should be combined size of 'a' and 'b' arrays:
 !PLS         if ( nb_pv /= size(hyb_a) + size(hyb_b) ) then
 !PLS            write (gol,'("numbers of hybride parameters do not match:")'); call goErr
 !PLS            write (gol,'("  expected : ",i6)') size(hyb_a)+size(hyb_b); call goErr
 !PLS            write (gol,'("  found    : ",i6)') nb_pv; call goErr
 !PLS            write (gol,'("  file     : ",a)') trim(F%fname); call goErr
 !PLS            TRACEBACK ; status = 1 ; return
 !PLS         end if
 !PLS         ! storage for all parameters:
 !PLS         allocate(pv(nb_pv))
 !PLS         ! read:
 !PLS         call grib_get(F%igrib,'pv',pv) 
 !PLS         ! split into 'a' and 'b' coeff:
 !PLS         hyb_a = pv(1:size(hyb_a))
 !PLS         hyb_b = pv(size(hyb_a)+1:nb_pv)
 !PLS         ! clear:
 !PLS         deallocate(pv)
 !PLS       else 
 !PLS          write (gol,'("there are no hybride parameters in the grib message")'); call goErr
 !PLS          write (gol,'("  file     : ",a)') trim(F%fname); call goErr
 !PLS          TRACEBACK ; status = 1 ; return
 !PLS       end if

    end if
    
    ! number of levels

    if ( present(nlev) ) then

       if (edition==1) then
          call grib_get(F%igrib,'indicatorOfTypeOfLevel', TypeOfLevel)
       
          select case (TypeOfLevel)
          case(levtype_hyb)
             call grib_get_size(F%igrib,'pv',nb_pv)   
             ! half level coeffs are hyba(0:nlev) and hybb(0:nlev), compute number of (full) levels:
             nlev = nb_pv/2 - 1
          case(levtype_sfc)
             !surface
             nlev=1
          case(levtype_land)
             nlev=1
          case default
             write (gol,'("do not know about that TypeOfLevel")'); call goErr
             TRACEBACK ; status = 1 ; return
          end select
       elseif (edition==2) then
          
          call grib_get(F%igrib,'PVPresent',PVPresent)
          if (PVPresent == 1) then 
             ! get total number of half level hyb. coeffs:
             call grib_get_size(F%igrib,'pv',nb_pv) 
             ! half level coeffs are hyba(0:nlev) and hybb(0:nlev), compute number of (full) levels:
             nlev = nb_pv/2 - 1
          else 
             write (gol,'("there are no hybrid parameters in the grib message")'); call goErr
             TRACEBACK ; status = 1 ; return
          end if
       end if
    endif

    ! --  reference time

    if ( present(reftime) ) then
       ! extract reference time:
       call grib_get(F%igrib,'dataDate',dataDate)
       call grib_get(F%igrib,'dataTime',dataTime)
       grib_reftime(1) = dataDate/10000           ! year
       grib_reftime(2) = Mod(dataDate/100,100)    ! month
       grib_reftime(3) = Mod(datadate,100)        ! day
       grib_reftime(4) = dataTime/100             ! hour
       grib_reftime(5) = Mod(dataTime,100)        ! minutes
       reftime = grib_reftime
    end if

    ! --  time range

    if ( present(timerange) ) then
       ! write (*,'(" ",a,"  WARNING timerange argument not fully supported ")') trim(name)
       ! Available from Grib-api: stepType, stepUnits, startStep, endStep, stepRange, step
       call grib_get(F%igrib,'step',step)
       grib_timerange(1:4) = (/ 1, step, 0, 0 /)
       timerange = grib_timerange
    end if

    ! ---  grid
    
    if ( present(gridtype)  ) then
       call grib_get(F%igrib,'gridType',gridtypestring)
       select case ( trim(gridtypestring) )
       case ( 'regular_ll' ) ; gridtype = gridtype_ll
       case ( 'regular_gg' ) ; gridtype = gridtype_gg     ! gribex gives that for reduced.
       case ( 'regular_sh' ) ; gridtype = gridtype_sh
       case ( 'sh'         ) ; gridtype = gridtype_sh     ! used by grid_api
       case ( 'reduced_gg' ) ; gridtype = gridtype_red_gg ! used by grid_api
       case default
          write (gol,'("unsupported grid type string : ",a)') trim(gridtypestring); call goErr
          TRACEBACK ; status = 1 ; return
       end select
    end if
    if ( present(lon_n)     ) call grib_get(F%igrib,'Ni',lon_n)
    if ( present(lat_n)     ) call grib_get(F%igrib,'Nj',lat_n)
    if ( present(lat_first) ) call grib_get(F%igrib,'latitudeOfFirstGridPointInDegrees',lat_first)
    if ( present(lon_first) ) call grib_get(F%igrib,'longitudeOfFirstGridPointInDegrees',lon_first)
    if ( present(lat_last)  ) call grib_get(F%igrib,'latitudeOfLastGridPointInDegrees',lat_last)
    if ( present(lon_last)  ) call grib_get(F%igrib,'longitudeOfLastGridPointInDegrees',lon_last)
    if ( present(lon_inc)   ) call grib_get(F%igrib,'iDirectionIncrementInDegrees',lon_inc)
    if ( present(lat_inc)   ) call grib_get(F%igrib,'jDirectionIncrementInDegrees',lat_inc)

    ! -- lat/lon grid
    
    if ( present(ll) ) then
       !call Check( F, gridtype=gridtype_ll )
       call grib_get(F%igrib,'gridType',gridtypestring)
       if ( trim(gridtypestring) /= 'regular_ll' ) then
          write (gol,'("grid types do not match:")'); call goErr
          write (gol,'("  expected : regular_ll")') ; call goErr
          write (gol,'("  found    : ",i4)') gridtype; call goErr
          TRACEBACK ; status = 1 ; return
       end if
       call grib_get(F%igrib,'Ni',ni)
       call grib_get(F%igrib,'Nj',nj)
       if ( (size(ll,1)/=ni) .or. (size(ll,2)/=nj) ) then
          write (gol,'("grid sizes do not match:")'); call goErr
          write (gol,'("  expected : ",i4," x ",i4)') size(ll,1), size(ll,2); call goErr
          write (gol,'("  found    : ",i4," x ",i4)') ni, nj; call goErr
          TRACEBACK ; status = 1 ; return
       end if
       !     get the size of the values array, and allocate storage
       call grib_get_size(F%igrib,'values',numberOfValues)
       allocate(values(numberOfValues), stat=iret)
       ! get the field
       call grib_get(F%igrib,'values',values)
       ll = reshape( values(1:size(ll)), shape(ll) )
       ! release storage
       deallocate(values)

    end if

    ! -- gaussian lat/lon grid
    
    if ( present(N) ) then
       call grib_get(F%igrib,'N',N)
    end if

    if ( present(gg) ) then
       ! Get and check size of field
       call grib_get_size(F%igrib,'values',numberOfValues)
       
       if ( size(gg) /= numberOfValues ) then
          write (gol,'("gg grid sizes do not match:")')     ; call goErr
          write (gol,'("  expected : ",i6)') size(gg)       ; call goErr
          write (gol,'("  found    : ",i6)') numberOfValues ; call goErr
          TRACEBACK; status=1; return
       end if

       ! Get the field
       call grib_get(F%igrib,'values',gg)
       
    end if
    
    ! -- spectral coef
    
    if ( present(T) ) then      ! truncation nb
       call grib_get(F%igrib,'M',T)
    end if

    if ( present(sh) ) then
       
       ! Get and check size of field
       call grib_get_size(F%igrib,'values',numberOfValues)
       
       if ( size(sh) /= numberOfValues/2 ) then
          write (gol,'("numbers of spectral coefficients do not match:")') ; call goErr
          write (gol,'("  expected : ",i6)') size(sh)                      ; call goErr
          write (gol,'("  found    : ",i6)') numberOfValues/2              ; call goErr
          TRACEBACK; status=1; return
       end if
       
       ! Get the field, convert from real to complex       
       allocate(values(numberOfValues), stat=iret)       
       call grib_get(F%igrib,'values',values)

       sh = transfer( values, (/(0.0,0.0)/) )
       deallocate( values )
       
    end if

    ! ok
    status = 0

  end subroutine gribfile_Get


  ! ===


  ! NOTE
  !   Do not check hybride coeff; different number of bits etc
  !   cause small differences ...
  
  subroutine gribfile_Check( F, status, debug, &
                             model_id, pid, &
                             levtype, level, &
                             reftime, timerange, &
                             gridtype, &
                             lon_first, lon_last, lon_inc, lon_n, &
                             lat_first, lat_last, lat_inc, lat_n, &
                             T  )
                
    use Grib_API,       only : grib_get
    use file_grib_base, only : gridtype_ll, gridtype_gg, gridtype_sh, gridtype_red_gg

    ! --- in/out -------------------------

    type(TGribFile), intent(in)      ::  F
    integer, intent(out)             ::  status

    integer, intent(in), optional    ::  debug


    integer, intent(in), optional    ::  model_id
    integer, intent(in), optional    ::  pid
    integer, intent(in), optional    ::  levtype, level
    !real, intent(in), optional       ::  hyb_a(:), hyb_b(:)      
    integer, intent(in), optional    ::  reftime(5)
    integer, intent(in), optional    ::  timerange(4)
    integer, intent(in), optional    ::  gridtype
    integer, intent(in), optional    ::  lon_first, lon_last, lon_inc, lon_n
    integer, intent(in), optional    ::  lat_first, lat_last, lat_inc, lat_n
    integer, intent(in), optional    ::  T

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/gribfile_Check'
    
    ! --- local --------------------------
    
    integer                ::  dataDate, dataTime
    integer                ::  grib_model_id
    integer                ::  grib_pid
    integer                ::  grib_levtype, grib_level
    integer                ::  grib_reftime(5)
    integer                ::  grib_timerange(4)
    integer                ::  grib_n
!    integer                ::  grib_T
!    
!    !real, allocatable      ::  grib_hyb_a(:), grib_hyb_b(:)

    logical                ::  verbose

    integer                ::  step, recGridType
    character(255)         ::  gridtypestring
    real                   ::  xdum

    ! --- begin --------------------------

    ! write error messages ?    
    verbose = .false.
    if ( present(debug) ) verbose = debug > 0

    ! ok by default
    status = 0

    ! -- check model id

    if ( present( model_id ) ) then
       call grib_get(F%igrib, 'generatingProcessIdentifier', grib_model_id )
       if ( grib_model_id /= model_id ) then
          if ( verbose ) then
             write (gol,'("model id''s do not match:")')      ; call goErr
             write (gol,'("  expected : ",i4)') model_id      ; call goErr
             write (gol,'("  found    : ",i4)') grib_model_id ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)          ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if


    ! -- check parameter

    if ( present( pid ) ) then
       call grib_get(F%igrib, 'paramId', grib_pid )
       if ( grib_pid /= pid ) then
          if ( verbose ) then
             write (gol,'("parameter id''s do not match:")') ; call goErr
             write (gol,'("  expected : ",i6)') pid          ; call goErr
             write (gol,'("  found    : ",i6)') grib_pid     ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)         ; call goErr
             TRACEBACK
          end if
          status = status+1
       end if
    end if


    ! -- check level

    if ( present(levtype) ) then
       call grib_get(F%igrib, 'levtype', grib_levtype )
       if ( grib_levtype /= levtype ) then
          if ( verbose ) then
             write (gol,'("level types do not match:")')     ; call goErr
             write (gol,'("  expected : ",i6)') levtype      ; call goErr
             write (gol,'("  found    : ",i6)') grib_levtype ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)         ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if
    
    if ( present(level) ) then
       call grib_get(F%igrib, 'level', grib_level )
       if ( grib_level /= level ) then
          if ( verbose ) then
             write (gol,'("levels do not match:")')        ; call goErr
             write (gol,'("  expected : ",i6)') level      ; call goErr
             write (gol,'("  found    : ",i6)') grib_level ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)       ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if
    
    !if ( present(hyb_a) .or. present(hyb_b) ) then
    !  if ( .not. (present(hyb_a) .and. present(hyb_b)) ) then
    !    print *, 'gribsections_Extract : error : extract both hybride params or none'
    !    !stop
    !  end if
    !  allocate( grib_hyb_a(size(hyb_a)) )
    !  allocate( grib_hyb_b(size(hyb_b)) )
    !  call Get( F, hyb_a=grib_hyb_a, hyb_b=grib_hyb_b )
    !  if ( any(grib_hyb_a/=hyb_a) .or. any(grib_hyb_b/=hyb_b) ) then
    !    if ( present(status) ) then
    !      status = 1
    !    else
    !      print *, 'gribsections_Check : hybride coeff in grib file do not match'
    !      stop
    !    end if
    !  end if
    !  deallocate( grib_hyb_a )
    !  deallocate( grib_hyb_b )
    !end if
    
    ! -- check reference time

    ! extract reference time:
    call grib_get(F%igrib,'dataDate',dataDate)
    call grib_get(F%igrib,'dataTime',dataTime)
    grib_reftime(1) = dataDate/10000           ! year
    grib_reftime(2) = Mod(dataDate/100,100)    ! month
    grib_reftime(3) = Mod(datadate,100)        ! day
    grib_reftime(4) = dataTime/100             ! hour
    grib_reftime(5) = Mod(dataTime,100)        ! minutes
    if ( present(reftime) ) then
       if ( any( grib_reftime /= reftime ) ) then
          if ( verbose ) then
             write (gol,'("reference times do not match:")')     ; call goErr
             write (gol,'("  expected : ",i4,4i3)') reftime      ; call goErr
             write (gol,'("  found    : ",i4,4i3)') grib_reftime ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)             ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if

    ! -- check time range

    if ( present(timerange) ) then
       write (*,'(a,"  WARNING timerange argument not fully supported ")') rname
       write(*,*) 
       ! Available from Grib-api: stepType, stepUnits, startStep, endStep, stepRange, step
       call grib_get(F%igrib,'step',step)
       grib_timerange(1:4) = (/ 1, step, 0, 0 /)
       if ( any( grib_timerange /= timerange ) ) then
          if ( verbose ) then
             write (gol,'("time ranges do not match:")')        ; call goErr
             write (gol,'("  expected : ",4i3)') timerange      ; call goErr
             write (gol,'("  found    : ",4i3)') grib_timerange ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)            ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if
  
    ! --- check grid
    
    if ( present(gridtype) ) then
       call grib_get(F%igrib,'gridType',gridtypestring)

       ! See https://software.ecmwf.int/wiki/display/GRIB/Grib+API+keys for availability
       
       select case ( trim(gridtypestring) )
       case ( 'regular_ll' ) ; recGridType = gridtype_ll
       case ( 'regular_gg' ) ; recGridType = gridtype_gg
       case ( 'regular_sh' ) ; recGridType = gridtype_sh
       case ( 'sh'         ) ; recGridType = gridtype_sh     ! used by grid_api
       case ( 'reduced_gg' ) ; recgridtype = gridtype_red_gg ! used by grid_api
       case default
          if ( verbose ) then
             write (gol,'("unsupported grid type string : ",a)') trim(gridtypestring); call goErr
          endif
          recGridType = -999
       end select

       if ( recGridType /= gridtype ) then
          if ( verbose ) then
             write (gol,'("gridtypes do not match:")')      ; call goErr
             write (gol,'("  expected : ",i6)') gridtype    ; call goErr
             write (gol,'("  found    : ",i6)') recGridType ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)        ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if

    if ( present(lon_n) ) then
       call grib_get(F%igrib,'Ni',grib_n)
       if ( grib_n /= lon_n ) then
          if ( verbose ) then
             write (gol,'("number of longitudes do not match:")') ; call goErr
             write (gol,'("  expected : ",i6)') lon_n             ; call goErr
             write (gol,'("  found    : ",i6)') grib_n            ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)              ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if

    if ( present(lat_n) ) then
       call grib_get(F%igrib,'Nj',grib_n)
       if ( grib_n /= lat_n ) then
          if ( verbose ) then
             write (gol,'("number of latitudes do not match:")') ; call goErr
             write (gol,'("  expected : ",i6)') lat_n            ; call goErr
             write (gol,'("  found    : ",i6)') grib_n           ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)             ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if
    
    if ( present(lat_first) ) then
       call grib_get(F%igrib,'latitudeOfFirstGridPointInDegrees',xdum)
       grib_n = nint(xdum*1000)
       if ( grib_n /= lat_first ) then
          if ( verbose ) then
             write (gol,'("first latitudes do not match:")') ; call goErr
             write (gol,'("  expected : ",i6)') lat_first    ; call goErr
             write (gol,'("  found    : ",i6)') grib_n       ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)         ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if

    if ( present(lon_first) ) then
       call grib_get(F%igrib,'longitudeOfFirstGridPointInDegrees',xdum)
       grib_n = nint(xdum*1000)
       if ( grib_n /= lon_first ) then
          if ( verbose ) then
             write (gol,'("first longitudes do not match:")') ; call goErr
             write (gol,'("  expected : ",i6)') lon_first     ; call goErr
             write (gol,'("  found    : ",i6)') grib_n        ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)          ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if

    if ( present(lat_last) ) then
       call grib_get(F%igrib,'latitudeOfLastGridPointInDegrees',xdum)
       grib_n = nint(xdum*1000)
       if ( grib_n /= lat_last ) then
          if ( verbose ) then
             write (gol,'("last latitudes do not match:")') ; call goErr
             write (gol,'("  expected : ",i6)') lat_last    ; call goErr
             write (gol,'("  found    : ",i6)') grib_n      ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)        ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if

    if ( present(lon_last) ) then
       call grib_get(F%igrib,'longitudeOfLastGridPointInDegrees',xdum)
       grib_n = nint(xdum*1000)
       if ( grib_n /= lon_last ) then
          if ( verbose ) then
             write (gol,'("last longitudes do not match:")') ; call goErr
             write (gol,'("  expected : ",i6)') lon_last     ; call goErr
             write (gol,'("  found    : ",i6)') grib_n       ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)         ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if

    if ( present(lon_inc) ) then
       call grib_get(F%igrib,'iDirectionIncrementInDegrees',xdum)
       grib_n = nint(xdum*1000)
       if ( grib_n /= lon_inc ) then
          if ( verbose ) then
             write (gol,'("longitude increments do not match:")') ; call goErr
             write (gol,'("  expected : ",i6)') lon_inc           ; call goErr
             write (gol,'("  found    : ",i6)') grib_n            ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)              ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if

    if ( present(lat_inc) ) then
       call grib_get(F%igrib,'jDirectionIncrementInDegrees',xdum)
       grib_n = nint(xdum*1000)
       if ( grib_n /= lat_inc ) then
          if ( verbose ) then
             write (gol,'("latitude increments do not match:")') ; call goErr
             write (gol,'("  expected : ",i6)') lat_inc          ; call goErr
             write (gol,'("  found    : ",i6)') grib_n           ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)             ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if

    ! --- check truncation
    
    if ( present(T) ) then
       write (*,'("ERROR, Spectral grid not yet supported for Grib-api")')
       stop
    end if

  end subroutine gribfile_Check


end module file_grib_api