ecearth3/sources/tm5mp/base/tmm_mf_grib.F90

!###############################################################################
!
! Input/output of meteofiles : grib version.
!
!### macro's ###################################################################
!
#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
!
#include "tmm.inc"
!
!###############################################################################

module tmm_mf_grib

  use GO       , only : gol, goErr, goPr
  use GO       , only : TDate
  use file_grib, only : TGribFile

  implicit none
  
  ! --- in/out ----------------------------
  
  private
  
  public  ::  TMeteoFile_grib
  public  ::  Init, Done
  public  ::  Get
  public  ::  ReadRecord
  
  ! --- const ------------------------------
  
  character(len=*), parameter  ::  mname = 'tmm_mf_grib'
  

  !--- type ---------------------------------
  
  type TMeteoFile_grib
    ! file name:
    character(len=256)         ::  fname
    ! time range covered by file:
    type(TDate)                ::  trange(2)
    logical                    ::  constant
    ! other time keys for this file:
    character(len=16)          ::  treskey
    type(TDate)                ::  tday
    ! current time range covered by grib record:
    type(TDate)                ::  tref, t1, t2
    !
    ! file description
    !
    character(len=16)          ::  ec_class, ec_type
    character(len=256)         ::  paramkeys
    !
    ! 3d field
    !
    type(TGribFile)            ::  grib
    !
    ! surface pressure field
    !
    logical                    ::  do_spm
    character(len=256)         ::  spm_fname
    type(TGribFile)            ::  spm_grib
    logical                    ::  spm_lnsp2sp
    !    
  end type TMeteoFile_grib
  
  
  ! --- interfaces -------------------------
  
  interface Init
    module procedure mf_Init
  end interface

  interface Done
    module procedure mf_Done
  end interface

  interface Get
    module procedure mf_Get
  end interface

  interface ReadRecord
    module procedure mf_ReadRecord
  end interface


contains


  ! ==============================================================
  ! ====
  ! ==== original
  ! ====
  ! ==============================================================


  subroutine mf_Init( mf, dir, archivekeys, paramkey, &
                               tday, t1, t2, status )
    
    use GO, only : TDate, Get, NewDate, operator(>)
    use GO, only : goVarValue, goWriteKeyNum
    
    ! --- in/out ----------------------------
    
    type(TMeteoFile_grib), intent(out)  ::  mf
    character(len=*), intent(in)        ::  dir
    character(len=*), intent(in)        ::  archivekeys
    character(len=*), intent(in)        ::  paramkey
    type(TDate), intent(in)             ::  tday, t1, t2
    integer, intent(out)                ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_Init'
    
    ! --- local --------------------------------
    
    character(len=16)     ::  form
    character(len=16)     ::  ec_class, ec_type, ec_tres, ec_levs, mfile
    integer               ::  ec_sh, ec_gg
    character(len=16)     ::  grid, gridN, gridT
    character(len=16)     ::  levs
    character(len=16)     ::  spm_tres
    
    type(TDate)           ::  tfile
    integer               ::  ccyy, mm, dd
    type(TDate)           ::  tc
    
    logical               ::  exist
    
    ! --- begin --------------------------------

    !
    ! extract fields from archivekey :
    !   form=tmpp;class=od;type=fg;nlev=60;sh=159;gg=80;tres=_fg006up4tr3
    !
    form   = 'tmpp'
    call goVarValue( archivekeys, ';', 'form', '=', form, status )
    IF_ERROR_RETURN(status=1)
    !
    ec_class   = 'od'
    call goVarValue( archivekeys, ';', 'class', '=', ec_class, status )
    IF_ERROR_RETURN(status=1)
    !
    ec_type    = 'fc'
    call goVarValue( archivekeys, ';', 'type', '=', ec_type, status )
    IF_ERROR_RETURN(status=1)
    !
    ec_levs    = 'ml60'
    call goVarValue( archivekeys, ';', 'levs', '=', ec_levs, status )
    IF_ERROR_RETURN(status=1)
    !
    ec_sh      = 159
    call goVarValue( archivekeys, ';', 'sh', '=', ec_sh, status )
    IF_ERROR_RETURN(status=1)
    !
    ec_gg      = 80
    call goVarValue( archivekeys, ';', 'gg', '=', ec_gg, status )
    IF_ERROR_RETURN(status=1)
    !
    ec_tres    = ''
    call goVarValue( archivekeys, ';', 'tres', '=', ec_tres, status )
    IF_ERROR_RETURN(status=1)
    
    ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
    ! fg data available up to 2000-09-12 00:00 ; after, use fc data:
    if ( ec_tres == '_fg006up4tr3' ) then
      tc = NewDate( 2000, 09, 12, 00, 00, 00 )
      if ( t1 > tc ) then
        write (gol,'("      WARNING - using fc data after 2000-09-12 00:00")'); call goPr
        ec_type = 'fc'
        ec_tres = '_fc012up2tr3'
      end if
    end if
    ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

    ! store some values ...
    mf%ec_class = ec_class
    mf%ec_type  = ec_type

    ! Date of the first forecast with cycle 37r3
    tc = NewDate( 2011, 11, 15, 12, 00, 00 )
    
    !
    ! meteo file
    !
    
    ! grid : T159, N80, etc 
    call goWriteKeyNum( gridT, 'T', ec_sh )
    call goWriteKeyNum( gridN, 'N', ec_gg )
    
    ! defaults
    levs           = ec_levs
    grid           = 'Xxx'
    mfile          = trim(paramkey)
    mf%paramkeys   = '-'//trim(paramkey)//'-'
    mf%constant    = .false.
    mf%do_spm      = .false.
    mf%spm_fname   = 'none'
    mf%spm_lnsp2sp = .false.
    spm_tres       = ec_tres
    
    ! specials
    select case ( paramkey )
      case ( 'LNSP' )
        levs           = 'ml1'
        grid           = gridT
        if ( form /= '3hr' ) ec_tres = '_03p03'
      case ( 'VO', 'D' )
        mfile = 'VOD'
        mf%paramkeys   = '-VO-D-'
        grid           = gridT
        mf%do_spm      = .true.
        spm_tres       = '_03p03'
        if ( form == '3hr' ) then
          mfile        = paramkey
          mf%paramkeys = '-'//trim(paramkey)//'-'
          spm_tres     = ec_tres
        end if
      case ( 'T', 'W' )
        grid           = gridT
        mf%do_spm      = .true.
        spm_tres       = '_03p03'
        if ( form == '3hr' ) then
          spm_tres     = ec_tres
        end if
      case ( 'Q' )
        grid           = gridN
        mf%do_spm      = .true.
        spm_tres       = '_03p03'
        mf%spm_lnsp2sp = .true.
        if ( form == '3hr' ) then
          spm_tres     = ec_tres
        end if
      case ( 'CLWC', 'CIWC', 'CC', 'CCO', 'CCU' )
        grid           = gridN
        mfile          = 'CLD'
        mf%paramkeys   = '-CLWC-CIWC-CC-CCO-CCU-'
        mf%do_spm      = .true.
        spm_tres       = '_03p03'
        mf%spm_lnsp2sp = .true.
        if ( form == '3hr' ) then
          spm_tres     = ec_tres
        end if
      case ( 'oro' )
        levs           = 'sfc'
        grid           = gridN
        mf%constant    = .true.
        ec_tres        = 'const'
      case ( 'lsm' )
        levs           = 'sfc'
        grid           = gridN
        mf%constant    = .true.
      case ( 'cvl', 'cvh', 'tvl', 'tvh', 'sr', 'al', 'lsrh' )
        ec_type        = 'an'
        levs           = 'sfc'
        grid           = gridN
        mfile          = 'S0'
        mf%paramkeys   = '-cvl-cvh-tvl-tvh-sr-al-lsrh-'
        ec_tres        = '_06p06'
! Before cycle 37r3
!      case ( 'ci', 'swvl1', 'swvl2', 'swvl3', 'swvl4', '10fg', 'sd', 'lsp', &
!             'cp', 'sf', 'sshf', 'slhf', 'blh', 'u10m', 'v10m', 't2m', 'd2m', &
!             'ssr', 'ewss', 'nsss', 'sstr' ,'src' )
!        levs           = 'sfc'
!        grid           = gridN
!        mfile          = 'S1'
!        mf%paramkeys   = '-'
!        mf%paramkeys   = trim(mf%paramkeys)//'ci-swvl1-swvl2-swvl3-swvl4-10fg-sd-lsp-'
!        mf%paramkeys   = trim(mf%paramkeys)//'cp-sf-sshf-slhf-blh-10u-10v-2t-2d-'
!        mf%paramkeys   = trim(mf%paramkeys)//'ssr-ewss-nsss-sstr-src-'
!        ec_tres        = '_03p03'
! After cycle 37r3
     case ( 'ci', 'swvl1', 'swvl2', 'swvl3', 'swvl4', '10fg','10fg3', 'sd', 'lsp', &
          'cp', 'sf', 'sshf', 'slhf', 'blh', 'u10m', 'v10m', 't2m', 'd2m', &
          'ssr', 'ewss', 'nsss', 'sstr' ,'src' )
        levs           = 'sfc'
        grid           = gridN
        mfile          = 'S1'
        mf%paramkeys   = '-'

        ! Before or after cycle 37r3
        if ( (t1 > tc).and.(trim(ec_class)=='od') ) then
           mf%paramkeys   = trim(mf%paramkeys)//'ci-swvl1-swvl2-swvl3-swvl4-10fg3-sd-lsp-'
        else
           mf%paramkeys   = trim(mf%paramkeys)//'ci-swvl1-swvl2-swvl3-swvl4-10fg-sd-lsp-'
        end if
        
        mf%paramkeys   = trim(mf%paramkeys)//'cp-sf-sshf-slhf-blh-10u-10v-2t-2d-'
        mf%paramkeys   = trim(mf%paramkeys)//'ssr-ewss-nsss-sstr-src-'
        ec_tres        = '_03p03'
    end select
    
    ! convert input times to file name times:
    call GetGribTime( ec_tres, tday, t1, t2, status, tfile=tfile, trange=mf%trange )
    IF_NOTOK_RETURN(status=1)
    
    ! store some time params
    mf%treskey = ec_tres
    if ( ec_type == 'an' ) mf%treskey = 'an'//ec_tres
    mf%tday    = tday
    
    ! create file name:
    !   dir/od-fc-2000-01-ml60-T159-T_20000101_fg006up4tr3.gb
    select case ( paramkey )
      case ( 'oro' )
        write (mf%fname,'(a,"/",a,"-",a,"-",i4.4,"-",i2.2,"-",a,"-",a,"-",a,a)') &
                  trim(dir), &
                  trim(ec_class), trim(ec_type), 0000, 00, trim(levs), trim(grid), &
                  'Z', '.gb'
      case ( 'lsm' )
        write (mf%fname,'(a,"/",a,"-",a,"-",i4.4,"-",i2.2,"-",a,"-",a,"-",a,a)') &
                  trim(dir), &
                  trim(ec_class), trim(ec_type), 0000, 00, trim(levs), trim(grid), &
                  'LSM', '.gb'
      case default
        ! extract time values:
        call Get( tfile, year=ccyy, month=mm, day=dd )
        ! path and file include date:
        write (mf%fname,'(a,"/",a,"-",a,"-",i4.4,"-",i2.2,"-",a,"-",a,"-",a,"_",i4.4,2i2.2,a,a)') &
                  trim(dir), &
                  trim(ec_class), trim(ec_type), ccyy, mm, trim(levs), trim(grid), &
                  trim(mfile), ccyy, mm, dd, trim(ec_tres), '.gb'
    end select
    
    ! exist ?
    inquire( file=mf%fname, exist=exist )
    if ( .not. exist ) then
      write (gol,'("grib file does not exist:")'); call goErr
      write (gol,'("  ",a)') trim(mf%fname); call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if

    
    !
    ! sp or lnsp ?
    !
    
    if ( mf%do_spm ) then

      ! grid : T159    
      grid = gridT
      
      ! levels: ml1
      levs = 'ml1'

      ! create file name:
      !   dir/od-fc-2000-01-ml60-T159-T_20000101_fg006up4tr3.gb
      write (mf%spm_fname,'(a,"/",a,"-",a,"-",i4.4,"-",i2.2,"-",a,"-",a,"-",a,"_",i4.4,2i2.2,a,a)') &
                trim(dir), &
                trim(ec_class), trim(ec_type), ccyy, mm, trim(levs), trim(grid), &
                'LNSP', ccyy, mm, dd, trim(spm_tres), '.gb'

      ! exist ?
      inquire( file=mf%spm_fname, exist=exist )
      if ( .not. exist ) then
        write (gol,'("grib file does not exist:")'); call goErr
        write (gol,'("  ",a)') trim(mf%spm_fname); call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
      end if

    end if

    ! ok
    status = 0
 
  end subroutine mf_Init
  
  
  ! ***
  
  
  subroutine mf_Done( mf, status )
  
    use file_grib, only : Done
    
    ! --- in/out ------------------------------------

    type(TMeteoFile_grib), intent(inout)  ::  mf
    integer, intent(out)                  ::  status

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

    ! nothing to be done ...

    ! ok
    status = 0
    
  end subroutine mf_Done
  

  ! ***

  
  subroutine mf_Get( mf, status, trange1, trange2, paramkeys )
    
    use GO, only : TDate
    
    ! --- in/out ----------------------------
    
    type(TMeteoFile_grib), intent(in)   ::  mf
    integer, intent(out)                ::  status

    type(TDate), intent(out), optional        ::  trange1, trange2
    character(len=*), intent(out), optional   ::  paramkeys
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_Get'
    
    ! --- local --------------------------------
    
    ! --- begin --------------------------------

    ! time range:
    if ( present(trange1) ) trange1 = mf%trange(1)
    if ( present(trange2) ) trange2 = mf%trange(2)

    ! contents:
    if ( present(paramkeys) ) paramkeys = trim(mf%paramkeys)
    
    ! ok
    status = 0
 
  end subroutine mf_Get
  

  ! ***

  
  subroutine mf_ReadRecord( mf, paramkey, t1, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll, sp_ll, &
                                ggi, gg, sp_gg, &
                                shi, sh, lnsp_sh, &
                                tmi, status )
                                
    use parray    , only : pa_Init, pa_Done
    use GO        , only : TDate
    use Grid      , only : TLevelInfo
    use Grid      , only : TllGridInfo, TggGridInfo, TshGridInfo
    use tmm_info  , only : TMeteoInfo

    ! --- in/out -------------------------------
    
    type(TMeteoFile_grib), intent(inout)  ::  mf
    character(len=*), intent(in)          ::  paramkey
    type(TDate), intent(in)               ::  t1, t2
    character(len=1), intent(in)          ::  nuv
    character(len=1), intent(in)          ::  nw
    character(len=2), intent(out)         ::  gridtype
    type(TLevelInfo), intent(out)         ::  levi
    type(TllGridInfo), intent(inout)      ::  lli
    real, pointer                         ::  ll(:,:,:)
    real, pointer                         ::  sp_ll(:,:)
    type(TggGridInfo), intent(inout)      ::  ggi
    real, pointer                         ::  gg(:,:)
    real, pointer                         ::  sp_gg(:)
    type(TshGridInfo), intent(inout)      ::  shi
    complex, pointer                      ::  sh(:,:)
    complex, pointer                      ::  lnsp_sh(:)
    type(TMeteoInfo), intent(out)         ::  tmi
    integer, intent(out)                  ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_ReadRecord'
    
    ! --- local -------------------------------
    
    real, pointer         ::  ll2(:,:,:)
    real, pointer         ::  gg2(:,:)
    complex, pointer      ::  sh2(:,:)
        
    ! --- begin ---------------------------------
    
    ! combined field ?
    select case ( paramkey )
    
      case ( 'sstr' )
      
        ! read first field:
        call mf_ReadRecord_1( mf, 'ewss', t1, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll, sp_ll, &
                                ggi, gg, sp_gg, &
                                shi, sh, lnsp_sh, &
                                tmi, status )
        IF_NOTOK_RETURN(status=1)  

        ! init pointer:
        call pa_Init( ll2 )
        call pa_Init( gg2 )
        call pa_Init( sh2 )

        ! read second field:
        call mf_ReadRecord_1( mf, 'nsss', t1, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll2, sp_ll, &
                                ggi, gg2, sp_gg, &
                                shi, sh2, lnsp_sh, &
                                tmi, status )
         IF_NOTOK_RETURN(status=1)      

        ! process:
        select case ( gridtype )
          case ( 'll' ) ; ll = sqrt( ll**2 + ll2**2 )
          case ( 'gg' ) ; gg = sqrt( gg**2 + gg2**2 )
          case ( 'sh' ) ; sh = sqrt( sh**2 + sh2**2 )
          case default
            write (gol,'("unsupported gridtype for substract :",a)') gridtype; call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select   

        ! clear pointers:
        call pa_Done( ll2 )
        call pa_Done( gg2 )
        call pa_Done( sh2 )

      case default
      
        call mf_ReadRecord_1( mf, paramkey, t1, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll, sp_ll, &
                                ggi, gg, sp_gg, &
                                shi, sh, lnsp_sh, &
                                tmi, status )
         IF_NOTOK_RETURN(status=1)
      
    end select
    
    ! ok
    status = 0

  end subroutine mf_ReadRecord
  


  ! ***

  
  subroutine mf_ReadRecord_1( mf, paramkey, t1, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll, sp_ll, &
                                ggi, gg, sp_gg, &
                                shi, sh, lnsp_sh, &
                                tmi, status )
                                
    use parray    , only : pa_Init, pa_Done
    use GO        , only : TDate, operator(<), operator(-), rTotal, wrtgol
    use Grid      , only : TLevelInfo
    use Grid      , only : TllGridInfo, TggGridInfo, TshGridInfo
    use tmm_info  , only : TMeteoInfo

    ! --- in/out -------------------------------
    
    type(TMeteoFile_grib), intent(inout)  ::  mf
    character(len=*), intent(in)          ::  paramkey
    type(TDate), intent(in)               ::  t1, t2
    character(len=1), intent(in)          ::  nuv
    character(len=1), intent(in)          ::  nw
    character(len=2), intent(out)         ::  gridtype
    type(TLevelInfo), intent(out)         ::  levi
    type(TllGridInfo), intent(inout)      ::  lli
    real, pointer                         ::  ll(:,:,:)
    real, pointer                         ::  sp_ll(:,:)
    type(TggGridInfo), intent(inout)      ::  ggi
    real, pointer                         ::  gg(:,:)
    real, pointer                         ::  sp_gg(:)
    type(TshGridInfo), intent(inout)      ::  shi
    complex, pointer                      ::  sh(:,:)
    complex, pointer                      ::  lnsp_sh(:)
    type(TMeteoInfo), intent(out)         ::  tmi
    integer, intent(out)                  ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_ReadRecord_1'
    
    ! --- local -------------------------------
    
    type(TDate)           ::  tref
    real, pointer         ::  ll1(:,:,:)
    real, pointer         ::  gg1(:,:)
    complex, pointer      ::  sh1(:,:)
    real                  ::  dt_sec
        
    ! --- begin ---------------------------------
    
    ! accumulated field ?
    select case ( paramkey )
    
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! accumulated fields 
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      case ( 'lsp', 'cp', 'sf', 'sshf', 'slhf', 'ssr', 'ewss', 'nsss' )
      
        ! get reference time for requested time interval:
        call GetGribTime( mf%treskey, mf%tday, t1, t2, status, tref=tref )
        IF_NOTOK_RETURN(status=1)

        ! should be a time interval ...
        if ( .not. (t1 < t2) ) then
          write (gol,'("accumulated fields requires time interval:")'); call goErr
          write (gol,'("  paramkey : ",a)') paramkey; call goErr
          call wrtgol( '  t1       : ', t1 ); call goErr
          call wrtgol( '  t2       : ', t2 ); call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
        end if
      
        ! read field accumulated over [tref,t2] :
        call mf_ReadRecord_2( mf, paramkey, tref, t2, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll, sp_ll, &
                                ggi, gg, sp_gg, &
                                shi, sh, lnsp_sh, &
                                tmi, status )
        IF_NOTOK_RETURN(status=1)  

        ! substract [tref,t1] if necessary:
        if ( tref < t1 ) then
        
          ! init pointer:
          call pa_Init( ll1 )
          call pa_Init( gg1 )
          call pa_Init( sh1 )

          ! read field accumulated over [tref,t1] :
          call mf_ReadRecord_2( mf, paramkey, tref, t1, t1, nuv, nw, &
                                  gridtype, levi, &
                                  lli, ll1, sp_ll, &
                                  ggi, gg1, sp_gg, &
                                  shi, sh1, lnsp_sh, &
                                  tmi, status )
           IF_NOTOK_RETURN(status=1)      

          ! substract:
          select case ( gridtype )
            case ( 'll' ) ; ll = ll - ll1
            case ( 'gg' ) ; gg = gg - gg1
            case ( 'sh' ) ; sh = sh - sh1
            case default
              write (gol,'("unsupported gridtype for substract :",a)') gridtype; call goErr
              write (gol,'("in ",a)') rname; call goErr; status=1; return
          end select   

          ! clear pointers:
          call pa_Done( ll1 )
          call pa_Done( gg1 )
          call pa_Done( sh1 )
          
        end if
        
        ! return time averages only:
        dt_sec = rTotal( t2 - t1, 'sec' )
        select case ( gridtype )
          case ( 'll' ) ; ll = ll / dt_sec
          case ( 'gg' ) ; gg = gg / dt_sec
          case ( 'sh' ) ; sh = sh / dt_sec
          case default
            write (gol,'("unsupported gridtype for time average :",a)') gridtype; call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select 
      
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! instantaneous fields
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      case default
      
        ! get reference time for requested time interval:
        call GetGribTime( mf%treskey, mf%tday, t1, t2, status, tref=tref )
        IF_NOTOK_RETURN(status=1)

        ! just read ..
        call mf_ReadRecord_2( mf, paramkey, tref, t1, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll, sp_ll, &
                                ggi, gg, sp_gg, &
                                shi, sh, lnsp_sh, &
                                tmi, status )
         IF_NOTOK_RETURN(status=1)
      
    end select
    
    ! ok
    status = 0

  end subroutine mf_ReadRecord_1
  

  ! ***

  
  subroutine mf_ReadRecord_2( mf, paramkey, tref, t1, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll, sp_ll, &
                                ggi, gg, sp_gg, &
                                shi, sh, lnsp_sh, &
                                tmi, status )

    use GO        , only : TDate, wrtgol
    use Grid      , only : Init, Done
    use Grid      , only : TLevelInfo
    use Grid      , only : TllGridInfo, TggGridInfo, TshGridInfo
    use Grid      , only : Interpol
    use file_grib , only : Init, Done, ReadRecord, Get, Check
    use file_grib , only : levtype_sfc, levtype_hyb
    use file_grib , only : gridtype_ll, gridtype_gg, gridtype_sh
    use grib_table, only : GetPid
    use PArray    , only : pa_Init, pa_Done, pa_SetShape
    use tmm_info  , only : TMeteoInfo, Init, AddHistory

    ! --- in/out -------------------------------
    
    type(TMeteoFile_grib), intent(inout)  ::  mf
    character(len=*), intent(in)          ::  paramkey
    type(TDate), intent(in)               ::  tref, t1, t2
    character(len=1), intent(in)          ::  nuv
    character(len=1), intent(in)          ::  nw
    character(len=2), intent(out)         ::  gridtype
    type(TLevelInfo), intent(out)         ::  levi
    type(TllGridInfo), intent(inout)      ::  lli
    real, pointer                         ::  ll(:,:,:)
    real, pointer                         ::  sp_ll(:,:)
    type(TggGridInfo), intent(inout)      ::  ggi
    real, pointer                         ::  gg(:,:)
    real, pointer                         ::  sp_gg(:)
    type(TshGridInfo), intent(inout)      ::  shi
    complex, pointer                      ::  sh(:,:)
    complex, pointer                      ::  lnsp_sh(:)
    type(TMeteoInfo), intent(out)         ::  tmi
    integer, intent(out)                  ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_ReadRecord_2'
    
    ! --- local -------------------------------
    
    logical                ::  isfirst
    logical                ::  reopened
    integer                ::  pid, tabid
    integer                ::  nlev, glevtype, glevel

    integer                ::  level

    integer                ::  ggridtype
    real                   ::  lon_first, lon_inc
    integer                ::  lon_n
    real                   ::  lat_first, lat_inc
    integer                ::  lat_n
    integer                ::  ggN
    integer                ::  shT
    
    integer                ::  greftime(5), gtimerange(4)
    character(len=64)      ::  key

    integer                ::  ilat
    real, pointer          ::  pat(:,:)

    type(TshGridInfo)      ::  tmp_shi
    complex, pointer       ::  tmp_sh(:)

    ! --- begin ---------------------------------
    
    ! no fluxes through boundaries ...
    if ( nuv /= 'n' ) then
      write (gol,'("unsupported nuv key : ",a)') nuv; call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if

    ! init pointer arrays:    
    call pa_Init( pat )

    !
    ! ~~~ 3d field or 2d stored in 3d array
    !
    
    ! open grib file
    call Init( mf%grib, mf%fname, 'r', status )
    IF_NOTOK_RETURN(status=1)

    ! arrays and grids not defined yet
    isfirst = .true.
    reopened = .false.
    
    ! loop over records
    level = 0
    do
    
      !
      ! read gribsection in file buffer
      !
      
      call ReadRecord( mf%grib, status )
      select case ( status )
        case ( 0 )
          ! no error
        case ( 1 )
          ! eof
          if ( .not. reopened ) then
            !write (*,'("grib read record: re-open ...")')
            ! close:
            call Done( mf%grib, status )
            IF_NOTOK_RETURN(status=1)
            ! reopen:
            call Init( mf%grib, mf%fname, 'r', status )
            IF_NOTOK_RETURN(status=1)
            reopened = .true.
            cycle
          else
            write (gol,'("reached eof before requested record was found")'); call goErr
            write (gol,'("  file     : ",a)') trim(mf%fname); call goErr
            write (gol,'("  paramkey : ",a)') trim(paramkey); call goErr
            call wrtgol( '  tref     : ', tref ); call goErr
            call wrtgol( '  t1       : ', t1 ); call goErr
            call wrtgol( '  t2       : ', t2 ); call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
          end if
        case default
          write (gol,'("error from grib ReadRecord; status=",i6)') status; call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select
      
      
      !
      ! checks ...
      !
      
      ! get param id for the requested field from grib table:
      select case ( paramkey )
        case ( 'spm' ) ; call GetPid( 'ec', 'SP'    , pid, tabid, status )
        case default   ; call GetPid( 'ec', paramkey, pid, tabid, status )
      end select
      IF_NOTOK_RETURN(status=1)

      ! check parameter; continue if not ok:
      call Check( mf%grib, pid=pid, debug=0, status=status )
      if (status/=0) cycle

      ! fill times ?
      if ( .not. mf%constant ) then
      
        ! extract time fields from grib, store in mf%tref/mf%t1/mf%t2
        call SetTime( mf, status )
        IF_NOTOK_RETURN(status=1)

        ! check time:
        call CheckTime( mf, tref, t1, t2, status )
        if (status/=0) then
          !write (*,'("grib read record: wrong time; skip ...")')
          cycle
          !write (gol,'("found unexpected times in grib file:")'); call goErr
          !write (gol,'("  paramkey     : ",a)') paramkey; call goErr
          !call wrtgol( '  req. t1 : ', t1 ); call goErr
          !call wrtgol( '       t2 : ', t2 ); call goErr
          !call wrtgol( '  grib t1 : ', mf%t1 ); call goErr
          !call wrtgol( '       t2 : ', mf%t2 ); call goErr
          !write (gol,'("  grib file    : ",a)') trim(mf%fname); call goErr
          !write (gol,'("in ",a)') rname; call goErr; status=1; return
        end if
        
      end if  ! time checking

      ! extract level stuff:
      call Get( mf%grib, nlev=nlev, levtype=glevtype, level=glevel, status=status )
      IF_NOTOK_RETURN(status=1)
      ! check level type:
      select case ( glevtype )
        case ( levtype_sfc )
          ! surface field
          nlev = 1
          glevel = 1
        case ( levtype_hyb )
          select case ( paramkey )
            case ( 'LNSP' )
              nlev = 1
            case default
              ! level in 3d field
          end select
        case default
          write (gol,'("found unexpected level type: ")'); call goErr
          write (gol,'("  leveltype : ",i3)') glevtype; call goErr
          write (gol,'("  paramkey  : ",a)') paramkey; call goErr
          write (gol,'("  grib file : ",a)') trim(mf%fname); call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select
      ! check level:
      level = level + 1
      if ( glevel /= level ) then
        write (gol,'("found unexpected level: ")'); call goErr
        write (gol,'("  paramkey   : ",a)') paramkey; call goErr
        write (gol,'("  req. level : ",i6)') level; call goErr
        write (gol,'("  grib level : ",i6)') glevel; call goErr
        write (gol,'("  grib nlev  : ",i6)') nlev; call goErr
        write (gol,'("  grib file  : ",a)') trim(mf%fname); call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
      end if
      
      !
      ! define grids and arrays
      !
      if ( isfirst ) then
        !
        ! * info
        !
        ! example of history:
        !   model=ecmwf;class=od;type=fc;tref=2000,12,31,12,00; ...
        !     trange=001,234,240,001;sh=159;nlev=60
        !
        call Init( tmi, paramkey, 'unknown', status )
        call AddHistory( tmi, 'model==ecmwf', status )
        call AddHistory( tmi, 'class=='//trim(mf%ec_class), status )
        call AddHistory( tmi, 'type=='//trim(mf%ec_type)  , status )
        !
        call Get( mf%grib, status, reftime=greftime, timerange=gtimerange )
        IF_NOTOK_RETURN(status=1)
        write (key,'("tref==",i4.4,4(",",i2.2))') greftime
        call AddHistory( tmi, trim(key), status )
        write (key,'("trange==",i3.3,3(",",i3.3))') gtimerange
        call AddHistory( tmi, trim(key), status )
        !
        write (key,'("nlev==",i3.3)') nlev
        call AddHistory( tmi, trim(key), status )

        !
        ! * define horizontal grid:
        !

        ! extract grid type:
        call Get( mf%grib, status, gridtype=ggridtype )
        IF_NOTOK_RETURN(status=1)

        ! setup:
        select case ( ggridtype )

          ! o lat/lon
          case ( gridtype_ll )

            ! routine returns lat/lon grid:
            gridtype = 'll'

            ! grib storage is north pole to south pole:
            call Get( mf%grib, status, &
                               lon_first=lon_first, lon_inc=lon_inc, lon_n=lon_n, &
                               lat_last =lat_first, lat_inc=lat_inc, lat_n=lat_n     )
            IF_NOTOK_RETURN(status=1)

            ! define grid structure:
            call Init( lli, lon_first, lon_inc, lon_n, &
                            lat_first, lat_inc, lat_n, status     )
            IF_NOTOK_RETURN(status=1)
            
            ! init array to store 2d field from grib file (north-south order):
            call pa_SetShape( pat, lon_n, lat_n )
            
            ! allocate output:
            call pa_SetShape( ll, lon_n, lat_n, nlev )

            ! add to history:
            write (key,'("longrid==",f7.2,",",f6.2,",",i4)') lon_first, lon_inc, lon_n
            call AddHistory( tmi, trim(key), status )
            write (key,'("latgrid==",f7.2,",",f6.2,",",i4)') lat_first, lat_inc, lat_n
            call AddHistory( tmi, trim(key), status )

          ! o gaussian grid
          case ( gridtype_gg )

            ! routine returns gg grid:
            gridtype = 'gg'

            ! extract grid number:
            call Get( mf%grib, status, N=ggN )
            IF_NOTOK_RETURN(status=1)

            ! define grid structure:
            call Init( ggi, ggN, .true., status )
            IF_NOTOK_RETURN(status=1)
            
            ! allocate output:
            call pa_SetShape( gg, ggi%np, nlev )

            ! add to history:
            write (key,'("gg==",i4.4)') ggN
            call AddHistory( tmi, trim(key), status )

          ! o spectral field:
          case ( gridtype_sh )

            ! routine returns sh grid:
            gridtype = 'sh'

            ! extract spectral truncation:
            call Get( mf%grib, status, T=shT )
            IF_NOTOK_RETURN(status=1)

            ! intialize spherical harmonic field info:
            call Init( shi, shT, status )
            IF_NOTOK_RETURN(status=1)

            ! allocate output:
            call pa_SetShape( sh, shi%np, nlev )

            ! add to history:
            write (key,'("sh==",i4.4)') shT
            call AddHistory( tmi, trim(key), status )

          case default
            write (gol,'("unsupported gridtype for setup : ",i6)') ggridtype; call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select

        !
        ! * levels
        !
        select case ( nlev )
          case ( 1 )
            call Init( levi, nlev, (/0.0,0.0/), (/0.0,0.0/), status )
            IF_NOTOK_RETURN(status=1)
          case ( 60 )
            call Init( levi, 'ec60', status )
            IF_NOTOK_RETURN(status=1)
        end select
        
        ! not again ...
        isfirst = .false.
      end if  ! isfirst (grid definition and allocation)
      
      !
      ! store
      !
      
      select case ( ggridtype )
      
        case ( gridtype_ll )
      
          ! read 2d pat from grib; storred from north to south
          call Get( mf%grib, status, ll=pat )
          IF_NOTOK_RETURN(status=1)

          ! store from south to north:
          do ilat = 1, lat_n
            ll(:,ilat,level) = pat(:,lat_n+1-ilat)
          end do
      
        case ( gridtype_gg )
      
          ! read 2d pat from grib:
          call Get( mf%grib, status, gg=gg(:,level) )
          IF_NOTOK_RETURN(status=1)

        case ( gridtype_sh )
      
          ! read 2d pat from grib:
          call Get( mf%grib, status, sh=sh(:,level) )
          IF_NOTOK_RETURN(status=1)

        case default
          write (gol,'("unsupported gridtype for 2d pat : ",i6)') gridtype; call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select

      ! last record for this field ?
      if ( glevel == nlev ) exit

    end do   ! records
    
    ! close grib file
    call Done( mf%grib, status )
    IF_NOTOK_RETURN(status=1)

    
    ! 
    ! ~~~ surface pressure
    !
    
    if ( mf%do_spm ) then
    
      ! open grib file
      call Init( mf%spm_grib, mf%spm_fname, 'r', status )
      IF_NOTOK_RETURN(status=1)
      
      ! loop over time records
      do

        ! read gribsection in file buffer
        call ReadRecord( mf%spm_grib, status )
        IF_NOTOK_RETURN(status=1)

        ! fill times
        call SetTime( mf, status )
        IF_NOTOK_RETURN(status=1)

        ! check time:
        call CheckTime( mf, tref, t1, t2, status )
        if (status/=0) then
          !write (*,'("grib read record: spm wrong time; skip ...")')
          cycle
          !write (gol,'("found unexpected times in grib file:")'); call goErr
          !write (gol,'("  paramkey     : ",a)') paramkey; call goErr
          !call wrtgol( '  req. t1 : ', t1 ); call goErr
          !call wrtgol( '       t2 : ', t2 ); call goErr
          !call wrtgol( '  grib t1 : ', mf%t1 ); call goErr
          !call wrtgol( '       t2 : ', mf%t2 ); call goErr
          !write (gol,'("  grib file    : ",a)') trim(mf%spm_fname); call goErr
          !write (gol,'("in ",a)') rname; call goErr; status=1; return
        end if
        
        ! time ok
        exit
        
      end do  ! time loop

      ! set param id:
      select case ( ggridtype )
        case ( gridtype_ll )
          call GetPid( 'ec', 'SP', pid, tabid, status )
          IF_NOTOK_RETURN(status=1)
        case ( gridtype_gg )
          if ( mf%spm_lnsp2sp ) then
            call GetPid( 'ec', 'LNSP', pid, tabid, status )
            IF_NOTOK_RETURN(status=1)
          else
            call GetPid( 'ec', 'SP', pid, tabid, status )
            IF_NOTOK_RETURN(status=1)
          end if
        case ( gridtype_sh )
          call GetPid( 'ec', 'LNSP', pid, tabid, status )
          IF_NOTOK_RETURN(status=1)
        case default
          write (gol,'("unsupported gridtype for setup sp/lnsp : ",i6)') ggridtype; call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select
      ! check parameter:
      call Check( mf%spm_grib, pid=pid, debug=1, status=status )
      IF_NOTOK_RETURN(status=1)

      ! check level:
      call Get( mf%spm_grib, levtype=glevtype, level=glevel, status=status )
      IF_NOTOK_RETURN(status=1)
      select case ( ggridtype )
        case ( gridtype_ll )
          if ( glevtype /= levtype_sfc ) then
            write (gol,'("found unexpected level type ")'); call goErr
            write (gol,'("  paramkey        : ",a)') paramkey; call goErr
            write (gol,'("  sfc level type  : ",i6)') levtype_sfc; call goErr
            write (gol,'("  grib level type : ",i6)') glevtype; call goErr
            write (gol,'("  grib file       : ",a)') trim(mf%spm_fname); call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
          end if
        case ( gridtype_gg )
          if ( mf%spm_lnsp2sp ) then
            if ( (glevtype /= levtype_hyb) .or. (glevel /= 1) ) then
              write (gol,'("found unexpected level type (lnsp for 3d gg)")'); call goErr
              write (gol,'("  paramkey        : ",a)') paramkey; call goErr
              write (gol,'("  hyb level type  : ",i6)') levtype_hyb; call goErr
              write (gol,'("  grib level type : ",i6)') glevtype; call goErr
              write (gol,'("  grib level      : ",i6)') glevel; call goErr
              write (gol,'("  grib file       : ",a)') trim(mf%spm_fname); call goErr
              write (gol,'("in ",a)') rname; call goErr; status=1; return
            end if
          else
            if ( glevtype /= levtype_sfc ) then
              write (gol,'("found unexpected level type ")'); call goErr
              write (gol,'("  paramkey        : ",a)') paramkey; call goErr
              write (gol,'("  sfc level type  : ",i6)') levtype_sfc; call goErr
              write (gol,'("  grib level type : ",i6)') glevtype; call goErr
              write (gol,'("  grib file       : ",a)') trim(mf%spm_fname); call goErr
              write (gol,'("in ",a)') rname; call goErr; status=1; return
            end if
          end if
        case ( gridtype_sh )
          if ( (glevtype /= levtype_hyb) .or. (glevel /= 1) ) then
            write (gol,'("found unexpected level type ")'); call goErr
            write (gol,'("  paramkey        : ",a)') paramkey; call goErr
            write (gol,'("  hyb level type  : ",i6)') levtype_hyb; call goErr
            write (gol,'("  grib level type : ",i6)') glevtype; call goErr
            write (gol,'("  grib level      : ",i6)') glevel; call goErr
            write (gol,'("  grib file       : ",a)') trim(mf%spm_fname); call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
          end if
        case default
          write (gol,'("unsupported gridtype for sp/lnsp levs : ",i6)') ggridtype; call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select
      
      ! read and store surface pressure field:
      select case ( ggridtype )

        case ( gridtype_ll )
          ! allocate storage
          call pa_SetShape( sp_ll, lon_n, lat_n )
          ! read 2d pat from grib; storred from north to south
          call Get( mf%spm_grib, status, ll=pat )
          IF_NOTOK_RETURN(status=1)
          ! store from south to north:
          do ilat = 1, lat_n
            sp_ll(:,ilat) = pat(:,lat_n+1-ilat)
          end do

        case ( gridtype_gg )
          ! allocate storage
          call pa_SetShape( sp_gg, ggi%np )
          ! convert from sh lnsp or read directly:
          if ( mf%spm_lnsp2sp ) then
            ! allocate output:
            call pa_SetShape( sh, shi%np, nlev )
            ! extract spectral truncation:
            call Get( mf%spm_grib, status, T=shT )
            IF_NOTOK_RETURN(status=1)
            ! intialize spherical harmonic field info:
            call Init( tmp_shi, shT, status )
            IF_NOTOK_RETURN(status=1)
            ! allocate storage for lnsp:
            call pa_Init( tmp_sh )
            call pa_SetShape( tmp_sh, tmp_shi%np )
            ! read lnsp on sh grid:
            call Get( mf%spm_grib, status, sh=tmp_sh )
            IF_NOTOK_RETURN(status=1)
            ! interpolate from sh to gg:
            call Interpol( tmp_shi, tmp_sh, ggi, sp_gg, status )  ! ln(Pa)
            IF_NOTOK_RETURN(status=1)
            ! clear
            call Done( tmp_shi )
            call pa_Done( tmp_sh )
            ! convert from lnsp to sp :
            sp_gg = exp(sp_gg)   ! Pa
          else
            ! read 2d pat from grib; storred from north to south
            call Get( mf%spm_grib, status, gg=sp_gg )
            IF_NOTOK_RETURN(status=1)
          end if
      
        case ( gridtype_sh )
          ! allocate storage
          call pa_SetShape( lnsp_sh, shi%np )
          ! read spectral coeff:
          call Get( mf%spm_grib, status, sh=lnsp_sh )
          IF_NOTOK_RETURN(status=1)
      
        case default
          write (gol,'("unsupported gridtype for reading sp/lnsp : ",i6)') ggridtype; call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select

      ! close grib file
      call Done( mf%spm_grib, status )
      IF_NOTOK_RETURN(status=1)

    end if
    
    ! 
    ! ~~~ end
    !
    
    ! deallocate arrays
    call pa_Done( pat )
    
    ! ok
    status = 0

  end subroutine mf_ReadRecord_2

  
  ! ****************************************************************************

  
  !
  ! In gribfile:
  !    reftime : for example time at which forecast is made
  !    timerange: increment or interval
  !
  ! arguments:          ok if:
  !
  !   time1 == time2      time1==time2 == reftime+timerange
  !
  !   time1 == 0          time2 == reftime+timerange
  ! 
  !   time2 == 0          time1 == reftime
  ! 
  !   time1 < time2       time1 == refitme, time2 == reftime+timerange
  !
  !
  !   grib [t1-----------t2]
  !                                                              status
  !    o                         time1/2        record too old      1
  !    o                time1/2                 ok                  0
  !    o       time1/2                          record too new      2
  !
  !
  ! SetTime( mf, status )
  !
  !   Extracts time values from current grib record, 
  !   store in mf%t1, mf%t2
  !
  !   return status:
  !     0       : ok
  !     other   : some error
  !
  ! CheckTime( mf, time1, time2, status )
  !
  !   return status:
  !     0   : times match
  !     1   : times do not match, try next record
  !     2   : current record is newer than requested (reopen ?)
  !     3   : some error
  !

  
  ! ***

  
  !
  ! Return time parameters for grib files:
  !  o tfile   :  date in filename
  !  o trange  :  time interval covered by fields in file
  !  o tref    :  reference time (forecast start?) for tday,[t1,t2]
  !
  
  subroutine GetGribTime( ec_treskey, tday, t1, t2, status, tfile, trange, tref )
  
    use GO, only : TDate, Get, Set, wrtgol, NewDate, IncrDate, operator(+), operator(-)
    
    ! --- in/out --------------------------------
    
    character(len=*), intent(in)            ::  ec_treskey
    type(TDate), intent(in)                 ::  tday, t1, t2
    integer, intent(out)                    ::  status
    
    type(TDate), intent(out), optional      ::  tfile
    type(TDate), intent(out), optional      ::  trange(2)
    type(TDate), intent(out), optional      ::  tref
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/GetGribTime'
    
    ! --- local --------------------------------

    integer          ::  hour2, time6(6)
    integer          ::  dd, hh, step

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

    ! set day shift, start hour, and step
    select case ( ec_treskey )

      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! constant field
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case ( 'const' )
      
        ! only tday is usefull ...
        if ( present(tfile ) ) tfile  = tday
        if ( present(trange) ) trange = (/t1,t2/)  ! any, any
        if ( present(tref  ) ) tref   = tday       ! dummy ...
      
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! fg, 3 hourly
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case ( '_fg006up4tr3' )
      
        ! extract end hour:
        call Get( t2, hour=hour2 )

        ! set forecast start time and step given hour:
        select case ( hour2 )
          case ( 00 ) ; dd =  0  ;  hh = 00  ;  step = 6
          case ( 03 ) ; dd =  0  ;  hh = 03  ;  step = 3
          case ( 06 ) ; dd =  0  ;  hh = 06  ;  step = 6
          case ( 09 ) ; dd =  0  ;  hh = 09  ;  step = 3
          case ( 12 ) ; dd =  0  ;  hh = 12  ;  step = 6
          case ( 15 ) ; dd =  0  ;  hh = 15  ;  step = 3
          case ( 18 ) ; dd =  0  ;  hh = 18  ;  step = 6
          case ( 21 ) ; dd =  0  ;  hh = 21  ;  step = 3
          case default
            write (gol,'("unsupported hour :")'); call goErr
            write (gol,'("  hour2     : ",i2)') hour2; call goErr
            write (gol,'("  timesteps : ",a )') ec_treskey; call goErr
            call wrtgol( '  time1     : ', t1 ); call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select
        
        ! file ccyymmdd contains fields for [00,24) :
        if ( present(tfile) ) then
          tfile = t1
          call Set( tfile, hour=00 )
        end if
        
        ! fields valid for [00,24) :
        if ( present(trange) ) then
          trange(1) = t1
          call Set( trange(1), hour=0, min=0, sec=0, mili=0 )
          trange(2) = trange(1) + IncrDate(hour=23,min=59,sec=59,mili=999)
          ! trap last fg time 2000-09-12 00:00 ; data valid forr 00:00 only:
          call Get( tfile, time6=time6 )
          if ( all(time6(1:3)==(/2000,09,12/)) ) trange(2) = trange(1)
        end if
        
        ! reference time = start of forecast
        if ( present(tref) ) then
          call Get( t1, time6=time6 )
          time6(4:6) = 0
          tref = NewDate( time6=time6 ) + IncrDate( day=dd, hour=hh )
        end if

      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! fc, 3 hourly
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case ( '_fc012up2tr3', '_03p03' )
      
        ! extract end hour:
        call Get( t2, hour=hour2 )

        ! set forecast start time and step given hour:
        select case ( hour2 )
          case ( 00 ) ; dd = -1  ;  hh = 12  ;  step = 12
          case ( 03 ) ; dd =  0  ;  hh = 00  ;  step =  3
          case ( 06 ) ; dd =  0  ;  hh = 00  ;  step =  6
          case ( 09 ) ; dd =  0  ;  hh = 00  ;  step =  9
          case ( 12 ) ; dd =  0  ;  hh = 00  ;  step = 12
          case ( 15 ) ; dd =  0  ;  hh = 12  ;  step =  3
          case ( 18 ) ; dd =  0  ;  hh = 12  ;  step =  6
          case ( 21 ) ; dd =  0  ;  hh = 12  ;  step =  9
          case default
            write (gol,'("unsupported hour :")'); call goErr
            write (gol,'("  hour2     : ",i2)') hour2; call goErr
            write (gol,'("  timesteps : ",a )') ec_treskey; call goErr
            call wrtgol( '  time1     : ', t1 ); call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select
        
        ! file ccyymmdd contains fields for (00,24] :
        if ( present(tfile) ) then
          ! current day by default:
          tfile = t2
          call Set( tfile, hour=0, min=0, sec=0, mili=0 ) ! 00:00
          ! trap (..,00:00], this should be previous day:
          call Get( t2, time6=time6 )
          if ( all(time6(4:6)==0) ) tfile = tfile - IncrDate(day=1)
        end if
        
        ! fields valid for (00,24] :
        if ( present(trange) ) then
          trange(1) = t2
          call Set( trange(1), hour=0, min=0, sec=0, mili=0 ) ! 00:00 
          ! trap 00:00, this should be previous day:
          call Get( t2, time6=time6 )
          if ( all(time6(4:6)==0) ) trange(1) = trange(1) - IncrDate(day=1)
          ! complete (00,24]
          trange(2) = trange(1) + IncrDate(day=1)  ! 24:00
          call Set( trange(1), mili=1 ) ! > 00:00 
        end if
        
        ! reference time = start of forecast
        if ( present(tref) ) then
          call Get( t2, time6=time6 )
          time6(4:6) = 0
          tref = NewDate( time6=time6 ) + IncrDate( day=dd, hour=hh )
        end if

      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! fc, 6 hourly
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case ( '_06p06' )
      
        ! extract end hour:
        call Get( t2, hour=hour2 )

        ! set forecast start time and step given hour:
        select case ( hour2 )
          case ( 00 ) ; dd = -1  ;  hh = 12  ;  step = 12
          case ( 06 ) ; dd =  0  ;  hh = 00  ;  step =  6
          case ( 12 ) ; dd =  0  ;  hh = 00  ;  step = 12
          case ( 18 ) ; dd =  0  ;  hh = 12  ;  step =  6
          case default
            write (gol,'("unsupported hour :")'); call goErr
            write (gol,'("  hour2     : ",i2)') hour2; call goErr
            write (gol,'("  timesteps : ",a )') ec_treskey; call goErr
            call wrtgol( '  time1     : ', t1 ); call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select
        
        ! file ccyymmdd contains fields for (00,24] :
        if ( present(tfile) ) then
          ! current day by default:
          tfile = t1
          call Set( tfile, hour=00 )
          ! trap 00:00, this should be previous day:
          call Get( t1, time6=time6 )
          if ( all(time6(4:6)==0) ) tfile = tfile - IncrDate(day=1)
        end if
        
        ! fields valid for (00,24] :
        if ( present(trange) ) then
          trange(1) = t1
          call Set( trange(1), hour=0, min=0, sec=0, mili=0 ) ! 00:00 
          ! trap 00:00, this should be previous day:
          call Get( t1, time6=time6 )
          if ( all(time6(4:6)==0) ) trange(1) = trange(1) - IncrDate(day=1)
          ! complete (00,24]
          trange(2) = trange(1) + IncrDate(day=1)  ! 24:00
          call Set( trange(1), mili=1 ) ! > 00:00 
        end if
        
        ! reference time = start of forecast
        if ( present(tref) ) then
          call Get( t1, time6=time6 )
          time6(4:6) = 0
          tref = NewDate( time6=time6 ) + IncrDate( day=dd, hour=hh )
        end if

      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! analysis, 6 hourly
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case ( 'an_06p06' )
      
        ! extract end hour:
        call Get( t2, hour=hour2 )

        ! set forecast start time and step given hour:
        select case ( hour2 )
          case ( 00 ) ; dd =  0  ;  hh = 00  ;  step = 00
          case ( 06 ) ; dd =  0  ;  hh = 06  ;  step = 00
          case ( 12 ) ; dd =  0  ;  hh = 12  ;  step = 00
          case ( 18 ) ; dd =  0  ;  hh = 18  ;  step = 00
          case default
            write (gol,'("unsupported hour :")'); call goErr
            write (gol,'("  hour2     : ",i2)') hour2; call goErr
            write (gol,'("  timesteps : ",a )') ec_treskey; call goErr
            call wrtgol( '  time1     : ', t1 ); call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select
        
        ! file ccyymmdd contains fields for (00,24] :
        if ( present(tfile) ) then
          ! current day by default:
          tfile = t1
          call Set( tfile, hour=00 )
          ! trap 00:00, this should be previous day:
          call Get( t1, time6=time6 )
          if ( all(time6(4:6)==0) ) tfile = tfile - IncrDate(day=1)
        end if
        
        ! fields valid for (00,24] :
        if ( present(trange) ) then
          trange(1) = t1
          call Set( trange(1), hour=0, min=0, sec=0, mili=0 ) ! 00:00 
          ! trap 00:00, this should be previous day:
          call Get( t1, time6=time6 )
          if ( all(time6(4:6)==0) ) trange(1) = trange(1) - IncrDate(day=1)
          ! complete (00,24]
          trange(2) = trange(1) + IncrDate(day=1)  ! 24:00
          call Set( trange(1), mili=1 ) ! > 00:00 
        end if
        
        ! reference time = start of forecast
        if ( present(tref) ) then
          call Get( t1, time6=time6 )
          time6(4:6) = 0
          tref = NewDate( time6=time6 ) + IncrDate( day=dd, hour=hh )
        end if

      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! ???
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case default
      
        write (gol,'("unsupported time resolution key:")'); call goErr
        write (gol,'("  ",a)') trim(ec_treskey); call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
      
    end select

    ! ok
    status = 0
 
  end subroutine GetGribTime


  ! ***
  
  
  !
  ! Extract time fields of current grib record,
  ! store in mf%tref, mf%t1, mf%t2 
  !
  

  subroutine SetTime( mf, status )

    use GO, only : TDate, NewDate, IncrDate, operator(+), wrtgol
    use file_grib, only : TGribFile, Check, Get
    
    ! --- const -------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/SetTime'
    
    ! --- in/out -------------------------------
    
    type(TMeteoFile_Grib), intent(inout)   ::  mf
    integer, intent(out)                   ::  status
    
    ! --- local -------------------------------
    
    integer          ::  reftime(5), timerange(4)

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

    ! extract time fields from grib record:
    call Get( mf%grib, status, reftime=reftime, timerange=timerange )
    IF_NOTOK_RETURN(status=1)
  
    ! Fill t1 and t2 with the time information; might be equal.
    ! Check time range indicator (WMO code table 5):
    select case ( timerange(4) )
    
      case ( 0, 1 )
        !
        ! 0 = Forecast product valid for reference time + P1 (P1>0),
        !     or uninitialized analysis product for reference time (P1=0)
        !
        ! 1 = Initialized analysis product for reference time (P1=0).
        !
        
        ! fill reference time:
        mf%tref = NewDate( time5=reftime )
        
        ! fill t1 with reftime+timerange;
        ! add P1 in hours; check time unit (WMO code table 4)
        mf%t1 = NewDate( time5=reftime )
        select case ( timerange(1) )
          case ( 1 )  ! hours
            mf%t1 = mf%t1 + IncrDate( hour=timerange(2) )
          case default
            write (gol,'("grib timerange units other than hours not supported yet")'); call goErr
            write (gol,'("  reftime   : ",i4,4i3)') reftime; call goErr
            write (gol,'("  timerange : ",4i3)') timerange; call goErr
            write (gol,'("  file      : ",a)') trim(mf%grib%fname); call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select
        
        ! instant time:
        mf%t2 = mf%t1
        
      case ( 2 )
        !
        ! 2 = Product with a valid time ranging between 
        !      reference time + P1 and reference time + P2
        !
        
        ! fill t1 with reftime+P1;
        ! add P1 in hours; check time unit (WMO code table 4)
        mf%t1 = NewDate( time5=reftime )
        select case ( timerange(1) )
          case ( 1 )  ! hours
            mf%t1 = mf%t1 + IncrDate( hour=timerange(2) )
          case default
            write (gol,'("grib timerange units other than hours not supported yet")'); call goErr
            write (gol,'("  file : ",a)') trim(mf%grib%fname); call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select
        
        ! fill t2 with reftime+P1;
        ! add P1 in hours; check time unit (WMO code table 4)
        mf%t2 = NewDate( time5=reftime )
        select case ( timerange(1) )
          case ( 1 )  ! hours
            mf%t2 = mf%t2 + IncrDate( hour=timerange(3) )
          case default
            write (gol,'("grib timerange units other than hours not supported yet")'); call goErr
            write (gol,'("  file : ",a)') trim(mf%grib%fname); call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select
        
      case default
        write (gol,'("unsupported time range indicator:")'); call goErr
        write (gol,'("  indicator : ",i6)') timerange(4); call goErr
        write (gol,'("  file      : ",a)') trim(mf%grib%fname); call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
        
    end select
    
    ! ok
    status = 0

  end subroutine SetTime
  

  ! ***
    
  
  subroutine CheckTime( mf, tref, t1, t2, status )

    use GO
    use file_grib, only : Check, Get
    
    ! --- const -------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/CheckTime'
    
    ! --- in/out -------------------------------
    
    type(TMeteoFile_Grib), intent(in)      ::  mf
    type(TDate), intent(in)                ::  tref, t1, t2
    integer, intent(out)                   ::  status
    
    ! --- local -------------------------------
    
    integer          ::  year1, year2
    
    ! --- begin -------------------------------

    ! requested year zero ? always ok
    call Get( t1, year=year1 )
    call Get( t2, year=year2 )
    if ( (year1 == 0) .and. (year2 == 0) ) then
      ! requested constant field, always ok
      status = 0; return    ! ok
    end if
    
    if ( year1 == 0 ) then
      ! do not test t1, only t2
      if ( t2 == mf%t2 ) then
        status = 0; return    ! ok
      else
        status = 1; return    ! not ok, try next
      end if
    else if ( year2 == 0 ) then
      ! do not test t2, only t1
      if ( t1 == mf%t1 ) then
        status = 0; return    ! ok
      else
        status = 1; return    ! not ok, try next
      end if
    end if

!    ! interval or instant time
!    if ( t1 < t2 ) then
!    
!      !! time interval: [t1,t2] should be inside [t1,t2]
!      !if ( (t1 >= mf%t1) .and. (t2 <= mf%t2) ) then
!      !  status = 0; return  ! ok
!      !else
!      !  status = 1; return  ! try next
!      !end if
!
!      ! time interval: [t1,t2] should be equal to [t1,t2]
!      if ( (t1 == mf%t1) .and. (t2 == mf%t2) ) then
!        status = 0; return  ! ok
!      else
!        status = 1; return  ! try next
!      end if
!
!    else if ( t1 == t2 ) then
!
!      ! instant time: t2 should match t2
!      if ( t2 == mf%t2 ) then
!        status = 0; return  ! ok
!      else
!        status = 1; return  ! try next
!      end if
!      
!    else
!
!      write (gol,'("t1 should not exceed t2:")'); call goErr
!      call wrtgol( '  t1 : ', t1 ); call goErr
!      call wrtgol( '  t2 : ', t2 ); call goErr
!      write (*,'("ERROR in ",a)') rname; status=3; return
!
!    end if


    ! compare all:
    if ( (tref == mf%tref) .and. (t1 == mf%t1) .and. (t2 == mf%t2) ) then
      status = 0; return  ! ok
    else
      status = 1; return  ! try next
    end if
    
    ! some error ...
    status = 1

  end subroutine CheckTime
  
  
end module tmm_mf_grib