ecearth3/sources/tm5mp/base/tmm_mf_tm5_hdf.F90

!###############################################################################
!
! Input/output of meteofiles produced by TM5 : HDF version.
!
!###############################################################################
!
#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_tm5_hdf

  use GO      , only : gol, goErr, goPr, goBug
  use GO      , only : TDate
  use file_hdf, only : THdfFile, TSds

  implicit none
  
  ! --- in/out ----------------------------
  
  private
  
  public  ::  TMeteoFile_tm5_hdf
  public  ::  TMM_MF_TM5_HDF_Init, TMM_MF_TM5_HDF_Done
  public  ::  Init, Done
  public  ::  Get
  public  ::  ReadRecord
  public  ::  WriteRecord
  
  ! --- const ------------------------------
  
  character(len=*), parameter  ::  mname = 'tmm_mf_tm5_hdf'
  
  ! ~~~ output keys and defaults
  
  ! current format version
  character(len=*), parameter  ::  output_format = 'tm5-hdf'
  
  ! integer/real kinds 
  integer, parameter   ::  iknd    = 4        ! integer
  integer, parameter   ::  rknd    = 8        ! real
  integer, parameter   ::  iknd_ds = 2        ! integer for big ds
  integer, parameter   ::  rknd_ds = 4        ! real    for big ds

  ! z compression
  character(len=*), parameter  ::  compression = 'deflate'
  integer, parameter           ::  deflate_level = 6


  !--- type ---------------------------------
  
  type TMeteoFile_tm5_hdf
    ! input/output ?
    character(len=1)       ::  io
    !
    ! field collection
    !
    character(len=256)     ::  fname
    type(THdfFile)         ::  hdf
    character(len=256)     ::  paramkeys
    type(TDate)            ::  trange(2)
    !
    ! file
    !
    type(TSds)             ::  sds
    logical                ::  selected
    character(len=16)      ::  paramkey
    !
    ! pw specials
    !
    logical                ::  mfw_redir
    !
    ! surface stress specials
    !
    logical                ::  sstr_to_ewss_nsss
    !
    ! surface pressure field
    !
    logical                ::  spm_load      ! load spm ?
    logical                ::  spm_selected  ! data set filled ?
    type(TSds)             ::  spm_sds       ! data set for sp
    logical                ::  spm_incl      ! record included in main file ?
    logical                ::  spm_extr      ! record in external spm file ?
    character(len=256)     ::  spm_fname
    type(THdfFile)         ::  spm_hdf
    logical                ::  spm_n_to_uv
    !
    ! output
    !
    logical                ::  output_initialised
    integer                ::  output_nrec
    integer                ::  output_ntrec
    character(len=20)      ::  output_names(20)
    integer                ::  output_nname
    !
    ! adhoc ...
    integer                ::  fixyear
    logical                ::  qad
  end type TMeteoFile_tm5_hdf
  
  
  ! --- 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

  interface WriteRecord
    module procedure mf_WriteRecord_2d
    module procedure mf_WriteRecord_3d
  end interface
  
  
  ! --- var --------------------------------------
  
  ! timer id's:
  integer           ::  itim_readrecord


contains


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


  subroutine TMM_MF_TM5_HDF_Init( rcf, status )
  
    use GO, only : TrcFile
    use GO, only : GO_Timer_Def

    ! --- in/out ---------------------------------
    
    type(TRcFile), intent(in)           ::  rcf
    integer, intent(out)                ::  status
    
    ! --- const ----------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/TMM_MF_TM5_HDF_Init'
    
    ! --- local ----------------------------------
    
    ! --- begin ----------------------------------
    
    ! define timers:
    call GO_Timer_Def( itim_readrecord, 'tmm tm5-hdf readrecord', status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0
    
  end subroutine TMM_MF_TM5_HDF_Init


  ! ***


  subroutine TMM_MF_TM5_HDF_Done( status )
  
    ! --- in/out ---------------------------------
    
    integer, intent(out)                ::  status
    
    ! --- const ----------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/TMM_MF_TM5_HDF_Done'
    
    ! --- local ----------------------------------
    
    ! --- begin ----------------------------------
    
    ! ok
    status = 0
    
  end subroutine TMM_MF_TM5_HDF_Done


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


  subroutine mf_Init( mf, io, dir, archivekeys, paramkey, &
                               tref, t1, t2, status )
    
    use GO, only : TDate, Set, Get, NewDate, AnyDate, IsAnyDate
    use GO, only : rTotal, operator(-), operator(>=)
    use GO, only : goVarValue, goWriteKeyNum, goReplace

    ! --- in/out ----------------------------
    
    type(TMeteoFile_tm5_hdf), intent(out)   ::  mf
    character(len=1), intent(in)        ::  io
    character(len=*), intent(in)        ::  dir
    character(len=*), intent(in)        ::  archivekeys
    character(len=*), intent(in)        ::  paramkey
    type(TDate), intent(in)             ::  tref, t1, t2
    integer, intent(out)                ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_Init'
    
    ! --- local --------------------------------
    
    character(len=64)     ::  mf_mdir, mf_tres
    character(len=64)     ::  mf_class, mf_type, mf_grid
    character(len=1)      ::  mf_psep, mf_nsep

    character(len=64)     ::  mf_filekey
    
    character(len=4)      ::  mf_fckey
    
    type(TDate)           ::  tfile
    integer               ::  ccyy, mm, dd, dh
    type(TDate)           ::  tc

    logical               ::  exist
    
    ! --- begin --------------------------------

    ! store i/o :
    mf%io = io
    
    ! default flags:
    mf%sstr_to_ewss_nsss = .false.
    mf%mfw_redir         = .false.

    !
    ! extract fields from archivekey :
    !   mdir=ec-fg_3h-ml60-glb3x2;tres=_21p06
    !
    mf%fixyear  = -1
    call goVarValue( archivekeys, ';', 'fixyear', '=', mf%fixyear, status )
    if (status>0) then; TRACEBACK; status=1; return; end if

    mf_mdir  = 'no_mdir'
    call goVarValue( archivekeys, ';', 'mdir', '=', mf_mdir, status )
    if (status>0) then; TRACEBACK; status=1; return; end if
    !
    mf_tres  = 'no_tres'
    call goVarValue( archivekeys, ';', 'tres', '=', mf_tres, status )
    if (status>0) then; TRACEBACK; status=1; return; end if
    !
    ! path seperation character:
    mf_psep  = '/'
    call goVarValue( archivekeys, ';', 'pathsep', '=', mf_psep, status )
    if (status>0) then; TRACEBACK; status=1; return; end if
    !
    ! name seperation character:
    mf_nsep  = '-'
    call goVarValue( archivekeys, ';', 'namesep', '=', mf_nsep, status )
    if (status>0) then; TRACEBACK; status=1; return; end if
    !
    ! adhoc flag
    mf%sstr_to_ewss_nsss = .false.
    call goVarValue( archivekeys, ';', 'sstr', '=', mf%sstr_to_ewss_nsss, status )
    if (status>0) then; TRACEBACK; status=1; return; end if


    !
    ! main file
    !

    ! by default, no surface pressure stuff ...
    mf%spm_load    = .false.
    mf%spm_incl    = .false.   ! no spm fields included in hdf file
    mf%spm_extr    = .false.   ! no extra spm file
    mf%spm_n_to_uv = .false.   ! no interpolation from 'n' to 'u' or 'v'

    ! * set mf_filekey (uvsp,t,etc) and parmeters:
    select case ( paramkey )
      case ( 'sp', 'tsp' )
        mf_filekey = paramkey
        mf%paramkeys = '-'//trim(paramkey)//'-'
        mf%spm_load  = .false.
      case ( 'spm' )
        mf_filekey = 'spm'
        mf_tres = '_00p06'
        mf%paramkeys = '-spm-'
      case ( 'mfu', 'mfv' )
        mf_filekey = 'mfuv'
        mf%paramkeys = '-mfu-mfv-'
        mf%spm_load  = .true.
      case ( 'mfw' )
        mf_filekey = 'mfw'
        mf%paramkeys = '-mfw-'
        mf%spm_load  = .true.
      case ( 'T' )
        mf_filekey   = 't'
        mf%paramkeys = '-T-'
        mf%spm_load  = .true.
      case ( 'Q' )
        mf_filekey   = 'q'
        mf%paramkeys = '-Q-'
        mf%spm_load  = .true.
      case ( 'CLWC', 'CIWC', 'CC', 'CCO', 'CCU', &
             'clwc', 'ciwc', 'cc', 'cco', 'ccu'  )
        mf_filekey   = 'cld'
        mf%paramkeys = '-CLWC-CIWC-CC-CCO-CCU-'
        mf%spm_load  = .true.
      case ( 'eu', 'ed', 'du', 'dd', 'cloud_base', 'cloud_top', 'cloud_lfs' )
        mf_filekey   = 'sub'
        mf%paramkeys = '-eu-ed-du-dd-cloud_base-cloud_top-cloud_lfs-'
        mf%spm_load  = .true.
      case ( 'K' )
        mf_filekey   = 'diffus'
        mf%paramkeys = '-K-'
        mf%spm_load  = .true.
      ! o constant fields
      case ( 'oro', 'lsm' )
        mf_tres      = 'constant'
        mf_filekey   = trim(paramkey)
        mf%paramkeys = '-'//trim(paramkey)//'-'
        mf%spm_load  = .false.
      ! o monthly fields:
      case ( 'srols' )
        mf_tres      = 'month'
        mf_filekey   = trim(paramkey)
        mf%paramkeys = '-'//trim(paramkey)//'-'
        mf%spm_load  = .false.
      ! o vegetation fields
      case ( 'cvl', 'cvh', &
             'tv01', 'tv02', 'tv03', 'tv04', 'tv05', &
             'tv06', 'tv07', 'tv08', 'tv09', 'tv10', &
             'tv11', 'tv12', 'tv13', 'tv14', 'tv15', &
             'tv16', 'tv17', 'tv18', 'tv19', 'tv20' )
        mf_filekey = 'veg'
        mf%paramkeys= '-'//&
             'cvl-cvh-'//&
             'tv01-tv02-tv03-tv04-tv05-'//&
             'tv06-tv07-tv08-tv09-tv10-'//&
             'tv11-tv12-tv13-tv14-tv15-'//&
             'tv16-tv17-tv18-tv19-tv20-'
        mf%spm_load  = .false.
      ! o each surface file in a seperate file; 
      !   might be collected in surf_*.tar files afterwards
      case ( 'sr', 'srmer', &
             'swvl1', &
             'albedo', 'lsrh', 'ci', 'g10m', 'fg10','u10m', 'v10m', 'sd', &
             'lsp', 'cp', 'sf', 'sshf', 'slhf', 'blh', &
             't2m', 'd2m', 'ssr', 'ssrd', 'str', 'strd', &
             'sstr', 'src', 'raero', 'ustar', &
             'sst', 'sps', 'skt', 'ch4fire' )
        mf_filekey = trim(paramkey)
        mf%paramkeys = '-'//trim(paramkey)//'-'
        mf%spm_load  = .false.
      case ( 'ewss', 'nsss' )
        if ( mf%sstr_to_ewss_nsss ) then
          mf_filekey   = 'sstr'
          mf%paramkeys = '-sstr-'
        else
          mf_filekey   = trim(paramkey)
          mf%paramkeys = '-'//trim(paramkey)//'-'
        end if
        mf%spm_load  = .false.
      case default
        write (gol,'("unsupported paramkey `",a,"`")') paramkey; call goErr
        TRACEBACK; status=1; return
    end select
    
    ! eventually replace by other value:
    call goVarValue( archivekeys, ';', 'filekey', '=', mf_filekey, status )
    if (status>0) then; TRACEBACK; status=1; return; end if

    ! convert input times to file name times:
    call GetTime( mf_filekey, mf_tres, tref, t1, t2, status, &
                       tfile=tfile, trange=mf%trange )
    IF_NOTOK_RETURN(status=1)

    ! adhoc: fixed year ?
    if ( mf%fixyear > 0 ) call Set( tfile, year=mf%fixyear )

    ! extract time values:
    call Get( tfile, year=ccyy, month=mm, day=dd )

    ! replace time values in directory name:
    call goReplace( mf_mdir, '<yyyy>', '(i4.4)', ccyy, status )
    IF_NOTOK_RETURN(status=1)
    call goReplace( mf_mdir, '<mm>', '(i2.2)', mm, status )
    IF_NOTOK_RETURN(status=1)

    ! special data set: trap change from fg to fc data:
    if ( mf_tres == '_fg006up4tr3' ) then
      tc = NewDate( 2000, 09, 12 )
      if ( tfile >= tc ) mf_tres = '_fc012up2tr3'
    end if

    ! forecast key: '', 'f1', .., 'f10' ;
    ! no key for constant files (t1 and t2 are anydate)
    ! or month files (t1 is begin of month thus probably < tref)
    mf_fckey = ''
    if ( (.not. IsAnyDate(t1)) .and. (t1 >= tref) ) then
      dh = floor( rTotal( t1 - tref, 'day' ) )
      if ( dh > 0 ) call goWriteKeyNum( mf_fckey, 'f', dh )
    end if

    ! create file name:
    !   dir / ec_od-ml60-T159 - oro.hdf
    !   dir / ec_od-ml60-T159 - T_20000101_fg006up4tr3.hdf
    select case ( mf_tres )
      case ( 'constant' )
        ! filename without date:
        write (mf%fname,'(a,a,a,a,a,".hdf")') trim(dir), mf_psep, trim(mf_mdir), mf_nsep, trim(mf_filekey)
      case ( 'month' )
        ! filename without day and forecast key:
        write (mf%fname,'(a,a,a,a,a,"_",i4.4,i2.2,a)') &
                trim(dir), mf_psep, trim(mf_mdir), mf_nsep, trim(mf_filekey), ccyy, mm, '.hdf'
      case default
        ! filename including date:
        write (mf%fname,'(a,a,a,a,a,"_",i4.4,2i2.2,a,a,a)') &
                trim(dir), mf_psep, trim(mf_mdir), mf_nsep, &
                trim(mf_filekey), ccyy, mm, dd, trim(mf_fckey), trim(mf_tres), '.hdf'
        ! trap request for sps ; might be in sp file ...
        if ( trim(mf_filekey) == 'sps' ) then
          inquire( file=trim(mf%fname), exist=exist )
          if ( .not. exist ) then
            write (mf%fname,'(a,a,a,a,a,"_",i4.4,2i2.2,a,a,a)') &
                trim(dir), mf_psep, trim(mf_mdir), mf_nsep, &
                'sp', ccyy, mm, dd, trim(mf_fckey), trim(mf_tres), '.hdf'
          end if
        end if
    end select

    ! spm is included in files ...
    mf%spm_incl  = mf%spm_load    ! used to be external, now included
    mf%spm_fname = mf%fname       ! same file
    mf%spm_extr  = .false.        ! not external for sure ...

    ! in case of output, not initialised yet ...
    mf%output_initialised = .false.
    mf%output_nname = 0
    ! number of expected time records in a file:
    call GetTime( mf_filekey, mf_tres, tref, t1, t2, status, nrec=mf%output_ntrec )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0
 
  end subroutine mf_Init
  
  
  ! ***

  
  subroutine mf_Done( mf, status )
  
    !use file_hdf, only : Done
    
    ! --- in/out ------------------------------------

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

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_Done'
    
    ! --- begin -------------------------------------
    
    ! files have been closed in ReadRecord/WriteRecord
    
    ! ok
    status = 0
    
  end subroutine mf_Done
  


  ! ***

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

    type(TDate), intent(out), optional        ::  trange1, trange2
    character(len=*), intent(out), optional   ::  paramkeys
    character(len=*), intent(out), optional   ::  filename
    
    ! --- 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)
    
    ! parameter names:
    if ( present(paramkeys) ) paramkeys = trim(mf%paramkeys)
    
    ! file name:
    if ( present(filename) ) filename = trim(mf%fname)

    ! ok
    status = 0
 
  end subroutine mf_Get
  


  ! ******************************************************************
  ! ***
  ! *** time range, parameters, file names
  ! ***
  ! ******************************************************************


  !
  ! Return time parameters:
  !  o tfile   :  date in filename
  !  o trange  :  time interval covered by fields in file
  !
  
  subroutine GetTime( filekey, tres, tref, t1, t2, status, tfile, trange, nrec )
  
    use GO, only : TDate, NewDate, AnyDate, Get, Set, wrtgol, IncrDate, IsAnyDate
    use GO, only : operator(<), operator(+), operator(-), rTotal
    
    ! --- in/out --------------------------------
    
    character(len=*), intent(in)            ::  filekey
    character(len=*), intent(in)            ::  tres
    type(TDate), intent(in)                 ::  tref, t1, t2
    integer, intent(out)                    ::  status
    
    type(TDate), intent(out), optional      ::  tfile
    type(TDate), intent(out), optional      ::  trange(2)
    integer, intent(out), optional          ::  nrec
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/GetTime'
    
    ! --- local --------------------------------

    integer          ::  year, month
    integer          ::  hour1, time6(6)
    integer          ::  dd, hh, step
    logical          ::  interval
    real             ::  dhr

    ! --- begin --------------------------------
    
    ! set day shift, start hour, and step
    select case ( tres )

      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! tmpp [21,21]
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case ( '_21p06', '_21p03', '_av21' )
      
        ! routine is called with  tref,t1,t2:
        !    t1,t1,t2
        !    t1,any,any    (oro and other constant fields)
        ! thus use tref to construct the file times
        if ( IsAnyDate(t1) ) then
          call Get( tref, hour=hour1 )
          interval = .false.
        else
          call Get( t1, hour=hour1 )
          interval = t1 < t2
        end if
      
        ! file ccyymmdd contains fields for (21,21];
        ! only uvsp is valid for [21,21] since it contains surface pressure for 21:00
        if ( present(tfile) ) then
          tfile = tref
          call Set( tfile, hour=0, min=0, sec=0 )
          if ( (hour1 > 21) .or. ((interval .or. filekey=='uvsp') .and. hour1==21) ) then
            tfile = tfile + IncrDate(day=1)
          end if
       end if
        
        ! fields by default valid for (21,21];
        ! only uvsp is valid for [21,21] since it contains surface pressure for 21:00
        if ( present(trange) ) then
          trange(1) = tref
          call Set( trange(1), hour=0, min=0, sec=0 )  !  00:00 today
          if ( (hour1 > 21) .or. ((interval .or. filekey=='uvsp') .and. hour1==21) ) then
            trange(1) = trange(1) + IncrDate(day=1)
          end if
          trange(1) = trange(1) - IncrDate(hour=3)     !  previous 21:00
          trange(2) = trange(1) + IncrDate(day=1)      !  next     21:00
          ! boundary not included in most cases:
          if ( filekey /= 'uvsp' ) trange(1) = trange(1) + IncrDate(mili=1)
        end if
        
        ! number of records in file:
        if ( present(nrec) ) then
          select case ( tres )
            case ( '_21p06' ) ; nrec = 24/6
            case ( '_21p03' ) ; nrec = 24/3
            case ( '_av21'  ) ; nrec = 24/24
            case default
              write (gol,'("unsupported tres for setting nrec : ",a)') tres; call goErr
              TRACEBACK; status=1; return
          end select
        end if
        
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! tm5 constant
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case ( 'constant' )
      
        ! no date in filename ...
        if ( present(tfile) ) tfile = AnyDate()
        
        ! fields always valid ...
        if ( present(trange) ) then
          trange(1) = AnyDate()
          trange(2) = AnyDate()
        end if
        
        ! only one output record in constant file:
        if ( present(nrec) ) nrec = 1
        
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! tm5 monthly file
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case ( 'month' )
      
        ! file ccyymmdd contains fields for this month:
        if ( present(tfile) ) then
          call Get( t1, year=year, month=month )
          tfile = NewDate( year=year, month=month, day=1 )
        end if
        
        ! field valid from begin to end of month:
        if ( present(trange) ) then
          call Get( t1, year=year, month=month )
          trange(1) = NewDate( year=year, month=month, day=1, hour=00  )
          month = month + 1
          if ( month > 12 ) then
            year = year + 1
            month = 1
          end if
          trange(2) = NewDate( year=year, month=month, day=1, hour=00  )
        end if
        
        ! only one output record in month file:
        if ( present(nrec) ) nrec = 1
        
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! tm5 [00,24]
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case ( '_00p06', '_00p03', '_an0tr6', '_fg006up4tr3', '_fc012up2tr3', '_00p01', '_00p24' )
      
        ! file ccyymmdd contains fields for [00,24) :
        if ( present(tfile) ) then
          tfile = t1
          call Set( tfile, hour=0, min=0, sec=0 )
        end if
        
        ! fields valid for [00,24) :
        if ( present(trange) ) then
          trange(1) = t1
          call Set( trange(1), hour=0, min=0, sec=0 )  !  00:00 today
          trange(2) = trange(1) + IncrDate(hour=24) - IncrDate(mili=1)
        end if
        
        ! number of records in file:
        if ( present(nrec) ) then
          select case ( tres )
            case ( '_00p06'  ) ; nrec = 24/6
            case ( '_an0tr6' ) ; nrec = 24/6
            case ( '_00p03', '_fg006up4tr3', '_fc012up2tr3' )
              ! by default: 3 hourly files
              ! for forecasts after 12+72, only 6 hourly available:
              !  f0  [ 00, 24)  : 00 03 06 09 12 15 18 21   :  nrec=8
              !  f1  [ 24, 48)  : 00 03 06 09 12 15 18 21   :  nrec=8
              !  f2  [ 48, 72)  : 00 03 06 09 12 15 18 21   :  nrec=8
              !  f3  [ 72, 96)  : 00 03 06 09 12    18      :  nrec=6
              !  f4  [ 96,120)  : 00    06    12    18      :  nrec=4
              !   :
              !  f9  [192,216)  : 00    06    12    18      :  nrec=4
              !  f10 [216,240)  : 00    06    12            :  nrec=3
              dhr = rTotal( t1 - tref, 'hour' )
              if ( dhr < 72.0 ) then
                nrec = 8
              else if ( dhr < 96 ) then
                nrec = 6
              else if ( dhr < 216 ) then
                nrec = 4
              else
                nrec = 3
              end if
            case ( '_00p01'       ) ; nrec = 24/1
            case ( '_00p24'       ) ; nrec = 24/24
            case default
              write (gol,'("unsupported tres for setting nrec : ",a)') tres; call goErr
              TRACEBACK; status=1; return
          end select
        end if
        
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      ! ???
      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      case default
      
        write (gol,'("unsupported time resolution key:")'); call goErr
        write (gol,'("  ",a)') trim(tres); call goErr
        TRACEBACK; status=1; return
      
    end select
      
  
    ! ok
    status = 0
 
  end subroutine GetTime


  ! ******************************************************************
  ! ***
  ! *** input
  ! ***
  ! ******************************************************************

  
  subroutine mf_SelectRecord( mf, paramkey, t1, t2, status )
  
    use GO, only : TDate, NewDate, IncrDate, wrtgol, Set, Get, rTotal
    use GO, only : operator(-), operator(+), operator(/), operator(==), operator(/=)
    use GO, only : operator(>=), operator(<=), operator(>)
    use file_hdf, only : GetInfo
    use file_hdf, only : TSds, Select, CheckInfo, CheckAttribute, ReadAttribute

    ! --- in/out --------------------------------
  
    type(TMeteoFile_tm5_hdf), intent(inout)     ::  mf
    character(len=*), intent(in)            ::  paramkey
    type(TDate), intent(in)                 ::  t1, t2
    integer, intent(out)                    ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_SelectRecord'
    
    ! --- local -------------------------------
    
    integer            ::  nsds, isds
    integer            ::  time1(6), time2(6)
    type(TDate)        ::  tmid, treq, thelp
    type(TDate)        ::  t_spm
    integer            ::  hour, dhour

    integer            ::  status1, status2
    integer            ::  time1s(6), time2s(6)
    type(TDate)        ::  tim1s, tim2s

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

    mf%selected = .false.

    ! initial no record is found ...    
    status = 1

    ! number of data sets:  
    call GetInfo( mf%hdf, status, num_datasets=nsds )
    IF_NOTOK_RETURN(status=1)

    ! loop over all data sets
    do isds = 1, nsds

      call Select( mf%sds, mf%hdf, isds-1, status )
      IF_NOTOK_RETURN(status=1)

      ! correct param ?
      status=-1; call CheckInfo( mf%sds, status, name=paramkey )
      if (status>0) then; TRACEBACK; status=1; return; end if
      
      ! adhoc: some data sets have other name ...
      if ( status < 0 ) then
        select case ( paramkey )
          case ( 'sp', 'spm' )
            status=-1; call CheckInfo( mf%sds, status, name='ps' )
          case ( 'sps' )
            status=-1; call CheckInfo( mf%sds, status, name='sp' )
          case ( 'mfu' )
            status=-1; call CheckInfo( mf%sds, status, name='pu' )
          case ( 'mfv' )
            status=-1; call CheckInfo( mf%sds, status, name='pv' )
          case ( 'mfw' )
            status=-1; call CheckInfo( mf%sds, status, name='pw' )
          case ( 'PVo' )
            status=-1; call CheckInfo( mf%sds, status, name='PV' )
          case ( 'sr' )
            status=-1; call CheckInfo( mf%sds, status, name='sr_ecm' )
          case ( 'srols' )
            status=-1; call CheckInfo( mf%sds, status, name='sr_ols' )
          case ( 'srmer' )
            status=-1; call CheckInfo( mf%sds, status, name='sr_mer' )
          case ( 'albedo' )
            status=-1; call CheckInfo( mf%sds, status, name='al' )
          case ( 'ewss', 'nsss' )
            if ( mf%sstr_to_ewss_nsss ) then
              status=-1; call CheckInfo( mf%sds, status, name='sstr' )
            end if
        end select
        if (status>0) then; TRACEBACK; status=1; return; end if
        ! try other data set ?
        if ( status < 0 ) cycle
      end if
      
      ! correct time ?
      select case ( paramkey )
        case ( 'oro' )
          ! constant field, skip time check
        case ( 'srols' )
          ! monthly field:  
          !   in TMPP files: daily samples, valid for [21,21]
          !   in TM5  files: monthly value
          ! read time2 :
          call ReadAttribute( mf%sds, 'time2', time2s, status )
          IF_NOTOK_RETURN(status=1)
          tim2s = NewDate( time6=time2s )
          ! [t1,t2] covers month; time2 should be in this month:
          if ( (t1 >= tim2s) .and. (tim2s <= t2) ) then
            ! ok
            status = 0
          end if
        case default
          ! extract time arrays:
          call Get( t1, time6=time1 )
          call Get( t2, time6=time2 )
          !
          ! replace year ?
          if ( mf%fixyear > 0 ) then
            ! year is valid for time interval:
            !   [ fixyear/01/01 00:00 , fixyear/12/31 24:00 )
            if ( t1 == t2 ) then
              ! instant time
              time1(1) = mf%fixyear
              time2    = time1
            else
              ! time interval
              time1(1) = mf%fixyear
              thelp = NewDate( time6=time1 ) + (t2-t1)
              call Get( thelp, time6=time2 )
            endif
          end if
          !
          if ( mf%qad ) then
            ! test if times match exactly:
            status1=-1 ; call CheckAttribute( mf%sds, 'time1', time1, status1 )
            if (status1>0) then; TRACEBACK; status=1; return; end if
            status2=-1 ; call CheckAttribute( mf%sds, 'time2', time2, status2 )
            if (status2>0) then; TRACEBACK; status=1; return; end if
            status = status1 + status2
            ! times do not match ?
            if ( status < 0 ) then
              ! read time1 and time2 :
              call ReadAttribute( mf%sds, 'time1', time1s, status )
              IF_NOTOK_RETURN(status=1)
              call ReadAttribute( mf%sds, 'time2', time2s, status )
              IF_NOTOK_RETURN(status=1)
              ! hours should match ...
              if ( (time1s(4) /= time1(4)) .or. (time2s(4) /= time2(4)) ) then
                status = -1
                cycle
              end if
              ! warning ...
              write (gol,'("WARNING - weak time check passed:")'); call goPr
              write (gol,'("WARNING -   parameter : ",a)') trim(paramkey); call goPr
              call wrtgol( 'WARNING -   t1        : ', t1 ); call goPr
              call wrtgol( 'WARNING -   t2        : ', t2 ); call goPr
              write (gol,'("WARNING -   file      : ",a)') trim(mf%fname); call goPr
              ! ok
              status = 0
            end if
            ! weak time check passed
          else
            !! tmpp fields valid for [21,21] should be read as valid for [00,24] ...
            !if ( rTotal(t2-t1,'hour') == 24.0 ) then
            !  tim1s = t1 + IncrDate(hour=3)
            !  tim2s = t2 + IncrDate(hour=3)
            !  call Get( tim1s, time6=time1 )
            !  call Get( tim2s, time6=time2 )
            !end if
            ! first try wether times exactely match;
            status1=-1 ; call CheckAttribute( mf%sds, 'time1', time1, status1 )
            if (status1>0) then; TRACEBACK; status=1; return; end if
            status2=-1 ; call CheckAttribute( mf%sds, 'time2', time2, status2 )
            if (status2>0) then; TRACEBACK; status=1; return; end if
            status = status1 + status2
            ! try mid time ?
            if ( status < 0 ) then
              if ( t1 == t2 ) then
                ! instant time requested; check wether this is mid of times in file:
                call ReadAttribute( mf%sds, 'time1', time1s, status )
                IF_NOTOK_RETURN(status=1)
                tim1s = NewDate( time6=time1s )
                call ReadAttribute( mf%sds, 'time2', time2s, status )
                IF_NOTOK_RETURN(status=1)
                tim2s = NewDate( time6=time2s )
                ! mid of times in file:
                tmid = tim1s + (tim2s-tim1s)/2
                ! requested time:
                treq = t1
                if ( mf%fixyear > 0 ) call Set( treq, year=mf%fixyear )
                ! not ok ? then try again
                if ( treq /= tmid ) then
                  status = -1
                  cycle
                end if
              else
                ! interval requested; check wether time in file is mid time:
                tmid = t1 + (t2-t1)/2
                if ( mf%fixyear > 0 ) call Set( tmid, year=mf%fixyear )
                call Get( tmid, time6=time1 )
                status1=-1 ; call CheckAttribute( mf%sds, 'time1', time1, status1 )
                if (status1>0) then; TRACEBACK; status=1; return; end if
                status2=-1 ; call CheckAttribute( mf%sds, 'time2', time1, status2 )
                if (status2>0) then; TRACEBACK; status=1; return; end if
                status = status1 + status2
                if ( status < 0 ) cycle
              end if
            end if
            ! time check passed
          end if
      end select
      
      ! found!
      exit
    end do

    ! not found ?
    if ( status /= 0 ) then
      write (gol,'("Unable to locate field in hdf file:")'); call goErr
      write (gol,'("  parameter : ",a)') trim(paramkey); call goErr
      call wrtgol( '  t1        : ', t1 ); call goErr
      call wrtgol( '  t2        : ', t2 ); call goErr
      call wrtgol( '  mf%tr(1)  : ', mf%trange(1) ); call goErr
      call wrtgol( '  mf%tr(2)  : ', mf%trange(2) ); call goErr
      write (gol,'("  file      : ",a)') trim(mf%fname); call goErr
      TRACEBACK; status = 1; return
    end if
    
    ! store paramkey
    mf%paramkey = paramkey
    
    ! ok
    mf%selected = .true.
    
    
    ! ~~~ surface pressure ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    if ( mf%spm_load ) then
    
      mf%spm_selected = .false.

      if ( paramkey /= 'sp' ) then

        ! only hours 00, 06, 12, 18
        dhour = nint(rTotal(t2-t1,'hour'))
        select case ( dhour )
          case ( 0 )
            t_spm = t1
          case ( 3 )
            call Get( t1, hour=hour )
            if ( modulo(hour,6) == 0 ) then
              t_spm = t1
            else
              t_spm = t2
            end if
          case ( 6 )
            t_spm = t1 + IncrDate( hour=3 )
          case default
            write (gol,'("do not know how to form time for mid surface pressure:")'); call goErr
            call wrtgol( '  t1 : ', t1 ); call goErr
            call wrtgol( '  t2 : ', t2 ); call goErr
            TRACEBACK; status = 1; return
        end select

        ! replace year ?
        if ( mf%fixyear > 0 ) then
          call Set( t_spm, year=mf%fixyear )
        end if

        ! initial no record is found ...    
        status = 1

        ! number of data sets:
        if ( mf%spm_extr ) then
          call GetInfo( mf%spm_hdf, status, num_datasets=nsds )
        else
          call GetInfo( mf%hdf, status, num_datasets=nsds )
        end if
        IF_NOTOK_RETURN(status=1)

        ! loop over all data sets:
        do isds = 1, nsds
        
          ! select data set:
          if ( mf%spm_extr ) then
            call Select( mf%spm_sds, mf%spm_hdf, isds-1, status )
          else
            call Select( mf%spm_sds, mf%hdf, isds-1, status )
          end if
          IF_NOTOK_RETURN(status=1)

          ! correct param ?
          status=-1; call CheckInfo( mf%spm_sds, status, name='ps' )
          if (status>0) then; TRACEBACK; status=1; return; end if
          
          ! not found ? also try other names:
          if ( status < 0) then
            select case ( paramkey )
              case ( 'mfu', 'pu' )
                status=-1; call CheckInfo( mf%spm_sds, status, name='spu' )
              case ( 'mfv', 'pv' )
                status=-1; call CheckInfo( mf%spm_sds, status, name='spv' )
              case default
                status=-1; call CheckInfo( mf%spm_sds, status, name='sp' )
            end select
            if (status>0) then; TRACEBACK; status=1; return; end if
            ! try next ?
            if ( status < 0 ) cycle
          end if

          ! correct time ?
          call Get( t1, time6=time1 )
          status=-1; call CheckAttribute( mf%spm_sds, 'time1', time1, status )
          if (status>0) then; TRACEBACK; status=1; return; end if
          if (status==0) then
            call Get( t2, time6=time1 )
            status=-1; call CheckAttribute( mf%spm_sds, 'time2', time2, status )
            if (status>0) then; TRACEBACK; status=1; return; end if
          else
            ! try special spm times
            call Get( t_spm, time6=time1 )
            status=-1; call CheckAttribute( mf%spm_sds, 'time1', time1, status )
            if (status<0) cycle
            if (status>0) then; TRACEBACK; status=1; return; end if
            ! try special spm times
            call Get( t_spm, time6=time2 )
            status=-1; call CheckAttribute( mf%spm_sds, 'time2', time2, status )
            if (status<0) cycle
            if (status>0) then; TRACEBACK; status=1; return; end if
          end if

          ! found!
          exit
        end do

        ! not found ?
        if ( status /= 0 ) then
          write (gol,'("Unable to locate surface pressure field in hdf file:")'); call goErr
          call wrtgol( '  t1 - t2 : ', t1, ' - ', t2 ); call goErr
          call wrtgol( '  t_spm   : ', t_spm ); call goErr
          write (gol,'("  file  : ",a)') trim(mf%spm_fname); call goErr
          TRACEBACK; status = 1; return
        end if

        ! ok
        mf%spm_selected = .true.
        
      end if
    
    end if
      
    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    ! ok
    status = 0

  end subroutine mf_SelectRecord

  
  ! ***
  

  !
  ! initialiase grid info from sds
  !
  
  subroutine lli_Init_mf( lli, nuv, mf, status )
  
    use file_hdf, only : ReadAttribute, CheckAttribute
    use Grid, only : TllGridInfo, Init

    ! --- in/out ----------------------------------
    
    type(TllGridInfo), intent(inout)       ::  lli
    character(len=1), intent(in)         ::  nuv
    type(TMeteoFile_tm5_hdf), intent(in)     ::  mf
    integer, intent(out)                 ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/lli_Init_mf'
    
    ! --- local -----------------------------------
    
    real           ::  lon_deg, dlon_deg
    integer        ::  nlon
    real           ::  lat_deg, dlat_deg
    integer        ::  nlat
    integer        ::  stat
    
    ! --- begin ------------------------------------

    ! record not selected ? then return
    if ( .not. mf%selected ) then
      write (gol,'("no record selected ...")'); call goErr
      TRACEBACK; status = 1; return
    end if
    
    ! check ...    
    call CheckAttribute( mf%sds, 'gridtype' , 'll', status )
    if ( status /= 0 ) then
      write (gol,'("sds does not seem to contain ll grid")'); call goErr
      TRACEBACK; status=1; return
    end if

    ! extract grid position parameters from sds:
    call ReadAttribute( mf%sds, 'lon_first', lon_deg , status )
    IF_NOTOK_RETURN(status=1)
    call ReadAttribute( mf%sds, 'lon_inc  ', dlon_deg, status )
    IF_NOTOK_RETURN(status=1)
    call ReadAttribute( mf%sds, 'lon_n    ', nlon    , status )
    IF_NOTOK_RETURN(status=1)
    call ReadAttribute( mf%sds, 'lat_first', lat_deg , status )
    IF_NOTOK_RETURN(status=1)
    call ReadAttribute( mf%sds, 'lat_inc  ', dlat_deg, status )
    IF_NOTOK_RETURN(status=1)
    call ReadAttribute( mf%sds, 'lat_n    ', nlat    , status )
    IF_NOTOK_RETURN(status=1)
    
    ! the just read values might define points on the cell boundaries,
    ! while lli should define the cell centers;
    ! fill correct values using nuv :
    select case ( nuv )
      case ( 'n' )
        ! sds attributes defined for lat/lon of center
        call Init( lli, lon_deg, dlon_deg, nlon, &
                        lat_deg, dlat_deg, nlat, status )
        IF_NOTOK_RETURN(status=1)
      case ( 'u' )
        ! sds attributes defined for lat/lon of east/west bound
        call Init( lli, lon_deg+0.5*dlon_deg, dlon_deg, nlon-1, &
                        lat_deg             , dlat_deg, nlat  , status )
        IF_NOTOK_RETURN(status=1)
      case ( 'v' )
        ! sds attributes defined for lat/lon of south/north bound
        call Init( lli, lon_deg             , dlon_deg, nlon  , &
                        lat_deg+0.5*dlat_deg, dlat_deg, nlat-1, status )
        IF_NOTOK_RETURN(status=1)
      case default
        write (gol,'("unsupported nuv `",a,"`")') nuv; call goErr
        TRACEBACK; status=1; return
    end select
    
    ! ok
    status = 0
    
  end subroutine lli_Init_mf
  
  
  ! ***
  

  !
  ! initialiase level info from sds
  !
  
  subroutine levi_Init_mf( levi, mf, status )
  
    use file_hdf, only : GetInfo, ReadAttribute
    use Grid, only : TLevelInfo, Init

    ! --- in/out ----------------------------------
    
    type(TLevelInfo), intent(out)        ::  levi
    type(TMeteoFile_tm5_hdf), intent(in)     ::  mf
    integer, intent(out)                 ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_Init_mf'
    
    ! --- local -----------------------------------
    
    integer                ::  data_rank
    integer                ::  lm
    real, allocatable      ::  a(:), b(:)
    
    ! --- begin ------------------------------------
    
    ! record not selected ? then return
    if ( .not. mf%selected ) then
      write (gol,'("no record selected ...")'); call goErr
      TRACEBACK; status = 1; return
    end if
    
    ! extract size of data array:
    call GetInfo( mf%sds, status, data_rank=data_rank )
    if ( status/= 0 ) then; TRACEBACK; status=1; return; end if
    
    ! 2D or 3D
    select case ( data_rank )
    
      case ( 2 )
      
        ! set dummy values ...
        call Init( levi, 1, (/0.0,0.0/), (/0.0,0.0/), status )
        if ( status/= 0 ) then; TRACEBACK; status=1; return; end if

      case ( 3 )
    
        ! read number of levels and hybride parameters:
        call ReadAttribute( mf%sds, 'lm', lm, status )

        ! extract hybride coeff
        allocate( a(lm+1), b(lm+1) )
        call ReadAttribute( mf%sds, 'at', a, status )
        IF_NOTOK_RETURN(status=1)
        call ReadAttribute( mf%sds, 'bt', b, status )
        IF_NOTOK_RETURN(status=1)

        ! fill ...
        call Init( levi, lm, a, b, status )
        IF_NOTOK_RETURN(status=1)

        ! done    
        deallocate( a, b )
        
      case default
      
        write (gol,'("unsupported data rank : ",i6)') data_rank; call goErr
        TRACEBACK; status=1; return

    end select
    
    ! ok
    status = 0
    
  end subroutine levi_Init_mf
    

  ! ***

  
  subroutine mf_GetField( mf, status, gridtype, ll, spm )
  
    use PArray, only : pa_SetShape
    use file_hdf, only : GetInfo, ReadData, ReadAttribute

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

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

    character(len=*), intent(out), optional   ::  gridtype
    real, pointer, optional                   ::  ll(:,:,:)
    real, pointer, optional                   ::  spm(:,:)

    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_GetField'
    
    ! --- local -------------------------------
    
    integer            ::  data_rank
    integer            ::  data_dims(3)

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

    ! record not selected ? then return
    if ( .not. mf%selected ) then
      write (gol,'("no record selected ...")'); call goErr
      TRACEBACK; status = 1; return
    end if
    
    ! initial data is not extracted ...
    status = 1
    
    ! return grid type ?
    if ( present(gridtype) ) then
      call ReadAttribute( mf%sds, 'gridtype', gridtype, status )
      IF_NOTOK_RETURN(status=1)
    end if
 
    ! return 3d data array ?
    if ( present(ll) ) then
      ! extract data rank and shape:
      call GetInfo( mf%sds, status, data_rank=data_rank, data_dims=data_dims )
      IF_NOTOK_RETURN(status=1)
      ! extract data array:
      select case ( data_rank )
        case ( 2 )
          data_dims(3) = 1
          call pa_SetShape( ll, data_dims )
          call ReadData( mf%sds, ll(:,:,1), status )
          IF_NOTOK_RETURN(status=1)
        case ( 3 )
          call pa_SetShape( ll, data_dims )
          call ReadData( mf%sds, ll, status )
          IF_NOTOK_RETURN(status=1)
        case default
          write (gol,'("unsupported data rank:",i6)') data_rank; call goErr
          TRACEBACK; status=1; return
      end select
      ! ReadData breaks on error ...
      status = 0
    end if

    ! return surface pressure array ?
    if ( present(spm) ) then
      if ( .not. mf%spm_selected ) then
        write (gol,'("no spm record selected ...")'); call goErr
        TRACEBACK; status = 1; return
      end if
      ! extract data rank and shape:
      call GetInfo( mf%spm_sds, status, data_rank=data_rank, data_dims=data_dims )
      IF_NOTOK_RETURN(status=1)
      ! extract data array:
      select case ( data_rank )
        case ( 2 )
          call pa_SetShape( spm, data_dims(1:2) )
          call ReadData( mf%spm_sds, spm, status )
          IF_NOTOK_RETURN(status=1)
        case default
          write (gol,'("unsupported data rank for spm:",i6)') data_rank; call goErr
          TRACEBACK; status=1; return
      end select
      ! ReadData breaks on error ...
      status = 0
    end if

    ! no arguments processed ? 
    ! probably subroutine is called incorrectly ...
    if ( status /= 0 ) then
      write (gol,'("no arguments processed; wrong call ?")'); call goErr
      TRACEBACK; status = 1; return
    end if
    
    ! ok
    status = 0

  end subroutine mf_GetField

  

  ! ***
  
  
  subroutine mf_ReadRecord( mf, paramkey, t1, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll, sp_ll, &
                                status )

    use PArray, only : pa_Done
    use GO, only : TDate
    use GO, only : GO_Timer_Start, GO_Timer_End
    use Grid, only : TllGridInfo, TLevelInfo
    use file_hdf_base, only : Init, Done   !needed on aster

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_hdf), intent(inout)  ::  mf
    character(len=*), intent(in)         ::  paramkey
    type(TDate), intent(in)              ::  t1, t2
    character(len=1), intent(in)         ::  nuv
    character(len=2), intent(out)        ::  gridtype
    type(TLevelInfo), intent(out)        ::  levi
    character(len=1), intent(in)         ::  nw
    type(TllGridInfo), intent(inout)     ::  lli
    real, pointer                        ::  ll(:,:,:)
    real, pointer                        ::  sp_ll(:,:)
    integer, intent(out)                 ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_ReadRecord'
    
    ! --- local -------------------------------
    
    logical               ::  exist
    
    real, allocatable     ::  mfw_n(:,:,:)
    
    real, allocatable     ::  sp_ll_n(:,:)
    integer               ::  i, j

    ! --- begin ---------------------------------
    
    ! start timing:
    call GO_Timer_Start( itim_readrecord, status )
    IF_NOTOK_RETURN(status=1)
    
    ! input ?
    if ( mf%io /= 'i' ) then
      write (gol,'("file should have been opened for input, but io=",a)') mf%io; call goErr
      TRACEBACK; status=1; return
    end if

    ! open for reading:
    ! if it is opened in Init, to many hdf files would be open:
    !
    inquire( file=trim(mf%fname), exist=exist )
    if ( .not. exist ) then
      write (gol,'("main file does not exist:")'); call goErr
      write (gol,'("  ",a)') trim(mf%fname); call goErr
      TRACEBACK; status=1; return
    end if
    !
    call Init( mf%hdf, trim(mf%fname), 'read', status )
    IF_NOTOK_RETURN(status=1)
    !
    call Init( mf%sds, status )
    IF_NOTOK_RETURN(status=1)
    !
    mf%selected = .false.

    ! open spm file if necessary;
    ! if it is opened in Init, to many hdf files would be open:
    if ( mf%spm_load ) then
      !
      if ( mf%spm_extr ) then
        !
        inquire( file=trim(mf%spm_fname), exist=exist )
        if ( .not. exist ) then
          write (gol,'("spm file does not exist:")'); call goErr
          write (gol,'("  ",a)') trim(mf%spm_fname); call goErr
          TRACEBACK; status=1; return
        end if
        !
        call Init( mf%spm_hdf, trim(mf%spm_fname), 'read', status )
        IF_NOTOK_RETURN(status=1)
        !
      end if
      !
      call Init( mf%spm_sds, status )
      IF_NOTOK_RETURN(status=1)
      !
      mf%spm_selected = .false.
    end if

    ! select records in hdf files:
    call mf_SelectRecord( mf, paramkey, t1, t2, status )
    if ( status/= 0 ) then; TRACEBACK; status=1; return; end if

    ! always regular lat/lon grid ..    
    gridtype = 'll'
    
    ! setup grid definition:
    call lli_Init_mf( lli, nuv, mf, status )
    if ( status/= 0 ) then; TRACEBACK; status=1; return; end if
    
    ! setup level definition:
    call levi_Init_mf( levi, mf, status )
    if ( status/= 0 ) then; TRACEBACK; status=1; return; end if

    ! fill data array:
    call mf_GetField( mf, status, ll=ll )
    if ( status/= 0 ) then; TRACEBACK; status=1; return; end if
    
    ! ***
    
    ! special treatment of pw fields ?
    if ( (paramkey == 'mfw') .or. (paramkey == 'pw') ) then
    
      ! add extra top level ?
      if ( size(ll,3) == levi%nlev ) then
        ! copy current pw ('n' levels)
        allocate( mfw_n(lli%im,lli%jm,levi%nlev) )
        mfw_n = ll
        ! renew target array
        deallocate( ll )
        allocate( ll(lli%nlon,lli%nlat,levi%nlev+1) )
        ! store new pw:
        ll(:,:,1:levi%nlev) = mfw_n
        ll(:,:,levi%nlev+1) = 0.0   ! kg/s
        ! clear
        deallocate( mfw_n )
      end if
      
      ! change flux direction ?
      if ( mf%mfw_redir ) ll = - ll     ! upwards (increasing level)
    
    end if

    ! ***

    ! special treatment of surface stress fields ?

    if ( ((paramkey=='ewss') .or. (paramkey=='nsss')) .and. mf%sstr_to_ewss_nsss ) then

      !   sstr**2  =  (sstr**2)/2 + (sstr**2)/2  =  nsss**2 + ewss**2
      !  thus  nsss = sqrt( (sstr**2)/2.0 )

      ll = sqrt( (ll**2) / 2.0 )

    end if
    
    ! ***
    
    ! fill surface pressure field ?
    if ( mf%spm_load .and. mf%spm_selected ) then

      ! read field:
      call mf_GetField( mf, status, spm=sp_ll )
      if ( status/= 0 ) then; TRACEBACK; status=1; return; end if
      
      ! convert from 'n' to 'u'/'v' ?
      if ( mf%spm_n_to_uv ) then
        select case ( nuv )
          case ( 'n' )
            ! no interpolation
          case ( 'u' )
            ! copy current sp ('n')
            allocate( sp_ll_n(lli%im,lli%jm) )
            sp_ll_n = sp_ll
            ! renew size of output array:
            deallocate( sp_ll )
            allocate( sp_ll(0:lli%im,lli%jm) )
            ! interpol:
            sp_ll(0,:) = 1.5*sp_ll_n(1,:) - 0.5*sp_ll_n(2,:)
            do i = 1, lli%im-1
              sp_ll(i,:) = 1.5*sp_ll_n(i,:) + 0.5*sp_ll_n(i+1,:)
            end do
            sp_ll(lli%im,:) = -0.5*sp_ll_n(lli%im-1,:) + 1.5*sp_ll_n(lli%im,:)
            ! clear
            deallocate( sp_ll_n )
          case ( 'v' )
            ! copy current sp ('n')
            allocate( sp_ll_n(lli%im,lli%jm) )
            sp_ll_n = sp_ll
            ! renew size of output array:
            deallocate( sp_ll )
            allocate( sp_ll(lli%im,0:lli%jm) )
            ! interpol:
            sp_ll(:,0) = 1.5*sp_ll_n(:,1) - 0.5*sp_ll_n(:,2)
            do j = 1, lli%jm-1
              sp_ll(:,j) = 1.5*sp_ll_n(:,j) + 0.5*sp_ll_n(:,j+1)
            end do
            sp_ll(:,lli%jm) = -0.5*sp_ll_n(:,lli%jm-1) + 1.5*sp_ll_n(:,lli%jm)
            ! clear
            deallocate( sp_ll_n )
          case default
            write (gol,'("unsupported nuv `",a,"`")') nuv; call goErr
            TRACEBACK; status=1; return
        end select
      end if
      
    else
      call pa_Done( sp_ll )
    end if
    
    ! close surface pressure file ?
    if ( mf%spm_load ) then
      call Done( mf%spm_sds, status )    
      IF_NOTOK_RETURN(status=1)
      !
      if ( mf%spm_extr ) then
        call Done( mf%spm_hdf, status )
        IF_NOTOK_RETURN(status=1)
      end if
    end if

    ! close field file:
    call Done( mf%sds, status )
    IF_NOTOK_RETURN(status=1)
    call Done( mf%hdf, status )
    IF_NOTOK_RETURN(status=1)
    
    ! end timing:
    call GO_Timer_End( itim_readrecord, status )
    IF_NOTOK_RETURN(status=1)
    
    ! ok
    status = 0

  end subroutine mf_ReadRecord
  

  ! ******************************************************************
  ! ***
  ! *** output
  ! ***
  ! ******************************************************************
  

  subroutine WriteHeader( mf, lli, status, levi )

    use Binas, only : grav, ae
    use Grid, only : TllGridInfo, TLevelInfo
    use file_hdf, only : WriteAttribute

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_hdf), intent(inout)      ::  mf
    type(TllGridInfo), intent(in)            ::  lli
    integer, intent(out)                     ::  status

    type(TLevelInfo), intent(in), optional   ::  levi

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/WriteHeader'
    
    ! --- begin ---------------------------------
    
    ! write header:
    call WriteAttribute( mf%hdf, 'fname'   , trim(mf%fname), status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute( mf%hdf, 'format'  , output_format, status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute( mf%hdf, 'gridtype', 'll'         , status )
    IF_NOTOK_RETURN(status=1)

    ! save first and last lon/lat (center) for use with HIPHOP
    call WriteAttribute( mf%hdf, 'lonmin', lli%lon_deg(1)       , status, knd=rknd )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute( mf%hdf, 'lonmax', lli%lon_deg(lli%nlon), status, knd=rknd )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute( mf%hdf, 'latmin', lli%lat_deg(1)       , status, knd=rknd )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute( mf%hdf, 'latmax', lli%lat_deg(lli%nlat), status, knd=rknd )
    IF_NOTOK_RETURN(status=1)

    ! other useful stuff ...
    call WriteAttribute( mf%hdf, 'grav'   , grav       , status, knd=rknd )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute( mf%hdf, 'ae'     , ae         , status, knd=rknd )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute( mf%hdf, 'area_m2', lli%area_m2, status, knd=rknd )
    IF_NOTOK_RETURN(status=1)

    ! level stuff
    if ( present(levi) ) then
      call WriteAttribute( mf%hdf, 'lm', levi%nlev, status, knd=iknd )
      IF_NOTOK_RETURN(status=1)
      call WriteAttribute( mf%hdf, 'at', levi%a   , status, knd=rknd )
      IF_NOTOK_RETURN(status=1)
      call WriteAttribute( mf%hdf, 'bt', levi%b   , status, knd=rknd )
      IF_NOTOK_RETURN(status=1)
    end if

    ! ok
    status = 0

  end subroutine WriteHeader


  ! ***
  
    
  subroutine WriteSdsHeader( sds, tmi, unit, tref, t1, t2, lli, nuv, status, &
                                  levi, nw, nlev )

    use Binas   , only : p_global
    use file_hdf, only : TSds
    use file_hdf, only : SetDim, WriteAttribute, Compress
    use GO      , only : TDate, Get, rTotal, operator(-), IsAnyDate
    use Grid    , only : TllGridInfo, TLevelInfo
    use tmm_info, only : TMeteoInfo

    ! --- in/out -----------------------------
    
    type(TSds), intent(inout)               ::  sds
    type(TMeteoInfo), intent(in)            ::  tmi
    character(len=*), intent(in)            ::  unit
    type(TDate), intent(in)                 ::  tref, t1, t2
    type(TllGridInfo), intent(in)           ::  lli
    character(len=1), intent(in)            ::  nuv
    integer, intent(out)                    ::  status
    type(TLevelInfo), intent(in), optional  ::  levi
    character(len=1), intent(in), optional  ::  nw
    integer, intent(in), optional           ::  nlev
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/WriteSdsHeader'
    
    ! --- local -----------------------------
    
    integer             ::  time6(6)
    integer             ::  dhour
    
    ! --- begin -----------------------------

    ! *** history
    
    ! write history of meteo field:
    if ( len_trim(tmi%history) < 1 ) then
      call WriteAttribute( sds, 'history', '-', status )
      IF_NOTOK_RETURN(status=1)
    else
      call WriteAttribute( sds, 'history', trim(tmi%history), status )
      IF_NOTOK_RETURN(status=1)
    end if
    
    ! *** unit
    
    ! write unit attribute    
    call WriteAttribute( sds, 'unit', unit, status )
    IF_NOTOK_RETURN(status=1)

    ! *** time
    
    ! time interval in hours
    if ( IsAnyDate(t1) .or. IsAnyDate(t2) ) then
      dhour = 0
    else
      dhour = nint(rTotal( t2 - t1, 'hour' ))
    end if
    
    ! write old time attributes
    call Get( t1, time6=time6 )
    call WriteAttribute( sds, 'idate', time6, status, knd=iknd )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute( sds, 'dthrs', dhour, status, knd=iknd )  ! hours
    IF_NOTOK_RETURN(status=1)
    
    ! write new time attributes
    call Get( tref, time6=time6 )
      call WriteAttribute( sds, 'tref', time6, status, knd=iknd )
      IF_NOTOK_RETURN(status=1)
    call Get( t1, time6=time6 )
      call WriteAttribute( sds, 'time1', time6, status, knd=iknd )
      IF_NOTOK_RETURN(status=1)
    call Get( t2, time6=time6 )
      call WriteAttribute( sds, 'time2', time6, status, knd=iknd )
      IF_NOTOK_RETURN(status=1)
    time6 = 0
    time6(4) = dhour
    call WriteAttribute( sds, 'dtime', time6, status, knd=iknd )
    IF_NOTOK_RETURN(status=1)
    
    ! *** grid
    
    select case ( nuv )
      case ( 'n' )
        ! sds dimensions:
        call SetDim( sds, 1-1, 'LON'   , 'deg', lli%lon_deg, status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call SetDim( sds, 2-1, 'LAT'   , 'deg', lli%lat_deg, status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        ! sds attributes:
        call WriteAttribute( sds, 'gridtype' , 'll'          , status )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lon_first', lli%lon_deg(1), status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lon_inc  ', lli%dlon_deg  , status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lon_n    ', lli%nlon      , status, knd=iknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lat_first', lli%lat_deg(1), status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lat_inc  ', lli%dlat_deg  , status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lat_n    ', lli%nlat      , status, knd=iknd )
        IF_NOTOK_RETURN(status=1)
      case ( 'u' )
        ! sds dimensions:
        call SetDim( sds, 1-1, 'LONP1' , 'deg', lli%blon_deg, status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call SetDim( sds, 2-1, 'LAT'   , 'deg', lli%lat_deg , status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        ! sds attributes:
        call WriteAttribute( sds, 'gridtype' , 'll', status )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lon_first', lli%blon_deg(0), status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lon_inc  ', lli%dlon_deg   , status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lon_n    ', lli%im+1       , status, knd=iknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lat_first', lli%lat_deg(1) , status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lat_inc  ', lli%dlat_deg   , status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lat_n    ', lli%nlat       , status, knd=iknd )
        IF_NOTOK_RETURN(status=1)
      case ( 'v' )
        ! sds dimensions:
        call SetDim( sds, 1-1, 'LON'   , 'deg', lli%lon_deg , status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call SetDim( sds, 2-1, 'LATP1' , 'deg', lli%blat_deg, status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        ! sds attributes:
        call WriteAttribute( sds, 'gridtype' , 'll', status )
        call WriteAttribute( sds, 'lon_first', lli%lon_deg(1) , status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lon_inc  ', lli%dlon_deg   , status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lon_n    ', lli%nlon       , status, knd=iknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lat_first', lli%blat_deg(0), status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lat_inc  ', lli%dlat_deg   , status, knd=rknd )
        IF_NOTOK_RETURN(status=1)
        call WriteAttribute( sds, 'lat_n    ', lli%jm+1       , status, knd=iknd )
        IF_NOTOK_RETURN(status=1)
      case default
        write (gol,'("unsupported nuv `",a,"`")') nuv; call goErr
        TRACEBACK; status=1; return
    end select
    
    ! *** levels
    
    if ( present(levi) ) then
      if ( .not. present(nw) ) then
        write (gol,'("optional levi requires nw")'); call goErr
        TRACEBACK; status=1; return
      end if
      ! sds dimensions:
      select case ( nw )
        case ( '*' )
          if ( present(nlev) ) then
            call SetDim( sds, 3-1, 'HYBRID_SELECTED', 'Pa' , levi%fp0(1:nlev), status, knd=rknd )
            IF_NOTOK_RETURN(status=1)
          end if
        case ( 'n' )
          call SetDim( sds, 3-1, 'HYBRID', 'Pa' , levi%fp0, status, knd=rknd )
          IF_NOTOK_RETURN(status=1)
        case ( 'w' )
          call SetDim( sds, 3-1, 'HYBRIDh', 'Pa', levi%p0, status, knd=rknd )
          IF_NOTOK_RETURN(status=1)
        case default
          write (gol,'("unsupported nw `",a,"` xxx")') nw; call goErr
          TRACEBACK; status=1; return
      end select
      ! sds attributes:
      call WriteAttribute( sds, 'lm', levi%nlev, status, knd=iknd )
      IF_NOTOK_RETURN(status=1)
      call WriteAttribute( sds, 'at', levi%a   , status, knd=rknd )
      IF_NOTOK_RETURN(status=1)
      call WriteAttribute( sds, 'bt', levi%b   , status, knd=rknd )
      IF_NOTOK_RETURN(status=1)
    end if
    
    ! *** data compression
    
    call Compress( sds, compression, status, deflate_level=deflate_level )
    IF_NOTOK_RETURN(status=1)
    
    ! *** end
    
    ! ok
    status = 0
    
  end subroutine WriteSdsHeader
  
  
  ! ***
  
  
  subroutine WriteStatus( mf, msg, status )

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_hdf), intent(inout)  ::  mf
    character(len=*), intent(in)         ::  msg
    integer, intent(out)                 ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/WriteStatus'
    
    ! --- local ------------------------------
    
    integer       ::  fu
    logical       ::  opened

    ! --- begin ---------------------------------
    
    ! select unused file unit:
    fu = 1234
    do
      inquire( unit=fu, opened=opened )
      if ( .not. opened ) exit
      fu = fu + 1
    end do

    ! open:
    open( fu, file=trim(mf%fname)//'.status', form='formatted', iostat=status )
    if (status/=0) then
      write (gol,'("opening status file:")'); call goErr
      write (gol,'("  file   : ",a)') trim(mf%fname); call goErr
      TRACEBACK; status=1; return
    end if

    ! write message:
    write (fu,'(a)',iostat=status) msg
    if (status/=0) then
      write (gol,'("writing status:")'); call goErr
      write (gol,'("  file   : ",a)') trim(mf%fname); call goErr
      write (gol,'("  msg    : ",a)') msg; call goErr
      TRACEBACK; status=1; return
    end if

    ! done:
    close( fu, iostat=status )
    if (status/=0) then
      write (gol,'("closing status file:")'); call goErr
      write (gol,'("  file   : ",a)') trim(mf%fname); call goErr
      TRACEBACK; status=1; return
    end if

    ! ok
    status = 0

  end subroutine WriteStatus
  

  ! ***
  
  
  subroutine AddOutputName( mf, name, status )

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_hdf), intent(inout)  ::  mf
    character(len=*), intent(in)         ::  name
    integer, intent(out)                 ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/AddOuputName'
    
    ! --- local ------------------------------
    
    integer       ::  iname

    ! --- begin ---------------------------------
    
    ! if not present yet, add name
    iname = 1
    do 
      ! add ?
      if ( iname > mf%output_nname ) then
        ! place to store ?
        if ( iname > size(mf%output_names) ) then
          write (gol,'("length of mf%output_names array too small:")'); call goErr
          do iname = 1, size(mf%output_names)
            write (gol,'("  ",i3," ",a)') iname, trim(mf%output_names(iname)); call goErr
          end do
          TRACEBACK; status=1; return
        end if
        ! long enough for name ?
        if ( len(mf%output_names(iname)) < len(name) ) then
          write (gol,'("length of mf%output_names too small:")'); call goErr
          write (gol,'("  len(mf%output_names(i)) : ",i4)') len(mf%output_names(iname)); call goErr
          write (gol,'("  len(mname)              : ",i4)') len(name); call goErr
          TRACEBACK; status=1; return
        end if
        ! store:
        mf%output_names(iname) = name
        ! increase counter:
        mf%output_nname = iname
        ! leave:
        exit
      end if
      ! found ? then leave loop
      if ( mf%output_names(iname) == name ) exit
      ! next index
      iname = iname + 1
    end do

    ! ok
    status = 0

  end subroutine AddOutputName
  

  ! ***
  
  
  subroutine mf_WriteRecord_2d( mf, tmi, paramkey, unit, tref, t1, t2, &
                                lli, nuv, ll, status )

    use GO      , only : TDate
    use Grid    , only : TllGridInfo
    use file_hdf, only : TSds, WriteData
    use file_hdf_base, only :  Init, Done   ! needed on aster
    use tmm_info, only : TMeteoInfo

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_hdf), intent(inout)  ::  mf
    type(TMeteoInfo), intent(in)         ::  tmi
    character(len=*), intent(in)         ::  paramkey, unit
    type(TDate), intent(in)              ::  tref, t1, t2
    type(TllGridInfo), intent(in)        ::  lli
    character(len=1), intent(in)         ::  nuv
    real, intent(in)                     ::  ll(:,:)
    integer, intent(out)                 ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_WriteRecord_2d'
    
    ! --- local ------------------------------
    
    type(TSds)    ::  sds

    ! --- begin ---------------------------------
    
    ! output ?
    if ( mf%io /= 'o' ) then
      write (gol,'("file should have been opened for output, but io=",a)') mf%io; call goErr
      TRACEBACK; status=1; return
    end if
    
    ! new or existing ?
    if ( .not. mf%output_initialised ) then
      ! open new file, destroy old:
      call Init( mf%hdf, trim(mf%fname), 'create', status )
      IF_NOTOK_RETURN(status=1)
      ! write file header:
      call WriteHeader( mf, lli, status )
      IF_NOTOK_RETURN(status=1)
      ! status new
      call WriteStatus( mf, 'in-progress', status )
      IF_NOTOK_RETURN(status=1)
      ! no records written yet:
      mf%output_nrec = 0
      ! now the file is initialised
      mf%output_initialised = .true.
    else
      ! re-open file:
      call Init( mf%hdf, trim(mf%fname), 'write', status )
      IF_NOTOK_RETURN(status=1)
    endif
    
      
    ! *** data set

    ! add record name:
    call AddOutputName( mf, paramkey, status )
    IF_NOTOK_RETURN(status=1)

    ! init data set:
    call Init( sds, mf%hdf, paramkey, shape(ll), 'real', status, knd=rknd_ds )
    IF_NOTOK_RETURN(status=1)

    ! write unit, time, and grid info:
    call WriteSdsHeader( sds, tmi, unit, tref, t1, t2, lli, nuv, status )
    IF_NOTOK_RETURN(status=1)
    
    ! write grid
    call WriteData( sds, ll, status )
    IF_NOTOK_RETURN(status=1)
    
    ! done:
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)
    
    ! *** completed ?
    
    ! next record has been written:
    mf%output_nrec = mf%output_nrec + 1
    
    ! completed ? then re-write status file:
    if ( mf%output_nrec == mf%output_ntrec*mf%output_nname ) then
      call WriteStatus( mf, 'completed', status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! *** close
    
    call Done( mf%hdf, status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0

  end subroutine mf_WriteRecord_2d
  

  ! ***
  
  
  subroutine mf_WriteRecord_3d( mf, tmi, spname, paramkey, unit, tref, t1, t2, &
                                lli, nuv, levi, nw, ps, ll, status, &
                                nlev )

    use GO      , only : TDate
    use Grid    , only : TllGridInfo, TLevelInfo
    use file_hdf, only : TSds, WriteData
    use file_hdf_base, only : Init, Done  !needed on aster
    use tmm_info, only : TMeteoInfo

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_hdf), intent(inout)  ::  mf
    type(TMeteoInfo), intent(in)         ::  tmi
    character(len=*), intent(in)         ::  spname, paramkey, unit
    type(TDate), intent(in)              ::  tref, t1, t2
    type(TllGridInfo), intent(in)        ::  lli
    character(len=1), intent(in)         ::  nuv
    type(TLevelInfo), intent(in)         ::  levi
    character(len=1), intent(in)         ::  nw
    real, intent(in)                     ::  ps(:,:)
    real, intent(in)                     ::  ll(:,:,:)
    integer, intent(out)                 ::  status
    
    integer, intent(in), optional        ::  nlev

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_WriteRecord_3d'
    
    ! --- local ------------------------------
    
    type(TSds)    ::  sds
    integer       ::  iname

    ! --- begin ---------------------------------
    
    ! output ?
    if ( mf%io /= 'o' ) then
      write (gol,'("file should have been opened for output, but io=",a)') mf%io; call goErr
      TRACEBACK; status=1; return
    end if
    
    ! new or existing ?
    if ( .not. mf%output_initialised ) then
      ! open new file, destroy old:
      call Init( mf%hdf, trim(mf%fname), 'create', status )
      IF_NOTOK_RETURN(status=1)
      ! write file header:
      call WriteHeader( mf, lli, status, levi )
      IF_NOTOK_RETURN(status=1)
      ! status new
      call WriteStatus( mf, 'in-progress', status )
      IF_NOTOK_RETURN(status=1)
      ! no records written yet:
      mf%output_nrec = 0
      ! now the file is initialised
      mf%output_initialised = .true.
    else
      ! re-open file:
      call Init( mf%hdf, trim(mf%fname), 'write', status )
      IF_NOTOK_RETURN(status=1)
    endif
    
      
    ! *** data set

    ! add record name:
    call AddOutputName( mf, paramkey, status )
    IF_NOTOK_RETURN(status=1)

    ! init data set:
    call Init( sds, mf%hdf, paramkey, shape(ll), 'real', status, knd=rknd_ds )
    IF_NOTOK_RETURN(status=1)

    ! write unit, time, and grid info:
    call WriteSdsHeader( sds, tmi, unit, tref, t1, t2, lli, nuv, status, levi, nw )
    IF_NOTOK_RETURN(status=1)
    
    ! write grid
    call WriteData( sds, ll, status )
    IF_NOTOK_RETURN(status=1)
    
    ! done:
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)
    
    ! *** surface pressure

    ! init data set:
    call Init( sds, mf%hdf, spname, shape(ps), 'real', status, knd=rknd_ds )
    IF_NOTOK_RETURN(status=1)

    ! write unit, time, and grid info:
    call WriteSdsHeader( sds, tmi, 'Pa', tref, t1, t2, lli, nuv, status )
    IF_NOTOK_RETURN(status=1)
    
    ! write grid
    call WriteData( sds, ps, status )
    IF_NOTOK_RETURN(status=1)
    
    ! done:
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)
    
    ! *** completed ?
    
    ! next record has been written:
    mf%output_nrec = mf%output_nrec + 1
    
    ! completed ? then re-write status file:
    if ( mf%output_nrec == mf%output_ntrec*mf%output_nname ) then
      call WriteStatus( mf, 'completed', status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! *** close

    call Done( mf%hdf, status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0

  end subroutine mf_WriteRecord_3d
  

end module tmm_mf_tm5_hdf