ecearth3/sources/tm5mp/base/file_grib_gribex.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
!
!#######################################################################

module file_grib_gribex

  use GO, only : gol, goErr, goPr
  
  implicit none

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

  private

  public   ::  TGribFile
  public   ::  Init, Done, Opened
  
  public   ::  Clear
  public   ::  CopySections, CopyAllSections
  public   ::  Get
  public   ::  Check
  public   ::  Set
  
  public   ::  ReadRecord
  
!  public   ::  pidname, pidmax
!  public   ::  pid_T
!  public   ::  pid_Q
!  public   ::  pid_U, pid_V, pid_W
!  public   ::  pid_SP, pid_LNSP
!  public   ::  pid_VO, pid_D
!  public   ::  pid_Z, pid_LSM
!  public   ::  pid_SR, pid_AL, pid_LSRH
!  public   ::  pid_SLHF
!  public   ::  pid_CLWC, pid_CIWC, pid_CC

  ! --- const ------------------------------
  
  character(len=*), parameter  ::  mname = 'file_grib_gribex'
  
  
!  ! *** Type of level (Code Table 3)
!  
!  integer, parameter  ::  levtype_sfc  = 1    ! ground or water surface 
!  integer, parameter  ::  levtype_hyb  = 109  ! hybrid level 
!  integer, parameter  ::  levtype_land = 112  ! layer below surface
!  
!  ! *** parameter identifier (Code Table 2)
!
!  ! o numbers
!
!  integer, parameter  ::  pid_Z    = 129  ! geopotential (orography)
!  integer, parameter  ::  pid_T    = 130  ! temperature
!  integer, parameter  ::  pid_U    = 131  ! u-velocity
!  integer, parameter  ::  pid_V    = 132  ! v-velocity
!  integer, parameter  ::  pid_Q    = 133  ! specific humidity
!  integer, parameter  ::  pid_SP   = 134  ! surface pressure
!  integer, parameter  ::  pid_W    = 135  ! vertical velocity
! 
!  integer, parameter  ::  pid_VO   = 138  ! vorticity
!
!  integer, parameter  ::  pid_SLHF = 147  ! surface latent heat flux  (W m**-2 s)
!
!  integer, parameter  ::  pid_LNSP = 152  ! ln surface pressure
! 
!  integer, parameter  ::  pid_D    = 155  ! divergence
!  
!  integer, parameter  ::  pid_LSM  = 172  ! land sea mask
!  integer, parameter  ::  pid_SR   = 173  ! Surface roughness                  m
!  integer, parameter  ::  pid_AL   = 174  ! Albedo                             (0-1)
!
!  integer, parameter  ::  pid_LSRH = 234  ! Logarithm of SR length for heat
!  
!  integer, parameter  ::  pid_CLWC = 246  ! cloud liquid water content
!  integer, parameter  ::  pid_CIWC = 247  ! cloud ice water content
!  integer, parameter  ::  pid_CC   = 248  ! cloud cover
!  
!  ! o maximum number
!  integer, parameter  ::  pidmax = 260
!  
!  ! o names
!  
!  character(len=4), parameter ::  pidname(pidmax) = (/ &
!    'a   ','b   ','c   ','d   ','e   ','f   ','g   ','p008','p009','p010', &
!    'p011','p012','p013','p014','p015','p016','p017','p018','p019','p020', &
!    'p021','p022','p023','p024','p025','p026','p027','p028','p029','p030', &
!    'p031','p032','p033','p034','p035','p036','p037','p038','p039','p040', &
!    'p041','p042','p043','p044','p045','p046','p047','p048','p049','p050', &
!    'p051','p052','p053','p054','p055','p056','p057','p058','p059','p060', &
!    'p061','p062','p063','p064','p065','p066','p067','p068','p069','p070', &
!    'p071','p072','p073','p074','p075','p076','p077','p078','p079','p080', &
!    'p081','p082','p083','p084','p085','p086','p087','p088','p089','p090', &
!    'p091','p092','p093','p094','p095','p096','p097','p098','p099','p100', &
!    'p101','p102','p103','p104','p105','p106','p107','p108','p109','p110', &
!    'p111','p112','p113','p114','p115','p116','p117','p118','p119','p120', &
!    'p121','p122','p123','p124','p125','p126','p127','p128','Z   ','T   ', &
!    'U   ','V   ','Q   ','SP  ','W   ','p136','p137','VO  ','p139','p140', &
!    'SD  ','LSP ','CP  ','SF  ','p145','SSHF','SLHF','p148','p149','p150', &
!    'p151','LNSP','p153','p154','D   ','zg  ','pw  ','pu  ','pv  ','sp  ', &
!    'eu  ','du  ','dk  ','ed  ','dd  ','p166','T2M ','D2M ','p169','p170', &
!    'p171','LSM ','SR  ','AL  ','p175','SSR ','p177','p178','p179','EWSS', &
!    'NSSS','p182','p183','p184','p185','p186','p187','p188','p189','clb ', &
!    'clt ','clfs','p193','p194','p195','p196','p197','SRC ','p199','p200', &
!    'p201','p202','p203','p204','p205','p206','p207','p208','p209','p210', &
!    'p211','p212','p213','p214','p215','p216','p217','p218','p219','p220', &
!    'p221','p222','p223','p224','p225','p226','p227','p228','p229','p230', &
!    'p231','p232','p233','LSRH','p235','p236','p237','p238','p239','p240', &
!    'p241','p242','p243','p244','p245','CLWC','CIWC','CC  ','cco ','ccu ', &
!    'p251','p252','p253','p254','p255','p256','p257','p258','p259','p260' /)
!
!  ! *** data representation type (Code Table 6)
!  
!  integer, parameter  ::  gridtype_ll = 0
!  integer, parameter  ::  gridtype_gg = 4
!  integer, parameter  ::  gridtype_sh = 50
  
  ! --- types -------------------------------

  type TGribFile
    integer                ::  fu
    character(len=256)     ::  fname

    ! *** buffer for grib message
    integer                ::  inbuff_bytes     ! length of input buffer
    ! Section 0 - Indicator Section.       
    integer                ::  isec0(2)
    ! Section 1 - Product Definition Section.
    integer                ::  isec1(512)
    ! Section 2 - Grid Description Section.
    integer                ::  isec2(512)
    real                   ::  rsec2(512)
    ! Section 3 - ??
    integer                ::  isec3(2)
    real                   ::  rsec3(2)
    ! Section 4 - Binary Data  Section.
    integer                ::  isec4(512)
    real, pointer          ::  rsec4(:)
    ! 
    integer                ::  kword

  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 Clear
    module procedure gribfile_Clear
  end interface
  
  interface CopySections
    module procedure gribfile_CopySections
  end interface
  
  interface CopyAllSections
    module procedure gribfile_CopyAllSections
  end interface
  
  interface Get
    module procedure gribfile_Get
  end interface

  interface Check
    module procedure gribfile_Check
  end interface
  
  interface Set
    module procedure gribfile_Set
  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 )

    ! --- 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'

    ! --- local ------------------------

    ! dummy buffer of insufficient length:
    integer        ::  inbuff(10)

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

    ! *** initialize grib sections

    ! Section 0 - Indicator Section.       
    F%isec0 = 0

    ! Section 1 - Product Definition Section.
    F%isec1 = 0

    ! Section 2 - Grid Description Section.
    F%isec2 = 0
    F%rsec2 = 0.0
    ! Section 3 - ??
    F%isec3 = 0
    F%rsec3 = 0.0

    ! Section 4 - Binary Data  Section.
    F%isec4 = 0
    nullify( F%rsec4 )

    F%kword = 0

    ! *** Open specified file with requested name.
    
    ! A free file unit seems to be assigned in pbOpen.
    call pbOpen( F%fu, file, mode, status )
    if ( status /= 0 ) then
      write (gol,'("from pbOpen:")'); call goErr
      select case ( status )
        case (-1) 
          write (gol,'("  could not open file")'); call goErr
        case (-2) 
          write (gol,'("  invalid file name")') ; call goErr
        case (-3) 
          write (gol,'("  invalid open mode")'); call goErr
        case default
          write (gol,'("  unknown return status: ",i6)') status; call goErr
      end select
      write (gol,'("  file: ",a)') trim(file); call goErr
      TRACEBACK; status=1; return
    end if

    ! save filename
    F%fname = file

    select case ( mode )
    
    case ( 'r' )
    
       ! Determine length of records ('messages' in grib context);
       ! open with insufficient buffer, and check error messages:
       
        status = 1   ! force return even if error occures
        call pbGrib( F%fu, inbuff, size(inbuff)*kind(inbuff), F%inbuff_bytes, status )
        select case ( status )
          case ( 0)
            ! no error ?
            write (gol,'("tried to read record into buffer with insuficient length,")'); call goErr
            write (gol,'("but no error occured ...")'); call goErr
            TRACEBACK; status=1; return
          case (-1)
            write (gol,'("end-of-file is hit before a GRIB product is read")'); call goErr
            TRACEBACK; status=1; return
          case (-3)
            ! size of input buffer was not sufficient (as expected)
          case default
            write (gol,'("unknown return status from pbGrib : ",i6)') status; call goErr
            TRACEBACK; status=1; return
        end select
        ! adhoc increment (sometimes first record is very different form others)
        F%inbuff_bytes = F%inbuff_bytes*2
        ! Close file:
        call pbClose( F%fu, status )
        if ( status /= 0 ) then
          write (gol,'("error in call to pbClose; status=",i6)') status; call goErr
          TRACEBACK; status=1; return
        end if

        ! Reopen file:
        call pbOpen( F%fu, file, mode, status )
        if ( status /= 0 ) then
          write (gol,'("from second open; status=",i6)') status; call goErr
          TRACEBACK; status=1; return
        end if
        
      case ( 'w' )
      
        ! nothing special
        
      case default
      
        write (gol,'(" unknown mode `",a,"`")') mode; call goErr
        TRACEBACK; status=1; return
        
    end select


  end subroutine gribfile_Init


  ! ===


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

  subroutine gribfile_Done( F, status )

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

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

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

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

    ! close file:
    call pbClose( 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
    end if
    
    ! clear section 4 :
    if (associated(F%rsec4)) deallocate( F%rsec4 )
    
  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 )

  end function gribfile_Opened
  

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


  subroutine gribfile_CopySections( gribIn, gribOut, status )
  
    ! ---in/out ------------------------------
    
    type(TGribFile), intent(in)     ::  gribIn
    type(TGribFile), intent(out)    ::  gribOut
    integer,         intent(inout)  ::  status

    ! --- begin -------------------------
    
    !gribOut%isec0 = gribIn%isec0

    gribOut%isec1 = gribIn%isec1 

    gribOut%isec2 = gribIn%isec2
    gribOut%rsec2 = gribIn%rsec2 

    gribOut%isec3 = gribIn%isec3

    gribOut%isec4 = gribIn%isec4
    
    ! ok
    status = 0

  end subroutine gribfile_CopySections


  ! ***


  subroutine gribfile_CopyAllSections( gribIn, gribOut )
  
    ! ---in/out ------------------------------
    
    type(TGribFile), intent(in)     ::  gribIn
    type(TGribFile), intent(out)    ::  gribOut
    integer,         intent(inout)  ::  status
    
    ! --- begin -------------------------
    
    !gribOut%isec0 = gribIn%isec0

    gribOut%isec1 = gribIn%isec1 

    gribOut%isec2 = gribIn%isec2
    gribOut%rsec2 = gribIn%rsec2 

    gribOut%isec3 = gribIn%isec3

    gribOut%isec4 = gribIn%isec4
    if ( associated(gribOut%rsec4) ) then
      if ( size(gribOut%rsec4) /= size(gribIn%rsec4) ) then
        deallocate( gribOut%rsec4 )
        allocate( gribOut%rsec4(size(gribIn%rsec4)) )
      end if
    else
      allocate( gribOut%rsec4(size(gribIn%rsec4)) )
    end if
    gribOut%rsec4 = gribIn%rsec4
    
    ! ok
    status = 0

  end subroutine gribfile_CopyAllSections


  ! ***


  subroutine gribfile_Clear( F )

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

    type(TGribFile), intent(out)    ::  F
    integer,         intent(inout)  ::  status
    
    ! --- begin ------------------------
    
    ! Section 0 - Indicator Section.       
    F%isec0 = 0

    ! Section 1 - Product Definition Section.
    F%isec1 = 0

    ! Section 2 - Grid Description Section.
    F%isec2 = 0
    F%rsec2 = 0.0

    ! Section 3 - ??
    F%isec3 = 0
    F%rsec3 = 0.0

    ! Section 4 - Binary Data  Section.
    F%isec4 = 0
    if ( associated(F%rsec4) ) F%rsec4 = 0.0

    ! default number of bits
    F%isec4(2) = 24

    F%kword = 0

    ! ok
    status = 0

  end subroutine gribfile_Clear


  ! ==================================================================
  
  ! USAGE
  !   call ReadRecord( gribfile [,debug=0|1] [,status=status] )
  !
  ! DESCRIPTION
  !   Reads next grib record 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 )

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

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

    ! --- local --------------------------
    
    integer, allocatable   ::  inbuff(:)
    integer                ::  length, len4
    
    ! --- begin --------------------------

    ! allocate array to store packed data:
    allocate( inbuff(ceiling(F%inbuff_bytes*1.0/kind(inbuff))) )

    ! read record:
    call pbGrib( F%fu, inbuff, size(inbuff)*kind(inbuff), length, status )
    select case ( status )
      case ( 0)
        ! no error !
      case (-1)
         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
      case (-3)
         ! size of input buffer is not sufficient for the GRIB product !
         write (gol,'("size of input buffer is not sufficient for the GRIB product:")') ; call goErr
         write (gol,'("  current size : ",i8)') size(inbuff)*kind(inbuff)               ; call goErr
         write (gol,'("  required     : ",i8)') length                                  ; call goErr
         write (gol,'("  grib file    : ",a)') F%fname                                  ; call goErr
         TRACEBACK; status=2; return
      case default
         write (gol,'("unknown return status from pbGrib : ",i6)') status ; call goErr
         write (gol,'("  grib file : ",a)') F%fname                       ; call goErr
         write (gol,'("in ",a)') rname                                    ; call goErr
         TRACEBACK; status=2; return
    end select

    ! Decode.
    ! The array to hold the unpacked data may need to be 4 times 
    ! as long as that for the packed data (or more!).
    len4 = ceiling(4*F%inbuff_bytes*1.0/kind(F%rsec4))
    if ( associated(F%rsec4) ) then
      if ( size(F%rsec4) < len4 ) then
        deallocate( F%rsec4 )
        allocate( F%rsec4(len4) )
      end if
    else
      allocate( F%rsec4(len4) )
    end if
    call GribEx(F%isec0,F%isec1,F%isec2,F%rsec2,F%isec3,F%rsec3,&
                F%isec4, F%rsec4, size(F%rsec4)*kind(F%rsec4), &
                inbuff, size(inbuff), F%kword, 'D' , status ) 
    if ( status /= 0 ) then
       select case ( status )
       case (-6)
          write (gol,'("in call to GribEx : ")'); call goErr
          write (gol,'("  found pseudo data in grib file: ",a)') F%fname; call goErr
          TRACEBACK; status=2; return
       case default
          write (gol,'("unknown return status from GribEx : ",i6)') status; call goErr
          TRACEBACK; status=2; return
       end select
    end if
    
    ! clear input buffer:
    deallocate( inbuff )
    
    ! ok
    status = 0
    
  end subroutine gribfile_ReadRecord


  ! ***
  
  
  ! encode and write grib sections

  subroutine gribfile_WriteRecord( F, status )
    
    ! --- in/out --------------------
    
    type(TGribFile), intent(in)        ::  F
    integer, intent(out)               ::  status

    ! --- const ---------------------
    
    character(len=*), parameter  ::  rname = mname//'/gribfile_WriteRecord'
    
    ! --- local ---------------------
    
    integer, allocatable   ::  outbuff(:)
        
    ! --- begin ---------------------

    ! *** encode
    
    status = 1
    allocate( outbuff(size(F%rsec4)) )
    call GribEx(F%isec0,F%isec1,F%isec2,F%rsec2,F%isec3,F%rsec3,&
                F%isec4, F%rsec4, size(F%rsec4)*kind(F%rsec4), &
                outbuff, size(outbuff), F%kword, 'C' , status ) 
    select case ( status )
    case (-2)
       write (gol,'("WARNING: bitmap was encountered with all bits set to 1")'); call goPr
    case (-3)
       write (gol,'("WARNING: predefined bitmap was encountered.")'); call goPr
    case (-4)
       write (gol,'("WARNING: bitmap was encountered.")'); call goPr
    case (-5)
       write (gol,'("WARNING: bitmap was encountered and the data has not been decoded.")'); call goPr
    case (-6)
       write (gol,'("WARNING: ECMWF pseudo-GRIB data (TIDE or BUDG) has been encountered.")'); call goPr
    case ( 0 )
       ! ok
    case default
       write (gol,'("unknown status from GribEx : ",i6)') status
       TRACEBACK; status=1; return
    end select

    ! *** write GRIB record

    ! PBWRITE: Writes a block of bytes to a file.
    !   status :    -1 = Could not write to file.
    !            >=  0 = Number of bytes written.
    !
    call pbWrite( F%fu, outbuff, F%isec0(1), status )
    if ( status == -1 ) then
       write (gol,'("could not write to file")'); call goErr
       write (gol,'("in ",a)') rname; status=1; return
    else if ( status >= 0 ) then
       ! bytes written correctly
    else
       write (gol,'("unknown status from GribEx : ",i6)') status; call goErr
       write (gol,'("in ",a)') rname; status=1; return
    end if
    
    ! done
    deallocate( outbuff )
    
    ! ok
    status = 0

  end subroutine gribfile_WriteRecord
  
  ! ============================================================

  ! 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 file_grib_base, only : pidName
    use file_grib_base, only : gridtype_ll, gridtype_gg, gridtype_sh

    ! --- 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  ::  name = mname//'/gribfile_Get'
    
    ! --- local --------------------------
    
    integer                ::  grib_reftime(5)
    integer                ::  grib_timerange(4)
    integer                ::  base
    integer                ::  grib_n
    real, allocatable      ::  rsec4_adhoc(:)

    ! --- begin --------------------------
    
    ! -- model id

    if ( present( model_id ) ) then
       model_id = F%isec1(3)
    end if

    ! --  parameter

    if ( present( pid ) ) then
       pid = F%isec1(6)
    end if

    ! --  parameter short name

    if ( present( pidShortName ) ) then
       pidShortName = pidName(F%isec1(6))
    end if

    ! --  level

    if ( present(levtype) ) levtype = F%isec1(7)
    if ( present(level  ) ) level   = F%isec1(8)

    ! -- 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 hybride params or none")'); call goErr
          TRACEBACK; status=1; return
       end if

       if ( F%isec2(12) /= 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)') F%isec2(12)                ; call goErr
          TRACEBACK; status=1; return
       end if
       base = 10            ;  hyb_a = F%rsec2(base+1:base+size(hyb_a))
       base = 10+size(hyb_a);  hyb_b = F%rsec2(base+1:base+size(hyb_b))

    end if
    
    ! number of levels

    if ( present(nlev) ) then
       nlev = F%isec2(12)/2 - 1
    end if

    ! --  reference time

    if ( present(reftime) ) then
       ! extract reference time:
       ! extract reference time:
       grib_reftime(1) = (F%isec1(21)-1)*100 + F%isec1(10)   ! year incl. century
       grib_reftime(2) = F%isec1(11)            ! month
       grib_reftime(3) = F%isec1(12)            ! day
       grib_reftime(4) = F%isec1(13)            ! hour
       grib_reftime(5) = F%isec1(14)            ! minutes
       reftime = grib_reftime
    end if

    ! --  time range

    if ( present(timerange) ) then
       grib_timerange(1:4) = F%isec1(15:18)
       timerange = grib_timerange
    end if

    ! ---  grid
    
    if ( present(gridtype)  ) gridtype  = F%isec2( 1)
    if ( present(lon_n)     ) lon_n     = F%isec2( 2)
    if ( present(lat_n)     ) lat_n     = F%isec2( 3)
    if ( present(lat_first) ) lat_first = F%isec2( 4) / 1000.0
    if ( present(lon_first) ) lon_first = F%isec2( 5) / 1000.0
    if ( present(lat_last)  ) lat_last  = F%isec2( 7) / 1000.0
    if ( present(lon_last)  ) lon_last  = F%isec2( 8) / 1000.0
    if ( present(lon_inc)   ) lon_inc   = F%isec2( 9) / 1000.0
    if ( present(lat_inc)   ) lat_inc   = F%isec2(10) / 1000.0

    ! -- lat/lon grid
    
    if ( present(ll) ) then
       !call Check( F, gridtype=gridtype_ll )
       ! >>> no recursive call >>>>>>>>>>>>>>
       grib_n = F%isec2(1)
       if ( grib_n /= gridtype_ll ) then
          write (gol,'("grid types do not match:")')             ; call goErr
          write (gol,'("  expected : ",i4," (ll)")') gridtype_ll ; call goErr
          write (gol,'("  found    : ",i4)') grib_n              ; call goErr
          TRACEBACK; status=1; return
       end if
       ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       if ( (size(ll,1)/=F%isec2(2)) .or. (size(ll,2)/=F%isec2(3)) ) 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)') F%isec2(2), F%isec2(3) ; call goErr
          TRACEBACK; status=1; return
       end if
       if ( .not. associated(F%rsec4) ) then
          write (gol,'("no grib section 4 read")'); call goErr
          TRACEBACK; status=1; return
       end if
       ll = reshape( F%rsec4(1:size(ll)), shape(ll) )
    end if

    ! -- gaussian lat/lon grid
    
    if ( present(N) ) then
       !call Check( F, gridtype=gridtype_gg )
       ! >>> no recursive call >>>>>>>>>>>>>>
       grib_n = F%isec2(1)
       if ( grib_n /= gridtype_gg ) then
          write (gol,'("grid types do not match:")')             ; call goErr
          write (gol,'("  expected : ",i4," (gg)")') gridtype_gg ; call goErr
          write (gol,'("  found    : ",i4)') grib_n              ; call goErr
          TRACEBACK; status=1; return
       end if
       ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       N = F%isec2(10)
    end if

    if ( present(gg) ) then
       !call Check( F, gridtype=gridtype_gg )
       ! >>> no recursive call >>>>>>>>>>>>>>
       grib_n = F%isec2(1)
       if ( grib_n /= gridtype_gg ) then
          write (gol,'("grid typess do not match:")')            ; call goErr
          write (gol,'("  expected : ",i4," (gg)")') gridtype_gg ; call goErr
          write (gol,'("  found    : ",i4)') grib_n              ; call goErr
          TRACEBACK
          status=1; return
       end if
       ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       if ( size(gg) /= F%isec4(1) ) then
          write (gol,'("gg grid sizes do not match:")') ; call goErr
          write (gol,'("  expected : ",i4)') size(gg)   ; call goErr
          write (gol,'("  found    : ",i4)') F%isec4(1) ; call goErr
          TRACEBACK; status=1; return
       end if
       ! data present ?
       if ( .not. associated(F%rsec4) ) then
          write (gol,'("no grib section 4 read")'); call goErr
          TRACEBACK; status=1; return
       end if
       gg = F%rsec4(1:size(gg))
    end if
    
    ! -- spectral coef
    
    if ( present(T) ) then
       !call Check( F, gridtype=gridtype_sh )
       ! >>> no recursive call >>>>>>>>>>>>>>
       grib_n = F%isec2(1)
       if ( grib_n /= gridtype_sh ) then
          write (gol,'("grid types do not match:")')             ; call goErr
          write (gol,'("  expected : ",i4," (sh)")') gridtype_sh ; call goErr
          write (gol,'("  found    : ",i4)') grib_n              ; call goErr
          TRACEBACK; status=1; return
       end if
       ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       T = F%isec2(4)
    end if

    if ( present(sh) ) then
       !call Check( F, gridtype=gridtype_sh )
       ! >>> no recursive call >>>>>>>>>>>>>>
       grib_n = F%isec2(1)
       if ( grib_n /= gridtype_sh ) then
          write (gol,'("grid typess do not match:")')            ; call goErr
          write (gol,'("  expected : ",i4," (sh)")') gridtype_sh ; call goErr
          write (gol,'("  found    : ",i4)') grib_n              ; call goErr
          TRACEBACK; status=1; return
       end if
       ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       ! check number of elements; complex is 2 elements !
       if ( size(sh) /= F%isec4(1)/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)') F%isec4(1)/2                  ; call goErr
          TRACEBACK; status=1; return
       end if
       ! convert from real to complex:
       ! >>> NOTE: this call lets 'f90 -C' fail on SUN compiler !
       !sh = transfer( F%rsec4(1:F%isec4(1)), (/(0.0,0.0)/) )
       ! <<< fix >>>
       allocate( rsec4_adhoc(F%isec4(1)) )
       rsec4_adhoc = F%rsec4(1:F%isec4(1))
       sh = transfer( rsec4_adhoc, (/(0.0,0.0)/) )
       deallocate( rsec4_adhoc )
       ! <<<
    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  )
                
    ! --- 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  ::  name = mname//'/gribfile_Check'
    
    ! --- local --------------------------

    logical                ::  match

    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
    
    logical                ::  verbose

    ! --- 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 Get( F, status, model_id=grib_model_id )
       if (status/=0) then
          TRACEBACK; status=1; return
       end if
       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 Get( F, status, pid=grib_pid )
       if (status/=0) then
          TRACEBACK ; status=1; return
       end if
       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 Get( F, status, levtype=grib_levtype )
       if (status/=0) then
          TRACEBACK ; status=1; return
       end if
       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 Get( F, status, level=grib_level )
       if (status/=0) then
          TRACEBACK ; status=1; return
       end if
       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

    ! -- check reference time

    ! extract reference time:
    grib_reftime(1) = (F%isec1(21)-1)*100 + F%isec1(10)   ! year incl. century
    grib_reftime(2) = F%isec1(11)            ! month
    grib_reftime(3) = F%isec1(12)            ! day
    grib_reftime(4) = F%isec1(13)            ! hour
    grib_reftime(5) = F%isec1(14)            ! 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

    grib_timerange(1:4) = F%isec1(15:18)
    if ( present(timerange) ) then
       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

    grib_n = F%isec2(1)
    if ( present(gridtype) ) then
       if ( grib_n /= gridtype ) then
          if ( verbose ) then
             write (gol,'("gridtypes do not match:")')            ; call goErr
             write (gol,'("  expected : ",i6)') gridtype          ; 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

    grib_n = F%isec2(2)
    if ( present(lon_n) ) then
       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

    grib_n = F%isec2(3)
    if ( present(lat_n) ) then
       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

    grib_n = F%isec2(4)
    if ( present(lat_first) ) then
       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

    grib_n = F%isec2(5)
    if ( present(lon_first) ) then
       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

    grib_n = F%isec2(7)
    if ( present(lat_last) ) then
       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

    grib_n = F%isec2(8)
    if ( present(lon_last) ) then
       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

    grib_n = F%isec2(9)
    if ( present(lon_inc) ) then
       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

    grib_n = F%isec2(10)
    if ( present(lat_inc) ) then
       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 ; status=1; return
          end if
          status=status+1
       end if
    end if

    ! --- check truncation

    if ( present(T) ) then
       grib_T = F%isec2(4)
       if ( grib_T /= T ) then
          if ( verbose ) then
             write (gol,'("spectral truncations do not match:")') ; call goErr
             write (gol,'("  expected : ",i6)') T                 ; call goErr
             write (gol,'("  found    : ",i6)') grib_T            ; call goErr
             write (gol,'("  in ",a)') trim(F%fname)              ; call goErr
             TRACEBACK
          end if
          status=status+1
       end if
    end if

    ! return with status equal to number of failures

  end subroutine gribfile_Check


  ! ===
  
  
  ! fill gribsections

  subroutine gribfile_Set( F, status, model_id, pid, &
                           levtype, hyb_a, hyb_b, level, &
                           reftime, timerange, &
                           gridtype, &
                           ll, &
                           lon_first, lon_inc, lon_last, lon_n, &
                           lat_first, lat_inc, lat_last, lat_n, &
                           scanning_mode, nbits )
  
    ! --- in/out ----------------------
    
    type(TGribFile), intent(inout)       ::  F
    integer, intent(out)                 ::  status

    integer, intent(in), optional        ::  model_id
    
    integer, intent(in), optional        ::  pid

    integer, intent(in), optional        ::  levtype
    integer, intent(in), optional        ::  level
    real, intent(in), optional           ::  hyb_a(:), hyb_b(:)

    integer, intent(in), optional        ::  reftime(5)   ! (/yy,mm,dd,hh,min/)
    integer, intent(in), optional        ::  timerange(4)

    integer, intent(in), optional        ::  gridtype

    real, intent(in), optional           ::  ll(:,:)
    integer, intent(in), optional        ::  lon_n, lat_n
    integer, intent(in), optional        ::  lon_first, lon_inc, lon_last  ! mili degree
    integer, intent(in), optional        ::  lat_first, lat_inc, lat_last  ! mili degree
    integer, intent(in), optional        ::  scanning_mode

    integer, intent(in), optional        ::  nbits

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

    ! --- local ------------------------
    
    integer           :: base
    
    ! --- begin ------------------------

    ! Section 1 - Product Definition Section.
    ! ---------------------------------------
    
    ! Code Table 2 Version Number.
    F%isec1(1) = 128

    ! Originating centre identifier.
    F%isec1(2) = 98    ! ECMWF
    
    ! Model identification.
    if ( present(model_id) ) F%isec1(3) = model_id
    
    ! Grid definition.
    F%isec1(4) = 255   ! ???
    
    ! Flag (Code Table 1)                   10000000
    F%isec1(5) = 128   ! ???
    
    ! Parameter identifier (Code Table 2, NO LOCAL TABLE FOR KNMI)
    if ( present(pid) ) F%isec1(6) = pid
    
    ! Type of level (Code Table 3).              109   ; hybride
    ! Value 1 of level (Code Table 3).             1
    ! Value 2 of level (Code Table 3).             0
    if ( present(levtype) ) F%isec1(7) = levtype
    if ( present(level    ) ) F%isec1(8) = level
    F%isec1(9) = 0

    ! 10. Year of reference time of data.   
    ! 11. Month of reference time of data.  
    ! 12. Day of reference time of data.    
    ! 13. Hour of reference time of data.   
    ! 14. Minute of reference time of data. 
    ! 21. Century of reference time of data.
    if ( present(reftime) ) then
      F%isec1(21)    = int(reftime(1)/100.0) + 1
      F%isec1(10)    = mod(reftime(1),100)
      F%isec1(11:14) = reftime(2:5)
    end if

    ! Time unit (Code Table 4).
    ! Time range one.
    ! Time range two.
    ! Time range indicator (Code Table 5)
    if ( present(timerange) ) F%isec1(15:18) = timerange

    ! Number averaged.
    F%isec1(19) = 0
    
    ! Number missing from average.
    F%isec1(20) = 0
    

    ! Section 2 - Grid Description Section.
    ! -------------------------------------
    ! (Southern latitudes and Western longitudes are negative.)
    
    !  1. Data represent type (Table 6)         0 = lat/long
    if ( present(gridtype) ) F%isec2(1) = gridtype
    
    !  2. Number of points along a parallel.                 144
    !  3. Number of points along a meridian.                  72
    !  4. Latitude of first grid point.                   -88750  mdeg
    !  5. Longitude of first grid point.                 -180000  mdeg
    !  6. Resolution and components flag.               10000000
    !  7. Latitude of last grid point.                     88750
    !  8. Longitude of last grid point.                   177500
    !  9. i direction (East-West) increment.                2500
    ! 10. j direction (North-South) increment.              2500
    
    if ( present(lon_n) .or. present(lon_first) .or. present(lon_last) .or. present(lon_inc) ) then
    
      if ( present(lon_n)     ) F%isec2(2) = lon_n
      if ( present(lon_first) ) F%isec2(5) = lon_first
      if ( present(lon_last)  ) F%isec2(8) = lon_last
      if ( present(lon_inc)   ) F%isec2(9) = lon_inc
       
      if ( present(lon_first) .and. present(lon_last) .and. present(lon_inc) ) then
        ! lon_n
        F%isec2(2) = (lon_last-lon_first)/lon_inc + 1
      else if ( present(lon_n) .and. present(lon_first) .and. present(lon_inc) ) then
        ! lon_last
        F%isec2(8) = lon_first + lon_inc*(lon_n-1)
      else if ( present(lon_n) .and. present(lon_first) .and. present(lon_last) ) then
        ! lon_inc
        F%isec2(9) = (lon_last-lon_first)/lon_n
      else
         write (gol,'("expected at least 3 of lon_n,lon_first,lon_last,lon_inc")'); call goErr
         TRACEBACK; status=1; return
      end if
    
    end if
    
    if ( present(lat_n) .or. present(lat_first) .or. present(lat_last) .or. present(lat_inc) ) then
    
      if ( present(lat_n)     ) F%isec2(3)  = lat_n
      if ( present(lat_first) ) F%isec2(4)  = lat_first
      if ( present(lat_last)  ) F%isec2(7)  = lat_last
      if ( present(lat_inc)   ) F%isec2(10) = lat_inc
       
      if ( present(lat_first) .and. present(lat_last) .and. present(lat_inc) ) then
        ! lat_n
        F%isec2(3) = (lat_last-lat_first)/lat_inc + 1
      else if ( present(lat_n) .and. present(lat_first) .and. present(lat_inc) ) then
        ! lat_last
        F%isec2(7) = lat_first + lat_inc*(lat_n-1)
      else if ( present(lat_n) .and. present(lat_first) .and. present(lat_last) ) then
        ! lat_inc
        F%isec2(10) = (lat_last-lat_first)/lat_n
      else
         write (gol,'("expected at least 3 of lon_n,lon_first,lon_last,lon_inc")'); call goErr
         TRACEBACK; status=1; return
      end if
    
    end if
    
    F%isec2(6) = 128
  
    ! 11. Scanning mode flags (Code Table 8)
    !       bit 1 : 0 = west to east   ,    1 = east to west
    !       bit 2 : 0 = north to south ,    1 = south to north
    !       bit 3 : 0 = store lon rows ,    1 = store lat columns
    ! default scanning mode 0
    !F%isec2(11) = 0
    !if ( F%isec2(8) < F%isec2(5) ) F%isec2(11) = F%isec2(11) + 4  ! east to west
    !if ( F%isec2(7) > F%isec2(4) ) F%isec2(11) = F%isec2(11) + 2  ! south to north
    ! PROBLEM >>> GribEx seems to accept scanning mode 0 only ..
    F%isec2(11) = 0
    ! <<<
    if ( present(scanning_mode) ) F%isec2(11) = scanning_mode
    
    ! 12. Number of vertical coordinate parameters.
    ! Vertical Coordinate Parameters strored in rsec2, start in element 11
    if ( present(hyb_a) .and. present(hyb_b) ) then
      F%isec2(12) = size(hyb_a) + size(hyb_b)
      base = 10            ;  F%rsec2(base+1:base+size(hyb_a)) = hyb_a
      base = 10+size(hyb_a);  F%rsec2(base+1:base+size(hyb_b)) = hyb_b
    end if
    
    ! set rest to zero for safety
    F%isec2(13:size(F%isec2)) = 0
    

    ! Section 4 - Binary Data  Section.
    ! -------------------------------------

    !  1. Number of data values coded/decoded.             11556
    !  2. Number of bits per data value.                      16
    !
    !  3. Type of data       (0=grid pt, 128=spectral).      128
    !  4. Type of packing    (0=simple, 64=complex).          64
    !  5. Type of data       (0=float, 32=integer).            0
    !  6. Additional flags   (0=none, 16=present).             0
    !  7. Reserved.                                            0
    !  8. Number of values   (0=single, 64=matrix).            0
    !  9. Secondary bit-maps (0=none, 32=present).             0
    ! 10. Values width       (0=constant, 16=variable).        0
    ! 11. Byte offset of start of packed data (N).          2214
    ! 12. Power (P * 1000).                                  500
    ! 13. Pentagonal resolution parameter J for subset.       20
    ! 14. Pentagonal resolution parameter K for subset.       20
    ! 15. Pentagonal resolution parameter M for subset.       20

    
    ! set lat/lon grid:
    !  3. Type of data       (0=grid pt, 128=spectral).        0
    !  4. Type of packing    (0=simple, 64=complex).           0
    !  5. Type of data       (0=float, 32=integer).            0
    !  6. Additional flags   (0=none, 16=present).             0
    !  7. Reserved.                                            0
    !  8. Number of values   (0=single, 64=matrix).            0
    !  9. Secondary bit-maps (0=none, 32=present).             0
    ! 10. Values width       (0=constant, 16=variable).        0
    if ( present(ll) ) then
    
      ! check grid size:
      if ( F%isec2(2) /= size(ll,1) .or. F%isec2(3) /= size(ll,2) ) 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)') F%isec2(2), F%isec2(3); call goErr
         TRACEBACK; status=1; return
      end if

      !  1. Number of data values coded/decoded.             11556
      F%isec4(1) = size(ll)

      !  3. Type of data       (0=grid pt, 128=spectral).        0
      !  4. Type of packing    (0=simple, 64=complex).           0
      !  5. Type of data       (0=float, 32=integer).            0
      !  6. Additional flags   (0=none, 16=present).             0
      !  7. Reserved.                                            0
      !  8. Number of values   (0=single, 64=matrix).            0
      !  9. Secondary bit-maps (0=none, 32=present).             0
      ! 10. Values width       (0=constant, 16=variable).        0
      F%isec4(3:10) = 0

      ! allocate space to store field
      ! NOTE: allocate extra, sometimes errors if fit exately!   
      if ( associated(F%rsec4) ) then
        if ( size(F%rsec4) < size(ll) ) then
          deallocate( F%rsec4 )
          allocate( F%rsec4(size(ll)*2) )
        end if
      else
        allocate( F%rsec4(size(ll)*2) )
      end if
     
      ! convert to 1D array:
      F%rsec4(1:size(ll)) = reshape( ll, (/ size(ll) /) )
      
    end if
    
    !  2. Number of bits per data value.                      16
    if ( present(nbits) ) F%isec4(2) = nbits
    
    ! ok
    status = 0

  end subroutine gribfile_Set



end module file_grib_gribex